亚像素相位相关法在获取汶川地震近场形变中的应用

刘云华 单新建 屈春燕 宋小刚 张桂芳

(中国地震局地质研究所，地震动力学国家重点实验室，北京 100029)

0 引言

地震是一种以力学运动为主体的自然现象。因此，测量并掌握地壳的力学状态与力学运动，是地震学与地震预报研究中最为重要的基础性工作(国家地震局科技监测司，1995)。获取大范围的地表连续覆盖的同震位移形变场图像，并利用这些地表同震形变场来揭示震源破裂面的断层位移和运动速率等的几何学和运动学特征，对深入了解地震重复间隔、潜在震源区等问题，从而认清地震应力应变的积累与释放过程，具有重要意义。地壳形变监测是InSAR技术应用较为成功的领域之一，自Massonnet等(1993)在《nature》上发表了利用ERS数据获取的美国Landers地震形变场以来，各国学者开展了大量的相关研究工作，成为研究地震形变场较为流行的一种方法。然而该方法也有其弱点:1)对相位失相关很敏感，因此对于震中失相关严重地区无法提供形变细节;2)为获得形变场必须进行相位解缠，需要进行大量计算;3)只能监测到雷达视线方向上的一维形变量，不能完全反映监测地区的形变(Van Puymbroeck et al.，2000)。

为了解决上述问题，一些学者开始探索使用光学影像偏移量法来作为InSAR技术的补充(Van Puymbroeck et al.，2000)。人们发现地表形变还可以通过计算震前与震后的两景卫星图像之间的偏移量来获得，这个偏移量反映了两幅图像上地表指定点的位置偏差。在该思想指导下，近年来利用亚像素级相位相关性获取地表形变逐渐成为研究的热点。并在一些震例(Massonnet et al.，2000;Michel et al.，2002;Dominguez，2003;Binet，2005)中获得了较好的结果，能够获取震中附近的形变量。

汶川地震后，一些学者利用ALOS数据获取了地表InSAR同震形变场(Linlin et al.，2008;孙建宝等，2008;单新建等，2009;屈春燕等，2010)，为理解此次地震的形变过程和发生机理，提供了重要的观测资料。但由于沿断层带形变梯度较大，在形变中心带出现了1条无干涉条纹的非相干条带，因此无法获得断层带近场形变量，使得对了解汶川地震总体同震位错量有一定难度。地震后，国内很多研究团队对汶川地震地表破裂带进行了野外实地考察和测量，取得了大量的第一手宝贵资料(李海兵等，2008;任俊杰等，2008)，但是人们在同震位移、地表破裂带长度等方面还存在认识上的差别，主要有3方面的原因:不同的测量方法、不同的观测地点以及很多观测点无法到达的限制、对同一现象的认识不同(李海兵等，2008)。其中，最主要的是由于汶川地震造成大量山体滑坡和崩塌，使道路不通，很多有地表破裂的地区无法到达，或被深埋，或破裂带穿过陡峭山坡无法测量。而遥感技术以其覆盖范围广、不受地面条件限制等独特优势，可以成为野外实地测量的有益补充。因此，本文拟利用光学影像相位相关估算偏移量法，来获取发震断裂带附近的形变量。

1 原理介绍

利用光学影像偏移量来获取断裂带附近的形变量需要图像配准达到亚像素级，这是由于地震形变的量级一般小于遥感影像的分辨率(少数大地震可以达到m级)。并且由于卫星的轨道运动、相机的扫描运动及数字高程模型的误差等都给精确配准带来了困难，需要用严格成像模型来修正原始遥感影像所存在的几何变形并投影到合适的参考系统中。成像几何模型反映了地面点3维空间坐标与相应像点在像平面坐标系的2维坐标之间的数学关系。另外，地面形变测量的精度，还依赖于图像重采样及相关性估计技术。在图像配准前，遥感影像需要投影并重采样到同一参考系下，一般有2种做法:方法1为选一主图像为参考基准，将从影像配准并投影到主影像的参考系中，雷达干涉测量中采用这种方法。方法2为将主、从影像均投影并重采样到独立于卫星成像几何的地面参考系中，之后对影像的变化检测均是基于地面参考系的。本文将采用方法2对SPOT数据进行处理来获得地表形变，整个过程涉及到SPOT正射纠正、影像重采样、相关性估计等步骤，下面从这3个方面来进行介绍。

