1 引言

强震破裂断层的位置、走向和长度等几何参数，对于震源机制的反演与理解具有重要作用.当前，对于地震断层地表破裂的识别仍主要依赖野外地质调查工作，但内陆型强震多发生于地形、地貌复杂的困难山区，存在大面积无法直接到达的区域，导致野外调查数据不完整，进而影响对断层位置、走向等重要几何参数的准确调查.

星载合成孔径雷达(SAR)像素偏移方法可观测沿卫星飞行方向(方位向)和雷达视线向的二维地表形变场，具有覆盖范围广、监测密度大、全天候观测等独特的技术优势(廖明生等, 2003; 曾琪明和解学通, 2004; Hamiel and Fialko, 2007; Prati et al., 2010; Hooper et al., 2012; Chen et al., 2014)，并在以冰川、沙丘、地震等为观测目标的研究中取得了有意义的科学进展(Strozzi et al., 2002; Fialko et al., 2005; Kobayashi et al., 2009; 许才军等, 2009; Casu et al., 2011; 刘云华等, 2012; Fielding et al., 2012; 杨莹辉等, 2014; 占文俊等2015).雷达影像像素偏移方法可用于测量震后断层近场地表形变，进而提取断层的位置、长度等几何参数信息，较常规的野外地质调查工作，其在提取精度、时间效率等方面具有明显优势.特别是近年来随着TerraSAR、COSMO-Skymed、Sentinel-1A/1B、ALOS-2以及GF-3等多颗雷达卫星的陆续发射升空，基于SAR像素偏移法测量同震地表形变场并提取同震断层破裂得到越来越广泛的应用.

但当前像素偏移技术提取地表形变场精度仍受到多种因素制约，其中电磁波传播受电离层干扰对像素偏移法的影响不容忽视(Meyer et al., 2006; Liu et al., 2010；Raucoules and De, 2010).电离层不仅引起InSAR干涉相位梯度和距离配准偏移量变化，电子密度扰动还造成SAR相位延迟，导致方位向位置偏移，已有研究已经发现采用offset-tracking技术，从SAR影像中提取的地震、冰川等方位向形变场中出现了典型的由电离层信号引起的“方位向条纹”(Joughin et al., 1999; Gray et al., 2000; Wegmuller et al., 2006).

电离层效应引起的SAR方位向形变场误差，导致地震地表形变信号被淹没于噪声中，使得基于断层近场形变的同震地表破裂线提取变得非常困难，极大地限制了对地震造成地表破坏和震源机制的理解.因此，如何精确地校正SAR偏移形变场中的电离层误差，并基于断层近场形变准确识别和提取同震地表破裂线显得尤为重要.本文在分析电离层对SAR像素偏移技术影响的基础上，采用了一种基于频率域的电离层误差校正方法，去除方位向形变场中的电离层误差，进而获得高精度的汶川地震SAR方位向形变场.在此基础上，提出一种基于SAR偏移形变场的断层破裂线自动提取算法，获得了精确的汶川地震同震断层破裂几何模型，并引入汶川地震野外地质调查地表破裂线数据，定量分析两者间的差异，验证了所提取地表破裂线的可靠性.

2 合成孔径雷达电离层影响

2000年，Gray等(Gray et al., 2000) 指出了电离层分布会影响卫星雷达数据且电离层电子密度扰动会导致SAR成像地面目标方位向位置的偏移，SAR方位向形变场中出现的非构造形变条纹即是受电离层活动影响的标志，且已有研究发现长波段雷达数据更易出现大的方位向偏移，目前合成孔径雷达成像中常用的L波段雷达数据在两极区域观测到的SAR方位向误差甚至能达到几个像素.电磁波通过电离层由电子分布导致的雷达波路径变化可由下式表达：

(1)

