地震裂缝是受地震活动影响，沿震源断层错动在地表形成的具有一定长度、宽度和深度的裂隙^[1]，是同震地表破裂的表现形式之一，发育于各类型地震同震地表破裂带内^[2]。除地震活动本身造成的灾害外，地面永久变形和破裂产生的地震裂缝同样会造成房屋等建筑物的开裂倒塌^[3]，严重威胁人民生命财产安全。同时，地震裂缝的空间分布特征及与其他破裂类型(地震鼓包、地震凹陷等)的组合形式是震源断层活动性质在地表最直接的体现，查明地震裂缝分布规律，有助于深入理解震源断层性质、构造运动状态和地震破裂过程等基础信息，在震后进行发震构造动力学、几何学等相关研究。

受制于地形和自然条件等因素，地表破裂野外调查工作往往难以开展且存在一定的盲目性^[4]。遥感技术具有宏观性强、时效性高、覆盖范围广等优势^[5]，在地表破裂调查研究中发挥了不可替代的作用。前人借助遥感手段，及时准确地获取了近20年7级以上地震同震地表破裂的空间位置并测量了同震位移量，为震灾评估、地质调查提供了重要的科学依据^[6-9]。然而，上述工作大都通过目视解译完成，工作量大且主观性较强，如何通过计算机结合专家知识实现自动、半自动提取，降低目视解译工作量，提高提取精度与效率，成为目前国内外亟需解决的重要问题。

诸多地表破裂类型中，地震裂缝结构简单，在遥感影像上通常表现为低亮度，纹理与光谱特征明显^[10-11]，最易于提取。杨进生等^[12]结合多种卷积滤波器以及纹理信息提取方法对影像进行增强，有效识别了华北平原地裂缝信息。王娅娟等^[13]根据地裂缝的方向性特征，通过方向滤波有效识别提取了矿区地裂缝。然而上述方法仅仅停留在图像增强处理后再进行目视解译的阶段，并未实现自动、半自动提取。汤伏全等^[14]利用最大似然法结合随机森林算法，通过二次分类实现了采动地裂缝的精细提取。张兴航等^[10]利用Canny边缘检测算子提取线性特征，通过分形维数区分地裂缝与其他线性要素，实现了地裂缝的精细提取。受文献[15]启发，Stumpf等^[16]发现地裂缝的灰度剖面图可近似拟合为高斯曲线，并通过高斯匹配滤波(matched filter with first-order derivative of Gaussian, MF-FDOG) 算法进行提取。韦博文等^[17]将该方法应用于西部某山体及矿区无人机影像，结果表明MF-FDOG算法对黄土地区地裂缝提取具有较强的适应性。在前人研究中，受岩土体性质及其强度影响，采动地裂缝、滑坡地裂缝的分布规律性不强，而震源断层错动形成的地震裂缝与发震断裂带位置基本重合，除光谱特征外，空间分布也具有明显的方向性特征。本文选取昆仑山口西地震同震地表破裂带为研究区，引入国产高分二号(GF-2)作为影像数据源，对比分析目前地裂缝提取常用方法，建立了一套基于地震裂缝光谱、空间特征的地震裂缝提取方法，以期在震后第一时间快速获取并分析地震裂缝空间分布及组合关系，指导地表破裂地面详细调查、灾后应急救援以及地震科学相关研究。

1 研究区概况以及数据源 1.1 东昆仑断裂带

2001年11月14日，昆仑山口西发生M_S8.1地震，震源断层为东昆仑断裂带。该断裂以左旋走滑为主要特征，同时兼有部分挤压分量，自第四纪以来活动强烈，曾多次发生7级以上地震^[18-21]。地震形成的地表破裂带西起库水浣，东至卡巴纽尔多湖，全长达426 km，地表同震最大左旋位移量为7.6 m, 最大垂直位移量为4 m^{[2, 22]}。地表破裂延伸之长，规模之大，为世界之罕见。其中，构造性地震裂缝在整个破裂带上发育最为广泛。此次地震由于地处藏北高海拔无人区，气候寒冷干燥，无植被覆盖，破裂形迹至今大都得以完整保留，为地震裂缝提取研究提供了一个得天独厚的试验场。