1.1 SPOT影像正射纠正

SPOT卫星为推扫式成像系统，即光学成像系统固定，随着卫星平台的向前移动，逐行以时序方式获取沿轨道的连续影像条带。SPOT 1A级产品是提供给用户的经过最低预处理级的影像，用户可在此基础上根据需要进行严格的集合校正，1A及数据各扫描行的像素与相应地面行的地面点之间保持严格的中心投影关系。依据SPOT提供的高精度卫星星历数据、姿态数据及时间数据，建立地面点坐标和像点坐标之间的几何关系进行正射纠正，这种几何关系可以通过正解法与反解法2种方式建立。

正解法最终的目的为求解像点P(c，r)在地心坐标系的视向矢量→u3与高程为h的地球椭球面的交点，期间要经过本体坐标系到轨道坐标系、再由轨道坐标系到地心坐标系的转换，用到的坐标系统如图1所示(Riazanoff，2002)。

正解法将像平面中规则的像素网格投影到不规则的地球表面上，在大范围尺度上，会由于地球自转导致畸变，小尺度范围内也会因为成像系统姿态变换及地形起伏而使图像变形。造成纠正影像上所得之像点非规则排列，有的像元素内可能“空白”(无像点)，有的可能重复(多个像点)，就难以实现灰度内插(张永生，2000)。由于受到图像分辨率的限制，要想测量地震形变，图像的重采样必须要保持来自原始数据中的亚像素级信息，因此以反解法为宜。

图1 参考坐标系Fig.1 Reference coordinate system.

图2 反解法示意图Fig.2 Schematic diagram of indirect image formation.

反解法是从空白的输出图像阵列(标准空间)出发，亦按行列的顺序依次对每个输出像元点位反求其原始图像坐标(畸变空间)中的位置。如图2，首先在地面建立规则的格网，因为要想对配准影像进行比对，则需要将它们都校正到相同格网上，这样2幅影像就具有相同的分辨率，不会在相关计算中出现影像像元错位。

为避免投影后造成像点不规则问题，Leprince等(2007)提出了一种反变换模型，以满足严格的重采样需求。下面对这一模型进行介绍。

地面点的位置仅与成像方向→u3和高程h有关，若把h确定下来，则成像方向与椭球的交点必定位于高程为h的椭球上，如果能通过像点坐标的初始值(c，r)初步确定→u3，然后通过迭代过程逐步修正(c，r)，即可获得实际的成像方向。

模型假设地面任一给定点都能在像平面上找到一个并且是惟一一个对应像点。如图3所示，M为地面一给定点，该点的高程由DEM数据经3次线性内插得到，并且为WGS84坐标系统。

考虑视线矢量U3(c，r)对所有的c，r=1，…，N，拓展到连续的二维实数空间，即u3(x，y)和(x，y)∈R2。那么在像空间中寻找像点坐标(x，y)即相当于求点(x，y)∈R2使得函数Ф(x，y)取得最小值:

式(1)中M'(x，y)为从视向矢量→u3所看到的地面点，另为视向矢量的卫星位矢，假设光线在大气层中为直线传播，那么的光线延长矢量，如果 M'位于 s→和地形面的交点，那么由于地形面的非线性特性而使得计算M'坐标变得相当繁琐，为了简化问题，为每一个点M构筑一投影平面P(M)，使得M'点实际位于该平面上，该投影平面P(M)经过点M并且垂直于(图3)。而M∈P(M)，使得方程(1)的解并没有变。所有位于投影平面P(M)的点M'(α，β，γ)必须满足，因此该投影平面的方程可以由下式决定:

图3 反变换模型(据Leprince et al.，2007)Fig.3 Inverse orthorectification model(after Leprince et al.，2007).

s和投影平面P(M)相交于

其中

则反变换的解可以通过使下式取得极小值得到

其中

