2021年认证杯SPSSPRO杯数学建模B题(第二阶段)依巴谷星表中的毕星团求解全过程文档及程序
2021年认证杯SPSSPRO杯数学建模
B题 依巴谷星表中的毕星团
原题再现：
依巴谷卫星（High Precision Parallax Collecting Satellite，缩写为 Hip-parcos），全称为“依巴谷高精度视差测量卫星”，是欧洲空间局发射的一颗天体测量卫星，用以精确测量恒星的视差和自行。通过视差可以推断出恒星距地球的距离。
毕星团位于金牛座，是离地球最近的疏散星团。其成员星在 300 个以上，有多颗肉眼可见的亮星。对毕星团的研究已经持续了许多年，包括确定它的距离，构建演化的模型，确认或排除成员，以及研究各成员星的特性等。依据依巴谷卫星的观测数据，我们可以以相当高的精度测量相关各星的距离和运动情况，以对毕星团进行更加精确的研究。
在依巴谷卫星的观测数据中，毕星团中的亮星平均视差在 22 毫角秒左右，意味着其平均距离在 45 秒差距左右。我们在依巴谷星表中选择了 2719颗恒星，选择的标准是视差在 20–25 毫角秒之间（也就是距离地球在 40–50秒差距之间），其中包括了许多毕星团的成员。
这个数据集有如下字段：
1. HIP：星体编号
2. Vmag：视星等
3. RA：赤经（度）
4. DE：赤纬（度）
5. Plx：视差角（毫角秒），1000/Plx 即为目标离观测点的距离（秒差距）
6. pmRA：恒星自行的 RA 分量（毫角秒/年）
7. pmDE：恒星自行的 DE 分量（毫角秒/年）
8. e_Plx：Plx 的测量误差（毫角秒）
9. B-V：恒星的色指数
第二阶段问题：在1869年，RichardA.Proctor观测到有一些距离毕星团相对较远的恒星，在空间中有着与毕星团相似的运动。此后的天文学家将这些恒星的集合称为毕宿星流。有人猜测这是一个更大的星团（被称为毕宿超级星团）在部分解体以后的遗迹，也有观点认为其中的大部分恒星来自不同的起源。
请你建立合理的数学模型，在依巴谷卫星的数据集中寻找毕宿星流的成员星。由于毕宿星流在空间中相对分散，要求对其成员星进行谨慎的界定。并请参考其赫–罗图来研究毕宿星流和毕星团的来源是否一致。
整体求解过程概述(摘要)
星流的研究一直是恒星的演化和银河系暗晕的形状等领域的重要研究对象。成员星辨认的准确性直接影响星团和星流基本物理参数的估计。本文主要建立的DBSCAN优化算法，基于依巴谷卫星的观测数据筛选出了毕宿星流成员星，并绘制了毕宿星流的赫罗图来研究毕宿星流和毕星团的来源是否一致。
针对问题一，需要在依巴谷卫星的数据集中寻找毕宿星流的成员星。首先，本文通过所提供的数据集中的e_Plx的测量误差对Plx视差角的数据进行误差处理，将不符合毕星团视差角的数据剔除。接着，通过恒星的自行的RA分量和DE分量建立二维坐标，运用DBSCAN 聚类算法初步筛选相似运动的恒星。接着，通过赤经、赤纬和视差角的相关数据并运用三维坐标转换公式计算出筛出的恒星在三维空间下的位置分布，在三维空间下运用DBSCAN算法进行第二次聚类。其中DBSCAN算法中的两个重要参数Eps和Minpts 通过对K-dist 图的一维聚类得出。最终筛选出具有相似运动的恒星为470颗。然后，通过赤经、赤纬和视差角的相关数据，对初步筛选出的470颗恒星进行二次筛选。利用DBSCAN聚类算法，在470颗恒星中筛选出了308颗毕星团成员星。最后，将所筛选出的毕星团成员星从初步筛选出的470颗恒星中剔除，最终得到162颗毕宿星流成员星。
为了保证本文筛选出的成员星的准确，本文通过计算所筛选出的毕星团成员星的平均视差数据与题干中所给的数据相比较，其误差不超过 0.4%。所以本文所筛选的毕星团成员星准确度极高，则本文所建立的算法有一定的可行性，可以间接证明所筛选出的毕宿星流成员星是可靠的。
针对问题二，需要参考赫罗图来研究毕宿星流和毕星团的来源是否一致。首先，将问题一中筛选出来的候选成员星的视星等通过普森公式求出各个侯选星成员星的绝对星等。接着，通过色指数与温度公式求出各个侯选星成员星的温度。最后以绝对星等为纵坐标，恒星表面温度为横坐标绘制出毕星团成员星的赫罗图和毕宿星流成员星的赫罗图。最后，将绘制的两张赫罗图绘制在同一坐标系内，通过观察发现两张赫罗图高度重合，所以可以证明毕宿星流和毕星团的来源一致。
本文在最后对模型的优缺点进行分析。
模型假设：
1、假设本题附件中所提供的观测数据具有准确性和可靠性；
2、假设中模型所涉及的成员星数量在短期内是不发生变化的，即不考虑星团内部成员的反复相撞造成星团成员缓慢地“蒸发”等非一般状况的发生；
3、恒星的自行或视向速度都满足高斯分布；
4、假设毕星团的每颗成员星的自行速度基本相同；
5、假设分析数据时的计算误差是可以忽略的，也包括小数的取舍。
问题分析
本题要求我们探索毕宿星流成员星的辨认，并参考其赫罗图来研究毕宿星流和毕星团的来源是否一致。由于恒星的位置分布不均，我们考虑对DBSCAN聚类算法进行改进。根据题干可知，毕宿星流为一些距离毕星团相对较远的恒星，在空间中有着与毕星团相似的运动。本文首先考虑运用恒星自行的RA分量和恒星自行的DE分量将恒星样本进行第一次聚类，得到运动相似的恒星。接着以毕星团为分布较为紧密的球状星团为依据，用赤经、赤纬和视差角三个量将已进行初步筛选的的恒星进行再次筛选，所筛选出的为毕星团成员星。最终剩余的恒星则为毕宿星流成员星。再依据候选成员星的Vmag视星等和B-V色指数通过普森公式和B-V色指数与温度的函数关系求出绝对星等与温度。最后以绝对星等为纵坐标，温度为横坐标绘制赫罗图。
论文缩略图：
程序代码：
function[IDX, isnoise]=DBSCAN(X, epsilon, MinPts) C=0; n=size(X, 1); IDX=zeros(n, 1); % 初始化全部为0，即全部为噪音点 D=pdist2(X, X); visited=false(n, 1); isnoise=false(n, 1); for i=1:n if ~visited(i) visited(i)=true; Neighbors=RegionQuery(i); if numel(Neighbors)<MinPts % X(i, :) is NOISE isnoise(i)=true; else C=C+1; ExpandCluster(i, Neighbors, C); end end end function ExpandCluster(i, Neighbors, C) IDX(i)=C; k=1; while true j=Neighbors(k); if ~visited(j) visited(j)=true; Neighbors2=RegionQuery(j); if numel(Neighbors2)>=MinPts Neighbors=[Neighbors Neighbors2]; %#ok end end if IDX(j)==0 IDX(j)=C; end k=k + 1; if k > numel(Neighbors) break; end end end function Neighbors=RegionQuery(i) Neighbors=find(D(i, :)<=epsilon); end end
function PlotClusterinResult(X, IDX) k=max(IDX); Colors=hsv(k); Legends= {
}
; for i=0:k Xi=X(IDX==i, :); if i~=0 Style='x'; MarkerSize=8; Color=Colors(i, :); Legends { end+1
}
=['Cluster #' num2str(i)]; else Style='o'; MarkerSize=6; Color=[0 0 0]; if ~isempty(Xi) Legends { end+1
}
='Noise'; end end if ~isempty(Xi) plot(Xi(:, 1), Xi(:, 2), Style, 'MarkerSize', MarkerSize, 'Color', Color); end hold on; end hold off; axis equal; grid on; legend(Legends); legend('Location', 'NorthEastOutside'); end
clc;
clear all pmra=xlsread(''); pmde=xlsread('); X=[pmra, pmde]; for t=1:1748 for k=1:1748 p(k, t)=sqrt((X(t, 1)-X(k, 1))^2+(X(t, 2)-X(k, 2))^2); end end P=sort(p);
%进行排序 A=1:1748;
figure('color', 'w'); subplot(4, 2, 1); plot(P(:, 1)/10, A, 'k');
ylabel('点数'); xlabel('距离'); title('1-dist') subplot(4, 2, 2); plot(P(:, 2)/10, A, 'k');
ylabel('点数'); xlabel('距离');
title('2-dist') subplot(4, 2, 3); plot(P(:, 3)/10, A, 'k');
ylabel('点数'); xlabel('距离');
title('3-dist') subplot(4, 2, 4); plot(P(:, 4)/10, A, 'k');
ylabel('点数'); xlabel('距离');
title('4-dist') subplot(4, 2, 5); plot(P(:, 5)/10, A, 'k');
ylabel('点数'); xlabel('距离');
title('5-dist') subplot(4, 2, 6); plot(P(:, 6)/10, A, 'k');
ylabel('点数'); xlabel('距离');
title('6-dist') subplot(4, 2, 7); plot(P(:, 7)/10, A, 'k');
ylabel('点数'); xlabel('距离');
title('7-dist') subplot(4, 2, 8); plot(P(:, 8)/10, A, 'k');
ylabel('点数'); xlabel('距离');