基于DS-InSAR的張雙樓煤礦長(zhǎng)時(shí)序地表形變監(jiān)測(cè)方法

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

（1.中國(guó)礦業(yè)大學(xué)環(huán)境與測(cè)繪學(xué)院，江蘇徐州221116；2.中國(guó)礦業(yè)大學(xué)自然資源部國(guó)土環(huán)境與災(zāi)害監(jiān)測(cè)重點(diǎn)實(shí)驗(yàn)室，江蘇徐州221116；3.徐州市生態(tài)文明建設(shè)研究院，江蘇徐州221009）

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

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

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

1 研究區(qū)概況和數(shù)據(jù)來(lái)源

1.1 研究區(qū)概況

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

1.2 SAR數(shù)據(jù)集

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

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

2 DS-InSAR方法原理

2.1 FaSHPS同質(zhì)點(diǎn)選取

區(qū)別于傳統(tǒng)的假設(shè)檢驗(yàn)方法通過(guò)判斷樣本數(shù)據(jù)的顯著性差異來(lái)識(shí)別同質(zhì)像元，F(xiàn)aSHPS（Fast SHP Selection）算法［18］利用置信區(qū)間估計(jì)，通過(guò)簡(jiǎn)單的邏輯運(yùn)算判斷像元的同質(zhì)點(diǎn)，大大提高了運(yùn)算效率。

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

確定平均強(qiáng)度的置信區(qū)間后，通過(guò)計(jì)算待估計(jì)像元時(shí)間維度上的平均振幅值是否落入目標(biāo)像元對(duì)應(yīng)的區(qū)間，判斷兩者是否屬于同質(zhì)點(diǎn)。

2.2 特征值分解相位優(yōu)化

時(shí)間序列SAR數(shù)據(jù)的像元相干矩陣的特征值對(duì)應(yīng)不同強(qiáng)度的散射信號(hào)，因此可通過(guò)特征值分解來(lái)分離多元散射機(jī)制對(duì)應(yīng)的特征信號(hào)［11，20］。

假設(shè)同質(zhì)區(qū)域Ω包含NP個(gè)具有相似散射特性的相鄰像素，則相干矩陣為

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

3 監(jiān)測(cè)結(jié)果與分析

3.1 時(shí)序監(jiān)測(cè)結(jié)果

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

3.2 監(jiān)測(cè)結(jié)果分析與驗(yàn)證

3.2.1 SBAS和DS-InSAR監(jiān)測(cè)結(jié)果對(duì)比分析

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

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

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

3.2.2 監(jiān)測(cè)結(jié)果驗(yàn)證

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

對(duì)得到的同名點(diǎn)對(duì)進(jìn)行統(tǒng)計(jì)分析可知，其整體形變結(jié)果誤差均值接近于0。其中，誤差絕對(duì)值小于10 mm的同名點(diǎn)數(shù)為62 414，占所選同名點(diǎn)總數(shù)的97.8%。因此，DS-InSAR在顯著提高礦區(qū)內(nèi)測(cè)量點(diǎn)密度的同時(shí)，其在小形變區(qū)域的監(jiān)測(cè)結(jié)果與SBAS相對(duì)吻合，在開(kāi)采沉陷區(qū)監(jiān)測(cè)能力明顯優(yōu)于SBAS，表明其在礦區(qū)時(shí)序地表沉降監(jiān)測(cè)中有很大的應(yīng)用潛力。

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

3.3 地表沉降與地下開(kāi)采關(guān)聯(lián)分析

研究時(shí)段內(nèi)，由于地下煤層的開(kāi)采，張雙樓煤礦共出現(xiàn)3個(gè)下沉盆地，如圖3所示。本研究重點(diǎn)分析II號(hào)下沉區(qū)工作面開(kāi)采引起的地表動(dòng)態(tài)變形特征。

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

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

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

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

4 結(jié) 論

為提高礦區(qū)自然地表環(huán)境下采煤沉陷的監(jiān)測(cè)精度，根據(jù)地面目標(biāo)散射的特性，采用融合分布式目標(biāo)的DS-InSAR方法，對(duì)2007—2011年采集的13景L波段ALOS影像進(jìn)行了處理和分析，獲得了研究區(qū)域的沉降時(shí)間序列，得到如下結(jié)論：

（1）與PS-InSAR、SBAS等傳統(tǒng)時(shí)序InSAR方法相比，DS-InSAR方法以分布式目標(biāo)為研究對(duì)象，通過(guò)相位優(yōu)化提高干涉質(zhì)量，經(jīng)時(shí)序處理獲取地表形變，可有效提高礦區(qū)地表形變監(jiān)測(cè)點(diǎn)的密度和監(jiān)測(cè)精度。

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

（3）本研究?jī)H利用L波段的ALOS-1數(shù)據(jù)驗(yàn)證和分析了DS-InSAR方法監(jiān)測(cè)礦區(qū)動(dòng)態(tài)地表形變的效果，數(shù)據(jù)獲取時(shí)間較早且數(shù)據(jù)量較少，將來(lái)還需深入研究多源、多波段SAR數(shù)據(jù)的DS-InSAR處理方法，以及聯(lián)合PS目標(biāo)和DS目標(biāo)的礦區(qū)地表時(shí)序形變解算方法，進(jìn)一步提高DS-InSAR監(jiān)測(cè)礦區(qū)地表形變的精度和適用性。