一种基于统计和多高斯的日面亮度分析方法与流程

本发明涉及一种基于统计和多高斯的日面亮度分析方法，属于天文技术和图像处理计算领域。

背景技术：

高斯混合模型(gaussianmixturemodel，gmm)是常用的拟合算法，现有技术中通过gmm来估计高斯分布的参数，从而直接在原始脏图中得到日面和背景的可能亮度值。在对脏图中的原始数据进行各像素点的统计结果中，会有一些异常值产生的峰值，对gmm产生影响，同时gmm对同一组数据同时拟合出多个高斯分布，不同高斯分布的数据互相影响，因此得到的对于日面亮度的分析不够准确。

技术实现要素：

本发明提供了一种基于统计和多高斯的日面亮度分析方法，将线性转换、最小二乘法结合起来，实现对muser原始脏图中日面和背景的可能亮度值分析。

本发明基于统计和多高斯的日面亮度分析方法是这样实现的：

首先对muser原始脏图进行一定组数的线性转换；接着对转换结果进行单高斯的最小二乘拟合，得到第一个高斯分布的均值μ1及标准差δ1，这个均值μ1即为背景亮度所在位置；然后根据三倍标准差原则，把介于(μ1-3δ1，μ1+3δ1)之间的数据全部剔除；接着对剩余数据进行单高斯的最小二乘拟合，得到第二个高斯分布的均值μ2及标准差δ2，这个均值μ2即为日面亮度所在位置；最后根据原始脏图中的最大值、最小值、组数、两个均值μ1和μ2计算得到日面和背景的可能亮度值。

所述基于统计和多高斯的日面亮度分析方法的具体方法如下：

参照日本野边山观测站数据处理程序中的源代码，对待处理的muser原始脏图进行一定组数下的线性转换，可以将组数的值定为2000。转换后，便于观察高斯分布情况，初步预测高斯分布数量，特别对于有太阳活动爆发的muser原始脏图，能够从转换结果中得到更多信息。

在完成第一次高斯拟合后，得到第一个高斯分布的均值μ1和标准差δ1，根据三倍标准差原则，把介于(μ1-3δ1,μ1+3δ1)之间的数据全部剔除，高斯分布中的3δ准则，又称为拉依达准则，是最常用的异常值判定与剔除准则。先假设一组数据只含有随机误差，对其进行计算处理得到标准偏差，按一定的概率确定一个区间，认为凡是超过这个区间的误差，就不属于随机误差而是粗大误差，含有该误差的数据应予以剔除。3δ原则的具体内容为：

(1)数值分布在(μ-δ,μ+δ)中的概率为0.6826；

(2)数值分布在(μ-2δ,μ+2δ)中的概率为0.9544；

(3)数值分布在(μ-3δ,μ+3δ)中的概率为0.9974。

可以认为，数据的取值几乎全部集中在(μ-3δ,μ+3δ)区间内，超出这个范围的取值概率不到0.3％，几乎不可能发生。根据此准则，可以在拟合出一个高斯分布之后，将该高斯分布所在区间的值全部剔除，进而做第二个高斯拟合，不同高斯分布之间互不影响，排除了第一次拟合出的高斯分布对第二次拟合的干扰。

本发明的有益效果是：

基于统计和多高斯的日面亮度分析方法，将线性转换、最小二乘法结合起来完成对脏图中日面和背景的可能亮度值分析，发挥了最小二乘拟合算法的优势，具备较高的扩展性与普适性。在进行线性转换时候，能够剔除部分异常数据对后续拟合效果的影响，同时能够更加方便地看出来高斯分布的情况；在拟合过程中，不同高斯分布间互不影响，相比同时拟合出多个高斯分布的拟合算法，如gmm，提高了拟合的准确性。

附图说明

图1是本发明基于统计和多高斯的日面亮度分析方法的总体流程图；

图2是实施例中的muser原始脏图；

图3(a)和(b)分别为经过500和1000组的线性转换的结果图；

图4(a)和(b)分别为经过2000和3000组的线性转换的结果图；

图5(a)和(b)分别为经过4000和5000组的线性转换的结果图；

图6是实施例中转换结果直方图；

图7是实施例中经过此方法的最终结果直方图及高斯分布曲线；

图8是利用依靠模型所得的结果进行太阳圆盘改正的结果图(muser最终脏图)。

具体实施方式

下面结合附图和具体实施例对本发明作进一步说明。