1.2 数据源

研究区共收集3景原始GF-2数据，覆盖布喀达坂峰冰舌至红水河口西破裂带(表 1)。以15 m分辨率Landsat ETM+ 波段8全色影像作为参考影像，结合有理多项式系数(rational polynomial coefficient，RPC)数据和90 m分辨率SRTM DEM，使用ENVI5.3流程化工具RPC Orthorectification Using Reference Image进行正射校正，该工具可自动从参考影像上寻找控制点并应用于正射校正，极大提高了校正精度。采用Pansharp算法进行图像融合，融合后影像分辨率为0.8 m，色彩信息丰富，影纹特征清晰，完全满足昆仑山口西地震同震地震裂缝的识别及提取(图 1)。

2 地震裂缝提取方法

基于地震裂缝与其他地物的光谱差异，先后采用自适应直方图均衡化和阈值分割对地震裂缝进行粗提取，再通过数学形态学算法对地震裂缝孔隙进行填充和连接，提取出地震裂缝完整形态；最后根据地震裂缝的空间特征，排除其他地表信息(如阴影、纹沟、河床)的干扰，得到地震裂缝的精确提取结果。

2.1 地震裂缝粗提取

光谱信息是遥感影像中对象的颜色及灰度或者波段间的亮度比等，是组成地物成分、结构等属性的直接反映，是遥感影像特征提取区别于另一种地物的本质特征^[23]。地震裂缝因其具有一定的深度而导致光线反射率低，在遥感影像中表现为低亮度值，可通过阈值分割将其完整形态提出(图 2)。同时也可以看出，地震裂缝两侧存在明显的阶跃型边缘，传统边缘检测算子如Sobel算子、Canny算子会对阶跃型边缘做出响应，检测出伪边缘，因此该方法在此并不适用^{[14, 16]}。

为避免光照等外界条件影响造成的明部或暗部细节丢失，捕捉更多的图像细节，首先采用限制对比度自适应直方图均衡化(contrast limited adaptive histgram equalization，CLAHE) 对影像进行增强处理^[24]。均衡化之后裂缝特征更为明显，在此基础上进行阈值分割，单阈值往往无法得到清晰有效的阈值分割图像，同样采用自适应阈值(adaptive threshold)分割^[25]，通过高斯方程计算得到每个像素点的权重值，并使用该阈值对当前像素点进行处理，提取出地震裂缝完整形态。

2.2 地震裂缝精提取

当地震裂缝较窄或由于部分被遮挡时，需要对提取结果进行形态学滤波，以填充裂缝之间的间隙及裂缝内的孔洞，使其形状趋于规则。本研究使用形态学闭运算对裂缝孔隙进行填充和连接。闭运算是先进行膨胀再进行腐蚀，具有平滑图像、连接间隙、消除孔洞的作用^[26]。

地震裂缝轮廓虽已被完整详尽的提取，但提取结果中包含大量背景干扰信息(如噪声点、纹沟、河床以及阴影等)，需根据地震裂缝的空间特征实现精确提取。经高分影像解译，昆仑山口西地震裂缝主要由走向100°~110°的剪切裂缝、走向80°~90°的张剪裂缝以及走向45°~60°的张裂缝构成(图 3)。不同性质地震裂缝几何形状均表现为直线线性特征，空间组合呈平行排列、雁列状斜列分布特征。因此，本文引入面积、长宽比和方向3个形状特征进行描述。

面积为对象区域内所有的像元个数之和乘以单个像元对应的地面面积，计算公式为

$ S=\sum\limits_{k=1}^n A\left(x_{k j}-x_{k i}\right). $

(1)

式中: A为单个像元对应的地面面积，x_kj和x_ki为第k行上对象边界上的x坐标，且x_kj大于x_ki；n为该对象的总行数。

长宽比为对象长轴与短轴的比值，计算公式为

$ r=\frac{l_1}{l_2}. $

(2)

