一种考虑地表形变的Landsat7条带影像修复方法与流程

本发明涉及条带影像修复领域，尤其是涉及一种考虑地表形变的Landsat 7条带影像修复方法。

背景技术：

2003年5月31日，Landsat 7ETM+机载扫描行校正器(ScanLinesCorrector，简称SLC)发生故障，导致获取的图像出现数据重叠和大约22％的数据丢失，因而在2003-2013年期间近10年的Landsat7影像存在异常，给这期间相关研究造成困难。大部分Landsat7影像条带修复研究都是针对没有地表形变的研究区域，包括 USGS(U.S.Geological Survey)在2004年提出的基于ENVI软件的多种条带修复方法，都无法在地表形变区域得到良好的修复结果。地表形变主要表现为地震、山体滑坡、冰川运动等现象，在影像成像时间间隔内同名地物点间产生了一定的位移量，如果直接采用以往的研究算法，则恢复出来的地物会产生错位现象。极地冰川一直处于不断的运动状态，其物质流失对全球海平面上升起着至关重要的影响。 Landsat影像是极地冰川研究的重要数据来源，将近10年的影像存在条带现象，极大地影响冰川研究。本文以南极冰川区域为例，选取该区域的Landsat 7影像进行条带修复，有助于冰川表面的流速计算以及冰架表面特征提取等进一步利用。

技术实现要素：

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种考虑地表形变的Landsat 7条带影像修复方法。

本发明的目的可以通过以下技术方案来实现：

一种考虑地表形变的Landsat 7条带影像修复方法，包括以下步骤：

1)获取周期间隔短且覆盖同一地区的Landsat 7ETM+的两景条带影像，其中，一景作为待修复影像，另一景作为填补影像；

2)获取待修复影像的灰度直方图，确定其灰度阈值，并根据此灰度阈值进行粗匹配；

3)对粗匹配后的待修复影像采用最小二乘影像匹配进行精匹配，并根据最小二乘影像匹配得到的辐射畸变改正参数，计算待修复影像条带处对应的像元值。

所述的步骤2)中，灰度阈值的确定方法如下：

在待修复影像上获取无明显特征点区域的灰度范围，取其灰度最小值为灰度阈值。

所述的步骤2)中，进行粗匹配的具体步骤为：

21)设定待修复影像的模板窗口尺寸以及填补影像的搜索窗口尺寸；

22)当模板窗口中所有像素的像元值均大于灰度阈值，则采用频率域互相关匹配法进行粗匹配；

23)当模板窗口中含有像元值小于灰度阈值的像素，则采用空间域归一化互相关匹配法进行粗匹配。

所述的步骤23)中，空间域归一化互相关匹配法中相似性测度的相关函数 CC(i,j)的计算式为：

其中，s为搜索窗口，r为模板，(i,j)为搜索窗口的位置，(k,l)为模板的位置，μr为模板窗口除条带区域外的所有像素的灰度平均值，μs为搜索窗口除条带区域外的所有像素的灰度平均值，

所述的步骤22)中，频率域互相关匹配法中的互相关系数CC(i,j)的计算式为：

CC(i,j)＝IFFT(F(u,v)G^*(u,v))

其中，F(u,v)为待修复影像匹配窗口由快速傅立叶变换计算表示，G^*(u,v)为填补影像匹配窗口由快速傅立叶变换的共轭复数，IFFT为快速傅立叶逆变换。

所述的步骤22)中，随机选取灰度直方图中的像元值赋值给条带像元，保证快速傅立叶变换的模板中不存在0值，这些随机选取的灰度值与匹配窗口没有一致的关联性，因此不会影响相关性计算。

所述的步骤3)具体包括以下步骤：

31)考虑线性灰度畸变，获取最小二乘影像匹配的误差v：

v＝c1dh0+c2dh1+c3da0+c4da1+c5da2+c6db0+c7db1+c8db2-Δg h0和h1的初值为：

其中，c1、c2、c3、c4、c5、c6、c7、c8为误差方程式系数，h0、h1为辐射畸变改正参数，a0、a1、a2、b0、b1、b2为仿射变换系数，Δg为灰度差，mean(p)为待修复影像搜索窗口内有效像元值的平均值，mean(f)为填补影像模板窗口内有效像元值的平均值；

