陆探一号（Lutan⁃1, LT⁃1）卫星是中国首个民用L波段合成孔径雷达（synthetic aperture radar, SAR）卫星星座,由A/B两颗卫星组成,分别于2022⁃01⁃26和2022⁃02⁃27成功发射,波长约23 cm,理论上单星重访周期为8 d,双星达到4 d^[1]。目前,LT⁃1卫星已正式投入使用,在地表形变监测领域开始发挥重要作用^[2-9]。因此,研究LT⁃1数据处理方法与策略,充分发挥其应用潜力,已成为合成孔径雷达干涉测量（interferometric synthetic aperture radar,InSAR）领域的重要研究内容。差分干涉测量技术需要进行高效的SAR影像配准运算,通过初配准获取可靠的初始偏移量是这一过程的关键步骤。目前,初配准方法主要包括目视法、基于轨道估算的几何配准法以及特征点法^[10-11]。目视法通过人工识别参考影像与待配准影像中的同名点,但受限于SAR影像纹理模糊、几何畸变和噪声明显,操作困难且效率不高,无法实现自动化批量处理^[12]。基于轨道的几何配准法在Sentinel⁃1、TerraSAR⁃X和ALOS⁃2等常用SAR卫星的数据处理中广泛应用^[13],但其精度极大受限于卫星轨道参数的准确性^[14-15]。一般的特征点法在处理光学影像时效果显著,但无法抵御SAR影像乘性斑点噪声的影响,不能直接适用于SAR影像^[16]。而目前LT⁃1卫星提供的是实时轨道数据,而非精密轨道数据,因此,采用几何配准法获得的初始配准偏移量往往误差较大,容易导致配准失败,进而影响数据处理的整体可靠性和效率。针对这一问题,本文研究基于特征点法的类尺度不变特征变换（scale invariant feature transform like , SIFT‐Like）的SAR影像初配准方法,旨在克服对精密轨道的依赖,提高LT‐1数据处理的效率,并为该数据的批量处理和大范围应用提供有效支持。

北京时间2023⁃12⁃18T23：59：30,甘肃省临夏回族自治州积石山县发生Ms 6.2地震,震中位于35.70°N,102.79°E,震源深度10 km^[17],地震最大烈度达到Ⅷ度。截至2023⁃12⁃31凌晨1点16分,共造成甘肃和青海两省151人遇难、979人受伤,导致大量房屋、电力、通信、交通等基础设施损毁^[18],同时诱发滑坡、崩塌、地裂缝、砂土液化等地质灾害^[19]。LT‐1 SAR卫星凭借其重返周期短和波长较长的优势,适用于此次地震的同震形变及其导致的地质灾害监测,可以为地震的应急响应提供支持。此外,该地震事件也为验证本文提出的LT‐1卫星初配准方法和策略的有效性以及应用提供了机会。

基于此,本文采用SIFT‐Like的初配准算法和处理策略,自动获取主、从SAR影像以及主SAR影像、数字高程模型（digital elevation model,DEM）模拟SAR影像之间特征点的初始偏移量,从而解决LT‐1卫星因轨道数据不精确而导致主、从SAR影像初配准以及去除地形相位时主SAR影像与DEM模拟SAR影像初配准失败的问题。将该算法应用于积石山地震震区的LT‐1卫星数据处理,获得同震形变场并进行分析,通过与Sentinel‐1结果比较验证其可靠性。此外,基于震前Sentinel‐1时序形变结果,利用滑坡自动识别算法识别震中附近的震前滑坡,并分析地震触发的草滩村液化滑坡⁃泥流事件的震前形变信号,讨论了LT‐1影像分幅不规则对算法有效性的影响,并展望了LT‐1卫星在未来地震和同震地质灾害监测中的应用潜力。

1 基于SIFT-Like的LT-1差分数据处理

本文基于SIFT⁃Like算法获取主、从SAR影像以及DEM模拟SAR影像、主SAR影像特征点的初始偏移量,来解决轨道不精确导致的配准初始偏移量偏差过大的问题,从而完成主、从SAR影像初配准以及DEM模拟SAR影像与主SAR影像初配准,并将该配准方法融入传统差分数据处理流程中,实现LT⁃1卫星差分干涉影像对的批量处理。