这个过程是一个反复迭代的过程，经测试两点梯度法(TPSS)在稳定性及效率方面均比准牛顿法、最速下降法要好(Barzilai et al.，1988)。

计算得到的最小值存储在2个与地面格网相同维数大小的矩阵中，如果格网左上角顶点坐标为(E0，N0)，格网的分辨率为r，在像平面中的像点［X(i，j)，Y(i，j)］即被投影到地面对应的(E0+i·r，N0－ j·r)点上。

1.2 影像重采样

由于计算后的点位置多半不在原图的像元中心处，即在原始影像中的相应位置的坐标值不为整数，因此必须利用在原始影像中该点附近的若干像元的灰度，利用重采样进行赋值。根据奈奎斯特采样定律，采样频率应＞2倍信号最高频率，理想的重采样函数为辛克(SINC)函数，在影像处理过程中经常用到的最近邻法、双线性内插及3次卷积等方法，在重采样过程中会造成一定程度的失真，会影响后续的相关性处理精度。为此，Leprince等(2007)构造了一个二维重采样核，具有理想的低通滤波形式:

d x和d y为重采样距离，代表分别在x和y方向上的邻采样的最大间隔。

1.3 相位相关法估算位错量

相位相关法的理论依据是傅立叶变换的相移定理，相移定理是指空间域内函数的位移将会引起频域内变换函数的相移，因此一对配准图像之间的相对位移可以从它们傅里叶变换后的相位差得到。

假设两幅影像i1和i2之间只存在位移关系，平移量为(x0，y0)，即则i1，i2对应的傅里叶变换I1，I2之间的关系为对应的频域中2个图像的互功率谱为

此函数在偏移位置处有明显的尖锐峰值，其他位置的值接近于零，估算峰值处的坐标就可找到两幅图像间的偏移量。

2 数据处理流程

2008年5月12日发生于龙门山断裂带的汶川地震举世震惊，造成了巨大损失，地表破裂带长约240km，汶川地震复杂的滑动分布可能表明几条活动断层同时参与了地震破裂过程(Xu et al.，2009)。震后不久葛林林教授便发表了用InSAR技术测量的地表形变初步结果(Linlin et al.，2008)，体现了遥感技术在时效性及覆盖范围方面的优势。然而由于四川地区多云天气较多，获得无云的可见光数据比较困难，因此一直未见有利用光学影像估算地表形变的文章发表。

本文选取了震前与震后各一景SPOT数据作为研究对象，级别为L1A级。基本参数见表1。

两幅影像时间上相距22个月，但均在冬季成像，这保证了地物差别不会很大。选取的1A级产品是SPOT数据经辐射校正处理后的产品，包含了用以进行后续的几何校正处理的辅助数据。正射纠正中用到的DEM数据为最新发布的ASTER GDEM version1，地面分辨率为1″(约30m)，在全球范围内满足垂直精度为20m的置信度为95%(NASA，2009)。下面介绍数据处理流程。

(1)选取1/10000的已经进行正射纠正的SPOT数据作为参考影像，对震前数据image1进行正射纠正，从影像上选取了12个GCP点，注意GCP点的选取尽量均匀分布在图像中间并避开断裂带，以避免后续的卫星视向矢量优化时将地震形变量抵消。

表1 SPOT数据相关参数Table 1 Relevant parameters for the SPOT scenes

(2)从1A级数据的头文件中提取辅助数据，并按照前面介绍的反变换模型对所选取的GCP点进行视向矢量优化。

(3)利用优化后的GCP点对image1进行正射纠正处理;得到经过精确纠正影像i1。

(4)重复上述步骤对震后数据image2进行处理，不同之处在于此时选取已处理好的i1为参考影像，得到精确纠正影像i2，且保证了i1和i2之间的配准。

(5)按照2.3节中的方法对i1和i2进行相位相关计算，选取的滑动窗口大小为32×32，步长为8。

