一种子像素偏移追踪算法的实例验证

杨士超，查显杰，郭晓丹，张志宏，张 琪

（1.辽宁省地震局，辽宁 沈阳 110034；2.中国科学技术大学，安徽 合肥 230026）

0 引言

大地震地表破裂的观测资料对于研究地震力学和断层区域结构性质都具有重要意义。断层几何结构和相关滑移分布也是决定地震动力学破裂过程的关键所在[1]。利用地表形变位移场可以反演断层的滑移分布，研究断层结构和断层破裂状态；同时也可以将地表位移场作为边界条件，对地震断裂运动进行约束。因此对于研究地震动力学和运动学问题来说，获取大地震同震位移场都是至关重要的科学问题。

目前，地表形变监测技术主要包括野外地质调查[2]、GPS 台站测量[3]、合成孔径雷达干涉（InSAR）测量[4]、子像素偏移追踪测量[5-6]。对于广泛区域的偏移场获取，野外地质调查需要大量的人力和物力，而偏远地区GPS 台站覆盖有限，GPS 技术很难有效发挥作用。InSAR 是一种理想的大地震近场形变监测技术。然而，大地震的近场形变变化梯往往较大，断层附近的InSAR 形变干涉图条纹常常非常密集（即出现信号饱和现象），因此InSAR 方法获取的断层附近区域的位移场精度会显著降低。子像素偏移追踪算法主要是通过处理地震前和地震后获取的两景光学遥感影像，采用互相关技术得到大地震的水平位移场[7-9]。该方法不存在形变信号饱和的问题，因此它非常适合监测大地震的地表形变。本研究小组在强度互相关、双线性插值和二维高斯回归三种技术基础上新发展了一套子像素偏移追踪算法[10]。与现有算法相比，该算法能以更高的精度从遥感影像对中获取位移场。为了验证算法的精度和可靠性，本文利用新发展的子像素偏移追踪算法对2014 年2 月12 日新疆于田地震和2013 年9 月24 日巴基斯坦地震观测资料进行处理，获取这两次地震的同震形变位移场。

1 算法原理与实现

1.1 算法原理

利用光学遥感影像求解形变位移场的精度取决于影像的分辨率和子像素偏移追踪算法的精度。由于遥感成像技术的限制，所获取的光学影像分辨率都是一定的（通常商用卫星数据空间分辨率为米级或亚米级），因此需要发展一种能获取高精度位移场的子像素偏移追踪算法。目前，子像素偏移追踪算法主要有两种：一种是频域互相关算法，例如COSI-Corr 程序包[11]；另一种是空间域互相关算法，例如ImCorr 程序包[12]。

1.1.1 频域互相关算法

频域互相关又称作相位互相关，这种算法主要基于傅里叶偏移原理[13]。空间域内两个函数存在的偏移值可以转变为频域内线性的相位差。可以用如下公式表达：

假设f1（x，y）和f2（x，y）分别表示两个图像的函数，（x0，y0）表示f1和f2之间的相对偏移值，那么

根据傅里叶偏移原理可得：

其中，F1（u，v）和F2（u，v）分别表示f1（x，y）和f2（x，y）的傅里叶变换。因此，标准互功率谱可以表示为：

其中*表示复共轭。

通过公式（3）求解偏移值（x0，y0）有两种方式：一种方式是直接在频率域内求解，拟合相位平面，相位平面的周期性使得这种方法很难实现；另一种方式是将标准互功率谱进行傅里叶反变换，便可以得到δ（x-x0，y-y0），通过找到狄拉克δ 函数的峰值位置来确定偏移值（x0，y0）[14-15]。后一种方法只能搜索到整数偏移值，要求得更高精度的偏移值还需要借助曲面拟合算法来求解。

1.1.2 空间域互相关算法

空间域互相关算法又称作强度互相关算法，通过查找主像与从像的相关系数矩阵的峰值位置来确定两幅影像的相对偏移值。相关系数的表达方式有很多种，因为零均值标准化相关系数（ZNCC）具有良好的抗噪性，对于偏移和亮度的线性变化不敏感[16]，本文选择ZNCC 算法计算相关系数。ZNCC 相关系数表示为：

其中|*|表示绝对值，m（i，j）和s（i，j）分别表示主像与从像第i 行j 列的像素灰度值，M和N 表示所选像对的大小，分别表示m（i，j）和s（i，j）的均值。