其中，δ_s表示由电离层引起的路径长度变化，a=40.3为常量，f为电磁波频率(单位为Hz)，N_t为每平米电子数量(Total Electron Content，简称TEC)，通常在描述单位面积的电子数目时，以TEC为单位度量，1 TEC=10¹⁶e·m^-2.而与之相关的由电离层引起的往返路径相位改变可根据下式计算：

(2)

其中，λ为电磁波波长.由此可以计算出C波段的雷达波经过1TEC电离层会导致约π的相位变化，而L波段则会产生约4π的相位变化.

对于单幅SAR影像，由SAR侧视成像原理可知，影响SAR数据的电离层不是位于成像区域的垂直上空而是地面目标点到卫星的雷达视线方向的电离层，因此同一地表位置上空的电离层会影响不同的雷达波路径变化，同时电离层电子密度分布差异还会引起合成孔径相位延迟变化，在方位向压缩中会导致方位向散焦，这种扭曲引起的方位向相位变化会引起方位向位置偏移(Wegmuller et al., 2006).

对于一对SAR影像，两次成像的卫星轨道不可能完全同轨，也会造成同一地物点由电离层造成的回波路径变化不一致，而且电离层自由电子本身也随着太阳活动、地球磁场和大气环境的改变而发生变化(李亮等, 2014)，因此基于星载SAR数据获得的雷达视线向和方位向形变场中包含了电离层影响分量.

3 电离层影响的频域校正方法

在采用像素偏移法(offset-tracking)提取的方位向形变场中，由地震引起的同震地表形变信号与电离层引起的路径变化信息在观测数据中存在显著差异.同震形变信号在方位向形变场中集中于断层近场，且形变量级表现为近场大，远场小；而电离层引起的路径变化在方位向形变场中呈现规则性条纹，条纹倾向一致，正负交替出现.考虑到这种呈规律正负交替变化的电离层条纹在频率域中表现为高频分量，为达到有效去除的目的，将方位向形变场转到频率域研究更为合理.本文利用电离层信号对应的高频信息，校正电离层自由电子在方位向形变场中引入的误差分量(Kobayashi et al., 2009)，校正流程如图 1所示，具体的“方位向条纹”校正方法与过程如下.

(1) 将受电离层影响显著的方位向形变场经傅里叶变换转到频率域.图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度.分析方位向形变场中电离层条纹的走势，可判断其沿条纹垂直方向灰度梯度变化最大，对应频率域中高能量值位置(椭圆区域).

(2) 提取电离层条纹对应的频率分量，根据两个条件联合约束，其一为位置范围信息，以对应两个角点为圆心，设置半径r，确定频率范围；其二为频域能量值，电离层条纹正负交替，对应的灰度梯度大，能量值大，设置能量阈值e.需要说明的是，不同的SAR影像因其轨道、获取时间和雷达波长的差异，半径r与能量阈值e可能差异较大，需根据实际具体情况设定.半径r可灵活设定，一般完整包络两角点附近的高能量值即可.能量阈值e需根据不同数据频谱能量值的分布确定，实验中需先求取全局频谱能量的均值M与标准差σ，考虑电离层信号在频率域呈现为高频信号分量，可设定能量阈值e>M+2σ，位于圈定半径范围内，所有能量值超过阈值e的均为电离层条纹高频信号分量.

(3) 将提取的电离层频域信息经傅里叶逆变换，即可获得电离层影响在空间域中对应的信号分量.

(4) 从原始方位向形变场中移除上述计算得到的电离层信号，进而获得校正后的方位向形变场，为进一步研究地震机制、地表破裂情况提供基础数据.

4 实验与结果

2008年汶川地震的矩震级为7.8级，造成断层近场发生显著的地表运动，地震近场的地表破裂极为严重(李志才等, 2009; Liu et al., 2009; Chen et al., 2015; 邹镇宇等, 2015)，本文以覆盖汶川地震主震区6个条带PALSAR影像为实验数据，开展SAR形变场电离层影响的校正和断层地表破裂线的提取，所使用的PALSAR影像的具体参数如表 1.

