CN105787901A - 一种用于太阳高分辨图像序列中相邻两帧间的多尺度速度场测量方法 - Google Patents

Abstract

本发明涉及一种用于太阳高分辨图像序列中相邻两帧间的多尺度速度场测量方法，属于天文技术和图像技术领域。本发明包括步骤：首先，分别对相邻两帧图像扩展图像边界对图像进行预处理；其次，分别将两帧图像分解成两帧金字塔图像；接着，对两帧金字塔图像逐层分块，测量两帧间的块位移，再通过加权平滑和插值计算两帧间的逐点位移量；然后，叠加图像所有层的两帧间的位移量；最终得到两帧间的综合位移量及速度场。本发明能更精确的计算两帧图像间的速度场，有效地提高了速度场的测量精度。

Description

一种用于太阳高分辨图像序列中相邻两帧间的多尺度速度场测量方法

技术领域

本发明涉及一种用于太阳高分辨图像序列中相邻两帧间的多尺度速度场测量方法，属于天文技术和图像技术领域。

背景技术

太阳图像中，蕴含着各种尺度的活动现象，可能既包含大尺度的运动(强耀斑、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\left(x\right)={\mathrm{ce}}^{-x^2/2{\mathrm{\δ}}^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)=\frac{\underset{i=1}{\overset{m}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}\{\[f(x+i,y+j)-{\stackrel{\‾}{f}}_{x,y}\]\×\[t(i,j)-\stackrel{\‾}{t}\]\}}{{\left\{\underset{i=1}{\overset{m}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}{\[f(x+i,y+j)-{\stackrel{\‾}{f}}_{x,y}\]}^2\underset{i=1}{\overset{m}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}{\[t(i,j)-\stackrel{\‾}{t}\]}^2\right\}}^{0.5}}$

其中，r(x,y)是子图像f_x,y与模板t的归一化互相关值；f是待匹配的图像，大小为M×N(M,N分别为待匹配图像的行数和列数)；t是模板图像大小为m×n(m,n分别为模板图像的行数和列数)；f_x,y是从图像f中任意选取一块像素大小为m×n的子图，(x,y)是图像f左上角的坐标，坐标范围是1≤x≤M-m，1≤y≤N-n；t(i,j)为模板图像t(m,n)中任一点图像值；f(x+i,y+j)为待搜索区域f(m,n)中任一点图像值；为待搜索区域图像的平均值 为模板图像的平均值

最后，对上述的粗匹配区域用归一化互相关函数再次获取相关平面最大值，以最大值为中心，设定阈值r为半径(在本发明中我们经过多次实验r取值为2)，取圆中的最小值为阈值t。将圆内区域减去阈值t后计算其重心。不断的重复以上操作，直至遍历金字塔的每层图像所有的块，获取每一层所有块的位移。

重心的函数表达式如下：

$\begin{array}{cc}{X}_0=\frac{\underset{i}{\overset{M}{\Σ}}\underset{j}{\overset{N}{\Σ}}{x}_i{I}_{ij}}{\underset{i}{\overset{M}{\Σ}}\underset{j}{\overset{N}{\Σ}}{I}_{ij}}& Y_0=\frac{\underset{i}{\overset{M}{\Σ}}\underset{j}{\overset{N}{\Σ}}{y}_i{I}_{ij}}{\underset{i}{\overset{M}{\Σ}}\underset{j}{\overset{N}{\Σ}}{I}_{ij}}\end{array}$

其中，X₀是重心的垂直坐标，Y₀是重心的水平坐标，M和N分别代表圆面的纵轴坐标和横轴坐标的范围；I_ij为(i,j)坐标处的灰度强度；x_iy_i为(i,j)处的坐标。

步骤5：逐点位移量计算：对两帧金字塔图像的每一层图像的所有块位移首先进行中值滤波，去除孤立点位移。

中值滤波是把局部区域的像素按灰度等级进行排序，取该领域中灰度的中值作为当前像素的灰度值。中值滤波不仅可以去除孤点，而且可以保持图像的边缘特性，不会使图像产生显著的模糊。

对于一幅图像的像素矩阵，取以目标像素为中心的一个子矩阵窗口，对窗口内的像素灰度排序，取中间一个值作为目标象素的新灰度值。本发明中，窗口的大小为3*3。

接着，用块匹配的最大互相关系数作为权值对块位移图像进行加权平滑。

一幅M×N的图像经过一个大小为m×n(m和n是奇数)的加权均值滤波器的公式如下：

$g(x,y)=\frac{\underset{s=-a}{\overset{a}{\Σ}}\underset{t=-b}{\overset{b}{\Σ}}w(s,t)f(x+s,y+t)}{\underset{s=-a}{\overset{a}{\Σ}}\underset{t=-b}{\overset{b}{\Σ}}w(s,t)}$

其中，g(x,y)为滤波器系数与由该滤波器包围的图像像素的乘积之和，w(s,t)为在点(s,t)处的滤波器的系数，本发明滤波器系数采用的是块图像匹配时的最大相关系数，f(x+s,y+t)为在点(x+s,y+t)图像的像素值，分母是滤波器的各系数之和；

最后，每一层图像所有块位移加权平滑之后需要进行三线性插值，因为在进行位移测量时并不是逐点进行，而是每隔4个点进行测量，所以需要进行三线性插值得到平滑后每一层逐点的位移量。

步骤6：逐层位移量叠加：首先通过后一帧金字塔图像最粗尺度层计算前一帧金字塔图像的最粗尺度层的运动矢量，得到矩阵A；然后通过后一帧金字塔图像的次高层计算前一帧次高层的运动矢量，对前一帧图像最高层的运动矢量A乘以2然后进行双线性插值得到前一帧次高层分辨率的图像矩阵B，利用矩阵B对前一帧次高层图像变形，即用前一帧金字塔图像的次高层减去前一帧金字塔图像最高层的逐点位移，利用加权局部相关跟踪以及单一尺度图像平滑计算残余的光流位移C，前一帧次高层最终的位移是B+C，利用这种方法反复循环，从金字塔的顶层一直到底层，保留前一帧原始图像每个像素对应的位移量得到前一帧图像的总位移量。

步骤7：通过总位移量计算两帧间的速度场：根据上述得出的总位移量除以两帧图像的时间间隔，最终求得两帧间的速度场，如图8所示。

图7是图3根据传统局部相关跟踪算法计算图像旋转5.12°得到的可视化速度场图像，图8是图3利用本方法(在局部相关跟踪技术的基础上利用多尺度算法)计算图像旋转5.12°得到的可视化速度场图像，通过主观的识别，很显然图8计算速度场的精度，要明显好于图7，图7除了图像中心速度场相对准确，中心外的速度计算的杂乱无章。接下来，我们通过误差测量来比较图7与图8的计算精度。我们主要采用平均角误差以及标准角误差对数据进行误差分析，其中平均角误差(AAE)表示整个速度场的平均角误差，反映了使用某种算法计算的速度场整体偏离标准速度场的程度，标准角误差反映了计算的速度场中平均角误差的波动情况。从表1中可以看出，我们采用多尺度速度场测量方法得出的计算精度要明显高于单一尺度速度场的测量结果。

表1

数据	AAE平均角误差	标准角误差
			图7(单一尺度)	0.97531°	0.90840°
图8(多尺度速度场测量)	0.07273°	0.23560°

实施例3：如图1-12所示，一种用于太阳高分辨图像序列中相邻两帧间的多尺度速度场测量方法，首先，分别对相邻两帧图像扩展图像边界对图像进行预处理；其次，分别将两帧图像分解成两帧金字塔图像；接着，对两帧金字塔图像逐层分块，测量两帧间的块位移，再通过加权平滑和插值计算两帧间的逐点位移量；然后，叠加图像所有层的两帧间的位移量；最终得到两帧间的综合位移量及速度场。

所述用于太阳高分辨图像序列中相邻两帧间的多尺度速度场测量方法的具体步骤如下：

步骤1：图像预处理：通过复制图像边界的像素得到边界被扩展的图像(本例中图像采用云南澄江抚仙湖1m新真空太阳望远镜在TiO波段的NOAA11598观测数据，图2为序列图像中的第一幅图像，图3为图2边界扩展的图像，图9是根据图3放大1.1倍得到的图像，图10为图3与图9得到的模拟标准速度场)采用如下公式对图像作边界扩展：

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：构建两帧金字塔图像：首先，对扩展后的图像采用高斯滤波对图像进行平滑处理，其中，经过边界扩展的图像f(x,y)构成第一层尺度(即最精细的尺度)，如图6(a)所示；

高斯滤波是用高斯核对图像进行卷积操作，二维高斯核函数卷积可以分两步进行，首先将图像与一维高斯核函数进行卷积，然后将卷积结果与方向垂直的相同一维高斯核函数卷积。

具体地，产生一维零均值高斯核函数做平滑滤波器是直接从高斯分布中计算模板值，函数表达式如下：

$g\left(x\right)={\mathrm{ce}}^{-x^2/2{\mathrm{\δ}}^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)=\frac{\underset{i=1}{\overset{m}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}\{\[f(x+i,y+j)-{\stackrel{\‾}{f}}_{x,y}\]\×\[t(i,j)-\stackrel{\‾}{t}\]\}}{{\left\{\underset{i=1}{\overset{m}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}{\[f(x+i,y+j)-{\stackrel{\‾}{f}}_{x,y}\]}^2\underset{i=1}{\overset{m}{\Σ}}\underset{j=1}{\overset{n}{\Σ}}{\[t(i,j)-\stackrel{\‾}{t}\]}^2\right\}}^{0.5}}$

其中，r(x,y)是子图像f_x,y与模板t的归一化互相关值；f是待匹配的图像，大小为M×N(M,N分别为待匹配图像的行数和列数)；t是模板图像大小为m×n(m,n分别为模板图像的行数和列数)；f_x,y是从图像f中任意选取一块像素大小为m×n的子图，(x,y)是图像f左上角的坐标，坐标范围是1≤x≤M-m，1≤y≤N-n；t(i,j)为模板图像t(m,n)中任一点图像值；f(x+i,y+j)为待搜索区域f(m,n)中任一点图像值；为待搜索区域图像的平均值 为模板图像的平均值

最后，对上述的粗匹配区域用归一化互相关函数再次获取相关平面最大值，以最大值为中心，设定阈值r为半径(在本发明中我们经过多次实验r取值为2)，取圆中的最小值为阈值t。将圆内区域减去阈值t后计算其重心。不断的重复以上操作，直至遍历金字塔的每层图像所有的块，获取每一层所有块的位移。

重心的函数表达式如下：

$\begin{array}{cc}{X}_0=\frac{\underset{i}{\overset{M}{\Σ}}\underset{j}{\overset{N}{\Σ}}{x}_i{I}_{ij}}{\underset{i}{\overset{M}{\Σ}}\underset{j}{\overset{N}{\Σ}}{I}_{ij}}& Y_0=\frac{\underset{i}{\overset{M}{\Σ}}\underset{j}{\overset{N}{\Σ}}{y}_i{I}_{ij}}{\underset{i}{\overset{M}{\Σ}}\underset{j}{\overset{N}{\Σ}}{I}_{ij}}\end{array}$

其中，X₀是重心的垂直坐标，Y₀是重心的水平坐标，M和N分别代表圆面的纵轴坐标和横轴坐标的范围；I_ij为(i,j)坐标处的灰度强度；x_iy_i为(i,j)处的坐标。

步骤5：逐点位移量计算：对两帧金字塔的每一层图像的所有块位移首先进行中值滤波，去除孤立点位移。

中值滤波是把局部区域的像素按灰度等级进行排序，取该领域中灰度的中值作为当前像素的灰度值。中值滤波不仅可以去除孤点，而且可以保持图像的边缘特性，不会使图像产生显著的模糊。

对于一幅图像的像素矩阵，取以目标像素为中心的一个子矩阵窗口，对窗口内的像素灰度排序，取中间一个值作为目标象素的新灰度值。本发明中，窗口的大小为3*3。

接着，用块匹配的最大互相关系数作为权值对块位移图像进行加权平滑。

一幅M×N的图像经过一个大小为m×n(m和n是奇数)的加权均值滤波器的公式如下：

$g(x,y)=\frac{\underset{s=-a}{\overset{a}{\Σ}}\underset{t=-b}{\overset{b}{\Σ}}w(s,t)f(x+s,y+t)}{\underset{s=-a}{\overset{a}{\Σ}}\underset{t=-b}{\overset{b}{\Σ}}w(s,t)}$

其中，g(x,y)为滤波器系数与由该滤波器包围的图像像素的乘积之和，w(s,t)为在点(s,t)处的滤波器的系数，本发明滤波器系数采用的是块图像匹配时的最大相关系数，f(x+s,y+t)为在点(x+s,y+t)图像的像素值，分母是滤波器的各系数之和。

最后，每一层图像所有块位移加权平滑之后需要进行三线性插值，因为在进行位移测量时并不是逐点进行，而是每隔4个点进行测量，所以需要进行三线性插值得到平滑后每一层逐点的位移量。

步骤6：逐层位移量叠加：首先通过后一帧金字塔图像最粗尺度层计算前一帧金字塔图像的最粗尺度层的运动矢量，得到矩阵A；然后通过后一帧金字塔图像的次高层计算前一帧次高层的运动矢量，对前一帧最高层的运动矢量A数值乘以2然后进行双线性插值得到前一帧次高层分辨率的图像矩阵B，利用矩阵B对前一帧次高层图像变形，即用前一帧金字塔图像的次高层减去前一帧金字塔图像最高层的逐点位移，利用加权局部相关跟踪以及单一尺度图像平滑计算残余的光流位移C，前一帧次高层最终的位移是B+C，利用这种方法反复循环，从金字塔的顶层一直到底层，保留前一帧原始图像每个像素对应的位移量得到前一帧图像的总位移量。

步骤7：通过总位移量计算两帧间的速度场：根据上述得出的总位移量除以两帧图像的时间间隔，最终求得两帧间的速度场，如图12所示。图11是图3根据传统局部相关跟踪算法计算图像放大1.1倍得到的可视化速度场图像，图12是图3利用本方法(在加权局部相关跟踪技术的基础上利用多尺度算法)计算图像放大1.1倍得到的可视化速度场图像，通过主观的识别，很显然图12计算速度场的精度，要明显好于图11，图11除了图像中心速度场相对准确，中心外的速度计算的杂乱无章。接下来，我们通过误差测量来比较图11与图12的计算精度。我们主要采用平均角误差以及标准角误差对数据进行误差分析，其中平均角误差(AAE)表示整个速度场的平均角误差，反映了使用某种算法计算的速度场整体偏离标准速度场的程度，标准角误差反映了计算的速度场中平均角误差的波动情况。从表2中可以看出，我们采用多尺度速度场测量方法得出的计算精度要明显高于单一尺度速度场的测量结果。

表2

实验方法对比	平均角误差(AAE)	标准角误差(STD)
			图11(单一尺度速度场测量)	0.84264°	0.84298°
图12(多尺度速度场测量)	0.19426°	0.47944°

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

Claims (4)

1.一种用于太阳高分辨图像序列中相邻两帧间的多尺度速度场测量方法，其特征在于：首先，分别对相邻两帧图像扩展图像边界对图像进行预处理；其次，分别将两帧图像分解成两帧金字塔图像；接着，对两帧金字塔图像逐层分块，测量两帧间的块位移，再通过加权平滑和插值计算两帧间的逐点位移量；然后，叠加图像所有层的两帧间的位移量；最终得到两帧间的综合位移量及速度场。

2.根据权利要求1所述的用于太阳高分辨图像序列中相邻两帧间的多尺度速度场测量方法，其特征在于：所述用于太阳高分辨图像序列中相邻两帧间的多尺度速度场测量方法的具体步骤如下：

步骤1：图像预处理：通过复制图像边界的像素得到边界被扩展的图像；

步骤2：构建两帧金字塔图像：首先，对扩展后的两帧图像采用高斯滤波对图像进行平滑处理；然后，进行降采样得到第二层尺度图像；接着，对第二层尺度图像再经过高斯滤波以及降采样得到第三层尺度图像，依此类推，从而形成一组尺度逐级模糊的两帧图像，最终构成基于原始图像的一组两帧金字塔图像，其中最顶层为最粗层；

步骤3：逐层图像分块：对上述构建的两帧金字塔图像逐层进行分块，将每一层图像划分成大小相同的子块；

步骤4：逐层块位移测量：首先，在后一帧图像同一尺度层中确定搜索区域，通过归一化互相关函数，获得一个互相关平面，以互相关平面的最大值作为位移量，得到粗匹配区域；然后，对上述的粗匹配区域用归一化互相关函数再次获取最大值，以最大值为中心，设定阈值r为半径，以阈值r为半径截取一个圆，取圆中像素的最小值为阈值t，将圆内区域各个像素减去阈值t后计算其重心；不断的重复以上操作，直至遍历金字塔的每层图像所有的块，获取每一层所有块的位移；

步骤5：逐点位移量计算：对两帧金字塔图像的每一层图像的所有块位移用块匹配的最大互相关系数作为权值进行加权平滑，然后进行三线性插值运算，得到平滑后每一层逐点的位移量；

步骤6：逐层位移量叠加：首先通过后一帧金字塔图像最粗尺度层计算前一帧金字塔图像的最粗尺度层的运动矢量，得到矩阵A；然后通过后一帧金字塔图像的次高层计算前一帧次高层的运动矢量，对前一帧最高层的运动矢量A数值乘以2然后进行双线性插值得到前一帧次高层分辨率的图像矩阵B，利用矩阵B对前一帧次高层图像变形，即用前一帧金字塔图像的次高层减去前一帧金字塔图像最高层的逐点位移，利用加权局部相关跟踪以及单一尺度图像平滑计算残余的光流位移C，前一帧次高层最终的位移是B+C，利用这种方法反复循环，从金字塔的顶层一直到底层，保留前一帧原始图像每个像素对应的位移量得到前一帧图像的总位移量；

步骤7：通过总位移量计算两帧间的速度场：根据上述得出的总位移量除以两帧图像的时间间隔，最终求得两帧间的速度场。

3.根据权利要求1所述的用于太阳高分辨图像序列中相邻两帧间的多尺度速度场测量方法，其特征在于：所述步骤2中，对平滑处理后图像进行降采样即删除图像的偶数行与偶数列。

4.根据权利要求1所述的用于太阳高分辨图像序列中相邻两帧间的多尺度速度场测量方法，其特征在于：所述步骤4中，选取阈值r时，通过高斯解析函数生成两幅偏移量已知的图像，通过半径r计算得出两幅图像的偏移量，当半径r取不同的值时，计算得出的两幅图像间的偏移量是不同的，由于高斯解析函数生成两幅偏移量是已知的称为标准量，将实际测量的偏移量与标准量相比，实际测量的偏移量越接近标准量，说明计算该实际测量的偏移量所选取的半径r是最佳的，就选此时的r值为阈值r。