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

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

0 引 言

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

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

以河北省邯鄲市某城鎮(zhèn)的異常損害為研究對(duì)象，利用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)對(duì)地表異常形變的有效識(shí)別和提取，可為異常損害規(guī)律研究和損害原因揭示提供數(shù)據(jù)基礎(chǔ)，對(duì)沉陷控制理論的發(fā)展完善亦具有借鑒意義。

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

1.1 研究區(qū)概況

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

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

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

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

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

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

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

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

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

2.1 原理與方法

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

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

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

（2）

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

（3）

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

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ù)如下：

成像時(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ù)影像的成像時(shí)間下載，歐空局提供了哨兵影像的Precise Orbit Ephemerides（POD精密定軌星歷數(shù)據(jù)）。此外，為消除地形相位，外部DEM采用的是美國(guó)宇航局SRTM3數(shù)據(jù)。

3 形變解譯與特征分析

3.1 InSAR形變解譯

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

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

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

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

3.2 InSAR解譯可靠性評(píng)定

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

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

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

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

3.3 InSAR解譯形變特征分析

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

圖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。此時(shí)，618工作面開采已接近結(jié)束，158工作面已進(jìn)入村莊煤柱進(jìn)行開采，但建筑物尚未出現(xiàn)明顯損害。表明周期性的InSAR形變監(jiān)測(cè)可實(shí)現(xiàn)對(duì)地表異常損害的預(yù)警。

圖4c—圖4f顯示，自2019-9-24日起（158工作面推進(jìn)距離約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ī)的地表沉陷以外，地表變形有向東擴(kuò)展的趨勢(shì)，且地表變形空間發(fā)育的終止位置近似成直線分布，而該終止位置恰為現(xiàn)場(chǎng)調(diào)研確定的城鎮(zhèn)異常損害的位置。

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

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

4 損害原因揭示

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

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

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

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

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

5 結(jié) 論

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

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

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

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