基于改進(jìn)MT-InSAR的日蘭高鐵巨野煤田段沉降監(jiān)測

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

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

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

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

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

Squee SAR采用相位三角算法(phase triangle algorithm,PTA)解算DS相位，理論上是真實(shí)值的最大似然估計(jì)，因而相位估計(jì)精度高。但是在低相干區(qū)域或者同質(zhì)像素較少時(shí)，協(xié)方差矩陣通常含有較大誤差，這會削弱PTA的估計(jì)精度。并且PTA算法需要迭代運(yùn)算，較為耗時(shí)。在SqueeSAR的基礎(chǔ)上又演化出幾種變體，如CAESAR (Component extrAction and Selection SAR)，PD-PSInSAR (Phase-Decomposition-based PS InSAR)以及EMI (Eigendecomposition-Based Maximum-likelihood-estimator of Interferometric phase)等。這幾種變體的原理與PTA算法在數(shù)學(xué)形式相互一致，只是加權(quán)策略不同。CAESAR、PD-PSIn SAR以及EMI的基本思想是通過特征分解(Eigenvalue Decomposition,EVD)技術(shù)解算DS相位：即對協(xié)方差矩陣(或加權(quán)后的協(xié)方差矩陣)進(jìn)行特征分解，最大或最小特征值所對應(yīng)的特征向量即為DS相位最優(yōu)估值，其他特征向量則視為噪聲而去除。因此，EVD算法中隱含了相位濾波處理。另外，EVD算法不需要迭代處理，因而計(jì)算速度快。相較于CAESAR和PD-PSInSAR，EMI算法利用相干矩陣調(diào)節(jié)各干涉對的權(quán)重，以期抑制低相干干涉對的影響，能夠在一定程度上降低協(xié)方差矩陣估計(jì)誤差的影響。然而，實(shí)際應(yīng)用中由于SAR數(shù)據(jù)量可能較少、同質(zhì)區(qū)域內(nèi)相位可能不平穩(wěn)、研究區(qū)域相干性低和/或同質(zhì)像元少等原因，協(xié)方差矩陣和相干估值通常是有偏的。

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

另外，標(biāo)準(zhǔn)SqueeSAR的處理流程是在識別出DS點(diǎn)集之后與PS點(diǎn)集合并，然后進(jìn)行單主影像干涉(PSI)處理?？紤]到研究區(qū)域的形變速度較快，筆者在合并PS和DS點(diǎn)集之后根據(jù)一定的時(shí)空基線閾值組合成小基線干涉序列，采用SBAS處理流程估計(jì)形變速率和時(shí)序形變。

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

1 DS相位估計(jì)方法

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

1.1 CAESAR和EMI算法

DS像元內(nèi)沒有強(qiáng)散射體，各基本散射單元的后向散射能力相似。因此，由中心極限定理可知DS的復(fù)反射率服從均值為零的復(fù)高斯分布。若有景SAR影像，則DS像元的時(shí)序復(fù)反射率=[,,…,]服從零均值多元復(fù)高斯分布。因?yàn)榫禐榱悖?span id="j5i0abt0b" class="emphasis_italic">的統(tǒng)計(jì)特征可用協(xié)方差矩陣描述。根據(jù)同質(zhì)區(qū)域內(nèi)的像元計(jì)算得到的樣本協(xié)方差矩陣是真實(shí)協(xié)方差矩陣的最大似然估計(jì)。樣本協(xié)方差矩陣的計(jì)算公式為

(1)

式中，為同質(zhì)區(qū)域；為同質(zhì)區(qū)域內(nèi)的像元數(shù)目；上標(biāo)H為共軛轉(zhuǎn)置。

為了提高DS的相干性，通常需要對其干涉相位進(jìn)行濾波處理(即空間自適應(yīng)多視處理)。濾波后干涉相位的一致性被破壞，導(dǎo)致協(xié)方差矩陣秩虧，因而無法根據(jù)式(1)直接解算DS的真實(shí)相位。因?yàn)?span id="j5i0abt0b" class="emphasis_italic">服從復(fù)Wishart分布，F(xiàn)erretti等提出了通過最大化的概率來估計(jì)的PTA算法。PTA算法需要迭代運(yùn)算，較為耗時(shí)。Fornaro等利用譜分解算法根據(jù)估計(jì)，即CAESAR算法：

(2)

(3)

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

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

(4)

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

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

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

=2(1-)

(5)

