楊嘉威，譚志祥，鄧喀中

(中國礦業(yè)大學 環(huán)境與測繪學院，江蘇 徐州 221116)

煤礦開采伴隨著上覆巖層及地表沉陷，嚴重時引發(fā)一系列地質災害，在監(jiān)測變形時，傳統(tǒng)方法通常在地表布設觀測站從而獲得實測數(shù)據(jù)，雖能獲得較高精度的形變，但成本較高且監(jiān)測點易受破壞，具有較大局限性。合成孔徑雷達干涉測量(Interferometric Synthetic Aperture Radar，InSAR)作為變形監(jiān)測領域的有效技術，憑借其全天時、全天候、大范圍、高精度等特點，已經(jīng)被廣泛應用于地震、冰川、滑坡、礦山等區(qū)域[1，2]。常規(guī)差分InSAR(Differential InSAR)對覆蓋同一地區(qū)兩景雷達SAR影像進行干涉處理，從而獲得時間間隔內(nèi)的地表微小形變。隨著SAR數(shù)據(jù)的積累，永久散射體(Persistent Scatterer，PS)InSAR、短基線集(Small Baseline Subset，SBAS)InSAR等時序InSAR技術[3-5]得到發(fā)展，時序技術使用多景覆蓋同一地區(qū)的SAR影像，先將影像對干涉處理，再篩選出時段內(nèi)具有高相干性表現(xiàn)的像元目標，重構其形變信息。

在城市地區(qū)，由于人工建筑等永久散射體的廣泛存在，高相干性像元數(shù)量適中，使時序技術能有效克服常規(guī)D-InSAR時空失相干的不足，得到良好的形變效果。但在人工建筑目標較少的礦區(qū)，地表通常被大量植被和農(nóng)田覆蓋，潛在永久散射體等像元目標不足，失相干影響依然存在，時序InSAR的應用受到一定局限性，因此引入了分布式散射體(Distributed Scatterer，DS)目標，它們在SAR影像中具有多個后向散射機制，如裸露的土壤、植被稀疏的土地和水體等[6，7]。關于DS的處理，2010年Ferretti等[8]首先提出了SqueeSAR技術，考慮DS像元目標的統(tǒng)計特征，將處理過程分為同質像元識別和相位優(yōu)化兩個過程，并重構出DS像元的形變信息。2015年，Cao等[9]揭示了相位優(yōu)化過程的數(shù)學原理，其中最大似然法、最小二乘法和特征值分解法等都具有相似數(shù)學形式，其不同點在于對權重矩陣的設定，并在此基礎上引出了等權和相干加權兩種優(yōu)化方法。同年，Jiang等[10]針對原有同質像元識別方法的局限性，將同質像元識別過程的假設檢驗問題轉換為置信區(qū)間求取問題，提出快速同質像元選取算法(FaSHPS)并驗證了其效果，算法在計算效率上相較于KS等檢驗方法具有較大提升。

在本文中，為了更高效監(jiān)測礦區(qū)地表和上覆建構筑物的形變，并獲得其變形規(guī)律，選用2018年11月至2019年11月間30景覆蓋研究區(qū)域的Sentinel-1A衛(wèi)星影像，使用最大似然估計優(yōu)化干涉圖相位抑制DS像元的噪聲，再聯(lián)合PS像元共同解算形變信息，從而獲得受采動影響的鐵路沿線移動變形情況與沉降規(guī)律。

1 DS-InSAR方法與優(yōu)化改進

本文利用時序InSAR技術，聯(lián)合PS和DS兩類像元提高相干點的密度與形變質量。主要過程概括為DS像元選取、PS像元選取和聯(lián)合像元形變解算三步，如圖1所示。

圖1 聯(lián)合PS和DS的時序InSAR處理流程圖

1.1 DS像元選取

