CN103236067A - 一种像素级sar影像时间序列构建的局部自适应配准方法 - Google Patents

Abstract

本发明涉及一种像素级SAR影像时间序列构建的局部自适应配准方法，方法如下：数据预处理之后从主从影像上提取同名特征点对，使用最小二乘法计算二次多项式参数并计算匹配总体误差，然后比较匹配总体误差与给定阈值的大小，若匹配总体误差小于或等于给定阈值，则从影像与主影像的位置关系由上述二次多项式确定，最后进行影像配准；反之，若总误差大于给定阈值，则进行误差点聚类获取畸变区域，将主、从影像的正常区域作为一对新主从影像，畸变区域作为另一对新主从影像，对两对新的主、从影像重复计算二次多项式参数及以后的步骤，直至所有新主、从影像的同名特征点对匹配总体误差小于给定阈值，然后进行影像配准。

Description

一种像素级SAR影像时间序列构建的局部自适应配准方法

技术领域

本发明涉及一种SAR影像配准方法，特别是涉及一种像素级SAR影像时间序列构建的局部自适应配准方法。

背景技术

遥感时间序列挖掘是时空数据挖掘的重要内容，研究如何有效地从海量、复杂、不同时间分辨率的遥感数据中挖掘潜在的有用知识，具有重要的理论价值和现实意义。遥感时间序列的构建是遥感时间序列挖掘的前提，因此，大数据背景下的遥感时间序列的建模方法研究变得十分迫切。目前应用比较广泛的遥感影像时间序列数据有MODIS、NOAA/AVHRR、SOPT/VEGETATI0N、Landsat等，它们在植被生长监测(Zhao et al.，2009；Melin et al.，2011；Maire etal，2011)、土地覆盖分类(L.Clark et al.，2010；Shao and Lunetta，2012；Klein et al.，2012)和地表变化监测(Verbesselt et al.，2010；Salmon et al.，2011；Ruiz et al.，2012)等方面发挥了重要作用。在微波遥感领域，SAR影像能够不受日照和天气条件的限制完成对地观测任务，并对某些地物具有一定的穿透能力，具有多极化、多入射角、多模式、短重访周期等特点，在全球变化和区域资源环境研究中具有独特的优势。SAR影像时间序列具有一定的稳定性，能够在一定程度上克服单幅影像的后散射值误差较大而难以应用于地物的信息提取、分类和变化监测的缺点，可以在海洋观测、环境监测、灾害监测、军事等方面的发挥独特的优势。

SAR影像时间序列主要包括三种类型：影像级SAR影像时间序列、区域级SAR影像时间序列和像素级SAR影像时间序列。相对于区域级和影像级的SAR影像时间序列，像素级SAR时间序列更多的保留了原始数据的信息量，尽可能多的保持了地物的细节层次信息，能够更准确地提取地物信息、揭示地物的变化规律。像素级SAR影像时间序列的构建的前提是不同时期影像的每个像元一一对应——即不同时期影像之间的高精度配准，而SAR影像中存在大量光斑噪声，特征点难以准确获取，高精度的影像匹配十分困难。现有的配准方法大多利用单一多项式方程约束主从影像的整体，在迭代的过程中，滤除误差大的匹配点对(如NEST、GAMMA等软件集成的方法)；然而，当多项式确定之后，误差较大的区域仍然要参与到内插计算完成整体匹配，配准之后误差仍会保留。因此若使用该方法进行配准，虽然整体上可以达到一定匹配要求，但局部畸变引起的误差却难以校正。

发明内容

本发明要解决技术问题是：克服现有技术缺点，提出一种像素级SAR影像时间序列构建的局部自适应配准方法，能够实现构建像素级SAR影像时间序列所需要的高精度配准。

为了解决上述技术问题，本发明提出的技术方案是：一种像素级SAR影像时间序列构建的局部自适应配准方法，包括以下步骤：

第一步、数据预处理——确定主影像和从影像，并对主、从影像进行辐射纠正、地形纠正和去噪处理；

第二步、提取同名特征点对——分别从主、从影像上提取同名特征点对；

第三步、计算二次多项式参数——根据主、从影像的同名特征点对的坐标，使用最小二乘法计算表征同名特征点对之间位置关系的二次多项式参数，所述二次多项式如下：

X＝c₀+c₁x+c₂y+c₃x²+c₄xy+c₅y²

Y＝d₀+d₁x+d₂y+d₃x²+d₄xy+d₅y²