需要补充说明的是，根据卫星自带的轨道数据和一些手动选取的初始控制点，得到1个由在参考影像和待配准影像上选取的一系列同名点对构成的改正矩阵，运用间接正射校正模型和改正矩阵，将原始影像即待配准影像投影至地面坐标系，采用相位相关技术对投影后的影像进行匹配运算，得到二者在EW向及SN向上的坐标偏移量及SNR值，根据该坐标差值在DEM的基础上就构成了一次迭代后的新的地面控制点;返回步骤2，再将上步得到的新的控制点作为初始控制点数据，生成新的改正矩阵。如此经过反复迭代直到方差达到稳定，即可得到精确的地面控制点，然后便可利用这些控制点进行正射校正和配准。如果轨道数据的精度已可满足初始正射校正的需要，则手动选取控制点的步骤可以省略。

3 结果分析

图4a为正射纠正后的震后影像image2，图4c，d分别为EW向及SN向的位错分布图，而图4b为相位相关计算过程中的SNR图，它代表了相位相关计算的质量，范围从0(去相关)到1(相关性很好)。从该图上可以看出去相关分布区域，去相关即失去相关性，表现为SNR较低或无值(Nan)，当相关性算子无法覆盖区域或位错量太大超过设定值(20m)会出现无值的情况。导致去相关的原因可归纳为3种:1)云层造成的去相关，在震后影像image2(Fig.4a)龙门山前缘上空漂浮着一些云朵，它们所在的位置和SNR图中去相关区域相对应;2)积雪覆盖及大量滑坡坍塌处，在龙门山中央海拔较高，覆盖有积雪(图像中央)，而在龙门山镇附近景区有大量坍塌(图像左下角)，导致了这些区域去相关面积较大;3)河流，在影像的右下方有2条河流经过，造成了2条线状去相关带。

图4 计算结果图Fig.4 Calculation result map.

SPOT4影像的扫描带宽度为60km，因此选取的这景影像仅覆盖茂县、什邡地区，看不到汶川地震破裂带的全貌，但在图幅覆盖范围内，地震造成的2条破裂带清晰可见。图中红线为根据图像不连续性画出的断层位置。在位错分布图上可以看出，地表破裂带主要分布于龙门山中央断裂带(北川-映秀断裂)及前山断裂带(彭灌断裂)上，在后山断裂带(茂汶断裂)未见明显的地表破裂带。横跨断层画了6条剖面，AA'、DD'、EE'、FF'横跨龙门山中央断裂带。AA'位于龙门山镇附近，该处水平缩短量最大，达到10m左右;FF'位于岳家山附近，水平缩短量最大可达6m;EE'位于清平乡附近，水平缩短量为2～4m;DD'位于高川乡附近，水平缩短量又有增大的趋势，最大可达5m，并且在该处断裂带有一转向。从EW向位错及剖面图上可看出，中央断裂带两端位错量大，中间略小，并且有分段的趋势，在高川附近有一横向转折，发震断层西北盘为抬升盘，南东盘断层附近仍然表现为隆起区，显示龙门山中央断裂带以逆冲为主、兼具走滑的断层性质，中央断裂带位错量一般达到4～6m，而龙门山镇附近的剖面AA'的最大位错达到10m左右。在龙门山镇附近显示有大片位错量较大的分布区，调查表明该地区为龙门山景区的回龙沟及银厂沟所在，地震造成自然景观损毁率达80%～90%。银厂沟的大龙潭、小龙潭受泥石流和崩塌山体阻隔而彻底破坏，银厂沟盖坪周围的居民点已被夷为平地，所有房屋都已垮塌(贾建中等，2008)。

图5 EW向位错及跨断层剖面图Fig.5 East-west displacement and profiles across fault.

BB'和CC'横跨龙门山前山断裂带，BB'位于蓥华镇附近，水平缩短量为2m左右，CC'位于汉旺附近，水平缩短量为2～4m，可见前山断裂带逆冲为主，在山前形成1条从八角镇至睢水镇长约40km的破裂带(因影像未能覆盖白鹿镇及以南地区，因此未能反映白鹿及小鱼洞地区的破裂情况)，位错量一般为1～2m，汉旺镇地区可达4m左右。综合以上几个剖面，可以看出中央断裂带位错量大于前山断裂带，且两端位错量大于中间。

