祝傳廣，張繼賢，鄧喀中，龍四春，張立亞，吳文豪

(1.湖南科技大學 地球科學與空間信息工程學院，湖南 湘潭 411201；2.湖南科技大學 測繪遙感信息工程湖南省重點實驗室，湖南 湘潭 411201；3.國家測繪產品質量檢驗測試中心，北京 100830；4.中國礦業(yè)大學 環(huán)境與測繪學院，江蘇 徐州 221116)

日照—蘭考高速鐵路(日蘭高鐵)自2017年開工建設，預計2021年底全線開通運營。該鐵路在菏澤市約有15 km路段穿越巨野煤田，地下煤礦開采引起的地表沉降對高鐵安全運營構成巨大隱患。已有研究顯示，高鐵沿線在2015—2019年間存在較大的地表沉降現象。為保證鐵路的安全運營，需要對高鐵軌道、路基及沿線地表進行持續(xù)的沉降監(jiān)測。

星載合成孔徑雷達干涉測量(InSAR)技術因其高分辨率、廣覆蓋、全天候等優(yōu)點成為常用的大地測量技術之一。其中，MT-InSAR技術借助特殊的相干目標能夠以mm/a～cm/a的精度實現長時間跨度的形變監(jiān)測。這些相干目標主要來自永久散射體(Persistent Scatterer,PS)和分布式散射體(Distributed Scatterer,DS)。PS集中分布在城區(qū)，以人工建/構筑物為主。在缺少人工建/構物的非城區(qū)，PS點過于稀少，難以詳實反映形變的空間分布情況。并且，低密度的空間點目標會增加相位解纏難度從而降低形變估計精度。

DS廣泛分布于自然界，利用DS能夠大幅度提高觀測點密度。DS的后向散射能力中等，能夠在短時期內保持一定的相干性。MT-InSAR中的小基線(SBAS-InSAR，下文簡稱SBAS)技術就是利用這一特點，通過限制干涉對的時空基線長度以便盡可能地保持干涉對的相干性，從而將DS納入MT-InSAR分析，達到增加觀測點密度的目的。相比于SBAS僅僅利用短時空基線干涉對，以Squee SAR為代表的DS-InSAR技術則是利用全組合干涉對解算DS相位最優(yōu)估值。采用全組合干涉對能夠進一步提高SAR信息的利用率，然后通過聯合DS和PS點集大幅度提高觀測點密度。

Squee SAR采用相位三角算法(phase triangle algorithm,PTA)解算DS相位，理論上是真實值的最大似然估計，因而相位估計精度高。但是在低相干區(qū)域或者同質像素較少時，協方差矩陣通常含有較大誤差，這會削弱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相位最優(yōu)估值，其他特征向量則視為噪聲而去除。因此，EVD算法中隱含了相位濾波處理。另外，EVD算法不需要迭代處理，因而計算速度快。相較于CAESAR和PD-PSInSAR，EMI算法利用相干矩陣調節(jié)各干涉對的權重，以期抑制低相干干涉對的影響，能夠在一定程度上降低協方差矩陣估計誤差的影響。然而，實際應用中由于SAR數據量可能較少、同質區(qū)域內相位可能不平穩(wěn)、研究區(qū)域相干性低和/或同質像元少等原因，協方差矩陣和相干估值通常是有偏的。

針對上述問題，筆者在EMI的基礎上提出了一種基于Fisher信息量的改進EMI算法(FEMI)。Fisher信息量是觀測到的變量(如干涉相位)中所攜帶的關于未知參數(如DS的真實相位)信息量的一種度量，通常用于計算最大似然估計中的協方差矩陣。利用Fisher信息量代替相干矩陣作為權陣能夠抑制相干值估計誤差對相位估計精度的影響。

