刁鑫鵬，孫全帥，白志輝，吳 侃，楊 靜

（1.中國礦業(yè)大學(xué) 環(huán)境與測繪學(xué)院，江蘇 徐州 221116；2.中國礦業(yè)大學(xué) 江蘇省資源環(huán)境信息工程重點實驗室，江蘇 徐州 221116；3.冀中能源峰峰集團有限公司，河北 邯鄲 056000；4.寧夏回族自治區(qū)地球物理地球化學(xué)調(diào)查院，寧夏 銀川 750000）

0 引 言

煤炭是我國的基礎(chǔ)能源，長期以來為經(jīng)濟社會發(fā)展和國家能源的安全穩(wěn)定供應(yīng)提供了有力保障[1-3]。煤炭資源大規(guī)模開采的同時，對礦區(qū)環(huán)境亦造成了極大損害，衍生出一系列的地質(zhì)災(zāi)害與社會問題。減輕或防治各類災(zāi)害發(fā)生的關(guān)鍵在于明確開采沉陷規(guī)律，當(dāng)前對沉陷規(guī)律的研究已取得了一系列顯著成果[4-8]。然而，因復(fù)雜地形及地質(zhì)構(gòu)造的存在，使得地表移動和變形分布的正常規(guī)律被破壞，地表異常變形現(xiàn)象持續(xù)出現(xiàn)，且因此而引起的地表損害較常規(guī)情況更為嚴(yán)重[9-13]。

目前，國內(nèi)外針對地表采動異常損害及其變形規(guī)律的研究正逐漸增多，亦取得部分成果[14-16]。如余學(xué)義等[17]從礦區(qū)地質(zhì)采礦條件入手，揭示了張坡村出現(xiàn)大量地裂縫和建筑物異常損害的原因為開采誘發(fā)的斷層滑移；郭迅等[18]分析了建筑物差異性沉降的根本原因是構(gòu)造和采煤沉陷相互作用下局部區(qū)域巖體的“多米諾骨牌”反應(yīng)；楊雋等[19]、DIAO等[20]多角度分析了昔陽縣某村房屋建筑異常損害的原因為采動引起的坡體滑移；DONNELLY等[21]研究發(fā)現(xiàn)已穩(wěn)沉區(qū)域的部分?jǐn)鄬尤源嬖凇盎罨爆F(xiàn)象，并會引起地表的嚴(yán)重?fù)p害，但因缺乏實測數(shù)據(jù)，未能確定異常形變的持續(xù)時間。準(zhǔn)確認(rèn)識和理解地表移動過程的前提是對研究對象進行高精度、高可靠性的監(jiān)測；然而，地表異常形變及損害相對而言具有特殊性和不可預(yù)知性，容易被忽視。同時，由于傳統(tǒng)測量手段存在著測點稀疏、測量范圍小、監(jiān)測周期長等不足[22]，使得多維度、長時間序列的地表實測數(shù)據(jù)相對缺乏，導(dǎo)致目前對地表采動異常形變時、空演變特征的認(rèn)識尚不充分。隨著SAR衛(wèi)星數(shù)據(jù)的不斷豐富及數(shù)據(jù)處理手段的發(fā)展，InSAR技術(shù)獲取地表形變的精度已達毫米級；特別是，該技術(shù)可利用研究區(qū)的存檔影像對歷史形變進行反演，成為研究采動區(qū)地表異常形變規(guī)律和揭示損害原因的有效手段。

以河北省邯鄲市某城鎮(zhèn)的異常損害為研究對象，利用2019-07-14—2020-04-27期間的25景Sentinel-1A SAR影像數(shù)據(jù)，通過短基線集技術(shù)（SBAS-InSAR）反演了異常損害區(qū)地表形變的空間分布形態(tài)及其演變過程。結(jié)合研究區(qū)地質(zhì)資料，揭示了前述損害是開采誘發(fā)F15斷層“活化”的結(jié)果，并得出地表變形存在跨斷層傳遞的結(jié)論。研究成果表明，SBAS能夠?qū)崿F(xiàn)對地表異常形變的有效識別和提取，可為異常損害規(guī)律研究和損害原因揭示提供數(shù)據(jù)基礎(chǔ)，對沉陷控制理論的發(fā)展完善亦具有借鑒意義。