式中，X、Y分别为主影像特征点的横坐标、纵坐标，x、y分别为其同名特征点在从影像中的横坐标、纵坐标，c₀、c₁、c₂、c₃、c₄、c₅及d₀、d₁、d₂、d₃、d₄、d₅为二次多项式各参数；

第四步、计算匹配总体误差——根据第三步得到的二次多项式参数对从影像的特征点作变换，并计算所有同名特征点对的匹配总体误差、以及每对同名特征点对的匹配误差；

第五步、匹配总体误差与给定阈值的比较——比较所述匹配总体误差与给定阈值的大小，若所述匹配总体误差小于或等于给定阈值，则从影像与主影像的位置关系由第三步得到的二次多项式确定，并转至第八步；若总误差大于给定阈值，则转至第六步，其中给定阈值的取值范围为：0.3-0.5个像素；

第六步、误差点聚类获取畸变区域——第四步获得的同名特征点对匹配误差超过所述给定阈值的特征点作为误差点，并对误差点进行局部聚类，误差点聚集区域为畸变区域，其余为正常区域；

第七步、影像划分——利用二分法分别从主、从影像上分离出正常区域与畸变区域，将主、从影像的正常区域作为一对新主从影像，畸变区域作为另一对新主从影像，对两对新的主、从影像重复第三步至第七步的操作，直至所有新主、从影像的同名特征点对匹配总体误差小于给定阈值；

第八步、影像配准——根据从影像各区域对应的二次多项式对从影像进行变换，完成配准。

本发明像素级SAR影像时间序列构建的局部自适应配准方法，相对于现有技术的区别是首先获取误差点，然后对误差点聚类，并据此将图像划分为正常区域和畸变区域，接着对得到的正常区域与畸变区域作为另一对新主从影像重新进行匹配。具体是：根据新主、从影像的同名特征点对的坐标得到二次多项式参数，利用二次多项式参数对新从影像的特征点作变换，并计算所有同名特征点对的匹配总体误差、以及每对同名特征点对的匹配误差，如此循环，直至所有新主、从影像的同名特征点对匹配总体误差小于给定阈值，最后根据从影像各区域对应的二次多项式对从影像进行变换，完成配准。

本发明像素级SAR影像时间序列构建的局部自适应配准方法，还具有如下改进：

1、第二步中，提取同名特征点对的方法如下：首先从主影像的行、列方向等距选取特征点，得到主影像特征点集；接着根据主影像特征点集内每个特征点的位置，确定从影像上同名特征点的初始位置，得到从影像初始特征点集；针对每一对特征点，分别以两特征点为中心在主、从影像上开设特征窗口，通过迭代计算主、从影像特征窗口的交叉相关影像，对确认为同名点的从影像特征点的位置进行修正，当一对特征点对无法确认为同名点对时，则删除该对特征点，从而得到准确的主从影像同名特征点对。

2、第二步计算主从两个特征窗口的交叉相关影像，其具体步骤如下：

1a)分别对成对的主影像特征窗口内的影像和从影像特征窗口内的影像进行二维离散傅里叶变换，分别得到主、从影像特征窗口内影像的波谱；

1b)将从影像特征窗口内影像的波谱的复共轭与主影像特征窗口内影像的波谱相乘，得到交叉波谱；

1c)对所述交叉波谱按2倍行比率和2倍列比率进行升采样；

1d)对交叉波谱的升采样结果进行傅里叶逆变换后取模，得到主、从影像特征窗口的交叉相关影像。

3、第六步中对匹配误差超过给定阈值的特征点进行局部聚类得到畸变区域，具体方法如下：

3a)对影像构建格网，格网的大小采用渐进式方法确定，即以一定步长从小到大取值，直到至少一半的格网内误差点对不少于30对为止；

3b)统计每个格网内的误差点对数目，并根据格网之间的距离计算格网之间的空间权重矩阵；

3c)计算每个格网的Getis-Ord统计量Z_i：

$Z_i=\frac{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{w}_{i,j}-\stackrel{\‾}{X}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{w}_{i,j}}{S*\sqrt{\frac{[n\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{w^2}_{i,j}-{\left(\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{w}_{i,j}\right)}^2]}{n-1}}}$

式中，Z_i为格网i的Getis-Ord统计量，

x_j是格网j中误差点对的数目，w_i，j是格网i和格网j的空间权重，n为格网总数；

3d)提取出Getis-Ord统计量处于前30％-40％的格网，并以Getis-Ord统计量最高的格网为基准，对提取出的格网进行连通性操作，将格网连接起来；

