基于改进MT-InSAR的日兰高铁巨野煤田段沉降监测

祝传广，张继贤，邓喀中，龙四春，张立亚，吴文豪

(1.湖南科技大学 地球科学与空间信息工程学院，湖南 湘潭 411201；2.湖南科技大学 测绘遥感信息工程湖南省重点实验室，湖南 湘潭 411201；3.国家测绘产品质量检验测试中心，北京 100830；4.中国矿业大学 环境与测绘学院，江苏 徐州 221116)

日照—兰考高速铁路(日兰高铁)自2017年开工建设，预计2021年底全线开通运营。该铁路在菏泽市约有15 km路段穿越巨野煤田，地下煤矿开采引起的地表沉降对高铁安全运营构成巨大隐患。已有研究显示，高铁沿线在2015—2019年间存在较大的地表沉降现象。为保证铁路的安全运营，需要对高铁轨道、路基及沿线地表进行持续的沉降监测。

星载合成孔径雷达干涉测量(InSAR)技术因其高分辨率、广覆盖、全天候等优点成为常用的大地测量技术之一。其中，MT-InSAR技术借助特殊的相干目标能够以mm/a～cm/a的精度实现长时间跨度的形变监测。这些相干目标主要来自永久散射体(Persistent Scatterer,PS)和分布式散射体(Distributed Scatterer,DS)。PS集中分布在城区，以人工建/构筑物为主。在缺少人工建/构物的非城区，PS点过于稀少，难以详实反映形变的空间分布情况。并且，低密度的空间点目标会增加相位解缠难度从而降低形变估计精度。

DS广泛分布于自然界，利用DS能够大幅度提高观测点密度。DS的后向散射能力中等，能够在短时期内保持一定的相干性。MT-InSAR中的小基线(SBAS-InSAR，下文简称SBAS)技术就是利用这一特点，通过限制干涉对的时空基线长度以便尽可能地保持干涉对的相干性，从而将DS纳入MT-InSAR分析，达到增加观测点密度的目的。相比于SBAS仅仅利用短时空基线干涉对，以Squee SAR为代表的DS-InSAR技术则是利用全组合干涉对解算DS相位最优估值。采用全组合干涉对能够进一步提高SAR信息的利用率，然后通过联合DS和PS点集大幅度提高观测点密度。

Squee SAR采用相位三角算法(phase triangle algorithm,PTA)解算DS相位，理论上是真实值的最大似然估计，因而相位估计精度高。但是在低相干区域或者同质像素较少时，协方差矩阵通常含有较大误差，这会削弱PTA的估计精度。并且PTA算法需要迭代运算，较为耗时。在SqueeSAR的基础上又演化出几种变体，如CAESAR (Component extrAction and Selection SAR)，PD-PSInSAR (Phase-Decomposition-based PS InSAR)以及EMI (Eigendecomposition-Based Maximum-likelihood-estimator of Interferometric phase)等。这几种变体的原理与PTA算法在数学形式相互一致，只是加权策略不同。CAESAR、PD-PSIn SAR以及EMI的基本思想是通过特征分解(Eigenvalue Decomposition,EVD)技术解算DS相位：即对协方差矩阵(或加权后的协方差矩阵)进行特征分解，最大或最小特征值所对应的特征向量即为DS相位最优估值，其他特征向量则视为噪声而去除。因此，EVD算法中隐含了相位滤波处理。另外，EVD算法不需要迭代处理，因而计算速度快。相较于CAESAR和PD-PSInSAR，EMI算法利用相干矩阵调节各干涉对的权重，以期抑制低相干干涉对的影响，能够在一定程度上降低协方差矩阵估计误差的影响。然而，实际应用中由于SAR数据量可能较少、同质区域内相位可能不平稳、研究区域相干性低和/或同质像元少等原因，协方差矩阵和相干估值通常是有偏的。

针对上述问题，笔者在EMI的基础上提出了一种基于Fisher信息量的改进EMI算法(FEMI)。Fisher信息量是观测到的变量(如干涉相位)中所携带的关于未知参数(如DS的真实相位)信息量的一种度量，通常用于计算最大似然估计中的协方差矩阵。利用Fisher信息量代替相干矩阵作为权阵能够抑制相干值估计误差对相位估计精度的影响。

