技术领域本发明涉及一种用于太阳高分辨图像序列中相邻两帧间的多尺度速度场测量方法，属于天文技术和图像技术领域。
背景技术：
太阳图像中，蕴含着各种尺度的活动现象，可能既包含大尺度的运动(强耀斑、CME、黑子、耀斑、暗条等)，又包含小尺度的运动(微耀斑、微暗条、埃勒曼炸弹以及存在于各种尺度的磁流体动力学波等)。叠加在大尺度上的小尺度变化会被大尺度运动所淹没，小尺度的速度场测量就变得非常困难，因而可以借助于多尺度方法分析速度场，这样就能够单独对每个不同分辨率的图像进行更为细致的分析。通过研究多尺度解决单一尺度无法精确测量速度场的问题，从而可以促进太阳磁场的研究。如果我们能从中定量提取更多的速度场有效信息，必将有利于研究各种太阳活动的演化规律，揭示物理本质从而提前预报太阳磁场活动，以更好地进行空间天气预报及对地球磁场的预报。目前，单一尺度下测量速度场计算大致可以分为两大类：基于相关性的方法以及基于梯度的方法。块匹配方法属于典型的基于相关性的方法，将图像序列的每一帧分成许多互不重叠的块，在某一给定的搜索区内根据一定的匹配准则找出最佳匹配块，对平移运动有较好的效果。Horn-Schunck(HS)以及Lucas-Kanade(LK)算法属于典型的基于梯度的方法，两种方法必须满足图像灰度一致性假设即相邻两幅图像上同一物体的灰度值保持不变，当对象运动速度较高或者大位移运动时，光流的基本等式会存在较大的误差；太阳图像各种现象都是非刚性运动，并且存在旋转以及大位移运动，传统的方法在解决小位移或者刚性运动中取得较好的结果，但是在大位移运动计算中存在缺陷。相邻两帧间的多尺度速度场测量方法，需要对两帧图像分别构建图像金字塔，然后分层对目标进行位移测量，目前常用的建立图像金字塔的方法主要有高斯卷积多尺度、小波多尺度以及数学形态学多尺度等，但一般都是结合传统的单一尺度的计算方法，与单一尺度的计算结果相比，效果有一定的改善，但是存在精度不高以及边界值测量不准确等问题。本发明正是为了解决这些问题而提出的一种用于太阳高分辨图像序列中相邻两帧间的多尺度速度场测量方法，该方法在加权局部相关跟踪技术(WLCT)的基础上，采用高斯卷积多尺度方法，精确地计算两帧图像间的速度场，解决了传统的单一尺度以及多尺度方法计算速度场不准确的问题。
技术实现要素：
本发明提供了一种用于太阳高分辨图像序列中相邻两帧间的多尺度速度场测量方法，该方法在加权局部相关跟踪技术(WLCT)的基础上，采用高斯卷积多尺度方法，更精确地计算两帧图像间的速度场，解决了传统的单一尺度以及传统的多尺度方法在计算大位移运动、目标旋转以及运动目标边界处速度场测量不准确的问题，本发明在很大程度上提高了速度场计算的精确度。本发明的技术方案是：一种用于太阳高分辨图像序列中相邻两帧间的多尺度速度场测量方法，首先，分别对相邻两帧图像扩展图像边界对图像进行预处理；其次，分别将两帧图像分解成两帧金字塔图像；接着，对两帧金字塔图像逐层分块，测量两帧间的块位移，再通过加权平滑和插值计算两帧间的逐点位移量；然后，叠加图像所有层的两帧间的位移量；最终得到两帧间的综合位移量及速度场。本发明在加权局部相关跟踪技术的基础上，采用高斯卷积多尺度方法，有效地提高了速度场的测量精度。所述用于太阳高分辨图像序列中相邻两帧间的多尺度速度场测量方法的具体步骤如下：步骤1：图像预处理：通过复制图像边界的像素得到边界被扩展的图像；步骤2：构建两帧金字塔图像：首先，对扩展后的两帧图像采用高斯滤波对图像进行平滑处理；然后，进行降采样得到第二层尺度图像；接着，对第二层尺度图像再经过高斯滤波以及降采样得到第三层尺度图像，依此类推，从而形成一组尺度逐级模糊的两帧图像，最终构成基于原始图像的一组两帧金字塔图像，其中最顶层为最粗层；步骤3：逐层图像分块：对上述构建的两帧金字塔图像逐层进行分块，将每一层图像划分成大小相同的子块；步骤4：逐层块位移测量：首先，在后一帧图像同一尺度层中确定搜索区域，通过归一化互相关函数，获得一个互相关平面，以互相关平面的最大值作为位移量，得到粗匹配区域；然后，对上述的粗匹配区域用归一化互相关函数再次获取最大值，以最大值为中心，设定阈值r为半径，以阈值r为半径截取一个圆，取圆中像素的最小值为阈值t，将圆内区域各个像素减去阈值t后计算其重心；不断的重复以上操作，直至遍历金字塔的每层图像所有的块，获取每一层所有块的位移；步骤5：逐点位移量计算：对两帧金字塔图像的每一层图像的所有块位移用块匹配的最大互相关系数作为权值进行加权平滑，然后进行三线性插值运算，得到平滑后每一层逐点的位移量；步骤6：逐层位移量叠加：首先通过后一帧金字塔图像最粗尺度层计算前一帧金字塔图像的最粗尺度层的运动矢量，得到矩阵A；然后通过后一帧金字塔图像的次高层计算前一帧次高层的运动矢量，对前一帧最高层的运动矢量A数值乘以2然后进行双线性插值得到前一帧次高层分辨率的图像矩阵B，利用矩阵B对前一帧次高层图像变形，即用前一帧金字塔图像的次高层减去前一帧金字塔图像最高层的逐点位移，利用加权局部相关跟踪以及单一尺度图像平滑计算残余的光流位移C，前一帧次高层最终的位移是B+C，利用这种方法反复循环，从金字塔的顶层一直到底层，保留前一帧原始图像每个像素对应的位移量得到前一帧图像的总位移量；步骤7：通过总位移量计算两帧间的速度场：根据上述得出的总位移量除以两帧图像的时间间隔，最终求得两帧间的速度场。优选地，所述步骤2中，对平滑处理后图像进行降采样即删除图像的偶数行与偶数列。优选地，所述步骤4中，选取阈值r时，通过高斯解析函数生成两幅偏移量已知的图像，通过半径r计算得出两幅图像的偏移量，当半径r取不同的值时，计算得出的两幅图像间的偏移量是不同的，由于高斯解析函数生成两幅偏移量是已知的称为标准量，将实际测量的偏移量与标准量相比，实际测量的偏移量越接近标准量，说明计算该实际测量的偏移量所选取的半径r是最佳的，就选此时的r值为阈值r。本发明的有益效果是：本发明能更加精确地计算两帧图像间的速度场，解决传统的单一尺度以及其他多尺度方法计算速度场不够准确的问题。附图说明图1是本发明太阳高分辨图像序列中相邻两帧间的多尺度速度场测量方法的总体流程图；图2是本发明中采用云南澄江抚仙湖1m新真空太阳望远镜(NewVacuumSolarTelescope,NVST)在TiO波段的NOAA11598观测数据高分辨图像序列中的第一张图；图3是本发明中对图2经预处理后边界扩展的图像；图4是本发明中对图3旋转5.12°的图像；图5是本发明中图3和图4模拟的标准速度场图；图6是本发明中对图3经高斯卷积多尺度方法构建的一组金字塔图像；图7是本发明中图3根据传统局部相关跟踪算法计算图像旋转5.12°得到的可视化速度场图像；图8是本发明中图3利用多尺度算法计算图像旋转5.12°得到的可视化速度场图像；图9是本发明中对图3放大1.1倍的图像；图10是本发明中图3和图9模拟的标准速度场图；图11是本发明中图3根据传统局部相关跟踪算法计算图像放大1.1倍得到的可视化速度场图像；图12是本发明中图3利用多尺度算法计算图像放大1.1倍得到的可视化速度场图像。具体实施方式实施例1：如图1-12所示，一种用于太阳高分辨图像序列中相邻两帧间的多尺度速度场测量方法，首先，分别对相邻两帧图像扩展图像边界对图像进行预处理；其次，分别将两帧图像分解成两帧金字塔图像；接着，对两帧金字塔图像逐层分块，测量两帧间的块位移，再通过加权平滑和插值计算两帧间的逐点位移量；然后，叠加图像所有层的两帧间的位移量；最终得到两帧间的综合位移量及速度场。所述用于太阳高分辨图像序列中相邻两帧间的多尺度速度场测量方法的具体步骤如下：步骤1：图像预处理：通过复制图像边界的像素得到边界被扩展的图像；步骤2：构建两帧金字塔图像：首先，对扩展后的两帧图像采用高斯滤波对图像进行平滑处理；然后，进行降采样得到第二层尺度图像；接着，对第二层尺度图像再经过高斯滤波以及降采样得到第三层尺度图像，依此类推，从而形成一组尺度逐级模糊的两帧图像，最终构成基于原始图像的一组两帧金字塔图像，其中最顶层为最粗层；步骤3：逐层图像分块：对上述构建的两帧金字塔图像逐层进行分块，将每一层图像划分成大小相同的子块；步骤4：逐层块位移测量：首先，在后一帧图像同一尺度层中确定搜索区域，通过归一化互相关函数，获得一个互相关平面，以互相关平面的最大值作为位移量，得到粗匹配区域；然后，对上述的粗匹配区域用归一化互相关函数再次获取最大值，以最大值为中心，设定阈值r为半径，以阈值r为半径截取一个圆，取圆中像素的最小值为阈值t，将圆内区域各个像素减去阈值t后计算其重心；不断的重复以上操作，直至遍历金字塔的每层图像所有的块，获取每一层所有块的位移；步骤5：逐点位移量计算：对两帧金字塔图像的每一层图像的所有块位移用块匹配的最大互相关系数作为权值进行加权平滑，然后进行三线性插值运算，得到平滑后每一层逐点的位移量；步骤6：逐层位移量叠加：首先通过后一帧金字塔图像最粗尺度层计算前一帧金字塔图像的最粗尺度层的运动矢量，得到矩阵A；然后通过后一帧金字塔图像的次高层计算前一帧次高层的运动矢量，对前一帧最高层的运动矢量A数值乘以2然后进行双线性插值得到前一帧次高层分辨率的图像矩阵B，利用矩阵B对前一帧次高层图像变形，即用前一帧金字塔图像的次高层减去前一帧金字塔图像最高层的逐点位移，利用加权局部相关跟踪以及单一尺度图像平滑计算残余的光流位移C，前一帧次高层最终的位移是B+C，利用这种方法反复循环，从金字塔的顶层一直到底层，保留前一帧原始图像每个像素对应的位移量得到前一帧图像的总位移量；步骤7：通过总位移量计算两帧间的速度场：根据上述得出的总位移量除以两帧图像的时间间隔，最终求得两帧间的速度场。优选地，所述步骤2中，对平滑处理后图像进行降采样即删除图像的偶数行与偶数列。优选地，所述步骤4中，选取阈值r时，通过高斯解析函数生成两幅偏移量已知的图像，通过半径r计算得出两幅图像的偏移量，当半径r取不同的值时，计算得出的两幅图像间的偏移量是不同的，由于高斯解析函数生成两幅偏移量是已知的称为标准量，将实际测量的偏移量与标准量相比，实际测量的偏移量越接近标准量，说明计算该实际测量的偏移量所选取的半径r是最佳的，就选此时的r值为阈值r。实施例2：如图1-12所示，一种用于太阳高分辨图像序列中相邻两帧间的多尺度速度场测量方法，首先，分别对相邻两帧图像扩展图像边界对图像进行预处理；其次，分别将两帧图像分解成两帧金字塔图像；接着，对两帧金字塔图像逐层分块，测量两帧间的块位移，再通过加权平滑和插值计算两帧间的逐点位移量；然后，叠加图像所有层的两帧间的位移量；最终得到两帧间的综合位移量及速度场。所述用于太阳高分辨图像序列中相邻两帧间的多尺度速度场测量方法的具体步骤如下：步骤1：图像预处理：通过复制图像边界的像素得到边界被扩展的图像(本例中图像采用云南澄江抚仙湖1m新真空太阳望远镜在TiO波段的NOAA11598观测数据，图2为序列图像中的第一幅图像，图3为图2边界扩展的图像，图4是根据图3旋转5.12°得到的图像，图5为图3与图4得到的模拟标准速度场)采用如下公式对图像作边界扩展：img2＝padarray(img1,padsize,'replicate'),padsize＝max_range*2(level-1)，其中，img2是扩展后的图像，img1是待处理图像，padsize为填充的行数与列数，replicate表示图像大小通过图像的边界值来扩展，max_range代表两张图像之间可能的最大位移，level是指金字塔层数。本发明中采用max_range＝10，level＝4，padarray是图像扩展函数。步骤2：构建两帧金字塔图像：首先，对扩展后的两帧图像采用高斯滤波对图像进行平滑处理，其中，经过边界扩展的图像构成第一层尺度(即最精细的尺度)，如图6(a)所示；高斯滤波是用高斯核对图像进行卷积操作，二维高斯核函数卷积可以分两步进行，首先将图像与一维高斯核函数进行卷积，然后将卷积结果与垂直方向的相同一维高斯核函数卷积。具体地，产生一维零均值高斯核函数做平滑滤波器是直接从高斯分布中计算模板值，函数表达式如下：g(x)=ce-x2/2δ2]]>其中，g(x)为需要产生的滤波器即整个模板的元素，c为规范化系数即对模板元素求和取倒数，x为滤波器水平方向上的点取值范围x＝-ceil(size/2):ceil(size/2)，size一维高斯核模板的大小，ceil为向上取整，δ为高斯函数的标准差，本发明中δ＝1.8。然后，对高斯滤波后的图像进行降采样得到第二层尺度图像，图像大小变为原来的四分之一，如图6(b)所示；接着，第二层尺度图像再经过高斯滤波以及降采样得到第三层尺度图像，如图6(c)所示，第三层尺度图像再经过高斯滤波以及降采样得到第四层尺度图像，如图6(d)所示，依此类推，从而形成一组尺度逐级模糊的两帧图像，最终构成基于原始图像的一组两帧金字塔图像，其中最顶层为最粗层；步骤3：逐层图像分块：对上述构建的两帧金字塔图像逐层进行分块，将每一层图像划分成大小相同的子块；具体的，从金字塔顶层图像f(x,y)(x＝1,2...M,y＝1,2...N,其中M代表图像的行数，N代表图像的列数)的起始坐标(x,y)开始分块，基本块的大小为m×n(m为块的行数，n为块的列数)，对图像f(x,y)按行优先的顺序进行重叠块划分，为了简化计算量，块与块之间的重叠行数与列数分别为d，可以划分为M/d×N/d个块，然后依次对金字塔其他图像层进行分块，块大小需要经实验后选择：若块过大，会影响估计的精度；过小则易受噪声的影响，而且运算量增加，经过多次实验在本发明中块的行数m为9，列数n为9，重叠的行数与列数d为4。步骤4：逐层块位移测量：首先，在后一帧图像同一尺度中确定搜索区域，假设每个图像块在X方向和Y方向的最大位移为max_range，在后一帧图像同一尺度层中的(M+2*max_range，N+2*max_range)的范围内进行初匹配，其中M，N为块图像的中心。接着，通过归一化互相关函数，获得一个互相关平面，以互相关平面的最大值作为位移量得到粗匹配区域。归一化互相关函数是在互相关函数的基础上借助施瓦茨不等式将得到的互相关矩阵归一化使结果介于[-1，1]之间。公式如下所示：r(x,y)=Σi=1mΣj=1n{[f(x+i,y+j)-f‾x,y]×[t(i,j)-t‾]