1.本发明属于天空背景亮度检测技术领域，尤其涉及一种日晕光度计检测系统及检测方法、系统、介质、设备。
背景技术：2.天空背景亮度是太阳周围区域散射光强度和日心强度的比值，是衡量地面日冕观测的重要参数。现代日晕光度计是在美国atst(advancedtechnology solar telescope)望远镜选址项目中设计和应用的，本发明将其引入中国西部太阳选址工作，用来测量天空背景亮度，因其易操作和便携性，尤其适合太阳选址工作，测量的精确度和数据采集的连续性也使其比传统的evans光度计更适合统计研究。
3.日晕光度计的数据处理比较复杂。从sbm的原始采集图像中，可以看到原始图像的一些特征，最里面的是经过nd4(中性减光片)减光之后的太阳像，中间的亮环是减光片边缘衍射光的离焦像，最外面的一个不完整的亮环是望远镜镜筒边缘衍射光的离焦像。三条阴影是用来固定nd4减光片的支架投影。天空背景亮度是天空区域的光强与日面光强之比，有效的天空区域是两个衍射环之间扣除减光片支架投影以后的区域。因为实际观测中缺乏精确的导星装置，不同观测时间太阳像的位置是有较大变化的。选址时的观测地点也经常变动，三个支架投影的位置也会随之改变。但是实际观测条件比较复杂，传统的方法不足以提取清晰的太阳边缘和nd4减光片边缘，导致有效天空区域和太阳区域的提取不够准确。因此，亟需一种新的日晕光度计检测方法。
4.通过上述分析，现有技术存在的问题及缺陷为：实际观测条件比较复杂，传统的方法不足以提取清晰的太阳边缘和nd4减光片边缘，导致有效天空区域和太阳区域的提取不够准确。
5.解决以上问题及缺陷的难度为：传统的任何基于强度变化的图像处理办法，均无法有效解决问题。通过调整图像处理参数，可以对短时间内的数据进行处理，然而这需要人为判断，不具备普适性，无法实现海量数据的自动化处理。故而我们采用了比较先进的相位一致性算法，并利用里斯变换的方法将其扩展应用到二维图像，以适应错综复杂的观测条件。
6.解决以上问题及缺陷的意义为：可以实现不同时间、不同观测地点、不同观测人员采集的海量数据的全自动化处理，无需依赖人为主观判断，极其适合选址工作中累积的大量数据的统计研究工作。
技术实现要素：7.针对现有技术存在的问题，本发明提供了一种日晕光度计检测系统及检测方法、系统、介质、设备。
8.本发明是这样实现的，一种日晕光度计检测方法，所述日晕光度计检测方法包括以下步骤：
9.步骤一，采用相位一致性方法从sbm图像上提取的特征结构边界；
10.步骤二，采用互相关的方法来分离出日面边缘和掩体边缘；
11.步骤三，根据太阳边缘点，采用最小二乘法拟合做圆方程的拟合，得到图像上的日心坐标、太阳半径、掩体中心坐标和掩体半径；
12.步骤四，将图像以掩体中心为原点进行极坐标展开，利用互相关方法确定掩体支架的位置；
13.步骤五，计算图像上扣除掩体边缘衍射环区域，和掩体支架投影区域之后的天空区域数值强度，以及日心附近区域的数值强度；
14.步骤六，计算不同太阳半径的天空背景亮度。
15.进一步，所述日晕光度计检测方法，还包括：采用互相关的办法从sbm图像特征结构边界中提取nd4减光片的边缘点，用最小二乘法拟合得到中心坐标和半径；以nd4减光片中心为原点，把图像从nd4边缘以外的一个环型区域进行极坐标展开，三个支架之间间隔120度；在极坐标中构造两个峰值间隔120度的高斯函数作为模板计算互相关，来确定支架投影的方向。至此，从sbm原始观测图像上把计算天空背景亮度所需要的日面区域和天空区域划分出来，分别计算天空区域的平均亮度和日面中心区域的亮度，二者之比就是衡量日冕观测条件的白昼天空背景亮度。
16.进一步，步骤一中，所述相位一致性方法，包括：
17.采用单演滤波的二维相位一致性方法进行图像的边界探测。
18.相位一致性原理认为，人眼感知的图像中的结构特征，集中在其傅立叶分量的相位一致性最高的地方，而非传统中认为的图像中强度梯度最大的地方。相位一致性的定义为：
[0019][0020]
一维信号的傅立叶分量，作为复矢量按照平行四边形法则加起来，得到一个局域能量，|e(x)|。如果所有傅立叶分量的相位角都相等，则pc(x)＝1，而若所有傅立叶分量的相位完全不相干，则pc(x)
→
0；沿用kovesi的方法，滤波函数用log
‑
gabor函数，该函数可以构造成任意带宽得到具有最小空间幅度的最优化的滤波函数，其频域中的传输函数为：
[0021][0022]
其中，ω0是中心频率，k/ω0控制滤波函数的带宽。由于实函数傅立叶变换的厄密共轭对称性，计算相位一致性时需要一对分别是偶对称和奇对称的滤波函数。对于二维图像而言，绕零频旋转对称的标量的奇对称滤波函数是不存在的，故采用单演滤波方法，利用里斯变换从矢量意义上得到了各项同性的二维奇对称滤波函数，使得整套局域相位振幅理论可以直接应用到二维图像上。
[0023]
里斯变换在频域的传输函数为：
[0024]
[0025][0026]
矢量满足作为希尔伯特变换在二维情况的推广，对应于解析信号的虚部，得到的三分量为：
[0027]
f
m
(x1，x2)＝(f(x1，x2)，(h1★
f)(x1，x2)，(h2★
f)(x1，x2))；
[0028]
就是一个二维的“解析信号”，称为单演信号，h1，h2是里斯变换在空间域中对应的滤波函数，
★
代表卷积。
[0029]
单演信号的三分量构成标准的球坐标，φ∈[0，2π)，θ∈[0，π)：
[0030]
f＝acos(φ)；
[0031]
(h1★
f)＝a sin(φ)cos(θ)；
[0032]
(h2★
f)＝a sin(φ)sin(θ)；
[0033]
则单演信号在(x1，x2)的局域振幅为：
[0034][0035]
局部相位为：
[0036][0037]
局部方位为：
[0038]
θ(x1，x2)＝tan
‑1(h2★
f，h1★
f)；
[0039]
采用log
‑
gabor函数作为多尺度小波变换的基函数。单演滤波里斯变换中，振幅和φ相位用来计算二维图像的相位一致性，θ相位代表图像的倾向性，参考经典的canny算法，利用所述倾向性做非局域极值抑制约束边界特征的提取。
[0040]
进一步，所述非局域极大抑制，包括：
[0041]
非局域极大抑制算法是图像边缘检测的后处理，对图像水平和垂直方向分别求一阶导数，通过两个方向的导数可以计算每个像素位置的“方向性”和“梯度强度”；然后按方向性，即两个方向导数的反正切划分8个主方向；在每个像素点沿着该处的主方向计算邻域两个位置的梯度强度，如果该像素的梯度强度值在这个邻域内不是一个局域极大值，就把它去掉，这称之为非局域极大抑制。相位一致性计算的是特征结构的方向性，而经典的边缘检测认为梯度强度最大的地方是特征结构，二者是互相垂直的，实际应用中需注意修正。
[0042]
滞后阈值是采用两个阈值，利用上限阈值找出所有可以确定的特征结构位置，寻找和这些位置相连通的点，如果这些连通的位置满足下限阈值，则同样被认为是特征结构，继续递归地寻找直到相连通的点处低于下限阈值为止。
[0043]
进一步，步骤二中，所述采用互相关的方法来分离出日面边缘，包括：
[0044]
通过观测时间计算当时的太阳角半径，根据ccd像元尺寸和望远镜焦距估算出图像的太阳像半径是多少像素；其中，所述互相关的模板根据太阳的位置参数和望远镜的仪器参数进行构造，所述标准的图像互相关方法还用于确定太阳边缘、nd4边缘和nd4支架位置。
[0045]
进一步，步骤三中，所述圆的最小二乘法拟合，包括：
[0046]
日冕图像的极坐标展开需要日心坐标和半径，利用非线性最小二乘法对圆方程做最佳拟合来得到这两个参数；其中，所述圆方程为：
[0047][0048]
其中，为待拟合的点。
[0049]
线性拟合为：
[0050][0051]
得到的参数称为代数拟合，它最小化数据点和代数方程的残差平方和，通常作为最佳拟合的初始参数。
[0052]
最佳拟合又称几何拟合，它最小化数据点和圆心距离的函数与半径的残差平方和。将圆方程写成如下形式：
[0053][0054]
几何距离函数为：
[0055][0056]
其中，是圆心坐标，是半径，是数据点的坐标。
[0057]
最佳拟合要求为：
[0058][0059]
最小，其中参数矢量为s是β
j
，j＝1，2，3的非线性函数，因此这一个非线性最小二乘法问题。
[0060]
进一步，所述日晕光度计检测方法，还包括：
[0061]
(1)坐标变换：直角坐标到极坐标；采用图像矩阵整体作为参数计算变换坐标，然后通过插值完成变换，计算效率较高；
[0062]
(2)坐标变换：极坐标到直角坐标；
[0063]
(3)瑞利分布
[0064]
相位一致性算法需要对噪声响应进行估计，令振幅的分布为|z|＝f(x；s)，假设噪声是高斯分布，频域中的实部和虚部为独立同的高斯分布，二者合成的振幅为瑞利分布：
[0065][0066]
其中，s>0是统计分布的展宽参数，又称为“mode”，x＝s刚好对应于分布的极大值：
[0067][0068]
在精度要求不高的情况下，可以通过中值来估计展宽参数，因为中值对应于概率密度累积分布的1/2，求解：
[0069][0070]
可得中值：
[0071][0072]
瑞利分布的平均值μ和标准偏差σ与展宽参数s之间的关系是：
[0073][0074][0075]
因此需要这个程序计算瑞利分布的展宽参数，从而得到平均值和标准偏差来估计噪声强度。
[0076]
本发明的另一目的在于提供一种应用所述的日晕光度计检测方法的日晕光度计检测方法系统，所述日晕光度计检测方法系统包括：
[0077]
特征结构边界提取模块，用于采用相位一致性方法从sbm图像上提取的特征结构边界；
[0078]
日面边缘与nd4掩体边缘分离模块，用于采用互相关的方法来分离出日面边缘和掩体边缘；
[0079]
圆方程拟合模块，用于根据太阳边缘点，采用最小二乘法拟合做圆方程的拟合，得到图像上的日心坐标和太阳半径；根据掩体边缘点，采用最小二乘法做圆方程的拟合，得到掩体边缘衍射环的中心坐标和半径。
[0080]
掩体支架定位模块，将图像以掩体中心为原点展开到极坐标，利用互相关方法确定与掩体刚性联结、相互间隔120度的三个掩体支架在图像上的投影位置。
[0081]
天空背景亮度计算模块，计算不同太阳半径处的有效天空区域的天空背景亮度参数。
[0082]
本发明的另一目的在于提供一种计算机设备，所述计算机设备包括存储器和处理器，所述存储器存储有计算机程序，所述计算机程序被所述处理器执行时，使得所述处理器执行如下步骤：
[0083]
采用相位一致性方法从sbm图像上提取的特征结构边界；
[0084]
采用互相关的方法来分离出日面边缘；
[0085]
根据太阳边缘点，采用最小二乘法拟合做圆方程的拟合，得到图像上的日心坐标和太阳半径。
[0086]
根据掩体边缘点，采用最小二乘法拟合做圆方程的拟合，得到图像上的掩体中心坐标和掩体半径。
[0087]
将图像以掩体中心为原点展开到极坐标，利用互相关方法确定与掩体刚性联结、相互间隔120度的三个掩体支架在图像上的投影位置。
[0088]
计算不同太阳半径处的有效天空区域的天空背景亮度参数。
[0089]
本发明的另一目的在于提供一种信息数据处理终端，所述信息数据处理终端用于实现所述的日晕光度计检测系统。
[0090]
结合上述的所有技术方案，本发明所具备的优点及积极效果为：本发明提供的日晕光度计检测方法，能够通过全自动的数据处理程序从原始观测图像中提取有效的区域，并从sbm原始观测图像上把计算天空背景亮度所需要的日面区域和天空区域划分出来，分别计算天空区域的平均亮度和日面中心区域的亮度，二者之比就是衡量日冕观测条件的白昼天空背景亮度。同时，本发明采用felsberg and sommer提出的单演滤波方法，利用里斯变换(riesz transform)从矢量意义上得到了各项同性的二维奇对称滤波函数，使得整套局域相位振幅理论可以直接应用到二维图像上。
附图说明
[0091]
为了更清楚地说明本发明实施例的技术方案，下面将对本发明实施例中所需要使用的附图做简单的介绍，显而易见地，下面所描述的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下还可以根据这些附图获得其他的附图。
[0092]
图1是本发明实施例提供的日晕光度计检测方法流程图。
[0093]
图2是本发明实施例提供的日晕光度计检测系统结构框图；
[0094]
图中：1、特征结构边界提取模块；2、日面边缘与nd4掩体边缘分离模块；3、圆方程拟合模块；4、掩体支架定位模块；5、天空背景亮度计算模块。
[0095]
图3是本发明实施例提供的日晕光度计采集图像示例示意图。
[0096]
图4是本发明实施例提供的用互相关确定nd4减光片支架位置的示意图。
[0097]
图4(a)是本发明实施例提供的太阳边缘和nd4减光片边缘的拟合结果示意图。
[0098]
图4(b)是本发明实施例提供的以nd4减光片中心做的极坐标展开示意图。
[0099]
图4(c)是本发明实施例提供的用两个高斯函数构造的互相关模板示意图。
[0100]
图4(d)是本发明实施例提供的互相关的结果示意图。
[0101]
图5是本发明实施例提供的边界探测示意图。
[0102]
图5(a)是本发明实施例提供的原始图像示意图。
[0103]
图5(b)是本发明实施例提供的传统方法(sobel 大津算法)示意图。
[0104]
图5(c)是本发明实施例提供的相位一致性示意图。
[0105]
图6是本发明实施例提供的单演滤波的意义示意图。
[0106]
图7是本发明实施例提供的2014年mnras文章中采用的传统的图像处理方法进行边界提取示意图，sobel 大津算法。
[0107]
图8是本发明实施例提供的互相关确定日面边缘示意图。
[0108]
图9是本发明实施例提供的用互相关确定nd4减光片支架位置示意图。
[0109]
图10是本发明实施例提供的采用相位一致性原理进行的边界提取示意图。
具体实施方式
[0110]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0111]
针对现有技术存在的问题，本发明提供了一种日晕光度计检测系统及检测方法，下面结合附图对本发明作详细的描述。
[0112]
如图1所示，本发明实施例提供的日晕光度计检测方法包括以下步骤：
[0113]
s101，采用相位一致性方法从sbm图像上提取的特征结构边界；
[0114]
s102，采用互相关方法提取日面边缘和掩体边缘；
[0115]
s103，得到图像上的日心坐标、太阳半径、掩体中心坐标、掩体半径；
[0116]
s104，将图像以掩体中心为原点进行极坐标展开，利用互相关方法确定掩体支架的位置；
[0117]
s105，计算图像上扣除掩体边缘衍射环区域，和掩体支架投影区域之后的天空区域数值强度，以及日心附近区域的数值强度；
[0118]
s106，计算不同太阳半径的天空背景亮度。
[0119]
掩体边缘是指nd4中性减光片在图像上减光区域边缘的衍射环
[0120]
如图2所示，本发明实施例提供的日晕光度计检测系统包括：
[0121]
特征结构边界提取模块1，用于采用相位一致性方法从sbm图像上提取的特征结构边界；
[0122]
日面边缘与nd4掩体边缘分离模块2，用于采用互相关的方法来分离出日面边缘和掩体边缘；
[0123]
圆方程拟合模块3，用于根据太阳边缘点，采用最小二乘法拟合做圆方程的拟合，得到图像上的日心坐标和太阳半径；根据掩体边缘点，采用最小二乘法做圆方程的拟合，得到掩体边缘衍射环的中心坐标和半径。
[0124]
掩体支架定位模块4，将图像以掩体中心为原点展开到极坐标，利用互相关方法确定与掩体刚性联结、相互间隔120度的三个掩体支架在图像上的投影位置。
[0125]
天空背景亮度计算模块5，计算不同太阳半径处的有效天空区域的天空背景亮度参数。
[0126]
下面结合实施例对本发明的技术方案作进一步描述。
[0127]
1、随着西部太阳选址工作的开展，逐渐累积了大量的sbm数据，为此本发明专门开
发了一种全自动的数据处理程序从原始观测图像中提取有效的区域。本发明使用的程序均采用idl(interactivedata language)编写，最初发表在2014年mnras杂志(https://doi.org/10.1093/mnras/stu1231)，后来对其中的边缘探测算法进行了重要改进，相关内容发表在2018年solar physics杂志上(https://doi.org/10.1007/s11207
‑
017
‑
1223
‑
z)。
[0128]
2、内容
[0129]
2.1主程序
[0130]
日晕光度计采集的图像例子见图3，首先采用相位一致性方法从sbm图像上提取的特征结构边界。日面边缘和nd4减光片边缘是在一起的，还包括一些可能的其他环型边界，来自镜筒衍射以及像差导致的鬼像。本发明采用互相关的方法来分离出日面边缘，互相关的模板根据太阳的位置参数和望远镜的仪器参数来进行构造。首先通过观测时间计算当时的太阳角半径，然后根据ccd像元尺寸和望远镜焦距估算出图像的太阳像半径是多少像素。因为sbm的图像上日面区域是经过nd4减光的，估算的太阳半径会比真实太阳半径略大，互相关并不能给出精确的太阳位置，但是正因为模板比真实太阳大，所以可以用该模板在相关最大的位置截取出太阳的边缘。然后本发明根据太阳边缘点，采用最小二乘法拟合做圆方程的拟合，得到图像上的日心坐标和太阳半径。
[0131]
日晕光度计图像上的亮环和支架投影在计算天空背景亮度时，是必须剔除的区域。亮环是减光片边缘衍射光的离焦像，因而和nd4减光片是同心的。而支架和nd4减光片刚性连接，因此这两种结构位置的确定，都需要知道nd4减光片的中心坐标和半径。太阳像在理想情况下应该在nd4减光片的中心，但实际观测中由于望远镜跟踪等问题，太阳中心和nd4减光片中心通常相差较大。本发明同样采用互相关的办法从sbm图像特征结构边界中提取nd4减光片的边缘点，用最小二乘法拟合得到中心坐标和半径。然后以nd4减光片中心为原点，把图像从nd4边缘以外的一个环型区域进行极坐标展开，三个支架之间近似间隔120度，所以本发明在极坐标中构造两个峰值间隔120度的高斯函数作为模板计算互相关，来确定支架投影的方向。上述过程参见图4。至此本发明就从sbm原始观测图像上把计算天空背景亮度所需要的日面区域和天空区域划分出来，分别计算天空区域的平均亮度和日面中心区域的亮度，二者之比就是衡量日冕观测条件的白昼天空背景亮度。
[0132]
图4中，为了方便，在计算互相关时本发明反转了图像的强度。
[0133]
2.2子程序
[0134]
2.2.1相位一致性
[0135]
图像的边界探测本发明采用了单演滤波的二维相位一致性方法，因为实际观测条件比较复杂，传统的方法不足以提取清晰的太阳边缘和nd4减光片边缘，导致有效天空区域和太阳区域的提取不够准确(见图5)。
[0136]
相位一致性原理认为，人眼感知的图像中的结构特征，集中在其傅立叶分量的相位一致性最高的地方，而非传统中认为的图像中强度梯度最大的地方。相位一致性(phase congruency)的定义为：
[0137][0138]
一维信号的傅立叶分量，作为复矢量按照平行四边形法则加起来，得到一个局域
能量(local energy)，|e(x)|。如果所有傅立叶分量的相位角都相等，则pc(x)＝1，而若所有傅立叶分量的相位完全不相干，则pc(x)
→
0。本发明沿用kovesi的方法，滤波函数用log
‑
gabor函数，该函数可以构造成任意带宽得到具有最小空间幅度的最优化的滤波函数，其频域中的传输函数为：
[0139][0140]
其中，ω0是中心频率，k/ω0控制滤波函数的带宽。由于实函数傅立叶变换的厄密共轭对称性，计算相位一致性时需要一对儿分别是偶对称和奇对称的滤波函数，对于二维图像而言，绕零频旋转对称的(标量的)奇对称滤波函数是不存在的，因此又采用felsberg and sommer提出的单演滤波方法，利用里斯变换(riesz transform)从矢量意义上得到了各项同性的二维奇对称滤波函数，使得整套局域相位振幅理论可以直接应用到二维图像上。
[0141]
里斯变换在频域的传输函数为：
[0142][0143][0144]
矢量满足作为希尔伯特变换在二维情况的推广，对应于解析信号的虚部。这样得到的三分量：
[0145]
f
m
(x1，x2)＝(f(x1，x2)，(h1★
f)(x1，x2)，(h2★
f)(x1，x2))；
[0146]
就是一个二维的“解析信号”，称为单演信号(monogenic signal)，这里h1，h2是里斯变换在空间域中对应的滤波函数，
★
代表卷积。单演信号的三分量构成了标准的球坐标，φ∈[0，2π)，θ∈[0，π)：
[0147]
f＝a cos(φ)
ꢀꢀ
(4)
[0148]
(h1★
f)＝a sin(φ)cos(θ)
ꢀꢀ
(5)
[0149]
(h2★
f)＝a sin(φ)sin(θ)
ꢀꢀ
(6)
[0150]
则单演信号在(x1，x2)的局域振幅(local amplitude)为：
[0151][0152]
局部相位(local phase)为：
[0153][0154]
局部方位(local orientation)为：
[0155]
θ(x1，x2)＝tan
‑1(h2★
f，h1★
f)
ꢀꢀ
(9)
[0156]
仍然采用log
‑
gabor函数作为多尺度小波变换的基函数，图6描述了单演滤波里斯变换的意义。其中的振幅和φ相位用来计算二维图像的相位一致性，而θ相位代表了图像的
倾向性，参考经典的canny算法，可以利用这种倾向性做非局域极值抑制(nonmaxima suppression)来进一步约束边界特征的提取。
[0157]
图6中，z方向对应于解析信号的实部，x
‑
y方向在一起对应于解析信号的虚部。滤波采用log
‑
gabor函数，图中的二维曲面为在一个特定的小波尺度上，里斯变换反演回空间域的滤波函数。
[0158]
2.2.2非局域极大抑制
[0159]
非局域极大抑制(nonmaxima suppression)算法是图像边缘检测的后处理，它对图像水平和垂直方向分别求一阶导数，通过两个方向的导数可以计算每个像素位置的“方向性”和“梯度强度”，然后按方向性(两个方向导数的反正切)划分8个主方向。在每个像素点沿着该处的主方向计算邻域两个位置的梯度强度，如果该像素的梯度强度值在这个邻域内不是一个局域极大值，就把它去掉，这称之为非局域极大抑制。但需要注意相位一致性计算的是特征结构的方向性，而经典的边缘检测认为梯度强度最大的地方是特征结构，二者是互相垂直的，实际应用中需注意修正。
[0160]
2.2.3滞后阈值
[0161]
滞后阈值(hysteresis thresholding)是采用两个阈值，首先用上限阈值找出所有可以确定的特征结构位置，然后寻找和这些位置相连通的点，如果这些连通的位置满足下限阈值，则同样被认为是特征结构，然后继续递归地寻找直到相连通的点处低于下限阈值为止。
[0162]
2.2.4图像互相关
[0163]
标准的图像互相关方法，用来确定太阳边缘，nd4边缘，nd4支架位置。
[0164]
2.2.5圆的最小二乘法拟合
[0165]
日冕图像的极坐标展开需要日心坐标和半径，我们利用非线性最小二乘法对圆方程做最佳拟合来得到这两个参数。圆方程为：
[0166][0167]
其中，为待拟合的点。线性拟合：
[0168][0169]
得到的参数称为代数拟合(algebraic fit)，它最小化数据点和代数方程的残差平方和，通常作为最佳拟合的初始参数。
[0170]
最佳拟合(best fit)又称几何拟合(geometric fit)，它最小化数据点和圆心距离的函数与半径的残差平方和。将圆方程写成如下形式：
[0171][0172]
几何距离函数为：
[0173][0174]
其中，是圆心坐标，是半径，是数据点的坐标。最佳拟合要求：
[0175][0176]
最小，其中参数矢量为显然s是β
j
，j＝1，2，3的非线性函数，因此这一个非线性最小二乘法问题。
[0177]
2.2.6坐标变换：直角坐标到极坐标
[0178]
采用图像矩阵整体作为参数计算变换坐标，然后通过插值完成变换，计算效率较高。
[0179]
2.2.7坐标变换：极坐标到直角坐标
[0180]
2.2.8瑞利分布
[0181]
降噪过程中需要对噪声响应进行估计，令振幅的分布为|z|＝f(x；s)，假设噪声是高斯分布，频域中的实部和虚部为独立同的高斯分布，二者合成的振幅为瑞利分布：
[0182][0183]
其中s>0是统计分布的展宽参数(scale parameter，决定着分布的宽度spread)，又称为“mode”，x＝s刚好对应于分布的极大值：
[0184][0185]
在精度要求不高的情况下，可以通过中值来估计展宽参数，因为中值对应于概率密度累积分布的1/2，求解：
[0186][0187]
可得中值：
[0188][0189]
瑞利分布的平均值μ和标准偏差σ与展宽参数s之间的关系是：
[0190][0191][0192]
因此需要这个程序计算瑞利分布的展宽参数，从而得到平均值和标准偏差来估计噪声强度。
[0193]
2.2.9索引数组
[0194]
2.2.10高斯函数
[0195]
创建互相关模板时使用的高斯函数。
[0196]
2.2.11画图程序
[0197]
2.2.12其他
[0198]
其他子程序headfits.pro，sxpar.pro，readfits.pro和get_rb0p.pro使用的是solar software(ssw)提供的库函数。前三个为fits格式数据的输入，最后一个是用来计算太阳位置参数。
[0199]
下面结合附图对本发明的技术效果作详细的描述。
[0200]
如图7所示，sbm图像的主要特征边界提取。(a)是对原始图像调整了动态范围使其特征边界更清晰。(b)是应用sobel算子计算的图像梯度。(c)是采用大津算法对梯度图像的直方图计算二值化阈值。(d)是根据这个阈值进行的二值化结果，并用2像素的结构元做了图像形态学腐蚀。
[0201]
如图8所示，用互相关从sbm图像主要特征边界中分离太阳边缘。白色实线为提取的特征边界，背景彩图是互相关结果，左下角是估算的互相关模板，箭头指向模板与图像边界相关最大的位置。
[0202]
如图9所示，用互相关确定nd4减光片支架位置。(a)给出了太阳边缘和nd4减光片边缘的拟合结果。(b)是以nd4减光片中心做的极坐标展开。(c)是用两个高斯函数构造的互相关模板。(d)给出了互相关的结果。为了方便，在计算互相关时反转了图像的强度。
[0203]
如图10所示，采用相位一致性原理进行的边界提取。边界提取的对象是(a)被受云影响的异常图像。图(a)是用单演滤波方法计算的相位一
‑
致性，(b)是对相位一
‑
致性进行非局域极值抑制后的结果。
[0204]
在上述实施例中，可以全部或部分地通过软件、硬件、固件或者其任意组合来实现。当使用全部或部分地以计算机程序产品的形式实现，所述计算机程序产品包括一个或多个计算机指令。在计算机上加载或执行所述计算机程序指令时，全部或部分地产生按照本发明实施例所述的流程或功能。所述计算机可以是通用计算机、专用计算机、计算机网络、或者其他可编程装置。所述计算机指令可以存储在计算机可读存储介质中，或者从一个计算机可读存储介质向另一个计算机可读存储介质传输，例如，所述计算机指令可以从一个网站站点、计算机、服务器或数据中心通过有线(例如同轴电缆、光纤、数字用户线(dsl)或无线(例如红外、无线、微波等)方式向另一个网站站点、计算机、服务器或数据中心进行传输)。所述计算机可读取存储介质可以是计算机能够存取的任何可用介质或者是包含一个或多个可用介质集成的服务器、数据中心等数据存储设备。所述可用介质可以是磁性介
质，(例如，软盘、硬盘、磁带)、光介质(例如，dvd)、或者半导体介质(例如固态硬盘solid state disk(ssd))等。
[0205]
以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，都应涵盖在本发明的保护范围之内。