另外，标准SqueeSAR的处理流程是在识别出DS点集之后与PS点集合并，然后进行单主影像干涉(PSI)处理。考虑到研究区域的形变速度较快，笔者在合并PS和DS点集之后根据一定的时空基线阈值组合成小基线干涉序列，采用SBAS处理流程估计形变速率和时序形变。

首先利用模拟数据和真实SAR数据验证了FEMI算法的可靠性，然后基于2020年10月至2021年4月间的Sentinel-1卫星数据，采用SBAS处理框架获取了日兰高铁巨野煤田段的地表沉降信息，并结合2015—2019年的沉降监测资料分析了地表沉降的成因及时空演化机理。

1 DS相位估计方法

首先介绍DS相位优化估计的基本思想，然后分析CAESAR和EMI算法原理及其特点，进而引出Fisher信息量并构建基于Fisher信息量的改进EMI算法，即FEMI。

1.1 CAESAR和EMI算法

DS像元内没有强散射体，各基本散射单元的后向散射能力相似。因此，由中心极限定理可知DS的复反射率服从均值为零的复高斯分布。若有景SAR影像，则DS像元的时序复反射率=[,,…,]服从零均值多元复高斯分布。因为均值为零，所以的统计特征可用协方差矩阵描述。根据同质区域内的像元计算得到的样本协方差矩阵是真实协方差矩阵的最大似然估计。样本协方差矩阵的计算公式为

(1)

式中，为同质区域；为同质区域内的像元数目；上标H为共轭转置。

为了提高DS的相干性，通常需要对其干涉相位进行滤波处理(即空间自适应多视处理)。滤波后干涉相位的一致性被破坏，导致协方差矩阵秩亏，因而无法根据式(1)直接解算DS的真实相位。因为服从复Wishart分布，Ferretti等提出了通过最大化的概率来估计的PTA算法。PTA算法需要迭代运算，较为耗时。Fornaro等利用谱分解算法根据估计，即CAESAR算法：

(2)

(3)

式中，,=[||],；′,=[||||],；,为DS的干涉相位的观测值，,=-。

因为协方差矩阵需要根据同质区域内的样本计算得到，而同质区域是通过假设检验得到，同质区域内可能存在不服从高斯分布的像元，这会导致存在偏差。另外，在低相干区域仅采用最大特征值所对应的特征向量代表协方差矩阵会引起较大偏差，因为此时可能需要多个甚至全部的谱分量才能有效反映。为此，ANSARI等提出了EMI算法，通过给各干涉对按其相干质量赋予不同的权重以抑制低相干的影响，即：

(4)

1.2 基于Fisher信息量的EMI算法：FEMI

由样本估计的相干矩阵是有偏的，尤其是在低相干区域，因而不能保证是正定-半正定矩阵。为此，应用式(4)求解DS相位之前需要预先判断是否满足正定-半正定。若不满足，还需要通过添加阻尼因子等方式调整。

Fisher信息量是一种度量观测量(如干涉相位的观测值,)中包含待求参数(如干涉相位的真值,)信息多寡的指标。观测到的干涉相位,是由同质区域内的样本估计得到，估计误差的方差渐进逼近Cramer-Rao下界，而Fisher信息量能够反映估计误差的方差与Cramer-Rao下界的逼近程度。干涉相位,的Fisher信息量可表示为

=2(1-)

(5)

可以看出，为单调递增函数。当=0时，表示干涉相位没有有效信息，此时=0；当=1时，说明干涉相位不含噪声，此时→∞。

相比于相干值只能在直观上反映观测量的好坏，Fisher信息量在数值上给出了观测量中包含待估参数信息的多寡。并且，在低相干区，利用加权比更稳健。另外，利用Fisher信息量的迭代优化算法不易陷入局部极值。因此，本文利用Fisher信息矩阵代替，实现基于Fisher信息量的DS相位估计算法，即FEMI算法。由于同一DS的同质像素数目在不同的干涉对上保持不变，因此可令=1。利用Fisher信息矩阵代替，式(4)可改写为

(6)

式中，′,=[],。

2 实验验证

2.1 模拟数据实验