32)根据待修复影像与填补影像的同名像元之间的线性关系，计算待修复影像中每一个条带覆盖的像元的真实灰度值pt，其计算式为：

pt＝aft+b

其中，pt为待修复影像中条带覆盖的像元的灰度值，ft为填补影像在同名像元位置的灰度值，a、b为待定系数，其值为最小二乘匹配中使得误差v最小时对应的辐射畸变改正参数和

所述的步骤21)中，待修复影像的模板窗口尺寸为三倍最大条带宽度，所述的填补影像的搜索窗口尺寸大于待修复影像的模板窗口。

与现有技术相比，本发明具有以下优点：

一、适用于地表形变区：现有影像修复技术大都应用在地表随时间无明显变化区域，而本发明加以匹配手段实现了在地表形变区的影像修复；

二、无需预先配准：现有影像修复技术需要先对待修复影像和填补影像做配准处理，而本发明在最小二乘匹配中考虑了影像的几何畸变，不需要先进行两景影像的配准预处理；

三、修复效果好：修复效果优于常用的基于ENVI软件平台的局部线性直方图匹配法。

附图说明

图1为本发明的方法流程图。

图2为Landsat 7ETM+条带影像，其中，图(2a)为待修复影像，图(2b) 为填补影像。

图3为影像修复结果图，其中，图(3a)为采用局部线性直方图匹配法的修复结果，图(3b)为采用本方法的修复结果。

图4为冰流区域局部细节特征，其中，图(4a)为修复前的条带影像，图(4b) 为局部线性直方图匹配法修复结果，图(4c)为采用本方法的修复结果。

图5为山地局部细节特征，其中，图(5a)为修复前的条带影像，图(5b) 为局部线性直方图匹配法修复结果，图(5c)为采用本方法的修复结果。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。

本方法主要有以下三个步骤：

(1)使用空间域的归一化互相关匹配和基于傅立叶变换的频率域互相关匹配方法相结合

考虑地表复杂的覆盖类型，不同的匹配方法适用于不同地表覆盖状况。空间域的归一化互相关匹配方法适用于有明显特征点的区域，基于傅立叶变换的频率域互相关匹配方法相较于其他匹配方法，在特征点不明显地区能获得较好的匹配结果。因此，首先获得待修复影像的灰度直方图，确定阈值大小，若匹配模板内所有像元值都大于该阈值，则采用频率域互相关匹配法；反之若有像元值小于该阈值，说明有明显的特征地物(其灰度值通常较小)，则采用空间域的归一化互相关匹配法。

归一化互相关匹配是通过计算模板和待匹配图像的互相关值来确定匹配的程度，互相关值最大时的搜索窗口位置决定了模板图像在待匹配图像中的位置。作为相似性测度的相关函数计算公式为：

其中，s为搜索窗口，r为模板，(i,j)为搜索窗口的位置，(k,l)模板的位置，μr为模板窗口除条带区域外的所有像素的灰度平均值，μs为搜索窗口除条带区域外的所有像素的灰度平均值。CC(i,j)的取值范围为[-1，1]，其绝对值越大，相似性越高。归一化互相关匹配算法有很高的准确性和适应性，并且不受图像灰度值线性变换的影响。

基于频率域互相关匹配是通过傅立叶变换将二维空间图像转换到频率域，根据卷积定理两景图像在频率域相乘，只做互相关计算，没有进行归一化，计算公式如下：

CC(i,j)＝IFFT(F(u,v)G^*(u,v)) (2)

其中，F(u,v)是待修复影像匹配窗口由快速傅立叶变换计算表示，G^*(u,v)是填补影像匹配窗口由快速傅立叶变换的共轭复数，IFFT则是快速傅立叶逆变换，将频率域图像转回到空间域。

由于模板中存在条带(像元值为0)，使用快速傅立叶变换要求模板中不能有0 值，实验中随机选取灰度直方图中的像元值赋值给条带像元，进而保证不受条带干扰。这些随机选取的灰度值与匹配窗口没有一致的关联性，因此不会影响相关性计算。