可以看出，為單調(diào)遞增函數(shù)。當(dāng)=0時(shí)，表示干涉相位沒有有效信息，此時(shí)=0；當(dāng)=1時(shí)，說明干涉相位不含噪聲，此時(shí)→∞。

相比于相干值只能在直觀上反映觀測量的好壞，F(xiàn)isher信息量在數(shù)值上給出了觀測量中包含待估參數(shù)信息的多寡。并且，在低相干區(qū)，利用加權(quán)比更穩(wěn)健。另外，利用Fisher信息量的迭代優(yōu)化算法不易陷入局部極值。因此，本文利用Fisher信息矩陣代替，實(shí)現(xiàn)基于Fisher信息量的DS相位估計(jì)算法，即FEMI算法。由于同一DS的同質(zhì)像素?cái)?shù)目在不同的干涉對上保持不變，因此可令=1。利用Fisher信息矩陣代替，式(4)可改寫為

(6)

式中，′,=[],。

2 實(shí)驗(yàn)驗(yàn)證

2.1 模擬數(shù)據(jù)實(shí)驗(yàn)

..數(shù)據(jù)模擬

模擬數(shù)據(jù)的相干值由下式計(jì)算：

=(-)e-Δ+

(7)

其中，Δ為干涉對的時(shí)間基線；為信號的時(shí)間相干長度，越大，干涉對的時(shí)間相干性越好；和為相干起始值和結(jié)束值,=0時(shí)表明相干值呈負(fù)指數(shù)快速衰減，干涉對只能在短期內(nèi)保持相干性，時(shí)間基線Δ≈3時(shí)相干性近乎于0。

采用式(7)模擬2種相干類型數(shù)據(jù)：① 短期相干型，即=0；②長期相干型，即=02。2種類型的設(shè)置為0.6，設(shè)為50 d。因?yàn)楹笪牟捎?6景真實(shí)Sentinel-1衛(wèi)星影像開展形變監(jiān)測研究，為此本節(jié)模擬16景SAR數(shù)據(jù)，并仿照Sentinel-1衛(wèi)星將時(shí)間間隔設(shè)為12 d、波長設(shè)為55.6 mm。據(jù)式(7)計(jì)算相干陣(圖1)，并由相干陣生成16景SAR數(shù)據(jù)，每景SAR數(shù)據(jù)包含300個(gè)同質(zhì)像元。為每景SAR數(shù)據(jù)添加變形速率為5 mm/a的形變相位。

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

..實(shí)驗(yàn)結(jié)果分析

根據(jù)2.1.1節(jié)的2組模擬數(shù)據(jù)，分別采用EMI和FEMI算法估計(jì)DS相位，重復(fù)進(jìn)行20萬次，然后計(jì)算均方根誤差():

(8)

圖2 EMI與FEMI算法精度評估Fig.2 Performance assessment of EMI and FEMI methods using simulated data

(1)由圖2(a)可知，對短期相干型數(shù)據(jù)，EMI的相位估計(jì)精度顯著降低。因?yàn)槎唐谙喔尚蛿?shù)據(jù)的時(shí)間失相干嚴(yán)重，干涉對的時(shí)間基線超過150 d時(shí)幾乎沒有相干性(圖1(a))，即干涉相位中幾乎全為噪聲，不包含任何有效信息。而通過樣本計(jì)算得到的相干值通常是有偏的，使得這些噪聲傳遞到估計(jì)的DS相位中。而FEMI通過使用Fisher信息作為權(quán)陣，大幅度降低低相干干涉對的權(quán)重，抑制了誤差的傳播，因而FEMI算法較EMI精度更高也更加穩(wěn)健。

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

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

2.2 真實(shí)數(shù)據(jù)實(shí)驗(yàn)

..實(shí)驗(yàn)數(shù)據(jù)

選擇Sentinel-1A衛(wèi)星獲取的SAR影像作為實(shí)驗(yàn)數(shù)據(jù)，波長約為5.6 cm，空間分辨率約為20 m × 5 m。Sentinel-1A衛(wèi)星的重訪周期為12 d，在6個(gè)月內(nèi)(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組數(shù)據(jù)的部分參數(shù)見表1。另外，采用空間分辨率為30 m的AW3D30 DSM數(shù)據(jù)補(bǔ)償?shù)匦蜗辔弧?/p>

表1 Sentinel-1衛(wèi)星數(shù)據(jù)主要參數(shù)

..實(shí)驗(yàn)結(jié)果分析

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