..数据模拟

模拟数据的相干值由下式计算：

=(-)e-Δ+

(7)

其中，Δ为干涉对的时间基线；为信号的时间相干长度，越大，干涉对的时间相干性越好；和为相干起始值和结束值,=0时表明相干值呈负指数快速衰减，干涉对只能在短期内保持相干性，时间基线Δ≈3时相干性近乎于0。

采用式(7)模拟2种相干类型数据：① 短期相干型，即=0；②长期相干型，即=02。2种类型的设置为0.6，设为50 d。因为后文采用16景真实Sentinel-1卫星影像开展形变监测研究，为此本节模拟16景SAR数据，并仿照Sentinel-1卫星将时间间隔设为12 d、波长设为55.6 mm。据式(7)计算相干阵(图1)，并由相干阵生成16景SAR数据，每景SAR数据包含300个同质像元。为每景SAR数据添加变形速率为5 mm/a的形变相位。

图1 模拟相干阵 Fig.1 Simulated coherence matrix

..实验结果分析

根据2.1.1节的2组模拟数据，分别采用EMI和FEMI算法估计DS相位，重复进行20万次，然后计算均方根误差():

(8)

图2 EMI与FEMI算法精度评估Fig.2 Performance assessment of EMI and FEMI methods using simulated data

(1)由图2(a)可知，对短期相干型数据，EMI的相位估计精度显著降低。因为短期相干型数据的时间失相干严重，干涉对的时间基线超过150 d时几乎没有相干性(图1(a))，即干涉相位中几乎全为噪声，不包含任何有效信息。而通过样本计算得到的相干值通常是有偏的，使得这些噪声传递到估计的DS相位中。而FEMI通过使用Fisher信息作为权阵，大幅度降低低相干干涉对的权重，抑制了误差的传播，因而FEMI算法较EMI精度更高也更加稳健。

(2)由图2(b)可知，对于长期相干型数据，FEMI和EMI算法的性能相近，并且接近CRLB。这是因为，对于长期相干型数据，时间基线最大的干涉对其相干性仍然大于0.2(图1(b))，所有的干涉对都能够提供有效信息。

综上可知，FEMI算法更加稳健。当研究区域的相干性较好时，FEMI和EMI算法性能相当。当研究区域的相干性较差时，则建议采用FEMI算法。虽然增大同质滤波的窗口理论上能够提高相干值的估计精度。但是，窗口过大时难以保证信号的空间平稳性。即使进行了地形相位补偿，大气效应和形变信息的空间变化也会导致信号空间失相关。

2.2 真实数据实验

..实验数据

选择Sentinel-1A卫星获取的SAR影像作为实验数据，波长约为5.6 cm，空间分辨率约为20 m × 5 m。Sentinel-1A卫星的重访周期为12 d，在6个月内(10月—次年4月)能够提供16景卫星影像。笔者收集了Sentinel-1A卫星在2020年10月—2021年4月间获取的日兰高铁巨野煤田段2组共32景SAR影像分别展开研究，其中1组(共16景影像)来自Path-40(下文简写为S1-40)，另一组(共16景影像)来自Path-142(下文简写为S1-142)。2组数据的部分参数见表1。另外，采用空间分辨率为30 m的AW3D30 DSM数据补偿地形相位。

表1 Sentinel-1卫星数据主要参数

..实验结果分析

通常，干涉对的时间基线越长相干性越差。检查干涉对的干涉相位图和相干系数图能够直观上评估算法的效果。对于本文数据而言，最大时间基线(180 d)对应的2组干涉对分别为：20201012—20210410(S1-40)和20201007—20210405(S1-142)。图3，4分别显示了这2组数据的干涉相位与相干系数的局部细节，为对比分析，同时显示了原始的以及EMI算法得到的干涉相位图和相干系数图。由图3，4可以看出，EMI和FEMI算法均能显著改善相位质量。在相干性高的区域，两种算法的效果近似一致。但是在低相干区域，FEMI算法的表现优于EMI。这与2.1节利用模拟数据得到的结果相一致。

图3 干涉对20201012—20210410的干涉相位和相干系数Fig.3 Interferograms and the coherence of 20201012-20210410