4.1 汶川地震方位向形变场的校正

采用像素偏移法(offset-tracking)对PALSAR影像进行精配准并获得高密度方位向偏移量，建模去除轨道系统偏移后，得到如图 2所示的方位向形变场.从图 2中可看出，由于形变场中存在显著的非地震形变趋势条纹(电离层对方位向形变场的影响表征)，掩盖了真实的地表方位向形变分量，无法正确判断出地表在卫星飞行方向上的形变情况.此外，六条带SAR影像的获取时间不同，而电离层自由电子密度随时间变化，也造成不同条带受电离层影响程度存在差异，如471-473和475-476条带受电离层影响严重，而474条带中电离层影响较弱.

为获得完整且准确的方位向形变场，本文采用基于频率域滤波的方法(上述第3节校正方法)去除这种条纹噪声.以472条带为例，图 3a为原始472条带的方位向形变场，图 3b为根据图 3a中电离层信号对应高频信息变换获得的空间域电离层分量，可以清晰看出图 3b与3a中电离层信号无论在形变量级还是分布走势都呈现高度相似性，从图 3a中移除电离层分量(3b)得到(3c)，从图 3c中可清晰分辨出汶川地震同震近场(黑色椭圆区域)沿方位向的形变趋势.

依次移除各条带方位向形变场中出现的电离层影响条纹，得到校正后的汶川地震方位向形变场(图 4).从图中可清晰判断出地表断裂线及同震位移的空间分布特征(黑色矩形框内).尤其是471-472条带，校正前电离层信号分量完全淹没了地表形变信息，该区域地表的运动模式和量级均无法有效地从原始数据中获取，但经过电离层校正后，该部分形变场精度得到极大改善，从中可清晰地判断出断层上盘北向和下盘南向的运动模式.校正后形变场显示在北川区域存在明显的断层走向变化，这与实际的野外调查一致；此外，474和475条带的形变场质量也得到改善，在绵竹附近(图中黑色箭头所示区域)发现，映秀—北川断层造成的地表形变在此处受到了另一条断层即彭灌断裂运动的影响.

4.2 地表破裂线的提取

基于同震形变场的空间分布特征提取断层的地表破裂痕迹是一种有效弥补野外地质调查数据不足的方法，有助于进一步分析震源机制.为得到准确的断层地表破裂线，本文采用形变梯度边界提取策率，基于校正后的SAR方位向形变场识别断层地表破裂线，具体提取过程如下.

(1) 采用中值滤波法对形变场进行去噪处理.中值滤波相比高斯滤波等其他去噪方法，在去除形变场噪声的同时能有效保留边缘线的清晰性，不会造成显著的边缘模糊化(Gonzalez and Woods, 2002; 赵高长等, 2011).

(2) 采用Sobel算子提取形变场的水平梯度和垂直梯度，进一步计算水平梯度与垂直梯度的平方和再开方，得到形变场的像素梯度(Kittler, 1983).

(3) 对上述(2) 提取的像素梯度进行非极大值抑制，提取形变场边缘点.根据水平梯度和垂直梯度值计算出可能的边缘线角度，取此边缘线两侧像素值(计算边缘线法线确定插值像素坐标)计算梯度M_left，M_right，判断中心像元的梯度M是否大于两侧M_left，M_right，如果都大于则认为中心像元处是极大值，标记作为边缘线备选点，如果至少小于一边，则剔除该中心像元.

(4) 为得到精确的边缘线点，选取高梯度阈值M_H和低梯度阈值M_L进行边缘像元的初步提取.其中备选点梯度大于M_H认为其为可靠的边缘点，进一步考虑到边缘线的连通性，在确定的边缘线点周围搜索梯度大于最低阈值的点，也标记为边缘线点，并基于上述边缘线点拟合粗略的边缘线.