1.1   基于特征点获取初始偏移量

1.1.1   特征点提取

获取主、从影像特征点是计算两者初始偏移量的关键,最常用的方法为1999年Lowe^[20-21]提出的尺度不变特征变换算法（scale invariant feature transform, SIFT）,但该方法主要适用于光学影像,在SAR影像上效果不佳。Dellinger等^[22]对SIFT算法进行改进,用比值法代替差值法计算像素点的梯度,并改进了特征描述向量的生成方法,提出了SIFT‐Like算法,该方法一定程度抑制了SAR影像的斑点噪声,因此具有更好的适用性。具体过程如图1所示,主要步骤如下：

图  1  基于SIFT⁃Like算法获取特征点初始偏移量的LT⁃1差分数据处理流程

Figure  1.  Differential Data Processing Flow of LT⁃1 Based on Initial Offset of Feature PointsObtained by SIFT⁃Like Algorithm

下载: 全尺寸图片 幻灯片

1）特征点探测。利用比值梯度法（gradient of ration, GR）构建一种新的多尺度SAR‐Harris矩阵和多尺度SAR‐Harris函数探测器。通过比较目标点在当前尺度、前一尺度和后一尺度共23个邻域点的极值,得到局部极大值点作为特征候选点,然后设定SAR‐Harris函数阈值,剔除低于阈值的边缘点和低对比度点。

比值梯度法把以某点$(a,b)$为中心的正方形窗口一分为二,分别对两部分进行指数加权局部均值计算,然后计算两部分和值的比值得到该方向的梯度值。指数加权局部均值计算公式如下：

式中,$I(x,y)$为像素的灰度值；${\mathrm{e}}^{-\frac{|x|+|y|}{\alpha }}$为指数加权函数；$\alpha$为指数权重参数和尺度因子；$a\mathrm{、}b$分别为中心点的横、纵坐标；$x\mathrm{、}y$分别为距离向和方位向；$i$为1表示沿垂直方向,$i$为3表示沿水平方向。

指数加权比值为：

沿水平方向和垂直方向的梯度为：

2）特征点的特征描述向量生成。首先,计算特征点的主方向。以特征点为圆心形成一个半径为$\mathrm{6}\alpha$（$\alpha$为尺度因子）的圆形邻域,并计算邻域内各像素高斯加权的梯度幅值和方向。建立直方图,峰值所对应的方向即为主方向,大于峰值80%的方向作为辅方向。其次,计算特征点的描述向量。将特征点的邻域图像旋转至主方向在$X$轴上,再将邻域范围划分为$d\times d(d=\mathrm{4})$个子区域,统计每个子区域内像素的梯度方向分布,将梯度方向分为8个范围,此时每个子区域有8个方向的梯度幅值,将这些梯度幅值按照规律排列并归一化,得到$\mathrm{128}(\mathrm{4}\times \mathrm{4}\times \mathrm{8})$维特征描述向量。

3）特征点匹配。两景不同影像的特征点根据它们各自的特征描述向量进行匹配,首先,计算从影像$I_{\mathrm{1}}$的特征描述向量$X_{\mathrm{1}}$与主影像的最近邻特征描述向量$X_{\mathrm{2}}^{\mathrm{\text{'}}}$和次近邻特征描述向量$X_{\mathrm{2}}^{\mathrm{\text{'}}\mathrm{\text{'}}}$的欧氏距离的比值,若该比值小于阈值$t$,则待匹配影像$I_{\mathrm{1}}$的特征描述向量$X_{\mathrm{1}}$与参考影像的最近邻特征描述向量$X_{\mathrm{2}}^{\mathrm{\text{'}}}$为匹配特征点。

1.1.2   特征点初始偏移量估算