另外，標準SqueeSAR的處理流程是在識別出DS點集之后與PS點集合并，然后進行單主影像干涉(PSI)處理。考慮到研究區(qū)域的形變速度較快，筆者在合并PS和DS點集之后根據一定的時空基線閾值組合成小基線干涉序列，采用SBAS處理流程估計形變速率和時序形變。

首先利用模擬數據和真實SAR數據驗證了FEMI算法的可靠性，然后基于2020年10月至2021年4月間的Sentinel-1衛(wèi)星數據，采用SBAS處理框架獲取了日蘭高鐵巨野煤田段的地表沉降信息，并結合2015—2019年的沉降監(jiān)測資料分析了地表沉降的成因及時空演化機理。

1 DS相位估計方法

首先介紹DS相位優(yōu)化估計的基本思想，然后分析CAESAR和EMI算法原理及其特點，進而引出Fisher信息量并構建基于Fisher信息量的改進EMI算法，即FEMI。

1.1 CAESAR和EMI算法

DS像元內沒有強散射體，各基本散射單元的后向散射能力相似。因此，由中心極限定理可知DS的復反射率服從均值為零的復高斯分布。若有景SAR影像，則DS像元的時序復反射率=[,,…,]服從零均值多元復高斯分布。因為均值為零，所以的統(tǒng)計特征可用協方差矩陣描述。根據同質區(qū)域內的像元計算得到的樣本協方差矩陣是真實協方差矩陣的最大似然估計。樣本協方差矩陣的計算公式為

(1)

式中，為同質區(qū)域；為同質區(qū)域內的像元數目；上標H為共軛轉置。

為了提高DS的相干性，通常需要對其干涉相位進行濾波處理(即空間自適應多視處理)。濾波后干涉相位的一致性被破壞，導致協方差矩陣秩虧，因而無法根據式(1)直接解算DS的真實相位。因為服從復Wishart分布，Ferretti等提出了通過最大化的概率來估計的PTA算法。PTA算法需要迭代運算，較為耗時。Fornaro等利用譜分解算法根據估計，即CAESAR算法：

(2)

(3)

式中，,=[||],；′,=[||||],；,為DS的干涉相位的觀測值，,=-。

因為協方差矩陣需要根據同質區(qū)域內的樣本計算得到，而同質區(qū)域是通過假設檢驗得到，同質區(qū)域內可能存在不服從高斯分布的像元，這會導致存在偏差。另外，在低相干區(qū)域僅采用最大特征值所對應的特征向量代表協方差矩陣會引起較大偏差，因為此時可能需要多個甚至全部的譜分量才能有效反映。為此，ANSARI等提出了EMI算法，通過給各干涉對按其相干質量賦予不同的權重以抑制低相干的影響，即：

(4)

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

由樣本估計的相干矩陣是有偏的，尤其是在低相干區(qū)域，因而不能保證是正定-半正定矩陣。為此，應用式(4)求解DS相位之前需要預先判斷是否滿足正定-半正定。若不滿足，還需要通過添加阻尼因子等方式調整。

Fisher信息量是一種度量觀測量(如干涉相位的觀測值,)中包含待求參數(如干涉相位的真值,)信息多寡的指標。觀測到的干涉相位,是由同質區(qū)域內的樣本估計得到，估計誤差的方差漸進逼近Cramer-Rao下界，而Fisher信息量能夠反映估計誤差的方差與Cramer-Rao下界的逼近程度。干涉相位,的Fisher信息量可表示為

=2(1-)

(5)

可以看出，為單調遞增函數。當=0時，表示干涉相位沒有有效信息，此時=0；當=1時，說明干涉相位不含噪聲，此時→∞。

相比于相干值只能在直觀上反映觀測量的好壞，Fisher信息量在數值上給出了觀測量中包含待估參數信息的多寡。并且，在低相干區(qū)，利用加權比更穩(wěn)健。另外，利用Fisher信息量的迭代優(yōu)化算法不易陷入局部極值。因此，本文利用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衛(wèi)星影像開展形變監(jiān)測研究，為此本節(jié)模擬16景SAR數據，并仿照Sentinel-1衛(wèi)星將時間間隔設為12 d、波長設為55.6 mm。據式(7)計算相干陣(圖1)，并由相干陣生成16景SAR數據，每景SAR數據包含300個同質像元。為每景SAR數據添加變形速率為5 mm/a的形變相位。