1 研究區(qū)及異常損害現(xiàn)象

1.1 研究區(qū)概況

研究區(qū)位于河北省邯鄲市，地理位置處于114°08′～114°11′E與36°36′～36°38′N，圖1展示了研究區(qū)域的井上下對照情況。其中，地面主要為農(nóng)田及城鎮(zhèn)建筑物；井下工作面開采運用綜采工藝，頂板管理方式為全部垮落法；煤層頂、底板巖性以粉砂巖為主。表1列出了研究區(qū)周邊所開采工作面的具體信息。

表1 研究區(qū)周邊各工作面開采信息

圖1 研究區(qū)井上下對照情況Fig.1 Surface-underground contrast plan of study area

1.2 損害特征及沉陷預(yù)計分析

1）損害特征。該城鎮(zhèn)于2019年10月開始出現(xiàn)損害現(xiàn)象，且隨時間推移損害加重。經(jīng)調(diào)查，至2020年5月，建筑墻體裂縫達30 mm，地表張開型裂縫及臺階裂縫寬度達10～20 mm。受損害區(qū)域呈條帶分布，具體分布位置如圖1中紅色曲線所示。

2）沉陷預(yù)計分析。為確定常規(guī)情況下工作面開采對該鎮(zhèn)的影響，采用概率積分法進行地表沉陷預(yù)計。預(yù)計參數(shù)基于該礦地表移動觀測站實測數(shù)據(jù)并參照相關(guān)規(guī)范指南確定，具體為：下沉系數(shù)q=0.78，主要影響角正切tanβ=1.70，開采影響傳播角θ=90°-0.4α，水平移動系數(shù)b=0.30，拐點偏移距S=0。圖2顯示了基于概率積分法預(yù)計的地表下沉等值線。

經(jīng)預(yù)計計算得到，該城鎮(zhèn)范圍內(nèi)地表下沉量、傾斜和水平變形的最大值分別為：200 mm、3.0 mm/m和2.5 mm/m；井下工作面開采會造成該區(qū)域房屋建筑一定程度上的損害。但基于《建筑物、水體、鐵路及主要井巷煤柱留設(shè)與壓煤開采規(guī)范》中的變形控制指標(biāo)，地下開采所致建筑物的損害程度不會超過II級，這與實際損害情況不符。此外，井下工作面開采的直接影響邊界（圖2中綠色曲線），并未波及至現(xiàn)場調(diào)研的嚴(yán)重?fù)p害區(qū)，兩者間水平投影距離超130 m。表明該城鎮(zhèn)范圍內(nèi)房屋建筑的損害屬異常損害。

2 SBAS原理方法與數(shù)據(jù)準(zhǔn)備

為進一步分析地表變形的演變過程，準(zhǔn)確揭示異常損害的原因，利用InSAR技術(shù)可通過存檔影像回溯地表歷史形變的特性，采用SBAS時序分析方法對研究區(qū)建筑受損害前、后的地表變形情況進行反演。

2.1 原理與方法

SBAS技術(shù)的提出是為削弱時空失相干影響，提高InSAR地表形變解譯精度。BERARDINO、LANARI等[23-24]先后探討了該技術(shù)在地表形變監(jiān)測中的應(yīng)用，而后HOOPER等[25]又對處理過程中高相干點的選取方法進行了優(yōu)化。

與常規(guī)差分干涉測量相比，SBAS技術(shù)的基本思想是：通過設(shè)置時、空基線條件，將同一研究區(qū)的多期SAR影像分組成若干個短基線集合；利用最小二乘準(zhǔn)則，獲取每個短基線集合的地表形變時間序列；通過奇異值分解方法再將分組后的短基線集進行聯(lián)合求解；進而得到研究區(qū)整個觀測時間段內(nèi)時序的地表形變信息[26]。具體[27-28]可表述為：假定同一研究區(qū)成像時間為t0，t1，t2，…，tn-1的N幅SAR影像，在短基線距的條件下可形成M幅差分干涉圖，且滿足N/2≤M≤N（N-1）/2。對于干涉圖i，在去除平地相位后，干涉圖中任意位置（x，y）的干涉相位δφ（x,y）可表示為

