基于DS-InSAR的張雙樓煤礦長時序地表形變監測方法

趙立峰 范洪冬，2 渠俊峰 李騰騰

（1.中國礦業大學環境與測繪學院，江蘇徐州221116；2.中國礦業大學自然資源部國土環境與災害監測重點實驗室，江蘇徐州221116；3.徐州市生態文明建設研究院，江蘇徐州221009）

長期高強度的地下開采導致礦區生態環境遭到破壞，引發了地表塌陷、山體滑坡和水土流失等各類地質災害，極大影響了礦區人民的生產生活［1-3］。傳統開采沉陷監測技術如全站儀、GNSS和水準測量等，雖然監測精度高，但存在勞動強度大、耗時長、工作效率低、點位易破壞、難以大范圍監測等不足。近年來，合成孔徑雷達差分干涉測量技術（D-InSAR）已在開采沉陷監測中得到諸多應用，然而該技術受時空失相干及大氣延遲等誤差的制約，監測精度受到很大限制。在此基礎上發展的諸多時序InSAR技術，如永久散射體干涉測量（PS-InSAR）［4］、小基線集（SBAS）［5］、相干目標分析方法［6］、斯坦福永久散射體方 法（StaMPS）［7］以 及 臨 時 相 干 點（TCP-InSAR）等［8-9］，克服了傳統單點變形測量技術和D-InSAR技術的不足，能夠實現低成本、高精度和大規模的地表監測，具有很大的優勢。

上述時序InSAR方法在具有較多強散射目標的城市區域應用較好，但在非城市區域、裸地、植被覆蓋區則難以達到有效的監測點位密度要求。為此，在常規時序InSAR技術的基礎上發展了融合分布式目標的時序地表形變監測方法，如SqueeSAR［10］、PDPSInSAR［11］，CAEInSAR［12］等，在考慮分布式散射體（Distributed Scatterer，DS）的高相干面狀目標的同時提高了其信噪比和測量點密度，已成為新的研究熱點。國內學者深入研究了DS-InSAR方法，其在填海區、山區滑坡、煤礦開采沉陷和水利設施等領域的變形監測中得到了成功應用。王明洲等［13］結合Kolmogorov-Smirnov（KS）檢驗和自適應非局域濾波，去除了相干目標的非平穩性信號，充分獲取了目標點的相位信息，得到了香港填海區域可靠的地表形變監測結果。但由于KS檢驗功效低，且在小樣本和高置信水平下容易引入許多異質點，使得難以實現同質點的高效提取，對后續時空濾波過程也存在影響。蔣彌等［14］提出了基于快速分布式目標探測的時序InSAR方法，該方法將同質點選取中的假設檢驗問題轉換為置信區間估計，大大提高了運算效率和檢驗功效，實現了油藏區的地表形變監測。王靜等［15］基于改進的干涉點目標分析（IPTA）方法，利用Anderson-Darling（AD）檢驗和相位三角算法，聯合PS和DS分析得到了貴州某鎮滑坡區的影響因素及變形機理，但由于植被密集和山區地形的復雜性，相位三角優化和相位解纏的準確性受到很大影響。張正佳等［16］利用分類信息和統計特性，基于區域增長策略識別和提取分布式目標，控制誤差傳播，增強了對煤礦區緩慢移動沉降的探測能力。劉友奉等［17］將DSInSAR方法應用到小浪底大壩的形變監測中，提高了InSAR在低相干區域的穩定性監測能力，與SBAS結果的交叉驗證表明了其監測結果的可靠性。但由于缺乏壩體結構、水位、蓄水量等資料，該研究未詳細分析壩體的變形機理。

與PS點相比，DS點相位容易受到失相干的影響且對垂直基線敏感，同時其在較大范圍非城市區域的計算效率也比較低。由于上述問題的存在且礦區地表形變范圍小、梯度大，失相干嚴重，沉降機理復雜［3］，目前DS-InSAR較少應用于礦區地表形變監測，尤其是在礦區長時間序列監測方面的研究涉及較少。為此，本研究提出了一種基于分布式目標InSAR的礦區長時序地表形變監測方法，以沛北礦區張雙樓煤礦為研究對象驗證了方法的可靠性。

1 研究區概況和數據來源