圖1 模擬相干陣 Fig.1 Simulated coherence matrix

..實驗結果分析

根據2.1.1節(jié)的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精度更高也更加穩(wěn)健。

(2)由圖2(b)可知，對于長期相干型數據，FEMI和EMI算法的性能相近，并且接近CRLB。這是因為，對于長期相干型數據，時間基線最大的干涉對其相干性仍然大于0.2(圖1(b))，所有的干涉對都能夠提供有效信息。

綜上可知，FEMI算法更加穩(wěn)健。當研究區(qū)域的相干性較好時，FEMI和EMI算法性能相當。當研究區(qū)域的相干性較差時，則建議采用FEMI算法。雖然增大同質濾波的窗口理論上能夠提高相干值的估計精度。但是，窗口過大時難以保證信號的空間平穩(wěn)性。即使進行了地形相位補償，大氣效應和形變信息的空間變化也會導致信號空間失相關。

2.2 真實數據實驗

..實驗數據

選擇Sentinel-1A衛(wèi)星獲取的SAR影像作為實驗數據，波長約為5.6 cm，空間分辨率約為20 m × 5 m。Sentinel-1A衛(wèi)星的重訪周期為12 d，在6個月內(10月—次年4月)能夠提供16景衛(wèi)星影像。筆者收集了Sentinel-1A衛(wèi)星在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衛(wèi)星數據主要參數

..實驗結果分析

通常，干涉對的時間基線越長相干性越差。檢查干涉對的干涉相位圖和相干系數圖能夠直觀上評估算法的效果。對于本文數據而言，最大時間基線(180 d)對應的2組干涉對分別為：20201012—20210410(S1-40)和20201007—20210405(S1-142)。圖3，4分別顯示了這2組數據的干涉相位與相干系數的局部細節(jié)，為對比分析，同時顯示了原始的以及EMI算法得到的干涉相位圖和相干系數圖。由圖3，4可以看出，EMI和FEMI算法均能顯著改善相位質量。在相干性高的區(qū)域，兩種算法的效果近似一致。但是在低相干區(qū)域，FEMI算法的表現優(yōu)于EMI。這與2.1節(jié)利用模擬數據得到的結果相一致。

圖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處理框架

礦區(qū)形變通常較為快速，并且高鐵沿線植被較多，時間失相干嚴重。為提高觀測點密度、增加觀測方程、保證每個干涉對的相干質量，采用聯合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 實驗區(qū)域概況及所用數據

選擇日蘭高鐵巨野煤田段作為實驗區(qū)域，具體位置如圖7所示。巨野煤田位于華北平原的菏澤市，主要分布在巨野、金鄉(xiāng)、梁寶寺等縣、鎮(zhèn)，總面積約960 km，是華東地區(qū)最大的煤田。該煤田已探明儲量達55億t，煤層平均厚度6.6 m，適合建設大型、特大型礦井。目前已建設7座大型礦井，規(guī)劃年采煤1 800萬t，于2007年3月正式投產。

圖7 研究區(qū)域與Sentinel-1衛(wèi)星影像位置， 底圖為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數據監(jiān)測了日蘭高鐵菏澤段沿線地表沉降情況，發(fā)現了2個沉降區(qū)域：① 一個沉降區(qū)域位于菏澤市東郊，沉降主要由于開采地下水引起;② 另一個沉降區(qū)域位于巨野煤田，高鐵沿線3 km范圍內的沉降速率在-8～-2 cm/a。目前，該區(qū)段已完成無砟軌道的鋪設，為保證高鐵的安全運營，有必要繼續(xù)對該區(qū)段進行持續(xù)的形變監(jiān)測。