图4 干涉对20201007—20210405的干涉相位和相干系数Fig.4 Interferograms and the coherence of 20201007-20210405

3 联合PS和DS的SBAS处理框架

矿区形变通常较为快速，并且高铁沿线植被较多，时间失相干严重。为提高观测点密度、增加观测方程、保证每个干涉对的相干质量，采用联合PS和DS的SBAS干涉处理策略估计形变。假设所有的SAR影像均已配准至公共主影像：即S1-40数据配准至2021-01-04的影像参考几何、S1-142数据配准至2021-01-11的影像参考几何。联合PS和DS的SBAS处理框架如图5所示，具体步骤：

(1)提取PS点集与相位。采用单主影像策略，利用PSI技术获取PS点集及相位，PSI原理与流程见文献[32]。PSI分析过程中，振幅离差设为0.4。

(2)提取DS点集与相位。首先根据SAR幅度序列采用K-S检验算法识别DS初选点，显著性水平设置为0.05。窗口大小设置为9×35，同质像元(SHP)阈值设为25。然后，根据配准后的SAR数据和AW3D30 DSM生成全组合(即16×152=120)差分干涉序列。根据DS初选点和全组合差分干涉序列，根据式(1)计算协方差矩阵，然后利用本文提出的FEMI算法(式(6))估计DS相位。

(3)筛选DS点集。按照式(9)计算DS相位估值的拟合度，选取>0.7的像素作为最终DS点。

(9)

图5 数据处理流程Fig.5 Framework of data processing

式中，Re为取复数的实部

(4)合并PS和DS点集，并生成小基线干涉序列。在组合小基线对时，时间基线阈值设置为36 d，空间基线阈值设置为200 m。另外，需要确保没有非连通子集出现以便增加解的可靠性。最终，2组数据集均组合得到42个小基线干涉对，其时空分布如图6所示。相比于PSI，采用SBAS时的观测方程数增加了1.6倍，有利于提高形变估计精度。

图6 小基线时空分布Fig.6 Conneted network of interferograms

(5)估计形变速率和时序形变。采用StaMPS的三维相位解缠算法进行相位解缠，并利用最小二乘算法将解缠后的多主影像解缠干涉相位转换到单主影像，根据时空域滤波算法分离出大气相位、残余高程相位，然后估计出形变速率和时序形变。

4 日兰高铁巨野煤田段地表形变

4.1 实验区域概况及所用数据

选择日兰高铁巨野煤田段作为实验区域，具体位置如图7所示。巨野煤田位于华北平原的菏泽市，主要分布在巨野、金乡、梁宝寺等县、镇，总面积约960 km，是华东地区最大的煤田。该煤田已探明储量达55亿t，煤层平均厚度6.6 m，适合建设大型、特大型矿井。目前已建设7座大型矿井，规划年采煤1 800万t，于2007年3月正式投产。

图7 研究区域与Sentinel-1卫星影像位置， 底图为Google Earth光学影像(© Google Maps) [1]Fig.7 Study area and outline of Sentinel-1 SAR data superimposed on a Google Earth optical image (© Google Maps)[1]

日兰高铁穿越巨野煤田(约15 km长)，地下开采引起的地表变形对高铁安全运营构成巨大隐患。ZHU等利用2015年7月至2019年11月的Sentinel-1数据监测了日兰高铁菏泽段沿线地表沉降情况，发现了2个沉降区域：① 一个沉降区域位于菏泽市东郊，沉降主要由于开采地下水引起;② 另一个沉降区域位于巨野煤田，高铁沿线3 km范围内的沉降速率在-8～-2 cm/a。目前，该区段已完成无砟轨道的铺设，为保证高铁的安全运营，有必要继续对该区段进行持续的形变监测。

所用SAR数据同第2.2.1节。由于缺乏地面实测数据，利用S1-40和S1-142两组数据分别采用第3节所述的SBAS干涉处理框架获取实验区域的地表形变信息，然后进行交叉验证。

4.2 PSI与SBAS结果比较

图8为根据S1-40和S1-142数据得到的平均形变速率。负号为地表发生了沉降。通过忽略水平移动，将LoS向形变按式(10)转为垂直向形变：

=cos

(10)