获取主、从影像的特征点后,首先计算每个匹配的特征点在距离向的偏移量$∆X_{\mathrm{1}},∆X_{\mathrm{2}},\cdots ,∆X_N$和方位向的偏移量$∆Y_{\mathrm{1}},∆Y_{\mathrm{2}},\cdots ,∆Y_N$（$N$为所匹配特征点的个数）,然后计算偏移量的均值并取整,作为主、从影像的初始偏移量$∆X\mathrm{、}∆Y$：

1.2   LT-1差分数据处理流程

相比传统差分干涉数据处理流程^[22],本文对LT-1卫星数据的差分干涉处理包含以下3个部分：（1）数据预处理。针对LT‐1卫星轨道数据不精确,采用SIFT‐Like算法获取主、从SAR影像特征点的初始偏移量,完成主、从影像的初配准,并进行前置滤波。（2）干涉图的生成。将主、从影像进行精配准,并裁剪至研究区范围,然后进行干涉处理得到原始干涉图,再去除平地相位得到去平后的干涉图。为解决轨道数据不精确导致的DEM模拟SAR影像与主SAR影像配准问题,再次应用SIFT‐Like算法获取DEM模拟SAR影像与主SAR影像特征点的初始偏移量,完成配准,从而去除地形相位。（3）形变结果生成。采用Goldstein滤波提高干涉图信噪比,利用最小费用流（minimum cost flow, MCF）算法进行相位解缠处理,将相位转换为地表形变,并通过多项式拟合去除轨道误差,建模去除地形相关大气影响,最终通过地理编码将形变结果转换至WGS84坐标系下的视线向（line of sight, LOS）形变。具体流程如图1所示。

2 地震形变监测

2.1   研究区概况及数据集

2.1.1   研究区地质构造背景

本次积石山Ms 6.2地震发生于青藏高原东北缘的甘东南活动构造区,震中位于积石山县柳沟乡。该区域因欧亚板块与印度板块的持续碰撞,形成了包括东昆仑断裂带、西秦岭北缘断裂带及拉脊山逆冲断裂带在内的多条走滑和逆冲断裂带^[23-25]。此次地震的发震断层为拉脊山北缘逆冲断裂的东段。自1949年以来,距震中200 km范围内共记录了14次5.0~5.9级和2次6.0~6.9级地震,均属于孤立型或主震⁃余震型地震^[26]。其中,2013年7月岷县发生的Ms 6.2地震距本次震中185 km,是该区域内震级最大且最强烈的一次地震^[17]。

2.1.2   研究区数据集

本文收集了积石山地震发生前、后的LT‐1和Sentinel‐1卫星数据。LT‐1数据包括震前2023‐12⁃18和震后2023⁃12⁃22的升轨影像,成像模式为条带模式1（Strip 1）,用于获取同震形变场。此外,还包括6个轨道覆盖的2023⁃07⁃19—2024⁃03⁃24的22景影像,用于验证基于SIFT⁃Like的SAR影像初配准方法的可靠性和自动化处理能力。Sentinel‐1数据包括震前2023⁃10⁃27和震后2023⁃12⁃26的升轨影像,以及震前2023⁃12⁃14和震后2023⁃12⁃26的降轨影像,也用于获取同震形变场并与LT‐1形变结果进行对比。由于震后的LT‐1和Sentinel‐1数据有限,无法达到时序数据处理的要求,因此,仅收集震前一段时间的Sentinel‐1升轨和降轨数据,用于获取震前的时间序列形变。影像覆盖见图2,具体参数见表1。从表1中可以看出,LT‐1卫星数据在震前和震后数据间隔只有4 d,能够很好地满足地震的应急响应要求,而Sentinel‐1卫星升、降轨数据在地震前后的数据间隔分别为61 d和12 d,没有LT‐1卫星响应及时。

图  2  LT‐1和Sentinel‐1卫星数据覆盖图

Figure  2.  Coverage Map of LT-1 and Sentinel-1 Satellite Data

下载: 全尺寸图片 幻灯片

表  1  LT-1和Sentinel-1卫星数据参数

Table  1.  Parameters of LT-1 and Sentinel-1 Satellite Data