所用SAR數據同第2.2.1節(jié)。由于缺乏地面實測數據，利用S1-40和S1-142兩組數據分別采用第3節(jié)所述的SBAS干涉處理框架獲取實驗區(qū)域的地表形變信息，然后進行交叉驗證。

4.2 PSI與SBAS結果比較

圖8為根據S1-40和S1-142數據得到的平均形變速率。負號為地表發(fā)生了沉降。通過忽略水平移動，將LoS向形變按式(10)轉為垂直向形變：

=cos

(10)

其中,為雷達波入射角。顯然，如果地表存在不可忽略的水平移動，則根據式(10)計算的垂直向形變會存在一定的誤差。

由圖8可以看出，PSI和SBAS兩種處理框架計算出的平均形變速率在空間分布上是一致的：2種方法獲得了相似的沉降輪廓和形變分布。2者在數值上的區(qū)別如圖9所示，相關系數分別達到0.959和0.944，說明2者在數值上也是一致的。但是，在形變速率較大時，如形變速率超過16 cm/a后，根據PSI處理流程得到的結果明顯小于SBAS結果，即發(fā)生了形變低估現象。因此，當形變速率較大時，推薦采用SBAS框架處理SAR數據。后續(xù)的分析均基于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處空洞位于郭屯煤礦主采區(qū))，由于地表沉降過于劇烈(對礦區(qū)地表沉降而言，通?？蛇_米級)，已經超出了InSAR技術的探測能力，SBAS和PSI兩種方法均未能探測到地表形變值。

4.3 形變監(jiān)測結果的可靠性與精度分析

由于沒有地面實測數據，筆者通過S1-40和S1-142兩組形變數據的交叉驗證來分析結果的可靠性與精度。實驗所用數據的空間位置以及成像幾何參數不同(表1)，導致S1-40和S1-142的相干點不能完全匹配。為此，首先建立一個規(guī)則格網，格網內的每個格子大小為50 m×50 m。然后分別將S1-40和S1-142的相干點根據其坐標值定位到相應的格子內，并將相干點的形變速率賦于該格子。若多個相干點落在同一個格子內，則取形變速率的平均值作為該格子的形變速率。

將S1-40和S1-142采樣到同一格網后，共有157 283個有效的同名格子，籍此即可采用交叉驗證的方法檢驗形變結果的可靠性并分析形變監(jiān)測精度，結果如圖10所示?？梢钥闯觯琒1-40和S1-142的結果總體上較為一致，皮爾遜相關系數達0.88，說明形變監(jiān)測結果是較為可靠的。2者不一致區(qū)域出現在沉降速度超過30 cm/a的快速形變區(qū)。由圖8可知形變快速區(qū)域位于采礦活動區(qū)。根據礦區(qū)地表移動規(guī)律可知，地下開采通常會引起較大的水平移動。一般而言，水平移動為沉降的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)所示?？梢钥闯觯?組形變速率差值的均值=0.1 cm/a，標準差=1.1 cm/a。這再次說明S1-40和S1-142形變結果是一致的，并且由差值標準差=1.1 cm/a推測本文的形變監(jiān)測結果的精度約為1.1 cm/a。這個精度略微低于PSI(一般為亞cm/a)。這是因為本文的形變監(jiān)測結果中包含了大量的DS目標，而DS的相干質量明顯不如PS，因而最終的形變監(jiān)測精度不如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)多處地表出現沉降加劇現象。菏澤地區(qū)2021年冬春季節(jié)較往年更為干旱，部分農田自2021年2月份以來已經澆灌2次及以上。農田灌溉用水主要抽自淺層地下水，從而導致地表沉降加劇。