δφi（x,y）=δφi（tB,x,y）-δφi（tA,x,y）≈

δφdef,i（x,y）+δφε,i（x,y）+δφα,i（x,y）+δφn,i（x,y）

（1）

其中，i為干涉圖序號；tA、tB為干涉圖i所對應(yīng)SAR影像的獲取時間；δφdef為tA～tB時間段內(nèi)地表

視線向形變相位；δφε為目標(biāo)區(qū)地形相位；δφα、δφn分別為大氣相位和噪聲相位。其中前3項可以表示為：

（2）

式中，di（tB,x,y）和di（tA,x,y）分別為tB、tA時刻地表沿雷達視線向的變形量；λ為雷達波長；R為斜距；B⊥為垂直基線；Δz為DEM高程差；θ為入射角。若不同干涉圖間的形變速率為vk,k+1，則tA～tB間的累積形變量可表達為：

（3）

對M幅干涉圖進行三維時空相位解纏和地理編碼后，便可求出不同SAR影像獲取時間的形變速率。

2.2 數(shù)據(jù)準(zhǔn)備

用于研究區(qū)形變信息提取的SAR數(shù)據(jù)為2019-07-14—2020-04-27期間的25景Sentinel-1A影像，C波段，地面分辨率為5 m×20 m，相關(guān)參數(shù)如下：

成像時間 2019-07-14、2019-07-26、2019-08-07、2019-08-19

2019-08-31、2019-09-12、2019-09-24、2019-10-06

2019-10-18、2019-10-30、2019-11-11、2019-11-23

2019-12-05、2019-12-17、2019-12-29、2020-01-10

2020-01-22、2020-02-03、2020-02-15、2020-02-27

2020-03-10、2020-03-22、2020-04-03、2020-04-15

2020-04-27

數(shù)據(jù)類型 IW

像元大?。╩×m） 2.3×14.1

Sentinel-1影像的軌道參數(shù)需根據(jù)影像的成像時間下載，歐空局提供了哨兵影像的Precise Orbit Ephemerides（POD精密定軌星歷數(shù)據(jù)）。此外，為消除地形相位，外部DEM采用的是美國宇航局SRTM3數(shù)據(jù)。

3 形變解譯與特征分析

3.1 InSAR形變解譯

SBAS時序分析主要包括：差分干涉對生成、點目標(biāo)選取、差分干涉圖解纏，以及時間形變序列獲取等幾個步驟，詳細(xì)的處理流程如圖3所示。

時序分析過程中相關(guān)參數(shù)設(shè)置如下：①連接圖生成時，為充分利用影像數(shù)據(jù)，時間基線設(shè)置為50 d，空間基線設(shè)置為臨界基線的10%；共連接成90個像對，其中最長時間基線為48 d，最長空間基線為175.56 m（平均59.78 m）。②干涉處理過程中距離向與方位向的視數(shù)比設(shè)置為4∶1，影像配準(zhǔn)利用精密軌道數(shù)據(jù)采用強度互相關(guān)算法，濾波方法選用自適應(yīng)濾波法，解纏方法采用最小費用流法（相干性閾值設(shè)置為0.3）。③軌道精煉和重去平處理過程中，地面控制點選擇相干性高、相位好的點，并且所選控制點位置的形變需為0。④第一次解譯反演選用較為穩(wěn)定的線性模型。其余參數(shù)均選用處理Sentinel-1數(shù)據(jù)的推薦參數(shù)。

圖3 SBAS地表形變解譯流程Fig.3 Flow of surface deformation with SBAS method

圖4為最終解譯的地表垂直位移的分布情況，顯示了不同時間段內(nèi)研究區(qū)地表變形的空間分布形態(tài)及其隨時間的演變過程。

3.2 InSAR解譯可靠性評定