利用以上算法可以求得一个相关系数矩阵，通过搜索峰值位置可以找到整数倍偏移值。接下来需要通过插值和曲面拟合来计算得到小数倍偏移值。一般插值算法可以提高影像的分辨率，但是高倍插值计算会耗费大量内存空间和运算时间。本文采用了一种基于双线性插值的相关系数插值方法[17]。这种算法可以将估算子像素偏移问题转化成一个关于偏移值的非线性反演问题，通过寻找最优解来确定偏移值。为了获取更高精度的偏移，本文算法进一步选取相关系数矩阵峰值周围的9 个点，对这9 个点进行9 点高斯回归求取微小偏移值[18]：

其中，C00，C01，C10，C11，C02和C20为通过回归分析推导出的6 个回归系数。

1.2 本文算法的实现

本文所验证的子像素偏移追踪算法主要基于上面提到的强度互相关算法、双线性插值算法和二维高斯回归算法，把这3 种算法巧妙的结合到一起来进行偏移计算。该算法的实现由以下5 步组成：

（1）利用两个大小一样的滑动窗在地震前和地震后获取的两景影像相应位置选取一组像片对，再用一个较小的查找窗分别放在主像片中心，用另一个一样大小的查找窗在从像片中滑动，这样可以得到像对矩阵。采用强度互相关算法计算出一个相关系数矩阵；

（2）在相关系数矩阵中查找峰值位置，即整数像素的偏移；

（3）再将从像片中查找窗移至峰值点位置，并利用相关系数插值算法在峰值周围对相关系数平面进行插值；

（4）在插值的相关系数矩阵中查找峰值位置，利用高斯回归方法获取子像素偏移量，即子像素偏移；

（5）将整数像素偏移和子像素偏移加在一起，构成位移场。

由于两景影像的成像时间存在差异，大气和光照的影响会使两景影像视觉效果上存在一定差异，降低影像对的相关性。并且在没有高精度的DEM 数据的情况下，地理坐标系的定位也可能存在差异，所以有必要在进行偏移场计算之前，对影像对进行预处理来消除这些影响。由于云、冰盖等因素的影响，部分区域的位移值可信度较低，需要利用峰值相关系数和相关系数比值来除去这些可靠性低的点。本文利用峰值相关系数和相关系数比值设置阈值，对云或冰盖覆盖区域进行掩膜处理。经过掩膜处理后的位移场仍然存在着由于成像过程和预处理产生的长波长畸变，长波长畸变需要采用适当的后处理算法来消除掉。经过以上处理便可得到平滑、正确的位移场。

2 应用实例

2.1 2014 年新疆于田地震

2014 年2 月12 日17 时19 分，我国新疆维吾尔自治区和田地区于田县发生强烈地震，据中国地震台网测定，震中位于36.1°N、82.5°E，震源深度12 km，震级Ms7.3（CENC）。此次于田地震震中位于阿尔金断裂西南段尾端的分支断裂的交汇处，为拉张与挤压变形的过渡带（图1），地质构造复杂，具有发生强震的背景[19]。经野外调查显示，这次地震在高海拔地区形成了由一系列张裂隙、张剪裂隙、剪切裂隙以及挤压鼓包和裂陷等斜列状组合而成的地表破裂带,整体呈NEE—SWW 走向，全长约25 km，分为不连续的两段，显示出左旋走滑伴随有正滑分量的特征，最大左旋位移约1 m[20]。

图1 2014 年2 月12 日于田地震构造背景及震中位置Fig.1 Construction background and epicenter location of Yutian earthquake on Feb.12，2014

图2 2014 年2 月12 日于田地震预处理后的光学遥感影像对Fig.2 Preprocessed optical remote sensing image pair of Yutian earthquake on Feb.12，2014

图3 影像对的峰值相关系数图Fig.3 Peak correlation coefficient map of image pairs

本文从美国地质调查局USGS 上选取了2013 年5 月31 日和2014 年5 月18 日的两景Landsat TM8 遥感数据分别作为监测2014 年2月12 日于田地震同震位移场的震前图像和震后图像，并选用图像分辨率为30 m 的第七波段。利用ENVI 软件，将两景图像进行影像增强、初配准、纹理增强、影像剪裁和灰度值匹配的预处理，获得一对7000×7000 像素的影像。图2a 和2b 分别为预处理之后的震前图像和震后图像。之后我们选用256×256、300×300、360×360、512×512、600×600 等像素大小的查找窗口，在行和列方向的采样间隔均为100 像素，利用上述子像素偏移追踪算法对预处理后的影像对进行计算得到位移场。图3 为峰值相关系数图。本文设置相关系数阈值和相关系数比值对偏移场掩膜，获取了2014 年于田地震同震偏移场（图4）。本文利用远场选取了部分位移场数据分析了算法的统计误差。经计算，东西向位移场误差为0.4170 m，南北向误差为0.4073 m，精度为百分之几个像素（图5）。

