光学影像匹配技术在地表形变监测中的应用研究

张丽

(上海市测绘院，上海 200129)

1 引 言

地表形变是在地壳表层发生的一种空间位移变化，由于地壳运动、火山运动、地震活动等自然因素，或者是地下水开采、矿山开采、城市建设等人为因素引发。近几年地震频发，有专家预测全球开始进入板块运动活跃期，因此，对地表进行变形监测显得尤为重要。

长期以来，研究者主要依靠传统的地面水准测量和GPS监测网来进行监测，其监测精度虽可高达亚毫米级，但不适宜大区域的长期重复监测，且要消耗大量人力物力。近年来，随着SAR技术的发展，其全天时、全天候工作，且可以大范围、高精度地从空间获取地形信息和地表微小的形变信息的优势弥补了这一缺陷；然而该技术容易受相位去相关影响，在形变过大区域常常无法监测，且只能获取沿雷达视线向即距离向的地表形变。

而光学影像匹配技术是借助遥感影像的幅度信息，采用互相关技术寻找两幅图像的同名点，从而获得同名点处对应像素在方位向和距离向的偏移量。由于避免了相位运算，该方法能很好地弥补了受相位去相关影响这一不足。因此，本文将该方法应用到具体地表形变监测案例中，并研究其可行性。

2 光学影像匹配

2.1 光学影像匹配技术原理

影像匹配技术的根本就是图像配准，在光学影像里，图像配准也称为几何精校正，它规避了成像的空间几何过程，直接利用地面控制点(GCP，Ground Control Points)数据，对卫星影像的几何畸变过程进行数学模拟，利用已知条件确定模型中的位置参数，建立原始的畸变图像空间与校正空间之间的相互对应关系，然后利用这种对应关系把畸变图像空间中的全部元素变换到校正图像空间中去，从而实现几何精校正。精校正后，光学遥感影像间的幅度信息仍存在差异，通过傅立叶变换，在频率域内采用高精度的相位相关技术即可提取偏移量。

2.2 光学影像匹配技术流程

透过原理可以看出，光学影像匹配技术主要分为正射校正和相位相关计算两步，本文所采用的数据处理流程大致包括：

①利用外部DEM数据及影像的太阳高度角，方位角模拟生成DEM阴影图；

②从震前影像与DEM阴影图上选取同名点，对震前影像进行正射纠正；

③从震后影像与已校正的震前影像上选取同名点，对震后影像进行正射纠正；

④对已校正后的震前震后影像做相位相关计算，获取有效形变值。

(1)选取同名点

在上述两次配准过程中均涉及同名点的选取，通常选取同名点分为人工选点和自动选点两种方法。本文采取人工选点法，即通过目视在两幅待配准影像上，手动选择若干同名点。同名点选取时需要遵循以下原则：

①至少选3个点，若待纠正影像和一幅色彩较暗的DEM配准时需要选15～30个点；

②同名点要连续选择，但需远离形变区域，若整幅图都有形变就在形变最小处选点；

③点选取时不能太靠近图像边界；

④同名点周围最好有容易识别的地物，最好选择无须校正的人工物。

(2)影像正射校正与重采样

正射校正是对图像空间和几何畸变进行校正生成多中心投影平面正射图像的处理过程。它除了能纠正一般系统因素产生的几何畸变外，还可以消除地形引起的几何畸变。采用少量的地面控制点与相机或卫星成像模型相结合，确立相机或传感器平台、图像和地面三者之间的空间位置关系，建立正确的校正公式，产生精确的正射影像。

影像重采样是影像数据重新组织过程中的灰度处理方法，影像采样是按一定的间隔采集影像灰度数值，当阈值不位于采样点上的原始函数的数值时，就需要利用已采样点进行内插，称为重采样。本文所采用的重采样法为三次卷积内插法(也叫辛克插值法)，使用内插点周围的16个像元值进行距离加权计算栅格值，先在Y方向内插四次(或X方向)，再在X方向(或Y方向)内插四次，最终得到该像元的栅格值。

(3)相位相关法估算偏移量

在经过正射校正和重采样后，两幅影像应该具有相同的分辨率，统一的坐标系，但因为存在地表形变，两幅影像并不能完全重合。因此，还需要利用COSI-CORR软件模块提供的高精度相位相关算法进行偏移量估计。该模块包含了两种算法，一种是基于频率域(Frequential)变换的方法，另一种是基于统计学(Statistical)理论的方法。相位相关法是用于配准图像进行平移变换的典型方法，它主要是依赖于傅立叶偏移理论。

假设有两幅影像i1和i2，用(△x，△y)来代表两幅影像上对应点的形变差异，则有公式：

i2(x，y)=i1(x-△x，y-△y)

(1)

用I1和I2分别代表两幅影像的傅立叶变换，则有以下关系式：

I2(wx，wy)=I1(wx，wy)e-j(wx△x+wy△y)

(2)

其中wx和wy分别代表行和列的频率变化。定义两幅影像的互能量谱为如下：

(3)

上式中，C(wx，wy)是i1(x，y)和i2(x，y)之间的互能量谱，*代表复数共轭相乘。

影像间的相对形变量可以通过互能量谱的二维斜率来计算，因此对式(3)进行傅立叶逆变换：

F-1{ej(wx△x+wy△y)}=δ(wx△x+wy△y)

(4)

而在应用中，C(wx，wy)仅为两幅影像间理论上的互能量谱，由实际计算得到的互能量谱为Q(wx，wy)。定义一个目标函数为φ(△x，△y)：

(5)