震后国家地震应急指挥部立即组织科技人员进行了科考，开展了地震地表破裂考察，通过对地形地貌、地质体和地面人工建(构)筑物的错动变形观测，完成了地震地表破裂带填图和沿地震断层同震位移的测量等科学考察(徐锡伟等，2008)，取得了大量宝贵的现场资料，本文从中选取一些现场调查点对SPOT影像获取结果进行验证。按照文献中的描述，在龙门山镇东林寺(31.285 50°N，103.818 30°E)抬升2.5m，右旋位移2.2m，按本文方法在该处附近得到的位错量约为2.7m。在龙门山镇附近(31.294 194 4°N，103.848 472 2°E)抬升5m，本文方法在该处附近为6.8m。在绵竹市九龙镇沙坝村(31.981 1°N，104.118 36°E)NE向陡坎上测量到的最大垂直位移(3.5±0.2)m，本文方法在该处附近为4.6m左右。在汉旺清水河东岸通往清平乡的公路上(31.461 56°N，104.166 39°E)，可见废弃公路在破裂带上地壳缩短量为(1.5±0.2)m，本文方法在该处在1.8m左右。虽然通过SPOT影像得到的位错量较野外测量值偏大20%左右，但是地表位错量的分布趋势达到了很好的吻合。遥感图像由于噪声的干扰，某一点的值和周边的值具有跳跃性，因此最好选取一定大小窗口的平均值作为该点的值。野外测量的值是3维的，而遥感影像计算得到的值只反映水平方向的位错量，因此在比较时我们更注重于整体形态的分布。图6为该文献中野外调查点与本文计算结果的叠合图，除岳家山至高川段因当时条件恶劣难以到达而没有取得调查点外，图中调查点都分布在根据遥感影像计算得到的地表破裂带上。表2为野外调查的滑动位错量与根据遥感影像计算得到的水平位错在对应区域的对照表。相比而言，野外调查是离散的各个观测点，而利用遥感影像获得的形变场则是连续的，更能反映破裂带的整体情况，包括整个破裂带分布的确切位置、沿走向的形态变化等。

表2 野外调查与遥感图像计算对比表Table 2 Comparison between field survey and computation from remote sensing image

4 结论

本文以覆盖茂县地区的SPOT影像为数据源，采用亚像素相位相关位错法计算了影像覆盖地区的水平位错量，取得了以下认识:

(1)在形变较大的近断层地区，In-SAR失相关严重，无法获得断层近场形变信息。本文的研究表明，利用光学影像相关法能够获得近断层位错量，可以成为InSAR手段的重要补充。但该方法对影像配准要求较高，需要加入卫星成像几何、姿态能辅助数据，利用反解法对SPOT数据进行了精确正射校正，对校正后影像再采用频域相关法计算。

图6 地表破裂与地面调查点对比图Fig.6 Surface rupture inferred from offset map compared with field survey points.

(2)汶川地震造成龙门山断裂带上至少2条断裂同时发生破裂，形成了主要地表破裂带(龙门山镇-高川破裂带)和次级地表破裂带(汉旺破裂带)，沿龙门山镇-高川破裂带平均位移量为4～6m，从跨断层剖面可看出，AA'水平缩短量最大，达到10m左右，向北至FF'、EE'开始减小，在北端DD'又开始增大，并且在高川附近有一转折;汉旺破裂带基本为纯逆冲性质，平均位移量一般为1～2m;与野外实际调查数据相比，位错量值偏大20%左右，但地表位错量的分布趋势与野外调查结果相吻合。

(3)该方法虽具有一定优势，能够监测水平方向运动，如对于走滑断层水平位错量，以及对于逆冲断层的水平缩短量都是可以监测的，但对于垂直方向位错很不敏感，且数据来源受到天气情况的制约，云量大会造成图像去相关而无法获得满意的结果。

本研究实现了亚像素级光学影像相位相关算法获取地表形变，并选取了汶川震区中的茂县、什邡为实验区，计算结果得到一些重要的地质信息，对于理解汶川地震的发震机理具有重要意义。研究结果表明，该方法具有较好的应用前景，为获取断层附近大形变分布提供了可能。