实施例1：如图1所示，一种基于统计和多高斯的日面亮度分析方法，线性转换和最小二乘的多高斯拟合的过程以此为开始，即首先在一定组数下对muser原始脏图进行线性转换；接着对转换结果进行单高斯的最小二乘拟合，得到第一个高斯分布的均值μ1及标准差δ1，这个均值μ1即为背景亮度所在位置；然后根据三倍标准差原则，把介于(μ1-3δ1,μ1+3δ1)之间的数据全部剔除；接着对剩余数据进行单高斯的最小二乘拟合，得到第二个高斯分布的均值μ2及标准差δ2，这个均值μ2即为日面亮度所在位置；最后根据原始脏图中的最大值、最小值、组数、两个均值μ1和μ2计算得到日面和背景的可能亮度值。

本实施例中的观测设备是明安图频谱射电日像仪(muser，mingantuspectralradioheliograph)，距离北京400km的国家天文台明安图观测站，是首次实现在厘米-分米波段上同时以高时间分辨率、高空间分辨率和高频率分辨率对太阳爆发活动进行观测的新一代太阳专用综合孔径射电成像望远镜。muser由低频阵和高频阵组成，其中包括40面4.5米口径的抛物面天线构成的观测频率范围为0.4ghz-2.0ghz的低频阵(muser-i)，60面2米口径的抛物面天线构成的观测频率范围为2.0ghz-15.0ghz的高频阵(muser-ii)。

此实施例中的实验数据是muser-i在2015年11月1日12时8分49秒354毫秒观测产生的原始观测数据经过数据格式转换后形成的fits数据文件20151101-120849_354161240.uvfits,其频率为1.7125ghz。通过muser数据处理系统执行相关操作生成muser原始脏图，如图2所示。

其中，组数的选择参考了日本野边山观测站数据处理程序中的源代码，我们将这一值带入muser进行多次试验之后，发现2000对于muser脏图处理同样适用，在本实施例中，分别对组数为500、1000、2000、3000、4000、5000对线性转换进行验证，如图3-5所示。在聚类数量超过2000以后，噪声增多；在聚类数量低于2000以前，峰数减少明显，不利于反映脏图中的像素点分布，故聚类数量确定为2000。

结合图1，所述基于统计和多高斯的日面亮度分析方法的步骤如下：

步骤1、根据组数对muser原始脏图进行线性转换，得到一个长度为2000的一维数组，画出直方图，如图6-7所示。转换后，数据保存在长度为2000的一维数组中。

步骤2、接着对转换结果进行单高斯的最小二乘拟合，得到第一个高斯分布的均值μ1及标准差δ1，这个均值μ1即为背景亮度所在位置。

此步骤结果为：μ1＝239.24468631，δ1＝120.621896786，峰值a1＝2557.23911985，d1＝0.0为渐近线到x轴的距离。

步骤3、接着根据三倍标准差原则，把介于(μ1-3δ1,μ1+3δ1)之间的数据全部剔除。

结果显示，剔除(μ1-3δ1,μ1+3δ1)后，剩余数据为1398个。

步骤4、对剩余数据进行单高斯的最小二乘拟合，得到第二个高斯分布的均值μ2及标准差δ2，这个均值μ2即为日面亮度所在位置。

此步骤结果为：μ2＝692.702436178，δ2＝206.655187192，峰值a2＝440.122883789，d2＝0.0为渐近线到x轴的距离。

步骤5、根据原始脏图中的最大值、最小值、组数bin、两个均值(μ1、μ2)计算得到日面和背景的可能亮度值。计算公式如下：

e＝(a-b)×μ1/d

f＝(a-b)×μ2/d

其中，e为背景亮度值，f为日面亮度值，a为原始脏图中的最大值即19586.0，b为原始脏图中的最小值即-3754，均值μ1为拟合得到的背景亮度值即239.24468631，均值μ2为拟合得到的日面亮度值即692.702436178，d为组数即2000。

结果显示，最终背景可能亮度值为：2791.98548924，日面可能亮度值为：8083.83743019。

在步骤2和步骤4中，都调用了scipy.optimize.leastsq函数(优化和拟合库中的最小二乘法)进行最小二乘的高斯拟合，其中高斯函数的标准形式为：a为峰值，e为常数，μ为均值，δ为标准差，d为渐近线。

最后，调用clean算法对模型进行准确性验证，结果如图8所示，经过改正后的太阳圆盘已经基本移到视场的中心，验证了拟合结果的准确性。

上面结合附图对本发明的具体实施方式作了详细说明，但是本发明并不限于上述实施方式，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下做出各种变化。