卫星	轨道方向	主影像时间/影像起始时间	从影像时间/影像终止时间	时间基线/d	空间基线/m	数量/个
LT⁃1（同震）	升轨	2023⁃12⁃18	2023⁃12⁃22	4	743.13	4
LT⁃1(算法实验)	升轨	2023⁃07⁃12	2024⁃03⁃24			22
Sentinel⁃1（同震）	升轨(T128)	2023⁃10⁃27	2023⁃12⁃26	60	64.28	2
	降轨(T135)	2023⁃12⁃14	2023⁃12⁃26	12	-116.22	2
Sentinel⁃1（震前）	升轨(T128)	2021⁃01⁃10	2023⁃06⁃05			67
	降轨(T135)	2021⁃01⁃05	2023⁃12⁃14			51

下载: 导出CSV 

| 显示表格

2.2   基于SIFT-Like的LT-1特征点初始偏移量获取实验

LT‐1卫星影像初配准的初始偏移量可通过目视法、基于轨道估算的几何配准法以及基于特征点的SIFT‐Like算法获取。目视法因其满足初配准所需的几十个像素的精度标准,可作为其他两种方法初配准结果的验证手段。比较几何配准法、SIFT‐Like算法分别与目视法在距离向和方位向的初始偏移量差值,若差值在几十个像素内,则初配准成功,若超过几百个像素,则初配准失败。进一步通过精配准后的两景SAR影像干涉图进行判断,若生成清晰可靠的干涉条纹,则精配准成功。分别利用上述3种方法对4景积石山地震前后的LT‐1升轨数据进行实验,获得的初始偏移量数据如表2所示。可以看出,在主、从SAR影像初配准阶段,几何配准法、SIFT‐Like配

表  2  LT-1卫星数据初始偏移量获取方法比较/像素

Table  2.  Comparison of LT-1 Satellite Data Initial Offset Acquisition Methods/Pixel

类别	几何配准法获取初始偏移量		目视法确定特征点初始偏移量		SIFT-Like自动确定特征点初始偏移量
	距离向	方位向	距离向	方位向	距离向	方位向
主、从SAR影像初配准	1 138	8 321	1 110	8 300	1 140	8 290
去除地形相位时初配准	0	0	-6	400	5	399

下载: 导出CSV 

| 显示表格

准法与目视配准法在距离向和方位向上的初始偏移量的差值均在几十个像素内,均满足初配准的精度要求,从而成功进入精配准阶段。将配准后的SAR影像进行干涉处理得到干涉图（如图3所示）,从整体干涉图（图3（a）~3（c））可以看出,3种方法的结果均产生了清晰的干涉条纹。从3个整体配准图中分别选取3处区域的细部配准图（图3（d）~3（f））,可以更加清楚地看到清晰的干涉条纹,表明均配准成功。然而,在去除地形相位时主SAR影像与DEM模拟SAR影像的初配准中,受限于轨道精度的几何配准法与目视法在方位向上的初始偏移量差值超过几百像素,无法有效配准和去除地形相位。在配准成功率方面,目视法和SIFT‐Like算法的成功率高于几何配准法。在配准结果可靠性方面,目视法确定1个特征点且受限于LT‐1 SAR影像的质量和目视解译人员的经验,几何配准法以参考影像的中心点作为特征点,而SIFT‐Like算法分别确定了51个（图4）和7个特征点（图5）。因此,SIFT‐Like算法的结果相对更为可靠。在配准自动化程度方面,目视法寻找特征点效率低下且无法自动化,几何配准法和SIFT‐Like算法均能实现批量自动化处理。综上所述,SIFT‐Like算法为LT‐1卫星影像提供了一种高成功率、高可靠性且高自动化程度的初配准方法。

图  3  配准后SAR影像的整体干涉图及细部干涉图

Figure  3.  Overall Interferograms and Detailed Interferograms of Registered SAR Images

下载: 全尺寸图片 幻灯片

图  4  初配准时主、从影像特征点匹配

Figure  4.  Feature Points Matching of Initial Registration for Master and Slave Images

下载: 全尺寸图片 幻灯片

图  5  去除地形相位时主SAR影像与DEM模拟SAR影像特征点匹配

Figure  5.  Matching of Feature Points Between Master SAR Image and SAR Image Simulated by DEMDuring Terrain Phase Removal