其中,为雷达波入射角。显然，如果地表存在不可忽略的水平移动，则根据式(10)计算的垂直向形变会存在一定的误差。

由图8可以看出，PSI和SBAS两种处理框架计算出的平均形变速率在空间分布上是一致的：2种方法获得了相似的沉降轮廓和形变分布。2者在数值上的区别如图9所示，相关系数分别达到0.959和0.944，说明2者在数值上也是一致的。但是，在形变速率较大时，如形变速率超过16 cm/a后，根据PSI处理流程得到的结果明显小于SBAS结果，即发生了形变低估现象。因此，当形变速率较大时，推荐采用SBAS框架处理SAR数据。后续的分析均基于SBAS算法得到的结果。

图8 采用PSI和SBAS算法得到的平均形变速率Fig.8 Average rates of land displacement results from PSI and SBAS framework

图9 PSI和SBAS算法结果对比Fig.9 Comparison between the PSI and SBAS results using

需要注意的是，SBAS的改善能力是有限度的：如图8中的沉降盆地处(图8中的2处空洞位于郭屯煤矿主采区)，由于地表沉降过于剧烈(对矿区地表沉降而言，通常可达米级)，已经超出了InSAR技术的探测能力，SBAS和PSI两种方法均未能探测到地表形变值。

4.3 形变监测结果的可靠性与精度分析

由于没有地面实测数据，笔者通过S1-40和S1-142两组形变数据的交叉验证来分析结果的可靠性与精度。实验所用数据的空间位置以及成像几何参数不同(表1)，导致S1-40和S1-142的相干点不能完全匹配。为此，首先建立一个规则格网，格网内的每个格子大小为50 m×50 m。然后分别将S1-40和S1-142的相干点根据其坐标值定位到相应的格子内，并将相干点的形变速率赋于该格子。若多个相干点落在同一个格子内，则取形变速率的平均值作为该格子的形变速率。

将S1-40和S1-142采样到同一格网后，共有157 283个有效的同名格子，籍此即可采用交叉验证的方法检验形变结果的可靠性并分析形变监测精度，结果如图10所示。可以看出，S1-40和S1-142的结果总体上较为一致，皮尔逊相关系数达0.88，说明形变监测结果是较为可靠的。2者不一致区域出现在沉降速度超过30 cm/a的快速形变区。由图8可知形变快速区域位于采矿活动区。根据矿区地表移动规律可知，地下开采通常会引起较大的水平移动。一般而言，水平移动为沉降的0.1～0.3倍。而图10中的结果是忽略水平移动后根据雷达波入射角将LoS向形变直接转换为垂直向。若地表存在不可忽略的水平移动，则入射角越大，转换误差也越大。从而导致了在形变快速的沉降盆地处S1-40和S1-142的结果之间存在一定的偏差。

图10 交叉验证结果 Fig.10 Cross-validation of InSAR results

为进一步分析S1-40和S1-142形变结果之间的差异，计算了S1-40和S1-142形变速率差值，根据速率差值绘制了差值直方图并计算了差值的均值和标准差，结果如图10(b)所示。可以看出，2组形变速率差值的均值=0.1 cm/a，标准差=1.1 cm/a。这再次说明S1-40和S1-142形变结果是一致的，并且由差值标准差=1.1 cm/a推测本文的形变监测结果的精度约为1.1 cm/a。这个精度略微低于PSI(一般为亚cm/a)。这是因为本文的形变监测结果中包含了大量的DS目标，而DS的相干质量明显不如PS，因而最终的形变监测精度不如PSI。

4.4 高铁沿线地表的形变分析

图11 高铁沿线平均形变速率Fig.11 Spatial distribution of average rates along RLHR

高铁沿线3 km范围内地表形变的空间分布如图11所示。图11(a)为文献[1]等利用2015年7月至2019年11月的Sentinel-1数据计算得到的平均形变速率；图11(b)为利用2020年10月至2021年4月的Sentinel-1数据计算得到的平均形变速率。由图11可知：

(1)基于短时期内的Sentinel-1数据，利用FEMI算法能够提取更多的相干点，甚至在农田处也提取出了大量相干点。利用密集的相干点可以更加细致的刻画形变的空间分布情况。