選取DS像元目標包括同質像元識別和相位優(yōu)化兩個過程[11，12]，本文先采用FaSHPS算法識別同質像元點，該算法的核心是將同質像元識別過程的假設檢驗問題轉換為求解置信區(qū)間[11]，其形式見式(1)：

式中，μ(P)=A(P)=∑Ai(P)/N為像元P的時序強度值的期望；Z1-α/2是標準正態(tài)分位點。將相鄰區(qū)域內(nèi)平均強度位于該區(qū)間的像元視為參考點P的同質像元，并對SAR影像中每個像元進行同質像元識別，同質像元數(shù)目大于給定閾值的像元為初始DS候選點。

采用最大似然法優(yōu)化同質像元相位進一步得候選DS點，優(yōu)化相位的最大似然估計形式見式(2)：

式中，θML=[θ1，θ2，θ3，…，θn]為優(yōu)化后相位，Λ=[ejθ1，ejθ2，…，ejθn]H，θ1可以被設為零，則僅需從N(N-1)/2個干涉相位中得到N-1個最佳相位，S為未知的真實相干矩陣，故用|T|來估計|S|，化簡得式(3)：

argmax {ΛH(W.*Φ)Λ}

(3)

式中，權重矩陣W為矩陣的絕對值部分。在求解式(3)時，每個像元都需對相干矩陣進行求逆運算從而獲得權重矩陣，導致其運算效率往往較低，為此采用相干矩陣加權法[13]有效避免矩陣求逆問題，見式(4)：

θML=argmax {ΛHTΛ}

(4)

用式(5)代替常規(guī)求解非線性系統(tǒng)的BFGS算法，迭代獲得優(yōu)化后的干涉相位：

(5)

式中，k為迭代步驟。最后對優(yōu)化后的相位進行質量指標評價，將時間相干值大于給定閾值的像元作為最終選取的DS像元，見式(6)：

1.2 PS像元選取

采用幅度離差指數(shù)法篩選PS像元[5]，見式(7)：

式中，σA和μA分別為每個像元的幅度序列平均值和標準差。幅度離差指數(shù)大于給定閾值的像元為初始PS候選點，再對初始候選點計算相干性，見式(8)：

式中，N為干涉對數(shù)量；φm，n為初始干涉相位；Δφm，n為誤差相位，包括視角誤差、大氣和軌道誤差等。最后通過迭代計算相干性，相干性較高的候選點為最終的PS像元。

1.3 像元形變解算

優(yōu)化后DS像元的相位與PS類似，包括形變、地形和大氣等相位[14]，故使用PS-InSAR解算方法，將兩類目標聯(lián)合后共同提取形變信息：先將像元構建Delaunay三角網(wǎng)分析相位，進行相位解纏，再在時間上進行高通濾波，空間上進行低通濾波去除大氣和軌道相位，得到最后的時序形變[15]。

2 研究區(qū)域與數(shù)據(jù)

研究區(qū)域地表存在植被、耕地、水域及少量建構筑物，還包括一條東西方向專用鐵路。開采的94101工作面距離鐵路最近約300m，工作面平均采厚約3.2m，煤層傾角20°，傾向方向向北，采寬180m，平均采深1020m，松散層厚度約260m。工作面開采時間2018年10月至2019年9月，走向方位角約270°，沿走向方向累計推進了約1100m。94101工作面北部為未采區(qū)域，上部為74101老采空區(qū)，故屬于重復采動，南部也為老采空區(qū)，其位置如圖2所示，老采空區(qū)停采時間較久，地表受其殘余沉降已不明顯。

圖2 94101工作面及周圍開采情況示意圖