3e)计算连通格网的最大外接矩形，得到初始畸变区域；

3f)对初始畸变区域按行列宽的2％外扩，计算外扩过程中增加的误差点对数与现有误差点对数的比值，若该比值小于10％，则停止外扩，以外扩后的区域为最终畸变区域；若该比值大于10％，则以外扩后的区域代替初始畸变区域并重复本步骤。

4、所述第二步中分别从主、从影像上提取对同名特征点对数量不低于200对。

本发明的有益成果是：对前一次匹配得到的正常区域与畸变区域作为另一对新主从影像重新进行匹配，如此循环，直至所有新主、从影像的同名特征点对匹配总体误差小于给定阈值，最后根据从影像各区域对应的二次多项式对从影像进行变换，完成配准，本发明客服了既有方法对图像畸变区域精度低的缺点，能够实现构建像素级SAR影像时间序列所需要的高精度配准。

附图说明

下面结合附图对本发明作进一步的说明。

图1为本发明实施例的流程图。

图2为本发明实施例的研究区及相应的SAR影像序列示意图。

图3-a为本发明实施例的主影像示意图。

图3-b为本发明实施例的从影像示意图。

图4-a为本发明实施例初次迭代误差大于给定阈值的误差点分布示意图。

图4-b为本发明实施例初次迭代误差点聚类的畸变区域示意图。

图5为本发明局部自适应高精度配准的迭代流程示意图。

具体实施方式

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

如图2所示，本实施例使用数据为构建时间序列所需的珠江流域广东省西江段2006-2009年不同时期的21景宽幅模式的ENVISAT ASAR影像，数据为标准格式的1B级产品，空间分辨率150m，像素间距75m，极化方式包含HH和VV两种。

如图1所示，本发明实施例像素级SAR影像时间序列构建的局部自适应配准方法，包括以下步骤：

第一步、数据预处理——确定主影像和从影像，并对主、从影像进行辐射纠正、地形纠正和去噪处理；

如图3-a所示，选择获取自2006年9月28日的一景SAR影像，在ENVI软件中手动裁出1000*1000像素大小的子影像作为主影像；如图3-b所示，选择获取自2008年2月2日的一景SAR影像，手动裁出与主影像区域重叠的1000*1000像素大小的子影像作为从影像；在ENVI软件中分别对主从影像进行辐射纠正、地形纠正和去噪处理；

第二步、提取同名特征点对——分别从主从影像上提取同名特征点对；具体方法如下：首先从主影像的行、列方向以40个像元为间隔选取525个特征点，得到主影像特征点集；接着根据主影像特征点集内每个特征点的位置，确定从影像上同名特征点的初始位置，得到从影像初始特征点集，将主影像特征点集和从影像特征点集中相同位置的点视为特征点对；针对每一对特征点，分别以两特征点为中心分别在主、从影像上开设32*32大小的特征窗口，通过迭代计算主、从影像特征窗口的交叉相关影像，对确认为同名点的从影像特征点的位置进行修正，当一对特征点对无法确认为同名点对时，则删除该对特征点，从而得到准确的主从影像同名特征点对；

其中，计算主从两个特征窗口的交叉相关影像，其具体步骤如下：

1a)令主影像特征窗口为mI，从影像特征窗口为sI，分别对成对的主影像特征窗口mI内的影像和从影像特征窗口sI内的影像进行二维离散傅里叶变换，分别得到主、从影像特征窗口内影像的波谱mSP和sSP；

$\mathrm{mSP}(u,v)=\frac{1}{10*10}\underset{x=1}{\overset{M}{\mathrm{\Σ}}}\underset{y=1}{\overset{N}{\mathrm{\Σ}}}\mathrm{mI}(x,y)\·e^{-j2\mathrm{\π}(\frac{\mathrm{ux}}{10}+\frac{\mathrm{vy}}{10})}$

$\mathrm{sSP}(u,v)=\frac{1}{10*10}\underset{x=1}{\overset{M}{\mathrm{\Σ}}}\underset{y=1}{\overset{N}{\mathrm{\Σ}}}\mathrm{sI}(x,y)\·e^{-j2\mathrm{\π}(\frac{\mathrm{ux}}{10}+\frac{\mathrm{vy}}{10})}$

式中，M、N分别代表特征窗口的高度和宽度；

1b)将从影像特征窗口内影像的波谱sSP的复共轭与主影像特征窗口内影像的波谱mSP相乘，得到两特征窗口的交叉波谱；