1.1 研究區概況

沛北礦區張雙樓煤礦（圖1）位于徐州市西北，在江蘇沛縣安國鎮境內，其地表屬黃泛沖積平原，地面標高+37～+39 m，地勢西高東低，地表水系密布，伴隨著采煤塌陷地的出現往往會形成積水坑，影響了該區域的生態環境。井田內含煤地層共3層，分別為石炭系上統太原組、二疊系下統山西組和二疊系下統石盒子組，目前開采煤層分別為山西組7、9煤。山西組煤層屬穩定中厚煤層，平均總厚為5.82 m。該礦井田邊界內地表大多被植被覆蓋，至今已有三十多年的開采歷史。長時間高強度的地下開采，導致礦區大面積地表沉降甚至塌陷現象時有發生，危及居民生命財產安全。

1.2 SAR數據集

本研究試驗數據為日本宇宙航空研究開發機構（JAXA）于2006年1月24日發射的L波段ALOS PALSAR衛星數據，其方位向像元尺寸約4.68 m，距離向像元尺寸約3.17 m，重訪周期46 d，入射角約38.72°。由于張雙樓煤礦植被覆蓋密集，且工作面一般布置在農田、荒地之下，相干性較低，而相控陣型L波段合成孔徑雷達（PALSAR）的波長為23.6 cm，較C波段雷達更容易穿透植被到達地表，有利于識別村莊和植被密集區域的采煤塌陷情況。

試驗采用2007年2月20日—2011年3月3日獲取的13景升軌影像，包括10景FBS HH極化圖像和3景FBD HH/HV極化圖像。干涉對垂直基線最長為4 307.24 m，最長時間基線為1 472 d。本研究采用90 m分辨率的SRTM（航天飛機雷達地形任務）數字高程模型（DEM），從干涉圖中去除地形的影響。

2 DS-InSAR方法原理

2.1 FaSHPS同質點選取

區別于傳統的假設檢驗方法通過判斷樣本數據的顯著性差異來識別同質像元，FaSHPS（Fast SHP Selection）算法［18］利用置信區間估計，通過簡單的邏輯運算判斷像元的同質點，大大提高了運算效率。

假設有N景SAR影像，根據中心極限定理，對于任一像元L，隨著樣本數N增加，振幅均值Aˉ()L逐漸靠近高斯分布，其區間估計可以表示為

確定平均強度的置信區間后，通過計算待估計像元時間維度上的平均振幅值是否落入目標像元對應的區間，判斷兩者是否屬于同質點。

2.2 特征值分解相位優化

時間序列SAR數據的像元相干矩陣的特征值對應不同強度的散射信號，因此可通過特征值分解來分離多元散射機制對應的特征信號［11，20］。

假設同質區域Ω包含NP個具有相似散射特性的相鄰像素，則相干矩陣為

相位優化之后，去除干涉圖中的地形相位并生成差分干涉圖，對解纏后的差分干涉圖進行濾波，在高相干點上建立觀測方程，采用奇異值分解估算地表形變相位，分離出大氣相位、DEM殘余相位和噪聲等，獲取研究區域視線向地表沉降時間序列。具體處理流程如圖2所示。

3 監測結果與分析

3.1 時序監測結果

根據圖2所示的處理流程，分別采用SBAS和DS-InSAR方法提取了研究區域地表時序沉降。SBAS和DS-InSAR兩種方法得到的張雙樓煤礦2007年2月—2011年3月監測時段內沿雷達視線方向的累計沉降值分布如圖3所示。SBAS方法共選取到88 728個相干點，累計最大沉降值為-451 mm；DSInSAR方法共選取到307 747個相干點，累計最大沉降值為-700 mm。分析圖3可知：兩種方法得到的監測結果較為準確地探測出了煤礦開采的沉陷位置，反映出了地表塌陷的影響范圍，具有較好的一致性；在沉降量級上小形變區域監測結果接近，差異較小，在開采沉陷區DS-InSAR監測結果明顯優于SBAS，能更好地反映張雙樓礦該時段煤礦開采引起的地表沉降情況。

3.2 監測結果分析與驗證

3.2.1 SBAS和DS-InSAR監測結果對比分析