下载: 全尺寸图片 幻灯片

为进一步评估SIFT‐Like算法的性能,本研究还利用收集的震中22景LT‐1卫星升轨影像数据（影像覆盖如图2（b）所示）,形成了29个影像对,并应用上述3种方法进行初配准实验。主、从SAR影像初配准以及去除地形相位时主SAR影像、DEM模拟SAR影像初配准时,几何配准法、SIFT‐Like算法与目视法初始偏移量的差值结果分别如图6（a）和图6（b）所示。可以看出,在主、从SAR影像初配准阶段,几何配准法仅对10个影像对实现了初配准（差值在几十个像素内）,最终也都实现了精配准,达到了34.5%的配准成功率。而SIFT‐Like算法则成功完成了所有29个影像对的初配准,最终也都完成了精配准,达到了100%的配准成功率。在去除地形相位时主SAR影像与DEM模拟SAR影像的初配准过程中,SIFT‐Like算法则再次展现了其优异的性能,29个影像对均完成了初配准,而几何配准法仅有7个影像对实现了初配准。最终SIFT‐Lik算法维持了100%的配准成功率,而几何配准法的成功率仅为24%。在配准成功率方面,目视法和SIFT‐Like算法均达到100%,而依赖于轨道精度的几何配准法表现较差；在配准可靠性方面,SIFT‐Like算法

图  6  几何配准法、SIFT⁃Like算法与目视法初始偏移量差值

Figure  6.  Difference in Initial Offset Between Geometric Registration Method, SIFT-Like Algorithm,and Visual Inspection Method

下载: 全尺寸图片 幻灯片

通过确定多个特征点表现出更高的可靠性；在配准自动化程度方面,几何配准法和SIFT‐Like算法均能实现自动化批量处理,而目视法效率低且不能实现自动化处理。综上所述,SIFT‐Like算法在实时轨道不精确的LT‐1卫星SAR影像初配准中表现出色,相比其他两种方法更具优势。

2.3   LT-1和Sentinel-1差分干涉数据处理

利用GAMMA软件处理LT‐1和Sentinel‐1卫星数据^[27-28],针对LT⁃1数据,在§2.2中基于SIFT‐Like算法获取了主、从SAR影像及DEM模拟SAR影像与主SAR影像特征点的初始偏移量后,完成主、从SAR影像初配准以及DEM模拟SAR影像与主SAR影像的初配准,再按照§1.2中的差分数据处理流程得到了积石山地震的LT‐1卫星同震形变场,结果如图7所示。

图  7  积石山地震LT-1卫星升轨同震形变场

Figure  7.  Coseismic Deformation Field of the Jishishan Earthquake from LT-1 Satellite Ascending Orbit

下载: 全尺寸图片 幻灯片

针对Sentinel⁃1卫星数据,除无需使用SIFT⁃Like算法来纠正轨道不精确导致的配准问题外,其余处理流程与LT⁃1类似,最终得到了Sentinel⁃1的同震形变场,如图8所示。

图  8  积石山地震Sentinel⁃1卫星数据的同震形变场

Figure  8.  Coseismic Deformation Field of the Jishishan Earthquake Using Sentinel-1 Satellite Data

下载: 全尺寸图片 幻灯片

2.4   震前Sentinel-1时序数据处理

利用GAMMA软件对震前67景Sentinel-1升轨数据和51景降轨数据进行小基线子集（small baseline subset, SBAS） InSAR处理^[29]。处理流程包括（1）配准与差分干涉：首先将所有影像与主影像进行精确配准,并设置时空基线阈值来选择合适的干涉对,对每个干涉对进行差分干涉处理,然后进行滤波、相位解缠和轨道误差去除处理；（2）奇异值分解法（singular value decomposition, SVD）及参数估计：在评估干涉对质量后,剔除质量较差的干涉对,再利用SVD法计算出初始形变序列和地形残差；（3）形变序列求解：基于初始形变序列进行时间序列的线性回归分析,得到线性形变,从初始形变序列中去除线性形变,得到残余相位。进一步对残余相位进行时空滤波处理和大气相位改正,得到非线性形变；将线性形变和非线性形变相加,得到完整的形变序列。升、降轨时序形变结果如图9所示。