研究數(shù)據(jù)采用了2018年11月至2019年11月間共30景Sentinel-1A衛(wèi)星數(shù)據(jù)，該衛(wèi)星由歐空局于2014年發(fā)射，搭載C波段合成孔徑雷達，影像的中心入射角為39.2°，衛(wèi)星軌道采用近極的太陽同步軌道，重訪周期為12d，干涉寬幅模式(Interferometric Wide Swath，IW)數(shù)據(jù)覆蓋范圍達20km×25km，方位向和距離向的分辨率分別為5m和20m，由于覆蓋范圍過大，將數(shù)據(jù)裁剪成約9km×6km覆蓋研究區(qū)域的范圍以便后續(xù)步驟，采用了美國宇航局航天飛機測圖任務(Shuttle Radar Topography Mission，SRTM)90m分辨率的DEM數(shù)據(jù)去除地形相位的影響。

3 實驗結果分析

以2019年4月28日的影像為主影像，按1.1和1.2節(jié)方法選取PS和DS像元，其中同質像元識別窗口采用了15×15，窗口內(nèi)同質像元數(shù)量大于20的像元為初始候選DS點，再采用最大似然法對原始干涉圖的相位迭代優(yōu)化，之后將時間相干閾值設為0.6篩選最終DS像元；PS像元的幅度離差閾值設為0.4。對實驗結果分別通過時間相干性和相干像元密度兩方面進行分析：

3.1 時間相干性

時間相干性可用來評價相位優(yōu)化的質量[8]，對初始候選的DS點按式(8)計算時間相干性Γ，如圖3所示，圖3(a)為最大似然估計相位優(yōu)化的時間相干性結果，圖3(b)為時間相干性頻率分布直方圖。

圖3 最大似然相位優(yōu)化時間相干性

由于建筑物密集地區(qū)存在大量穩(wěn)定的散射體，時間相干性接近于1；此外，農(nóng)田、裸露的土地等像元的相位優(yōu)化后，時間相干性也有良好的表現(xiàn)，其中大于閾值0.6的高相干像元數(shù)量達43358，約占總像元的14%。

3.2 相干像元密度

僅用幅度離差法選取PS像元和采用最大似然優(yōu)化干涉圖后選取DS像元的數(shù)量分別為10560和43358，說明與PS-InSAR方法相比，使用最大似然法優(yōu)化干涉圖相位使相干像元選取的數(shù)量大大提升，尤其在非人工建筑區(qū)域，由于SAR影像分辨率的限制，DS像元存在廣泛，提升效果明顯。去除重復選取的像元后將PS和DS聯(lián)合，計算具體像元密度見表1，其增幅高達374.9%。

表1 相干像元密度

3.3 變形區(qū)域地表平均形變速率

最后，對篩選出的高相干像元進行形變解算，獲得了研究區(qū)域衛(wèi)星視線方向的年平均形變速率。解算結果顯示，2018年11月至2019年11月間，由于94101工作面的開采，衛(wèi)星視線方向的地表最大年平均沉降速率大于100mm/a。工作面北部為裸地和農(nóng)田區(qū)域，具有較多分布式散射體目標，但永久散射體目標不足，故在PS-InSAR的解算結果中形變信息較少，但聯(lián)合PS與DS像元后，則能成功獲得該區(qū)域的形變信息；工作面的南部為老采空區(qū)積水區(qū)域，永久散射體與分布式散射體數(shù)量均不足，獲得的形變信息較少。

4 礦區(qū)鐵路動態(tài)變形規(guī)律分析

工作面向北約300m的鐵路受采動影響產(chǎn)生了變形，為了得到詳細的鐵路變形情況，確保鐵路日常安全運營，在形變解算結果中按鐵路沿線提取相干像元作進一步分析。基于影像分辨率等因素，點位間隔約為80m左右，共得到點30個，點號自東向西(與工作面回采方向一致)分別命為1～30。以2018年11月11日為起始時間，根據(jù)時間基線計算相對該起始時間的時序累積沉降值，并將衛(wèi)星視線方向的沉降轉換為垂直方向，鐵路沿線動態(tài)沉降情況如圖4所示，隨著工作面的推進，鐵路沉降量和沉降范圍逐漸增大，至最后一期，最大沉降量累積達271mm。