由于張雙樓煤礦地物覆蓋類型多以裸地和植被為主，地物的后向散射性較弱，傳統SBAS方法選取的相干點的分布不足以反映出礦區具體沉降情況。與SBAS相比，融合了分布式目標的DS-InSAR方法在非城市區域能選取到更多的測量點目標，探測出沉降區的形變信息。由圖3可以看出，隨著工作面推進，研究時間段內，張雙樓煤礦共出現3處明顯的下沉盆地，分別為礦區內東側靠近邊界的沉降區Ⅰ、中部的沉降區Ⅱ以及西側邊緣的沉降區Ⅲ，各沉降區有明顯蔓延的趨勢。進一步分析可知：Ⅰ號沉降中心的沉降速率為-147 mm/a，Ⅱ號和Ⅲ號沉降中心的沉降速率分別為-174 mm/a和-51 mm/a，最大下沉點位于下沉區域Ⅱ。

為分析監測時段內SBAS和DS-InSAR兩種方法解算得到的時序累計沉降的差異，以SBAS方法的監測結果作為參照，分別選取了3個下沉盆地的沉降中心區和小形變區域的共同點進行分析，繪制了各點的時序沉降下沉曲線，如圖4所示。通過圖中參考點時序累計沉降對比可以看出，兩種方法得到的結果時間變化趨勢基本吻合。在小形變區域，SBAS和DS-InSAR的監測結果有高度的一致性，最大差值僅為6 mm，一定程度上說明了監測結果的可靠性；但在開采沉陷區，SBAS解算得到的沉降量級明顯小于DS-InSAR方法，最大差值達184 mm。原因可能是：①張雙樓礦采煤工作面位于農田和荒地之下，研究區中SBAS選擇的相干目標分布不均勻，多集中在建筑物和裸露的道路、陡坎區域，工作面上方有效監測點數量稀少，無法準確反映開采沉陷區變形情況；②稀疏且離散分布的SBAS相干目標點使得相位解纏時相鄰兩點間的解纏距離較大，增加了解纏誤差，而DS-InSAR在非城市區域仍能獲取到充足的測量點，且由于相位優化后噪聲得到了有效抑制，提高了形變區域干涉條紋的質量，故而使得相位解纏和形變解算更加精確。

為進一步對比分析兩種方法得到的各下沉盆地的沉降特征變化情況，分別于3個下沉盆地處構建了一條剖面線，并沿剖面線方向繪制了100 m緩沖區內的累計沉降分布圖。其中，由于SBAS方法在下沉盆地Ⅲ號沉降區中心所在區域未選取到相干點，在其下方偏離沉降中心一定距離處繪制了剖面線。如圖5所示，Ⅰ號、Ⅱ號和Ⅲ號沉降區的累計沉降量分別為-591、-700、-207 mm，沉降趨勢與地下采礦活動基本一致。與SBAS方法相比，融合分布式目標的DSInSAR方法有效提高了測量點空間分布范圍和密度，更能夠直觀地描述地下開采引起的礦區地表下沉盆地的變化特征。

3.2.2 監測結果驗證

由于缺乏研究時段內的水準數據，為分析分布式目標礦區地表沉降監測的可靠性，以SBAS監測結果為參照，選取與SBAS完全重合的分布式目標測量點，形成SBAS-DS點對，繪制了如圖6所示的形變結果相關圖及誤差分布圖。共選取63 867對同名點，由Pearson相關系數計算方式得到同名點對沉降值的整體相關性系數為0.97，但由于SBAS方法選取的相干點較少且分布較為集中，使得相位解纏時相鄰兩點間的距離較大，引入較大的解纏誤差，導致沉陷區SBAS結果明顯小于DS-InSAR，但在小形變區域兩種方法的監測結果相近，相關性高，表明兩種方法得到的監測結果有較好的一致性。

對得到的同名點對進行統計分析可知，其整體形變結果誤差均值接近于0。其中，誤差絕對值小于10 mm的同名點數為62 414，占所選同名點總數的97.8%。因此，DS-InSAR在顯著提高礦區內測量點密度的同時，其在小形變區域的監測結果與SBAS相對吻合，在開采沉陷區監測能力明顯優于SBAS，表明其在礦區時序地表沉降監測中有很大的應用潛力。