(5) 为防止由于过高的阈值引起的边缘点遗漏，降低M_H和M_L，获得更高密度的边缘目标点，但由于阈值的减小，新获取的数据中存在一定比例的非边缘目标.

(6) 根据上述(4) 确定的粗略边缘线，约束边缘线点的范围，剔除(5) 中的非边缘目标点，根据剩余的边缘线，寻找误差最小的一、二阶多项式拟合出地表破裂线.

根据上述的地震地表破裂线自动提取流程，获得了如图 5所示的汶川地震地表破裂线，其中表 2列出了图 5中1-8号点的经纬度位置信息.提取结果显示，主破裂带北川—映秀破裂带从汶川县映秀镇向东北方向单侧破裂至青川县境内，断层平均走向角约46°，在北川区域出现明显的走向变化，破裂长度约300 km，破裂形态与已有研究结论一致(药晓东等, 2015; 朱守彪和袁杰, 2016).灌县—江油破裂带长约72 km，其破裂方向与北川—映秀断裂带近似平行，与北方向夹角约47°，并在破裂末端约5 km长度处断层走向偏向东南.

5 结果与分析

为验证上述SAR形变场提取地表破裂线的可靠性，引入野外地质调查获取的地表破裂数据进行对比.汶川地震发生两周后，为确定此次地震造成的地表破裂情况，Liu等人在北川—映秀破裂带设站156个，灌县—江油破裂带设站57个开展了实地野外调查，并公布了研究成果(Liu et al., 2009).图 6中的黑色实线即为以Liu等人的调查数据为依据绘制的汶川地震地表破裂线，其中红色实线为基于本文方法自动提取的汶川地震断层地表破裂线.对比两者发现，两者具有一致的断层走向，但自动提取的断层破裂线略长于野外地质调查结果，这是因为野外地质调查是基于断层的地表出露开展的，难以获得无明显地表破裂断层的具体情况，特别是无法正确识别出破裂线起始位置与结束位置，而断层的自动提取是基于近断层地表运动开展的，可有效提取获得整个断层，因此自动提取的地表破裂线长度超过野外地质调查结果是合理的.

为进一步对比两者的差异，我们绘制了断层末端和转折位置处的细节图像，发现在震中以及断层末端等位置(椭圆区域)，两条破裂线存在一定的差异.为定量地进行对比与分析，我们获取了每个椭圆内自动提取破裂线和野外地质调查获取破裂线的地理坐标，并以最左端为起算点，等间距计算两条破裂线的间距，最终获得两条破裂线间的距离变化曲线，如图 6所示.

由于整条破裂线长约300 km，从“a”到“e”地表破裂的情况有较大区别，影响因素也各有不同.首先在断层破裂起始位置的a处，提取的破裂线与野外调查破裂线之间的差距在端点处最大，为2.5 km，然后差距逐渐变小，最后两者完全重合(图 6a).造成此差异的主要原因为，一方面该位置地处震中，同震形变量级较大，地表破裂严重，使得地震前后SAR影像相干系数低，造成自动提取的破裂线位置受影响；此外，野外地质调查主要参考断层的地表出露位置，无法获得断层末端地表形变量级较小且无明显出露的断层分布情况，这些综合因素导致两种方法获得结果存在一定偏差.

在断层中段的b处，两条破裂线从完全重合开始，距离逐渐变大至2.6 km，后又迅速变小直至重合(图 6b)，这是因为b位置处于汶川地震地表运动最为剧烈的北川县城区域，自然地貌与人工建筑混杂分布，而剧烈的同震运动造成严重的人工建筑损坏与自然地貌变化，这些因素导致InSAR相关系数较低，自动提取断层破裂线位置与真实情况存在一定的差异；其次，在北川区域，汶川地震的断层破裂方向出现明显的变化，而野外地质调查相对离散的数据不利于精确反应断层的真实位置，上述综合因素最终导致了两者的差异.