其中，W(wx，wy)表示加权矩阵。目标函数φ(△x，△y)在其他位置都为零，只有在平移位置处不为零，因此这个位置就是影像间的相对偏移，当目标函数φ(△x，△y)达到最小值时的(△x，△y)就是所求的偏移量。

3 实验数据处理

北京时间2010年4月14日，在我国青海省玉树县发生了震级为里氏7.1级的大地震。此次地震的发震构造是青藏高原东部规模巨大的甘孜-玉树断裂带，该断裂带长度约 500 km，总体走向呈285°～315°弧形展布，是一条大型左旋走滑断裂带。

3.1 实验数据参数

本实验选用了两幅SPOT5卫星影像。影像数据采集的时间分别为2009年11月5日和2010年4月15日，影像数据具体参数如表1所示。选用时间相距较近的卫星影像可以最大限度排除人工建筑等因素的影响，协调时间、太阳方位角、入射角的一致，确保了不会由于角度问题导致正射校正的偏差。

实验选用的SPOT5卫星影像参数 表1

续表1

选用的SRTM数据绝对高程精度为 16 m，相对高程精度为 10 m左右，坐标系为WGS84地理坐标系，投影方式为UTM投影。

3.2 实验结果分析

本文利用基于ENVI平台开发的COSI-CORR软件模块，按照上述技术路线进行数据处理。

(1)生成DEM阴影图

由于已有SRTM数据与SPOT5影像分辨率存在较大差异，数据处理的第一步是先将DEM插值到空间分辨率 10 m，然后根据太阳高度角和太阳方位角生成一幅DEM阴影图，将其作为参考影像，用于震前影像的正射校正中，如图1所示。

(2)SPOT5影像匹配

两次配准时，在SPOT影像和DEM阴影图上，通过人工选点法分别选取了18个和19个均匀分布的同名点。在选取过4个及以上同名点后，软件模块会自动预测后续同名点位置，为提高配准精度，在经过重新选取和删除误差较大的点后，保证两次的均方根误为分别为0.48和0.43，均小于1个像元。

图1 生成的DEM阴影图

选取同名点后，需将其转换为GCP点并进行优化，才能用于后续校正与重采样。而COSI-CORR软件模块提供两种优化方法(相位相关法和统计相关法)，可以使同名点间配准误差值达到分米级甚至是厘米级。第一次正射校正时，由于参考影像为DEM阴影图，其影像特点与SPOT影像存在很大差异，因而选用统计相关法。而第二次正射校正时，由于参考影像为校正过后的震前影像，两幅影像特点相一致，因此选用相位相关法来提高计算精度。两次优化的结果表明，选取5次迭代计算能够使误配准值即残差偏移量达到收敛。点优化过程的迭代结果如图2所示。

图2 5次迭代后GCP点处的残差偏移量均值和标准差

(3)影像校正与重采样

当GCP点优化结果达到收敛且满足精度要求后，影像进行正射校正与重采样。此模块采用的正射校正模型为严格轨道模型(Pushroom Sensor)，校正过程利用头文件中卫星状态参数提高几何校正的精度。在选用精确度最高的辛克插值法进行重采样后，得到校正后的震前、震后影像如图3和图4所示。

图3 校正后的震前影像

图4 校正后的震后影像

(4)相位相关计算偏移量

本文经过反复对比，进行相位相关法计算时，选用滑动窗口为64×64的单一型窗口，即X和Y方向的Initial和final都选择64；滑动窗口在影像X和Y方向移动的像素数都采用默认值16；迭代次数选择2，掩模阈值选择0.9。获得了此次地震东西向和南北向的同震位移量，结果如图5所示。此一步可输出偏移量文件，共包含三个波段数据，分别为地表东西向形变值，南北向形变值及信噪比SNR值。

图5玉树地震同震位移图(左为东西向形变，右为南北向形变)

图6 去噪后玉树地震局部同震位移分布图

(5)形变结果优化

通过对结果的SNR值(0效果差，1效果优)分析，得到了同震位移场的质量评价。图6中，上图为震前SPOT5卫星影像，红线为同震地表破裂带，中图为东西向的同震位移分布图，下图为南北向的同震位移分布图。从图中可以看到，虽有一些效果比较差的地方，但主要分布在远离地表破裂带处，在破裂带展布区域分布很少。这是由于玉树地处高原地区，地形起伏比较大，在一些地区得到的数据结果不是很理想，而且由于震后SPOT5卫星影像为第二天应急拍摄，没有考虑云量，导致震后影像云量比较多，但幸运的是云量基本没有覆盖地表破裂带的范围，而且地形起伏导致的不理想数据结果也很少在地表破裂带上，这保证了得到的同震位移场至少在近地表破裂带处是比较可信的。

从图中可以明显地看出地表破裂带的展布以及同震位移量的分布。此次玉树地震是左旋走滑型地震，产生的同震位移主要为水平位移，垂直位移分量比较小。其最大位移不超过 2 m，地表裂长度为 30 km左右。

4 结 语

本文对光学影像匹配技术进行了初步研究，利用光学影像空间分辨率高的特点，获取了玉树地震水平方向的形变场。根据野外地表破裂带调查，此次地震较清晰的地表破裂带由3条主破裂左阶组成，总长约 31 km，左旋走滑性质，最大左旋位移量约为 1.8 m。这与本文得到的破裂带长度和最大同震位移量基本一致，这也从另一个角度说明了实验得到的同震位移场的合理性和正确性，显示了光学影像匹配技术在监测地表形变方面的巨大应用价值。