圖3 干涉對20201012—20210410的干涉相位和相干系數(shù)Fig.3 Interferograms and the coherence of 20201012-20210410

圖4 干涉對20201007—20210405的干涉相位和相干系數(shù)Fig.4 Interferograms and the coherence of 20201007-20210405

3 聯(lián)合PS和DS的SBAS處理框架

礦區(qū)形變通常較為快速，并且高鐵沿線植被較多，時(shí)間失相干嚴(yán)重。為提高觀測點(diǎn)密度、增加觀測方程、保證每個(gè)干涉對的相干質(zhì)量，采用聯(lián)合PS和DS的SBAS干涉處理策略估計(jì)形變。假設(shè)所有的SAR影像均已配準(zhǔn)至公共主影像：即S1-40數(shù)據(jù)配準(zhǔn)至2021-01-04的影像參考幾何、S1-142數(shù)據(jù)配準(zhǔn)至2021-01-11的影像參考幾何。聯(lián)合PS和DS的SBAS處理框架如圖5所示，具體步驟：

(1)提取PS點(diǎn)集與相位。采用單主影像策略，利用PSI技術(shù)獲取PS點(diǎn)集及相位，PSI原理與流程見文獻(xiàn)[32]。PSI分析過程中，振幅離差設(shè)為0.4。

(2)提取DS點(diǎn)集與相位。首先根據(jù)SAR幅度序列采用K-S檢驗(yàn)算法識別DS初選點(diǎn)，顯著性水平設(shè)置為0.05。窗口大小設(shè)置為9×35，同質(zhì)像元(SHP)閾值設(shè)為25。然后，根據(jù)配準(zhǔn)后的SAR數(shù)據(jù)和AW3D30 DSM生成全組合(即16×152=120)差分干涉序列。根據(jù)DS初選點(diǎn)和全組合差分干涉序列，根據(jù)式(1)計(jì)算協(xié)方差矩陣，然后利用本文提出的FEMI算法(式(6))估計(jì)DS相位。

(3)篩選DS點(diǎn)集。按照式(9)計(jì)算DS相位估值的擬合度，選取>0.7的像素作為最終DS點(diǎn)。

(9)

圖5 數(shù)據(jù)處理流程Fig.5 Framework of data processing

式中，Re為取復(fù)數(shù)的實(shí)部

(4)合并PS和DS點(diǎn)集，并生成小基線干涉序列。在組合小基線對時(shí)，時(shí)間基線閾值設(shè)置為36 d，空間基線閾值設(shè)置為200 m。另外，需要確保沒有非連通子集出現(xiàn)以便增加解的可靠性。最終，2組數(shù)據(jù)集均組合得到42個(gè)小基線干涉對，其時(shí)空分布如圖6所示。相比于PSI，采用SBAS時(shí)的觀測方程數(shù)增加了1.6倍，有利于提高形變估計(jì)精度。

圖6 小基線時(shí)空分布Fig.6 Conneted network of interferograms

(5)估計(jì)形變速率和時(shí)序形變。采用StaMPS的三維相位解纏算法進(jìn)行相位解纏，并利用最小二乘算法將解纏后的多主影像解纏干涉相位轉(zhuǎn)換到單主影像，根據(jù)時(shí)空域?yàn)V波算法分離出大氣相位、殘余高程相位，然后估計(jì)出形變速率和時(shí)序形變。

4 日蘭高鐵巨野煤田段地表形變

4.1 實(shí)驗(yàn)區(qū)域概況及所用數(shù)據(jù)

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

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

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

4.2 PSI與SBAS結(jié)果比較

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

=cos

(10)

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

由圖8可以看出，PSI和SBAS兩種處理框架計(jì)算出的平均形變速率在空間分布上是一致的：2種方法獲得了相似的沉降輪廓和形變分布。2者在數(shù)值上的區(qū)別如圖9所示，相關(guān)系數(shù)分別達(dá)到0.959和0.944，說明2者在數(shù)值上也是一致的。但是，在形變速率較大時(shí)，如形變速率超過16 cm/a后，根據(jù)PSI處理流程得到的結(jié)果明顯小于SBAS結(jié)果，即發(fā)生了形變低估現(xiàn)象。因此，當(dāng)形變速率較大時(shí)，推薦采用SBAS框架處理SAR數(shù)據(jù)。后續(xù)的分析均基于SBAS算法得到的結(jié)果。

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

圖9 PSI和SBAS算法結(jié)果對比Fig.9 Comparison between the PSI and SBAS results using