為驗證InSAR解譯結(jié)果的可靠程度，收集了48個地面觀測點在2020-01-12—2020-05-03期間的水準(zhǔn)測量數(shù)據(jù)，水準(zhǔn)點位的具體分布位置如圖2所示。因水準(zhǔn)測量時間與InSAR解譯時間段并不完全一致，為降低因時間差異而造成的對比誤差，用于參與對比的InSAR數(shù)據(jù)選用時間段最為接近的2020-01-10—2020-04-27期間的解譯結(jié)果。圖5為2類數(shù)據(jù)各監(jiān)測點位的變形曲線及對比情況。

基于圖5中地表水準(zhǔn)實測數(shù)據(jù)，可分析得出：①對于Q測線，Q13點以東地表變形不明顯，Q13與Q14點間地表變形存在突變現(xiàn)象，Q14點以西的地面沉降先減小后又逐漸增加。②對于P測線，P5點以東地表變形不明顯，P5與P7點間地表變形存在突變現(xiàn)象，P7至P18測點地面沉降有減小趨勢。③地表變形突變點位連線，與城鎮(zhèn)嚴(yán)重受損區(qū)位置具有一致性。

經(jīng)統(tǒng)計得到，對于上述48個監(jiān)測點位，水準(zhǔn)測量與InSAR解譯的平均差值為+5.32 mm，均方根差為11.49 mm。兩者監(jiān)測結(jié)果存在差異的主要原因有：①觀測時間段不同（水準(zhǔn)測量較InSAR解譯的時間間隔長），導(dǎo)致上述對比數(shù)據(jù)中水準(zhǔn)測量變形值較InSAR解譯結(jié)果偏大；②InSAR解譯過程中，因控制點選取及大氣效應(yīng)去除等因素導(dǎo)致解譯結(jié)果存在誤差。

雖然2類數(shù)據(jù)的監(jiān)測結(jié)果間存在一定偏差，但平均差值占最大下沉值的比例僅為5%；同時，2類監(jiān)測數(shù)據(jù)的變形趨勢具有一致性。因此，我們可以基于InSAR解譯結(jié)果分析地表變形的時、空演變特征。

3.3 InSAR解譯形變特征分析

圖5中的水準(zhǔn)實測曲線一定程度上展示了采動地表的變形情況，但因時間和空間分辨率低的缺陷，對于地表突變出現(xiàn)的時機、變形的動態(tài)演變過程和空間分布特征難以準(zhǔn)確表達，而基于SBAS形變解譯結(jié)果可對上述特征做進一步分析。

圖4a顯示，2019-07-14—2019-09-12期間，周邊區(qū)域地下工作面的開采并未引起城鎮(zhèn)相關(guān)位置出現(xiàn)明顯的變形現(xiàn)象。圖4b顯示，至2019-09-24，受損區(qū)域地表的零星位置開始出現(xiàn)變形，量值約為10 mm。此時，618工作面開采已接近結(jié)束，158工作面已進入村莊煤柱進行開采，但建筑物尚未出現(xiàn)明顯損害。表明周期性的InSAR形變監(jiān)測可實現(xiàn)對地表異常損害的預(yù)警。

圖4c—圖4f顯示，自2019-9-24日起（158工作面推進距離約130 m），城鎮(zhèn)相關(guān)位置已經(jīng)開始出現(xiàn)明顯的地表變形現(xiàn)象，且地表變形發(fā)育的停止位置與鎮(zhèn)內(nèi)受損嚴(yán)重房屋的分布位置具有很好的一致性。圖4g—圖4i顯示，在158、618、619工作面作面開采引起常規(guī)的地表沉陷以外，地表變形有向東擴展的趨勢，且地表變形空間發(fā)育的終止位置近似成直線分布，而該終止位置恰為現(xiàn)場調(diào)研確定的城鎮(zhèn)異常損害的位置。

自圖4b顯示2019-9-24受損區(qū)域地表的零星位置開始出現(xiàn)變形后，圖4g、圖4h及圖4i進一步顯示上述位置地表的變形量逐漸增大、變形范圍逐漸增加，但與工作面開采引起的主要變形區(qū)域在空間上相對獨立。圖4k、圖4l亦表明，地下工作面開采在引起地表的常規(guī)變形以外，在沉陷盆地以東存在較為明顯的地表變形，且該異常形變與常規(guī)沉陷盆地在空間分布上具有較為明顯的獨立性。圖4m及圖4n則展示出，隨地下開采范圍的不斷擴展，上述兩獨立的沉陷區(qū)域逐漸發(fā)育為一體，發(fā)育的終止位置城鎮(zhèn)異常損害的分布具有一致性。