断层中北段的两处差异区c和d呈现了相似的现象，即基于本文算法自动提取的断层线与地表正负形变的分界线完全贴合，但野外调查结果与此存在差异，考虑到自c位置开始，断层地表浅层滑动已迅速衰减，断层的地表出露情况相对减弱，这极可能影响了野外地质调查结果的可靠性，进而导致了两者之间的差异.因此，我们认为在这两处，自动提取的断层破裂线比野外地质调查结果具有更高的可信度.

汶川地震前山断裂的彭灌断层也存在一个明显的差异区(图 6e处)，此处两条破裂线的距离呈现先增大后减小的趋势，最大距离差为2.55 km，该断层浅层滑动量级相对较小，造成的地表形变特别是在断层末端并不显著，因此我们认为依据地表形变自动提取的断层破裂线在e处精度不高，造成两者不一致；但同时，由于无明显的地表出露，在野外地质调查工作中无法查明的断层破裂段(图 6中自动提取彭灌断层的末端)，在基于地表形变提取的断层数据中具有很好的体现，这也很好地凸显出了本文方法的有效性.

综上所述，从SAR形变场提取的破裂线与野外地质调查获取的破裂线整体走向一致，破裂线位置基本吻合，在几处存在差异的地方，经分析发现，其中映秀、北川区域由于地表破碎严重影响了自动提取地表破裂线的可靠性，北川至青川段，自动提取获得的地表破裂线相比野外地质调查数据与断层上下盘正负形变分界位置具有更好的一致性.此外，研究也显示基于断层近场地表形变的断层破裂线自动提取技术，具有发现破裂末端不具有明显地表出露情况的隐伏破裂的能力，这可很好地弥补野外地质调查在此方面的缺陷.

尽管校正后方位向形变场的精度对提取断层破裂线具有一定的影响，但考虑汶川地震为典型的浅源逆冲型强震，已有同震滑动反演与野外调查显示，断层近场上盘相对下盘存在显著的抬升运动，形变梯度差显著大于方位向形变场的观测误差，因此对破裂线提取的负面影响不显著.同时也需指出，针对中强型地震，特别是可造成显著地表破裂的地震，本文提出的方法具有较好的适用性；而对于小型地震或未造成显著地表破裂的事件，本文方法的应用会受到一定的局限.

6 结论

针对汶川地震方位向形变场电离层条纹噪声严重的问题，本文使用基于频域校正方法，从方位向形变场的频域中提取电离层信号对应的频域信息，通过建模获得空间域的电离层分量，进而从原始空间域形变场中扣除了电离层影响分量，校正并还原了近震区沿卫星飞行方向的真实地表形变场.此外，本文构建了一套实用的地表破裂线自动提取算法，采用边缘算子计算像素梯度并进行非极大值抑制得到边缘点，再根据高阈值提取可靠的边缘线，同时为保证边缘连通性根据较低阈值获得大量边缘线点，并以粗略边缘线约束其范围，获得具有较好连通性且可靠性较高的边缘点，最终提取得到汶川地震映秀—北川断层与彭灌断层的地表破裂线.

通过引入野外地质调查地表破裂数据验证本文方法的可靠性，发现从校正后方位向形变场中自动提取的地表破裂与野外地质调查获取的地表破裂线整体走向一致，震中和北川区域由于地表破碎严重导致自动提取的地表破裂线与野外调查结果相比存在一定的差异，其他地区自动提取的地表破裂线具有较高的可靠性.研究结果表明，基于电离层校正后的SAR方位向形变场能提取完整的断层地表破裂线，可作为野外调查地震断层地表破裂的重要补充手段，以弥补野外地质调查方法在完整确定断层破裂方面的不足.此外，利用本文方法提取的断层地表破裂数据，可为震源模型反演提供准确的断层几何参数，有助于提高反演断层滑动分布的可靠性，并可进一步联合用于计算同震应力传输、余震能量与分布估计等研究.