1c)对交叉波谱按2倍行比率和2倍列比率进行升采样；

1d)对交叉波谱的升采样结果进行傅里叶逆变换后取模，得到主、从影像特征窗口的交叉相关影像。

第三步、计算二次多项式参数——根据主、从影像的同名特征点对的坐标，使用最小二乘法计算表征同名特征点对之间位置关系的二次多项式参数，二次多项式如下：

X＝c₀+c₁x+c₂y+c₃x²+c₄xy+c₅y²

Y＝d₀+d₁x+d₂y+d₃x²+d₄xy+d₅y²

式中，X、Y分别为主影像特征点的横坐标、纵坐标，x、y分别为其同名特征点在从影像中的横坐标、纵坐标，c₀、c₁、c₂、c₃、c₄、c₅及d₀、d₁、d₂、d₃、d₄、d₅为二次多项式各参数；

第四步、计算匹配总体误差——根据第三步得到的二次多项式参数对从影像的特征点作变换，并计算所有同名特征点对的匹配总体误差、以及每对同名特征点对的匹配误差；

第五步、匹配总体误差与给定阈值的比较——比较匹配总体误差与给定阈值的大小，若匹配总体误差小于或等于给定阈值0.5个像素，则从影像与主影像的位置关系由第三步得到的二次多项式确定，并转至第八步，计算结束；若总误差大于给定阈值0.5个像素，则转至第六步；

第六步、误差点聚类获取畸变区域——第四步获得的同名特征点对匹配误差超过给定阈值0.5个像素的特征点作为误差点，对误差点进行局部聚类，误差点聚集区域为畸变区域，其余为正常区域，图4-a为误差点，图4-b虚线框区域所示部分为畸变区域；具体方法如下：

3a)对影像构建格网，格网的大小采用渐进式方法确定，即以一定步长从小到大取值，直到至少一半的格网内误差点对不少于30对为止；

3b)统计每个格网内的误差点对数目，并根据格网之间的距离计算格网之间的空间权重矩阵；

3c)计算每个格网的Getis-Ord统计量Z_i：

$Z_i=\frac{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{w}_{i,j}-\stackrel{\‾}{X}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{w}_{i,j}}{S*\sqrt{\frac{[n\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{w^2}_{i,j}-{\left(\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{w}_{i,j}\right)}^2]}{n-1}}}$

式中，Z_i为格网i的Getis-Ord统计量，

x_j是格网j中误差点对的数目，w_i，j是格网i和格网j的空间权重，n为格网总数；

3d)提取出Getis-Ord统计量处于前30％的格网，并以Getis-Ord统计量最高的格网为基准，对提取出的格网进行连通性操作，将格网连接起来；

3e)计算连通格网的最大外接矩形，得到初始畸变区域；

3f)对初始畸变区域按行列宽的2％外扩，计算外扩过程中增加的误差点对数与现有误差点对数的比值，若该比值小于10％，则停止外扩，以外扩后的区域为最终畸变区域；若该比值大于10％，则以外扩后的区域代替初始畸变区域并重复本步骤。

第七步、影像划分——如图5所示，利用二分法分别从主、从影像上分离出正常区域与畸变区域，将主、从影像的正常区域作为一对新主从影像，畸变区域作为另一对新主从影像，对两对新的主、从影像重复第三步至第七步的操作，直至所有新主、从影像的同名特征点对匹配总体误差小于给定阈值0.5个像素时停止迭代；

第八步、影像配准——根据从影像各区域对应的二次多项式使用三次卷积法对从影像进行变换，完成配准。

除上述实施例外，本发明还可以有其他实施方式。凡采用等同替换或等效变换形成的技术方案，均落在本发明要求的保护范围。

Claims (5)

1.一种像素级SAR影像时间序列构建的局部自适应配准方法，步骤包括：

第一步、数据预处理——确定主影像和从影像，并对主、从影像进行辐射纠正、地形纠正和去噪处理；

第二步、提取同名特征点对——分别从主、从影像上提取同名特征点对；

第三步、计算二次多项式参数——根据主、从影像的同名特征点对的坐标，使用最小二乘法计算表征同名特征点对之间位置关系的二次多项式参数，所述二次多项式如下：

X＝c₀+c₁x+c₂y+c₃x²+c₄xy+c₅y²

Y＝d₀+d₁x+d₂y+d₃x²+d₄xy+d₅y²

式中，X、Y分别为主影像特征点的横坐标、纵坐标，x、y分别为其同名特征点在从影像中的横坐标、纵坐标，c₀、c₁、c₂、c₃、c₄、c₅及d₀、d₁、d₂、d₃、d₄、d₅为二次多项式各参数；