图  9  积石山地震震前Sentinel-1卫星数据的时序形变

Figure  9.  Time Series Deformation of Sentinel-1 Satellite Data Before the Jishishan Earthquake

下载: 全尺寸图片 幻灯片

3 地震形变分析

3.1   同震形变分析

从LT⁃1卫星升轨数据获得的积石山地震LOS向同震形变场（图7）可以看出,形变场大致呈椭圆形,长轴方向为北西-南东走向,与多家机构公布的震源机制走向相符^[24]。震中形变场位于拉脊山北缘断裂和南缘断裂带之间,说明此次地震可能与该断裂有关。以地表形变靠近卫星飞行方向为正（抬升）,远离卫星飞行方向为负（沉降）,可以看出震中附近的同震形变场以抬升为主,基本没有沉降区域,最大抬升量约6.2 cm。从Sentinel‐1卫星升、降轨LOS向同震形变场（图8）可以看出,升、降轨最大抬升量分别为5.3 cm和6.3 cm。震中附近的形变场符号相同,说明同震形变为朝向卫星飞行方向的抬升,且以垂直方向为主,这与逆冲型地震特征相符,与各机构公布的震源机制解一致^[24]。

为验证LT‐1卫星同震形变场的准确性,通过比较图7中LT‐1升轨数据和图8中Sentinel‐1升、降轨数据所得的同震形变场,发现三者形变结果基本一致。在形变场图中选取一块基本无

形变的区域（如图7、图8（a）中棕色矩形框所示）,做形变量直方图,如图10所示,可以看出3种数据的形变量均集中在0 mm附近,且均值接近0 mm。从标准差可以看出,L波段的LT‐1数据精度水平略低于C波段的Sentinel‐1数据,这主要归因于卫星设计波长的差异。

图  10  LT⁃1升轨、Sentinel⁃1升、降轨形变量统计直方图

Figure  10.  Statistical Histogram of Displacement for LT-1 Ascending Orbit, Sentinel-1 Ascendingand Descending Orbits

下载: 全尺寸图片 幻灯片

为进一步比较LT‐1与Sentinel‐1卫星同震形变场的一致性,本文在两种升轨同震形变场中选取了1条跨拉脊山南、北缘断裂的剖线AB,并绘制其剖线图,如图11所示。结果显示,从A到B,两种卫星数据均呈现由缓慢形变到显著抬升再到沉降的趋势,有明显的抬升中心,与逆冲型地震特征相符。两种形变趋势较为相近,无大范围形变位移跳变,表明此次地震的活动断层未破裂至地表,同震滑动未传播到地表。虽然LT‐1卫星数据获得的形变场中存在一定的空值,形变曲线存在一定的波动,而Sentinel‐1的形变结果较为平滑,但两者总体形变趋势一致,表明两种数据的形变结果具有一定的一致性。

图  11  AB剖线LT⁃1升轨、Sentinel⁃1升轨同震形变曲线图

Figure  11.  Coseismic Deformation Curves of Ascending Orbit and Sentinel-1 Ascending and Descending Orbits for AB Profile Line

下载: 全尺寸图片 幻灯片

3.2   震前时序形变分析

研究区内Sentinel‐1升轨和降轨震前的时序形变速率结果（图9）显示,研究区正北部存在一个形变速率显著的区域,研究区西北部零散分布着多个形变速率较大的区域,既有沉降也有抬升。同时,研究区南部也分布着多处形变速率较大的区域,以沉降为主。鉴于地震可能引发同震地质

灾害,尤其是滑坡灾害,但目前震后的LT‐1和Sentinel‐1卫星数据相对有限,不足以获得震后的时序形变信息,因此,本研究充分利用震前的时序形变信息探测震前的潜在滑坡区域。基于图9的震前Sentinel‐1升、降轨时序地表平均形变速率,采用一种基于时序地表平均形变速率的滑坡自动识别算法^[30],在积石山地震震中研究区内自动识别潜在滑坡灾害。升轨数据识别出195个潜在滑坡区域（图9（a））,总面积达83.19 km²,降轨数据识别出179个潜在滑坡区域（图9（b））,总面积达76.95 km²。两种数据识别的滑坡区域有93处重叠,具体滑坡信息如表3所示。由于升轨和降轨数据的观测几何特性不同,识别的潜在滑坡区域存在一定的差异。