图4 2014 年2 月12 日于田地震同震位移场Fig.4 Displacement fields of 2014 Yutian earthquake

从我们得到的同震偏移场（图4）可以看出，两条黑色线段上部呈北西向运动，黑色线段下部呈南东向运动，两个条带存在左旋走滑趋势的断层，分别为北硝尔库勒断裂和阿什库勒—硝尔库勒断裂，断裂以左旋走滑为主，带有一定量的正滑分量。北硝尔库勒断裂和阿什库勒—硝尔库勒断裂全长约20 km，整体走向约N60°E，最大偏移量可达到约1.5 m，与李海兵等（2014）野外地质调查结果一致。震源南方地貌复杂，不易调查，本文得到的位移场表明该区域存在明显的地震地表破裂。

图5 2014 年2 月12 日于田地震同震位移场误差图东西向误差Fig.5 Error map of 2014 Yutian earthquake displacement fields

随着印度板块的向北俯冲，青藏高原沿着郭扎错断裂、阿尔金断裂中段和东段向北东方向推挤。阿尔金断裂在西南尾端分为两支：一支为康西瓦断裂呈西南走向延伸；另一支为郭扎错断裂整体呈现北东向走滑形式，并带有张裂特征[21]。2008 年于田Ms7.3 级地震的破裂带主要呈NS—NNE 向，位于2014 年于田地震的SWW 方向约120 km（图1），该断裂为左旋张裂破裂带[22]。2014 年于田地震正处于三条断裂所夹的三角区域内，可能是在这些影响下发生了2014 年于田地震。

2.2 2013 年巴基斯坦地震

2013 年9 月24 日的巴基斯坦西南地区发生Mw7.7 地震，经Global CMT 定位地震震源位于26.70°N、65.04°E，震源深度12 km。地震位于欧亚板块、印度板块和阿拉伯板块的交界处附近，该次地震也是板块间的相互作用引起的。

图6 2013 年9 月24 日巴基斯坦地震预处理后的光学遥感影像对Fig.6 Preprocessed optical remote sensing image pair of Pakistan earthquake on Sep.24，2013

与上面的例子相同，从美国地质调查局USGS 上选取了2013 年9 月10 日和2013 年9月26 日的两景Landsat TM8 遥感数据分别作为监测2013 年9 月24 日于田地震同震偏移场的震前图像和震后图像（图6），并选用图像分辨率为15 m 的第八波段。由于单景影像无法覆盖断裂区域，我们选取了两个连续景的光学遥感数据构成了两个影像对。经预处理，再分别求出位移场，最后将两个位移场绘制到一起构成一幅覆盖2013 年巴基斯坦地震震中区域的位移图。本文得到的2013 年9 月巴基斯坦地震同震偏移场如图7 所示。同样，我们在远场处选取了部分位移数据，统计分析位移计算误差。结果表明东西向偏移场误差为0.5318 m，南北向误差为0.3931 m，精度为百分之三个像素（图8）。

图7 2013 年9 月26 日巴基斯坦地震同震形变位移场Fig.7 Displacement fields of Pakistan earthquake on Sep.26，2013

所得同震偏移场给出，该地震断裂为左旋走滑型的圆弧形断裂，基本无正断分量，平均滑移量6 m，最大可达到10 m，并且断层西北侧的滑移量明显大于东南侧，结果与Barnhart（2014）、Avouac（2014）的研究结果一致[23-24]。

图8 求取的2013 年9 月26 日巴基斯坦地震同震位移场的误差Fig.8 Error map of Pakistan earthquake displacement fields on Sep.26，2013

3 讨论

本文利用一种新子像素偏移追踪算法计算了2013 年巴基斯坦和2014 年新疆于田两次地震的同震形变场。该子像素偏移追踪算法将强度互相关算法、双线性插值算法和二维高斯回归算法结合到一起。应用实例表明，该子像素偏移追踪算法测量精度可以在百分之三个像素水平上，可以很好的运用到求取中强震的同震形变场中。该算法精度略高于目前被广泛使用的ImCorr 和COSI-Corr 算法。由于本文算法内部使用强度互相关，因此该算法可以使用任意大小的查找窗口进行计算。该算法在后处理、遥感影像的条带状残差、两幅影像合并过程等方面仍需要改进。在中强震的偏移场监测中，可以与InSAR、GPS 相结合，进行互补，对该算法进行优化。