式中: l₁、l₂分别为对象长轴和短轴，该参数反映对象的延伸性，可进一步将裂缝与其他对象区分开来。

方向定义为对象长轴与水平方向的夹角，取值范围为0°~180°。

不同类型裂缝的描述和阈值范围如表 2所示。地震裂缝表现为连续延伸，具有一定面积的8连通区域，因此可将小于面积阈值的区域判定为噪点去除。剪切裂缝平直延伸，其长宽比最大；张裂缝波折状断续延伸，长宽比最小，张剪裂缝受剪切和拉张共同作用；其长宽比介于上述二者之间。受构造运动控制，地震裂缝延伸具有明显的方向性特征，可以有效与纹沟、河床等线性特征进行区分，河床虽可能与地震裂缝方向一致，但其长度应大于地震裂缝，且具有地震裂缝所不具备的延伸性，可通过设定地震裂缝长宽比这一空间特征的阈值去除其干扰，提取结果如图 4所示。

在本研究中，CLCHE处理、自适应阈值分割、形态学处理均基于Python语言开发环境实现，基于空间特征的误分类剔除通过ImageJ软件实现^[27]。

3 地震裂缝提取精度评价

本研究的提取结果体现了地震裂缝的空间分布与延伸特征。受制于地理位置和恶劣自然条件，地震裂缝精度评价很难开展野外实地验证工作。为对提取精度进行定量评估，将专家目视解译图作为真值与提取结果进行比较，并与改进的边缘检测^[28]以及MF-FDOG算法^[16]对比分析。根据检测和参考图中是否存在裂缝，将每个像素分配为裂缝或非裂缝，从提取结果和专家目视解译图两组数据中计算真阳性率(true positive rate，R_TP)和假阳性率(false positive rate，R_FP)。

$ R_{\mathrm{TP}}=\frac{T_{\mathrm{P}}}{T_{\mathrm{P}}+F_{\mathrm{N}}}, $

(3)

$ R_{\mathrm{FP}}=\frac{F_{\mathrm{P}}}{F_{\mathrm{P}}+T_{\mathrm{N}}}. $

(4)

式中: T_P-检测为裂缝，实际为裂缝; F_N-检测为非裂缝，实际为裂缝; F_P-检测为裂缝，实际为非裂缝; T_N-检测为非裂缝，实际为非裂缝。

各方法提取结果如图 5所示，提取精度评价如表 3所示，由于LoG-Canny边缘检测算法检测出的是裂缝两侧边缘，得到的结果中含有大量孔洞，采用闭运算时，较大尺寸的滤波核虽能将其完全填充，但会造成相邻地震裂缝堆叠。同样，该算法无法区别地裂缝与其他线性特征，导致部分河床的线性结构与地裂缝相连，从而造成R_FP较大。MF-FDOG算法则因为步骤复杂，且需要设置参数阈值较多，对于不同分辨率、不同类型的地震裂缝提取缺乏适应性，提取结果未能反映地震裂缝的延伸性特征，从而造成R_TP较小。而本研究方法提取精度最优，R_TP最高，R_FP最低，较好地反映了地震裂缝的形状特征。

4 地震裂缝空间分布特征及组合形式

查明地震裂缝的空间分布特征及组合形式，有助于掌握分析震源断层的最新几何学、运动学和动力学特征，同时也能为理解里德尔剪切模型提供良好示例。

4.1 分布特征

几何上，里德尔剪切破裂^[29]包括R和共轭剪切破裂R′、P破裂、Y破裂和共轭剪切X破裂、局部张性破裂(T)和局部挤压结构(图 6)。在走滑型地震中，当地表岩性松散、强度较低时，如洪积扇(台地)、湖积台地，震源断层错动常常发育纯剪切裂缝(Y破裂)；当位于山前陡倾位置，受拉张力和剪切力共同作用，往往发育张剪裂缝(R破裂)，其小角度相交于主剪切带(交角一般小于15°)，走滑方向与主剪切方向一致；张裂缝(T破裂)与主剪切带大致呈45°交角，开口较大，平面上主要呈楔形或平行四边形，呈雁列状排列(图 3)。在调查同震地震裂缝的力学性质、空间分布及延伸特征之后，将裂缝几何形状与里德尔剪切结构相关联，可确定震源断层的主剪切带方向和区域应力场^[30]，有助于进行发震构造动力学、几何学等相关研究。