表  3  识别的震前滑坡信息表

Table  3.  Identified Pre-Seismic Landslide Information

轨道方向	数量/个	总面积/km2	最小面积/km2	最大面积/km2	重叠个数
升轨	195	83.19	0.085	6.28	93
降轨	176	76.95	0.085	8.78

下载: 导出CSV 

| 显示表格

4 讨论

4.1   LT-1卫星数据配准方法

针对LT‐1卫星未提供精密轨道数据导致的配准失败问题,本文提出一种基于SIFT‐Like的

特征点初始偏移量自动获取方法,用于LT‐1卫星主、从SAR影像的初配准以及DEM模拟SAR影像与主SAR影像间的初配准。为比较3种方法的性能,本研究进行了多组实验,从配准成功率、可靠性、自动化程度3个方面,对比了基于实时轨道的几何配准法、目视法以及SIFT‐Like算法。实验结果显示,在配准成功率方面,几何配准法受限于轨道精度,两类初配准的成功率分别为34.5%和24%,而目视法和SIFT‐Like算法均达到100%,显著提升了配准成功率。在配准可靠性方面,目视法通常每次只确定1组特征点,且受限于LT‐1 SAR影像的质量和目视解译人员的经验,虽然可行但效率低且易出错；几何配准法每次也仅选择参考影像的一个中心点作为特征点；而SIFT‐Like算法通过确定多个特征点,展现了更高的可靠性。在配准的自动化程度方面,目视法无法实现数据的批量自动化处理,通常用于结果验证；而几何配准法和SIFT‐Like算法均能实现配准的批量化和自动化,显著提升工作效率。综上所述,SIFT‐Like算法相比其他两种方法能更有效地解决LT‐1卫星配准问题。

随着LT‐1卫星数据的持续积累,海量影像需要自动化配准处理。通过本研究的SIFT‐Like算法可以解决配准问题,但在实际应用中仍会存在失败的情况。主要由于LT‐1卫星数据分幅的不规则性,导致部分同一轨道在不同时间获取的同一覆盖数据间存在较大偏差。为此,需要按照重叠率对同一轨道不同时间获取的同一覆盖数据进行分类。但由于目前LT‐1卫星的存档数据相对有限,数据间的重叠率不高,存在将不同覆盖数据错分为一类的情况。即使计算了特征点的初始偏移量,也可能因可用特征点数量过少而导致配准失败。随着LT‐1卫星数据的增多,数据间的重叠率将逐渐增大,有望在一定程度上缓解上述问题。

4.2   积石山地震触发滑坡分析

一般地震发生后会触发各种地质灾害,此次积石山地震触发了一处位于震中东北方向的青海省海东市民和县中川乡草滩村的大型液化滑坡-泥流灾害,造成20人遇难、95间房屋损毁以及19条道路阻塞,是本次地震中人员伤亡和财产损失最为严重的地质灾害^[31-33]。该滑坡⁃泥流分为滑源区和流通堆积区,滑源区位于黄河三级阶地台塬上,是地势平坦的耕地,流通堆积区位于黄河二级阶地,呈条带状分布^[33]。根据李亚军等^[34]的研究,地震引发的台塬面下饱水黄土层液化是该灾害发生的主要原因。根据李为乐等^[32]的研究,该液化滑坡‐泥流在滑源区正南部形成的积水池堵塞地下水通道,导致地下水位上升,下部黄土饱和,在地震的震动下,饱水黄土液化,导致此次灾害的发生。尽管震前（2021‐01—2023‐12）的形变速率结果（图12（a）和图12（b））显示该区域整体上较为稳定,且未探测到滑坡灾害,但滑源区出现了局部形变。通过选取滑源区的一个时序点（如图12（a）~图12（d）红色三角形所示）,从时序形变图（图12（e）、图12（f））可以看出,升、降轨影像的时序点在监测时间内持续沉降,特别在震前一段时间内出现加速沉降趋势,升、降轨年均沉降速率均超过9 mm/a,可能与农作物灌溉活动有关。震前农作物灌溉活动导致地下水位上升,地下黄土处于饱水状态,地震的震动使得饱水黄土液化,导致此次灾害的发生。