需要注意的是，SBAS的改善能力是有限度的：如圖8中的沉降盆地處(圖8中的2處空洞位于郭屯煤礦主采區(qū))，由于地表沉降過于劇烈(對礦區(qū)地表沉降而言，通常可達(dá)米級)，已經(jīng)超出了InSAR技術(shù)的探測能力，SBAS和PSI兩種方法均未能探測到地表形變值。

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

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

將S1-40和S1-142采樣到同一格網(wǎng)后，共有157 283個(gè)有效的同名格子，籍此即可采用交叉驗(yàn)證的方法檢驗(yàn)形變結(jié)果的可靠性并分析形變監(jiān)測精度，結(jié)果如圖10所示?？梢钥闯?，S1-40和S1-142的結(jié)果總體上較為一致，皮爾遜相關(guān)系數(shù)達(dá)0.88，說明形變監(jiān)測結(jié)果是較為可靠的。2者不一致區(qū)域出現(xiàn)在沉降速度超過30 cm/a的快速形變區(qū)。由圖8可知形變快速區(qū)域位于采礦活動區(qū)。根據(jù)礦區(qū)地表移動規(guī)律可知，地下開采通常會引起較大的水平移動。一般而言，水平移動為沉降的0.1～0.3倍。而圖10中的結(jié)果是忽略水平移動后根據(jù)雷達(dá)波入射角將LoS向形變直接轉(zhuǎn)換為垂直向。若地表存在不可忽略的水平移動，則入射角越大，轉(zhuǎn)換誤差也越大。從而導(dǎo)致了在形變快速的沉降盆地處S1-40和S1-142的結(jié)果之間存在一定的偏差。

圖10 交叉驗(yàn)證結(jié)果 Fig.10 Cross-validation of InSAR results

為進(jìn)一步分析S1-40和S1-142形變結(jié)果之間的差異，計(jì)算了S1-40和S1-142形變速率差值，根據(jù)速率差值繪制了差值直方圖并計(jì)算了差值的均值和標(biāo)準(zhǔn)差，結(jié)果如圖10(b)所示?？梢钥闯觯?組形變速率差值的均值=0.1 cm/a，標(biāo)準(zhǔn)差=1.1 cm/a。這再次說明S1-40和S1-142形變結(jié)果是一致的，并且由差值標(biāo)準(zhǔn)差=1.1 cm/a推測本文的形變監(jiān)測結(jié)果的精度約為1.1 cm/a。這個(gè)精度略微低于PSI(一般為亞cm/a)。這是因?yàn)楸疚牡男巫儽O(jiān)測結(jié)果中包含了大量的DS目標(biāo)，而DS的相干質(zhì)量明顯不如PS，因而最終的形變監(jiān)測精度不如PSI。

4.4 高鐵沿線地表的形變分析

圖11 高鐵沿線平均形變速率Fig.11 Spatial distribution of average rates along RLHR

高鐵沿線3 km范圍內(nèi)地表形變的空間分布如圖11所示。圖11(a)為文獻(xiàn)[1]等利用2015年7月至2019年11月的Sentinel-1數(shù)據(jù)計(jì)算得到的平均形變速率；圖11(b)為利用2020年10月至2021年4月的Sentinel-1數(shù)據(jù)計(jì)算得到的平均形變速率。由圖11可知：

(1)基于短時(shí)期內(nèi)的Sentinel-1數(shù)據(jù)，利用FEMI算法能夠提取更多的相干點(diǎn)，甚至在農(nóng)田處也提取出了大量相干點(diǎn)。利用密集的相干點(diǎn)可以更加細(xì)致的刻畫形變的空間分布情況。

(2)煤田西側(cè)(圖11中的Zone 1)多處地表出現(xiàn)沉降加劇現(xiàn)象。菏澤地區(qū)2021年冬春季節(jié)較往年更為干旱，部分農(nóng)田自2021年2月份以來已經(jīng)澆灌2次及以上。農(nóng)田灌溉用水主要抽自淺層地下水，從而導(dǎo)致地表沉降加劇。

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

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

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

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

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

5 結(jié) 論

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

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

(3)煤田西側(cè)的農(nóng)田區(qū)域出現(xiàn)沉降加劇現(xiàn)象，平均沉降速率在 -3.0～ -0.5 cm/a。考慮到當(dāng)?shù)卦趯?shí)驗(yàn)期間農(nóng)田灌溉次數(shù)增多，推測該處形變加劇現(xiàn)象是由大量抽取淺層地下水所致。

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