(2)最小二乘影像匹配

由于影像匹配窗口尺寸均很小，一般只考虑一次畸变，可以用仿射变换模型描述：

其中，(a0,a1,a2,b0,b1,b2)为仿射变换系数，(x1,y1)为左方影像坐标，(x2,y2) 为右方影像坐标。

同时，考虑右方影像相对于左方影像的线性灰度畸变，可得：

g1(x,y)+n1(x,y)＝h0+h1g2(a0+a1x+a2y,b0+b1x+b2y)+n2(x,y) (4)

经线性化后，可得最小二乘影像匹配的误差方程式：

v＝c1dh0+c2dh1+c3da0+c4da1+c5da2+c6db0+c7db1+c8db2-Δg (5)

h0和h1的初值计算如下：

h0＝0 (6)

其中，mean(p)为待修复影像模板窗口内有效像元值的平均值，mean(f)为填补影像模板窗口内有效像元值的平均值。

(3)条带像元值计算

待修复影像与填补影像的同名像元之间存在线性关系，每一个条带覆盖的像元都可以通过该线性转换关系计算得到其真实灰度值，计算公式如下：

pt＝aft+b (8)

其中，pt为待修复影像在某个像元处的灰度值，ft为填补影像在同名位置的灰度值，a、b为待定系数，系数的最优值对应于最小二乘匹配中的辐射畸变改正参数和

实施例1：

本实施例中采用南极松岛冰川区域的Landsat 7ETM+全色波段影像，拍摄时间为2010年12月8号和2010年12月1号，两景影像都已正射纠正，影像分辨率为15米。两景条带影像如图2和表1所示。

表1实验数据

实验结果

南极冰川区域地物类型虽不复杂，但是除了流速较快的冰流区有明显的特征点，其它地区则是大面积冰雪覆盖，纹理特征不明显，考虑上述两个不同特征，首先在待修复影像上确定无明显特征点的冰雪覆盖区域的灰度范围，取灰度最小值为阈值确定粗匹配方法。由匹配计算得到同名点对应的位置，再根据线性转换关系求得条带处的真实像元值。此外，为了比较，采用USGS开发的基于ENVI平台的局部线性直方图匹配算法修复上述影像。

图3(a)是本发明提出的算法修复结果，3(b)是采用局部线性直方图匹配算法修复的结果。整体上看，两个修复结果与原条带影像(图2a)相比，被条带遮挡区域已经恢复，被修复的区域与周围像元间没有明显的灰度色差。图4和图5 分别展示了快速冰流区域和山地区域的局部细节特征图，所在处为图3方框标示处。根据图4的三幅图比较，局部线性直方图匹配法虽然能够恢复出条带遮挡处的像元值，但由于不同时期的两景影像间存在位移偏差，恢复出来的像元值不对应真实地物样貌，并且可以明显看出条带覆盖处与周围区域存在错位现象。然而，本发明运用的算法，不仅修复区域与整体影像没有灰度色差，而且没有出现错位现象，呈现出真实的地物样貌。同样，图5在有山的地区，由于待修复影像与填补影像之间存在几何偏差的原因，局部线性直方图匹配法在没有配准的基础上无法改正这一畸变差，导致恢复区域产生错位。而在南极地区，由于特征点(山等岩石)分布稀疏且不均匀，大部分地区只能采用简单的平移方法进行配准，对影像之间的旋转和缩放无法进行改正。而本发明的算法，不仅在最小二乘匹配中对辐射畸变进行改正，而且同时考虑了几何畸变的改正，因此修复后的影像质量很好，恢复出了山的轮廓样貌。

本发明提供一种考虑地表形变的Landsat 7影像条带修复方法，首先根据匹配算法确定同名地物在不同时期两影像上的位置(精确到亚像素级别)，然后由最小二乘匹配方法计算得到的辐射畸变改正参数，计算条带遮挡处的像元灰度值。将此方法应用到南极有快速冰流的区域，成功修复了Landsat 7条带影像，提高了影像的可利用性。