图  12  滑坡-泥流震前形变速率图、滑源区形变速率图以及时序形变图

Figure  12.  Landslide-Debris Flow Pre-Seismic Deformation Rate Map, Deformation Rate Map of the Slip Source Area, and Temporal Deformation Map

下载: 全尺寸图片 幻灯片

此外,基于Sentinel‐1卫星的升、降轨时序形变速率,成功探测出震中附近多处面积大、形变量显著的潜在滑坡灾害（如图9所示）。这些形变量级较大的区域在地震影响下极可能引发严重的次生灾害,需加强关注和监测。当前,一些研究采用光学遥感影像来识别积石山地震震后潜在的地质灾害^[31-32],但不可避免受到云、雾等因素的干扰而造成漏判或误判。尽管LT‐1和Sentiel‐1卫星的震后数据有限,无法提供震后的时序形变信息,但随着SAR数据的积累,未来可以更深入研究积石山地震触发的地质灾害。

4.3   LT-1卫星对地震防灾减灾的展望

LT‐1卫星因其重访周期短、波长长及穿透能力强,能较好地抵御时空失相干,可以有效地捕获地震引起的大梯度形变信号,因此广泛应用于地震同震形变监测和同震地灾探测。目前有研究利用LT‐1数据对泸定、花莲和积石山等中强地震进行同震形变监测,与同期其他SAR卫星和GNSS结果一致性较高^[9]。此外,还利用LT‐1的纹理特征提取震后滑坡,并评估对道路造成的损毁情况,为应急救灾提供路线规划,凸显其在地震应急救灾领域的巨大潜力。总的来说,目前LT‐1卫星除了可以提供地震的形变信息外,还能提供重点防震区地壳形变产品、震灾基础设施破坏评估产品和活动断层解译产品等^[2]。这些产品为地震防灾减灾提供有力支持,具有广泛应用前景。此外,地震是触发滑坡、崩塌、泥石流等同震地质灾害的重要因素,历史上地震造成的人员伤亡和财产损失中很大一部分是由于同震地质灾害导致的^[17],LT‐1数据在此类同震地灾监测中具有显著优势。目前地震、同震地灾监测主要以C波段中等分辨率的Sentinel‐1卫星数据为主,未来联合L波段高分辨率的LT‐1数据,可提高形变监测的时间分辨率,更及时、有效地响应复杂地区的地震灾害和同震地质灾害。

5 结语

本文针对LT‐1卫星实时轨道数据的不精确性,提出了一种基于SITF⁃Like算法的LT‐1卫星数据配准方法。该方法自动获取LT‐1卫星SAR影像特征点的初始偏移量,完成了主、从SAR影像初配准和去除地形相位时主SAR影像与DEM模拟SAR影像的初配准,成功获取了积石山地震的LT‐1升轨同震形变场。研究结果如下：

1）此次积石山地震以抬升的垂直形变为主,最大形变量达6.3 cm,属于逆冲型地震,且活动断裂未破裂至地表。LT‐1与Sentinel‐1形变结果有较好的一致性。

2）SIFT‐Like算法在主、从SAR影像初配准和去除地形相位时主SAR影像与DEM模拟SAR影像的初配准方面,比基于实时轨道的几何配准法和目视方法更准确、高效。

3）基于Sentinel‐1升、降轨时序数据,成功识别了多处震前滑坡。草滩村液化滑坡‐泥流在震前已有持续沉降,滑源区正南部的积水池堵塞地下水通道,以及震前农作物的灌溉活动,使得地下水位抬升,导致下部黄土饱和,最终在地震的作用下引发液化滑坡-泥流。