第四步、计算匹配总体误差——根据第三步得到的二次多项式参数对从影像的特征点作变换，并计算所有同名特征点对的匹配总体误差、以及每对同名特征点对的匹配误差；

第五步、匹配总体误差与给定阈值的比较——比较所述匹配总体误差与给定阈值的大小，若所述匹配总体误差小于或等于给定阈值，则从影像与主影像的位置关系由第三步得到的二次多项式确定，并转至第八步；若总误差大于给定阈值，则转至第六步，其中给定阈值的取值范围为：0.3-0.5个像素；

第六步、误差点聚类获取畸变区域——第四步获得的同名特征点对匹配误差超过所述给定阈值的特征点作为误差点，并对误差点进行局部聚类，误差点聚集区域为畸变区域，其余为正常区域；

第七步、影像划分——利用二分法分别从主、从影像上分离出正常区域与畸变区域，将主、从影像的正常区域作为一对新主从影像，畸变区域作为另一对新主从影像，对两对新的主、从影像重复第三步至第七步的操作，直至所有新主、从影像的同名特征点对匹配总体误差小于给定阈值；

第八步、影像配准——根据从影像各区域对应的二次多项式对从影像进行变换，完成配准。

2.根据权利要求1所述的像素级SAR影像时间序列构建的局部自适应配准方法，其特征在于：第二步中，提取同名特征点对的方法如下：首先从主影像的行、列方向等距选取特征点，得到主影像特征点集；接着根据主影像特征点集内每个特征点的位置，确定从影像上同名特征点的初始位置，得到从影像初始特征点集；针对每一对特征点，分别以两特征点为中心在主、从影像上开设特征窗口，通过迭代计算主、从影像特征窗口的交叉相关影像，对确认为同名点的从影像特征点的位置进行修正，当一对特征点对无法确认为同名点对时，则删除该对特征点，从而得到准确的主从影像同名特征点对。

3.根据权利要求2所述的像素级SAR影像时间序列构建的局部自适应配准方法，其特征在于：第二步计算主从两个特征窗口的交叉相关影像，其具体步骤如下：

1a)分别对成对的主影像特征窗口内的影像和从影像特征窗口内的影像进行二维离散傅里叶变换，分别得到主、从影像特征窗口内影像的波谱；

1b)将从影像特征窗口内影像的波谱的复共轭与主影像特征窗口内影像的波谱相乘，得到交叉波谱；

1c)对所述交叉波谱按2倍行比率和2倍列比率进行升采样；

1d)对交叉波谱的升采样结果进行傅里叶逆变换后取模，得到主、从影像特征窗口的交叉相关影像。

4.根据权利要求1所述的像素级SAR影像时间序列构建的局部自适应配准方法，其特征在于：第六步中对匹配误差超过给定阈值的特征点进行局部聚类得到畸变区域，具体方法如下：

3a)对影像构建格网，格网的大小采用渐进式方法确定，即以一定步长从小到大取值，直到至少一半的格网内误差点对不少于30对为止；

3b)统计每个格网内的误差点对数目，并根据格网之间的距离计算格网之间的空间权重矩阵；

3c)计算每个格网的Getis-Ord统计量Z_i：

$Z_i=\frac{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{w}_{i,j}-\stackrel{\‾}{X}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{w}_{i,j}}{S*\sqrt{\frac{[n\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{w^2}_{i,j}-{\left(\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{w}_{i,j}\right)}^2]}{n-1}}}$

式中，Z_i为格网i的Getis-Ord统计量，

x_j是格网j中误差点对的数目，w_i，j是格网i和格网j的空间权重，n为格网总数；

3d)提取出Getis-Ord统计量处于前30％-40％的格网，并以Getis-Ord统计量最高的格网为基准，对提取出的格网进行连通性操作，将格网连接起来；

3e)计算连通格网的最大外接矩形，得到初始畸变区域；

3f)对初始畸变区域按行列宽的2％外扩，计算外扩过程中增加的误差点对数与现有误差点对数的比值，若该比值小于10％，则停止外扩，以外扩后的区域为最终畸变区域；若该比值大于10％，则以外扩后的区域代替初始畸变区域并重复本步骤。

5.根据权利要求1所述的像素级SAR影像时间序列构建的局部自适应配准方法，其特征在于：所述第二步中分别从主、从影像上提取对同名特征点对数量不低于200对。

Images