為進一步分析DS-InSAR監測的可靠性，對13景ALOS影像采用累積D-InSAR方法計算了時序地表形變。限于篇幅，以Ⅰ號沉降區為例，沿剖面線AA'方向繪制了5 m緩沖區內測點的累計沉降分布圖（圖7）。由圖7可知：DS-InSAR方法和D-InSAR方法沿剖面線方向沉降趨勢基本吻合，具有很好的一致性。但由于D-InSAR方法受到植被和噪聲的影響較大，離群值較多，其監測結果分布不如DS-InSAR方法穩定。

3.3 地表沉降與地下開采關聯分析

研究時段內，由于地下煤層的開采，張雙樓煤礦共出現3個下沉盆地，如圖3所示。本研究重點分析II號下沉區工作面開采引起的地表動態變形特征。

A工作面開采時間為2007年1月—2008年1月，其開采方向、范圍以及累計沉降分布如圖8（a）所示。雖然工作面上方測量點稀少，但仍可以發現隨著工作面推進，該時間段內地表下沉量沿工作面推進方向明顯增加，且最大下沉區域逐漸向下山方向偏移，最大下沉值為-209 mm。由于該區域采深約800 m，工作面寬度126 m，明顯未達到充分采動，且屬于深部開采，因此，導致地表形變量較小是合理的。

2008—2009年期間，有B和C兩個工作面在開采，開采時間分別為2008年1月—2008年12月和2008年6月—2008年12月，如圖8（b）所示。隨時間推移，B、C工作面開采的影響逐漸傳遞到地表，加劇了地表下沉。其中，由于B工作面煤層傾角較大，為29°，南部淺北部深，最大下沉點并未在工作面正上方，而是位于B工作面中部北側，其下沉值相對于2008年增加了約288 mm。此時，A、B兩工作面開采后，導致地下采空區范圍增大，地表沉降影響面積和量級也隨之增大。C工作面沿走向方向地表下沉逐漸增大，此期間相對下沉值達116 mm，但此時該工作面周圍缺少其它工作面開采，四周煤柱有支護作用，導致地表變形量較小。因此，總體上本研究獲取的地表沉降分布與期間工作面開采較為吻合。

2009年2—12月，隨著下組煤D工作面的開采，地表沉陷加劇，如圖8（c）所示。根據井下資料，D工作面煤層傾角約30°，且其位于原B工作面所在位置，受重復采動的影響，最大下沉區域逐漸向外蔓延，累計最大下沉值達-547 mm。

2010—2011年，E和F兩個工作面正在開采，開采時間分別為2010年1—12月和2011年1—6月，如圖8（d）所示，其中D與E工作面開采條件一致。隨著工作面E開采，D工作面處下沉逐漸趨緩，E工作面下沉速度逐漸增大，同時，到2011年F工作面的開采加劇了C工作面處殘余沉降的影響，使得下沉范圍擴大，最大下沉量達700 mm。

4 結 論

為提高礦區自然地表環境下采煤沉陷的監測精度，根據地面目標散射的特性，采用融合分布式目標的DS-InSAR方法，對2007—2011年采集的13景L波段ALOS影像進行了處理和分析，獲得了研究區域的沉降時間序列，得到如下結論：

（1）與PS-InSAR、SBAS等傳統時序InSAR方法相比，DS-InSAR方法以分布式目標為研究對象，通過相位優化提高干涉質量，經時序處理獲取地表形變，可有效提高礦區地表形變監測點的密度和監測精度。

（2）研究時間段內，利用DS-InSAR方法成功識別了張雙樓煤礦3處明顯的開采沉陷區，最大沉降約700 mm。結合井下開采資料分析表明，DS-InSAR方法探測的開采沉陷位置、影響范圍和發展趨勢較為準確，可為采煤塌陷地邊界確定、開采沉陷規律研究等提供參考。

（3）本研究僅利用L波段的ALOS-1數據驗證和分析了DS-InSAR方法監測礦區動態地表形變的效果，數據獲取時間較早且數據量較少，將來還需深入研究多源、多波段SAR數據的DS-InSAR處理方法，以及聯合PS目標和DS目標的礦區地表時序形變解算方法，進一步提高DS-InSAR監測礦區地表形變的精度和適用性。