上述地表變形的空間分布形態(tài)及特征，一定程度上符合斷層賦存環(huán)境下開采沉陷的規(guī)律。

4 損害原因揭示

綜合地表實測數(shù)據(jù)與InSAR動態(tài)解譯結(jié)果的分析情況，研究區(qū)地表異常形變與損害大概率與某斷層的“活化”有關(guān)。此外，現(xiàn)場調(diào)研發(fā)現(xiàn)損害區(qū)地裂縫兩側(cè)存在較為明顯落差（圖1b），符合斷層露頭處地表損壞特征。

井上下對照情況（圖2）顯示，研究區(qū)周邊發(fā)育有3條較為明顯的斷層，分別為：F49（H=140 m，∠65°）、F13（H=100 m，∠70°）和F15（H>600 m，∠70°），圖6為8號勘探線的地質(zhì)剖面圖。

圖6 勘探線地質(zhì)剖面Fig.6 Geological section of exploration line

由圖6可得，地下工作面的開采容易擾動到F13斷層，引起斷層的“活化”；以45°作為松散層移動角所劃定的F13斷層露頭位置并未出現(xiàn)在城鎮(zhèn)屋建筑的受損區(qū)域，距離損害區(qū)約306 m。若F13斷層露頭出現(xiàn)在建筑物受損區(qū)域，則松散層移動角約為26°，這與常規(guī)情況下松散層移動角的取值范圍不符；即工作面開采引起F13斷層“活化”，而直接導(dǎo)致建筑物受損的可能性較低。

當(dāng)以56°作為松散層移動角時，F(xiàn)15斷層的露頭位于城鎮(zhèn)的異常受損區(qū)，且移動角的選擇滿足常規(guī)取值要求；此外，F(xiàn)15斷層落差大于600 m，屬特大型斷層，相對而言更易受到擾動。因此，結(jié)合前節(jié)地表實測數(shù)據(jù)及InSAR形變反演的分析結(jié)果，認(rèn)為F15斷層受到擾動而發(fā)生“活化”是研究區(qū)城鎮(zhèn)建筑發(fā)生異常損害的主要原因?；谏鲜鼋Y(jié)論，可進一步得到斷層賦存環(huán)境下，采動影響范圍會遠(yuǎn)超常規(guī)沉陷理論劃定的影響邊界，地表變形存在跨斷層傳遞現(xiàn)象。

5 結(jié) 論

1）現(xiàn)場調(diào)研及水準(zhǔn)實測顯示，研究區(qū)受損害位置及程度超出常規(guī)沉陷理論劃定的采動影響邊界和級別，地表變形存在突變現(xiàn)象，建（構(gòu)）筑物損害屬異常損害。

2）通過InSAR解釋結(jié)果與水準(zhǔn)數(shù)據(jù)的對比得到，兩者監(jiān)測地表變形的趨勢具有一致性，監(jiān)測變形的平均差值為5.32 mm，證明SBAS形變解譯結(jié)果具有可靠性，時序分析技術(shù)能夠有效識別并提取采動區(qū)地表的異常變形。

3）InSAR解譯地表變形空間分布的最終形態(tài)特征，一定程度上符合斷層存在時開采沉陷的規(guī)律；但從地表變形的動態(tài)發(fā)育過程看，異常形變初始發(fā)育位置與常規(guī)沉陷盆地在空間上相對獨立，而后才逐漸發(fā)育為一體；進一步揭示了斷層受開采擾動的過程。

4）地表異常損害為開采誘發(fā)F15斷層“活化”的結(jié)果，斷層賦存環(huán)境下采動影響范圍會遠(yuǎn)超以移動角劃定的影響邊界，且地表變形存在跨斷層傳遞現(xiàn)象。研究成果對沉陷控制理論的發(fā)展完善具有借鑒意義。