4.2 组合形式

经遥感解译和实地调查，地表破裂带虽在宏观上表现为一条直线，但往往由地震裂缝与其他破裂单元组合而成，地震裂缝呈羽列状排列，右阶阶区发育表征地表缩短构造，左阶阶区发育表征地表伸展构造。

4.2.1 地震裂缝与伸展构造

在地震破裂过程中，受两条左行左阶破裂相对运动，地表松散堆积或风化残坡积层上因局部拉张而发育伸展构造，如张裂缝、地震凹陷和拉分盆地(图 7)。松散覆盖层发育的拉分盆地从构造上与黏土物质中的雁列状张裂缝类似，其规模与剪切位移呈正相关的关系^[31-32]。演化模式如图 7(d)，随着变形的继续，雁列状斜列的张裂缝发生叠覆，单个张裂缝或小型拉分盆地逐渐变宽、变深，最后联合起来形成规模更大的盆地。图 7(c)中拉分盆地规模及形态发育程度高，可能对应多次地震事件累计而成。

4.2.2 地震裂缝与缩短构造

在地震破裂过程中，受两条左行右阶破裂相对运动，在地表松散堆积或风化残坡积层上因局部压扭形成缩短构造，如地震鼓包、鼓梁(图 8)。前者为宽度略小于长度的小幅度隆起, 一般在平面上呈椭圆形状，后者隆起长度远大于宽度，平面上呈长条状。地震鼓梁长轴方向一般在N45 °~ 70 °W，在地表破裂带各段落均有不同程度发育。张(剪)性裂缝段与压扭性鼓包段相间排列，是整个破裂带分布最为普遍的破裂形态组合^[33]。相邻鼓包呈左阶斜列展布, 其间以右阶斜列展布的张(剪)性裂缝首尾相连。前述图 3(a)中张剪裂缝呈波折状不连续延伸，其阶区位置实质上也为地震鼓包(梁)相连。

地震裂缝的空间分布特征及组合形式反映地表破裂的几何学、运动学特征，揭示了构造活动下地震裂缝的动态发育过程。需要指出，地震裂缝的空间分布对于其他伸展、缩短型地表破裂的圈定具有重要的指示意义，后者虽然形态不一，无法通过遥感手段直接提取，但可根据地震裂缝阶区间接圈划出其可能发育位置。

5 结论与展望

震源断层错动形成的地震裂缝与发震断裂带位置基本重合，除光谱特征外，空间分布也具有明显的方向性特征。本文选取昆仑山口西地震同震地表破裂带为研究区，建立了一套基于光谱、空间特征的地震裂缝提取方法。对比目前常用地裂缝提取方法，结果表明该方法结果最优，MF-FDOG算法次之，改进的边缘检测算法效果较差。该方法适用于重点区域地震裂缝的精细刻画描述，本文选取试验区各类型地震裂缝形态完整，具备其典型的光谱特征和空间特征，在地表破裂提取方法研究和应用中能起到一定的参考价值。如何在更为复杂的背景中提取地震裂缝，将是今后的研究重点。

运用里德尔剪切理论在遥感影像中识别地震裂缝，分析走滑型地震同震地表破裂的运动特征与空间展布关系。将裂缝几何形状及延伸方向与里德尔剪切结构相关联，能快速确定主剪切带方向和区域应力场，有助于进行发震构造动力学、几何学等相关研究。同时指出，地震裂缝的空间分布对于其他伸展、缩短型地表破裂的圈定具有重要的指示意义。

本研究引入国产高分卫星GF-2作为影像数据源，显示了国产高分辨率遥感数据在地表破裂考察等地震科学研究方面的应用价值。提取结果证明，GF-2卫星完全满足大地震后地表破裂考察的需求，可为灾后应急与灾情评估提供有效的空间信息依据。