(2)煤田西侧(图11中的Zone 1)多处地表出现沉降加剧现象。菏泽地区2021年冬春季节较往年更为干旱，部分农田自2021年2月份以来已经浇灌2次及以上。农田灌溉用水主要抽自浅层地下水，从而导致地表沉降加剧。

(3)煤田中部(图11中的Zone 2)地表沉降仍在持续。该处的形变速率高于农田区域，但又明显小于采煤活动区(图8中白色方框中上部)的形变速率。文献[1]根据形变的空间分布形态以及地下煤炭开采易导致断层“活化”的特点，推测该处下方存在断层，并与采煤活动区的断层连通在一起，使得该处的深层地下水通过断层流入采煤活动区后被抽排走，最终导致地表发生沉降。另外，借助于密集的相干点，图11(b)显示出在日兰高铁南侧也存在明显的地表沉降。图11(a)中由于相干点稀少，沉降的空间分布情况不明显。由Google earth卫星影像可看出该区域基本为农田，由此推测该处地表沉降可能是由于抽取浅层地下水所致，与高铁北侧地表沉降诱因不同。

为进一步研究形变的时间演化过程，提取了4个沉降点(图11中红色三角形)的时序形变量。另外，因为本次实验数据的采集时间为晚秋—早春季节(10月至次年4月初)。为了同前期InSAR观测结果对比，笔者从文献[1]中提取了这4个点在2015—2020年间共4个晚秋—早春时间段的时序形变量，结果如图12所示。每一个时间段的时序形变量均是相对于该时间段中第1景SAR数据的累计形变。由图12可以看出：沉降速度没有发生明显变化，地表在以较为稳定的速度沉降。而日兰高铁菏泽段即将开通运营，可能需要采取必要的控沉措施。

(4)煤田东侧地表总体上较为稳定，巨野高铁站附近没有明显的形变。

图12 M1～M4点在2015—2021年晚秋—早春季节的时序形变(2015—2020年数据来自文献[1])Fig.12 Time-series displacement of M1-M4 between late autumn and early spring from 2015 to 2021 (The time-series displacements from 2015 to 2020 are from Reference[1])

图13(a)为高铁沿线3 km内的地表平均形变速率，结果取自文献[1]。图13(b),(c)显示的是本次实验结果，区别是图13(b),(c)分别为高铁沿线3 km和1 km范围内的地表平均形变速率。可以看出：本次实验结果与文献[1]的结果基本一致，地表沉降仍在继续，在煤田边界以内没有出现加剧现象，但是在边界西侧的农田区域沉降略有增加。高铁沿线的形变速率基本分布在-3.5～-0.5 cm/a。在煤田中心区域有2个沉降漏斗，形变速率集中分布在-3.0～-1.0 cm/a，最大形变速率超过-5.0 cm/a。

图13 日兰高铁沿线形变速率Fig.13 Profile of average rates along RLHR

5 结 论

(1)利用FEMI算法和SBAS干涉处理准确获取了农田区域高密度的相干点及高铁沿线地表沉降的空间分布信息。说明在我国北方地区可以根据10月至次年4月初的少量SAR数据采用本文方法监测非城区沉降。

(2)高铁沿线地表以沉降为主，平均形变速率集中在-3.5～-0.5 cm/a，与前期观测结果(2015—2019年)基本一致，地表仍在持续沉降，但未发现明显加剧现象。

(3)煤田西侧的农田区域出现沉降加剧现象，平均沉降速率在 -3.0～ -0.5 cm/a。考虑到当地在实验期间农田灌溉次数增多，推测该处形变加剧现象是由大量抽取浅层地下水所致。

(4)煤田中部仍然存在两处沉降中心，平均沉降速率集中在 -3.0～ -1.0 cm/a。但是该处沉降原因与西侧不同：根据沉降的空间分布形态与速率大小以及地下煤炭开采容易造成断层等特点，推测该处沉降应是远处地下煤炭开采时的抽排深层地下水所致。在该区域，其影响范围远离工作面2～4 km，并且偏向高铁一侧。因此，在高铁沿线附近开采地下煤炭时应当格外慎重，需要考虑开采导致断层活化以及断层连通造成深层地下水流失等可能性。