(3)煤田中部(圖11中的Zone 2)地表沉降仍在持續(xù)。該處的形變速率高于農田區(qū)域，但又明顯小于采煤活動區(qū)(圖8中白色方框中上部)的形變速率。文獻[1]根據形變的空間分布形態(tài)以及地下煤炭開采易導致斷層“活化”的特點，推測該處下方存在斷層，并與采煤活動區(qū)的斷層連通在一起，使得該處的深層地下水通過斷層流入采煤活動區(qū)后被抽排走，最終導致地表發(fā)生沉降。另外，借助于密集的相干點，圖11(b)顯示出在日蘭高鐵南側也存在明顯的地表沉降。圖11(a)中由于相干點稀少，沉降的空間分布情況不明顯。由Google earth衛(wèi)星影像可看出該區(qū)域基本為農田，由此推測該處地表沉降可能是由于抽取淺層地下水所致，與高鐵北側地表沉降誘因不同。

為進一步研究形變的時間演化過程，提取了4個沉降點(圖11中紅色三角形)的時序形變量。另外，因為本次實驗數據的采集時間為晚秋—早春季節(jié)(10月至次年4月初)。為了同前期InSAR觀測結果對比，筆者從文獻[1]中提取了這4個點在2015—2020年間共4個晚秋—早春時間段的時序形變量，結果如圖12所示。每一個時間段的時序形變量均是相對于該時間段中第1景SAR數據的累計形變。由圖12可以看出：沉降速度沒有發(fā)生明顯變化，地表在以較為穩(wěn)定的速度沉降。而日蘭高鐵菏澤段即將開通運營，可能需要采取必要的控沉措施。

(4)煤田東側地表總體上較為穩(wěn)定，巨野高鐵站附近沒有明顯的形變。

圖12 M1～M4點在2015—2021年晚秋—早春季節(jié)的時序形變(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)顯示的是本次實驗結果，區(qū)別是圖13(b),(c)分別為高鐵沿線3 km和1 km范圍內的地表平均形變速率?？梢钥闯觯罕敬螌嶒灲Y果與文獻[1]的結果基本一致，地表沉降仍在繼續(xù)，在煤田邊界以內沒有出現加劇現象，但是在邊界西側的農田區(qū)域沉降略有增加。高鐵沿線的形變速率基本分布在-3.5～-0.5 cm/a。在煤田中心區(qū)域有2個沉降漏斗，形變速率集中分布在-3.0～-1.0 cm/a，最大形變速率超過-5.0 cm/a。

圖13 日蘭高鐵沿線形變速率Fig.13 Profile of average rates along RLHR

5 結 論

(1)利用FEMI算法和SBAS干涉處理準確獲取了農田區(qū)域高密度的相干點及高鐵沿線地表沉降的空間分布信息。說明在我國北方地區(qū)可以根據10月至次年4月初的少量SAR數據采用本文方法監(jiān)測非城區(qū)沉降。

(2)高鐵沿線地表以沉降為主，平均形變速率集中在-3.5～-0.5 cm/a，與前期觀測結果(2015—2019年)基本一致，地表仍在持續(xù)沉降，但未發(fā)現明顯加劇現象。

(3)煤田西側的農田區(qū)域出現沉降加劇現象，平均沉降速率在 -3.0～ -0.5 cm/a。考慮到當地在實驗期間農田灌溉次數增多，推測該處形變加劇現象是由大量抽取淺層地下水所致。

(4)煤田中部仍然存在兩處沉降中心，平均沉降速率集中在 -3.0～ -1.0 cm/a。但是該處沉降原因與西側不同：根據沉降的空間分布形態(tài)與速率大小以及地下煤炭開采容易造成斷層等特點，推測該處沉降應是遠處地下煤炭開采時的抽排深層地下水所致。在該區(qū)域，其影響范圍遠離工作面2～4 km，并且偏向高鐵一側。因此，在高鐵沿線附近開采地下煤炭時應當格外慎重，需要考慮開采導致斷層活化以及斷層連通造成深層地下水流失等可能性。