繪制最大沉降點的動態(tài)沉降曲線和沉降速度如圖5所示，在工作面開采期間，沉降速度不斷增加，最大為1.66mm/d，直到2018年9月工作面停采后，每月沉降量逐漸減小，沉降速度降低，趨勢減緩，其中10月1日至11月18日間沉降量小于10mm，沉降速度減少至0.4mm/d，鐵路沉降基本穩(wěn)定。

圖4 鐵路沿線動態(tài)沉降曲線

圖5 最大點沉降值及速度曲線

對于鐵路等延伸性構筑物，地表的不均勻下沉將進一步造成鐵路的傾斜與曲率變形[16，17]，影響鐵路行車安全，鐵路沿線的動態(tài)傾斜和曲率情況計算分別如圖6和圖7所示。從圖中可以看出，在工作面推進過程中，鐵路沿線傾斜和曲率變形逐漸增大，至2019年10月底逐漸趨于穩(wěn)定，傾斜在點12至點17附近(采空區(qū)中心)近似為零，在點6和點21附近(采空區(qū)邊界)達到最大，為0.82mm/m；曲率變形較小，均小于0.01mm/m2。

圖6 鐵路沿線動態(tài)傾斜曲線

圖7 鐵路沿線動態(tài)曲率曲線

為獲取一組鐵路移動變形預計的可靠參數(shù)，利用鐵路沿線11月18日的沉降數(shù)據(jù)，采用“礦區(qū)沉陷預測預報系統(tǒng)”軟件[18]求得概率積分預計參數(shù)，其中重復采動條件下的下沉系數(shù)q=0.79，主要影響角正切tanβ=1.63，拐點偏移距s=0，開采影響傳播角θ0=79°。

對于重復采動而言，一次重采下沉系數(shù)與初采下沉系數(shù)的關系如式(9)：

q復1=(1+a)q初

(9)

式中，a=0.15(堅硬巖層)或a=0.20(中硬巖層)。本區(qū)覆巖巖性中硬，按式(9)計算進一步獲得的初次采動條件下沉系數(shù)q初=0.66。

礦區(qū)鐵路行車速度較慢，因此對鐵路稍加維修可消除變形，使其安全運營，鐵路變形雖趨于穩(wěn)定，但未來還需開采鐵路附近的其它煤炭資源，因此仍必要持續(xù)監(jiān)測鐵路沿線及附近的變形情況，同時根據(jù)變形及時按照鐵路管理的相關規(guī)程提出鐵路綜合整治方案，建立完整的礦區(qū)鐵路沉降監(jiān)測系統(tǒng)，使地下煤炭資源合理生產(chǎn)的同時，地表環(huán)境和上覆建構筑物不受破壞。

5 結 論

本文基于DS-InSAR方法，對PS和DS像元的聯(lián)合求解形變展開研究，提高了時序InSAR技術在PS點數(shù)量不足的礦區(qū)的形變監(jiān)測效果，并分析了礦區(qū)鐵路的形變規(guī)律，主要成果如下：

1)基于最大似然估計的相位優(yōu)化方法能有效優(yōu)化原始干涉圖的相位，抑制相位噪聲，克服失相干影響，使高相干性像元數(shù)量與形變解算的質量相比PS-InSAR具有顯著提升，研究區(qū)域內(nèi)的相干像元密度增幅達374.9%。

2)受地下開采影響，工作面附近地表受持續(xù)形變，鐵路沿線最大沉降值達271mm，最大沉降速度1.66mm/d，橫向最大傾斜0.8mm/m。

3)根據(jù)監(jiān)測的鐵路沿線沉降結果求得最優(yōu)概率積分預計參數(shù)：q=0.79，tanβ=1.63，s=0，θ0=79°。本區(qū)屬一次重復采動，覆巖巖性中硬，按經(jīng)驗公式修正獲得初次采動條件下沉系數(shù)q初=0.66。研究所得的參數(shù)可以為后續(xù)開采沉陷預計提供數(shù)據(jù)支持。