熊 威，孫志杰，郝曉陽(yáng)

（1.山西省交通科技研發(fā)有限公司 黃土地區(qū)公路建設(shè)與養(yǎng)護(hù)技術(shù)交通行業(yè)重點(diǎn)實(shí)驗(yàn)室，山西 太原 030032；2.山西省煤炭地質(zhì)物探測(cè)繪院 資源環(huán)境與災(zāi)害監(jiān)測(cè)山西省重點(diǎn)實(shí)驗(yàn)室，山西 晉中 030053）

山西省作為傳統(tǒng)的煤炭資源大省，為國(guó)家的發(fā)展建設(shè)提供了源源不斷的能源供應(yīng)，同時(shí)也形成了近5 000 km2的采空區(qū)，其中沉陷區(qū)約3 000 km2，給采空區(qū)周邊居民生活、公路運(yùn)營(yíng)、工程建設(shè)等造成了較大的安全隱患。為了對(duì)采空區(qū)的影響范圍和發(fā)育程度進(jìn)行有效評(píng)估，有必要采取一種高效便捷的觀測(cè)手段對(duì)采空區(qū)進(jìn)行識(shí)別和監(jiān)測(cè)。傳統(tǒng)地表沉陷監(jiān)測(cè)方法主要有GNSS、水準(zhǔn)儀、全站儀等[1-3]，它們均是以點(diǎn)為單元進(jìn)行測(cè)量，效率低下，無(wú)法覆蓋全區(qū)，需要耗費(fèi)大量的人力財(cái)力。合成孔徑雷達(dá)干涉測(cè)量（interferometric synthetic aperture radar，InSAR）技術(shù)作為一項(xiàng)新型空間對(duì)地觀測(cè)技術(shù)，具有大面積、全天候、連續(xù)、高精度、高空間分辨率等優(yōu)點(diǎn)，由此衍生而來(lái)的差分干涉測(cè)量（D-InSAR）技術(shù)在地震同震形變[4]、城市地面沉降[5]、礦區(qū)地表沉陷[6]等方面做了廣泛的研究，取得了一定的成效。由于D-InSAR 容易受到大氣效應(yīng)、DEM誤差、基線誤差等因素的影響，其中，大氣效應(yīng)對(duì)雷達(dá)信號(hào)傳播產(chǎn)生的延遲，在中緯度地區(qū)一年內(nèi)的影響區(qū)間是20 cm，是最主要的誤差源[7]。為了消減各項(xiàng)誤差對(duì)監(jiān)測(cè)結(jié)果的影響，在此基礎(chǔ)上發(fā)展了多時(shí)相InSAR 技術(shù)（MT-InSAR）來(lái)進(jìn)行高精度形變監(jiān)測(cè)，永久散射體（PS-InSAR）和小基線集（SBAS-In-SAR）是其中最具代表性且廣泛使用的兩種方法，二者均具有毫米級(jí)精度，且能回溯地表歷史形變過(guò)程，在采空區(qū)沉陷監(jiān)測(cè)與成因分析[8-9]和形變預(yù)測(cè)[10-11]方面發(fā)揮了重要作用。PS-InSAR 方法基于一景主影像與其余輔影像配準(zhǔn)后差分，獲取PS點(diǎn)的位移時(shí)間序列，適用于相干性較強(qiáng)的城鎮(zhèn)地區(qū)，且對(duì)數(shù)據(jù)量要求較高，SBAS-InSAR技術(shù)克服了PS-InSAR因選取一幅影像作為公共主影像而引起的部分干涉圖相干性較差的不足，同時(shí)降低了對(duì)SAR數(shù)據(jù)的需求量，運(yùn)算效率較高，適用性更廣。

本文針對(duì)山西省陽(yáng)城地區(qū)煤礦采空區(qū)現(xiàn)狀，利用SBAS-InSAR技術(shù)對(duì)2019年1月至2020年7月間共25景Sentinel-1A影像進(jìn)行處理，根據(jù)時(shí)間序列形變結(jié)果分析黃土地區(qū)采空區(qū)沉陷的動(dòng)態(tài)變化過(guò)程，以及引起沉陷變化的直接因素，為黃土地區(qū)采空區(qū)的特征識(shí)別及危險(xiǎn)性分析提供技術(shù)支持。

1 SBAS-InSAR技術(shù)

當(dāng)M≥N時(shí)，通過(guò)最小二乘法即可求解ν；當(dāng)M＜N時(shí)，B為秩虧矩陣，可利用SVD方法對(duì)矩陣B進(jìn)行分解，進(jìn)而求得ν的最小范數(shù)二乘解。最后對(duì)ν在時(shí)間維度上積分即可得到相應(yīng)時(shí)刻的累積形變量。

2 測(cè)區(qū)概況與數(shù)據(jù)處理

2.1 測(cè)區(qū)概況

測(cè)區(qū)位于山西陽(yáng)城縣北部，晉運(yùn)高速、濟(jì)陽(yáng)高速在這里交匯（圖1），其中在晉運(yùn)高速陽(yáng)城段兩側(cè)分布有大大小小幾十個(gè)采煤區(qū)，開(kāi)采時(shí)間從20 世紀(jì)80 年代至今，高速南側(cè)大部分開(kāi)采時(shí)間早于2010年，北側(cè)大多為近年來(lái)形成的采空區(qū)，部分礦區(qū)甚至下穿高速。該地區(qū)主要為黃土丘陵地貌，地層從上往下主要為第四系全新統(tǒng)人工堆積層（Q4 ml）素填土，全新統(tǒng)沖洪積層（Q4 al+pl）粉質(zhì)黏土、碎石土，第四系上更新統(tǒng)坡洪積層（Q3 dl+pl）黏質(zhì)黃土和二疊系上統(tǒng)（P2）泥巖、砂巖。

圖1 測(cè)區(qū)地理位置圖

2.2 數(shù)據(jù)處理

本文利用2019-01-08—2020-07-01 的歐洲空間局Sentinel-1A 衛(wèi)星數(shù)據(jù)對(duì)采空區(qū)的地表形變進(jìn)行監(jiān)測(cè)，數(shù)據(jù)詳情如表1 所示，并配合使用與影像數(shù)據(jù)相匹配的精密軌道數(shù)據(jù)，其位置精度優(yōu)于5 cm。外部DEM 使用的是美國(guó)宇航局提供的SRTM-1，其分辨率為30 m，高程精度約為10 m。上述數(shù)據(jù)的使用能最大限度的減少DEM誤差和軌道誤差對(duì)干涉圖及解纏結(jié)果的影響。

表1 影像信息表

本文利用GMTSAR 軟件中SBAS 模塊對(duì)溝東隧道上覆邊坡的SAR數(shù)據(jù)進(jìn)行處理，該軟件具有開(kāi)源、自動(dòng)處理、易于繪圖等優(yōu)點(diǎn)[13]，在震后大范圍地表形變監(jiān)測(cè)方面取得了較好的效果[14]。由于測(cè)區(qū)位于丘陵地帶，季節(jié)性的灌木叢生，為了減少時(shí)空失相干的影響，增加相干點(diǎn)數(shù)，設(shè)置時(shí)間基線閾值為180 d，空間基線閾值為最大基線長(zhǎng)度的30%，組合形成的影像連接圖如圖2 所示，共形成163 對(duì)有效干涉像對(duì)，每幅影像至少與3 幅以上的影像產(chǎn)生干涉，網(wǎng)型健壯性較好。相干性閾值設(shè)置為0.1。在經(jīng)過(guò)影像裁剪、配準(zhǔn)、干涉圖生成、去地形相位、相位解纏、形變反演、地理編碼等一系列步驟之后最終得到累積形變結(jié)果如圖3 所示。其中距離向與方位向多視比為5∶1，成圖分辨率為20 m；解纏方法為經(jīng)典的最小費(fèi)用流法，解纏分解等級(jí)設(shè)置為1；最終選取的有效控制點(diǎn)數(shù)量為30個(gè)，均遠(yuǎn)離采空區(qū)且均勻分布于全圖；初始反演模型選擇較為穩(wěn)定的線性模型；軌道精煉選擇3 次多項(xiàng)式擬合方法。

圖2 影像連接圖

圖3 累積形變量圖

3 結(jié)果驗(yàn)證與精度評(píng)價(jià)

為了驗(yàn)證SBAS-InSAR技術(shù)監(jiān)測(cè)結(jié)果的精度，利用前期布設(shè)在6#采空區(qū)周邊的水準(zhǔn)線路進(jìn)行對(duì)比分析。該條水準(zhǔn)線路沿東西向的鄉(xiāng)村公路布設(shè)，包含1 個(gè)基準(zhǔn)點(diǎn)KZ1 和12 個(gè)監(jiān)測(cè)點(diǎn)（JC1～JC12），基準(zhǔn)點(diǎn)位于曲堤村，點(diǎn)位穩(wěn)定。線路總長(zhǎng)3.3 km，平均點(diǎn)位間距0.3 km，按照二等精密水準(zhǔn)的精度要求測(cè)設(shè)，監(jiān)測(cè)頻率1次∕月，均在每月上旬完成。由于InSAR技術(shù)得到的是衛(wèi)星視線向上的形變，水準(zhǔn)監(jiān)測(cè)得到的是垂直方向形變結(jié)果，兩者之間有如下轉(zhuǎn)換關(guān)系：

式中，θ為衛(wèi)星視線向（LOS）與地表法向的夾角；dH為精密水準(zhǔn)得到的地表垂向監(jiān)測(cè)結(jié)果。將水準(zhǔn)監(jiān)測(cè)結(jié)果轉(zhuǎn)換到衛(wèi)星視線向，并根據(jù)基準(zhǔn)點(diǎn)沉降位移值做校正，再與InSAR結(jié)果做對(duì)比（圖4），并根據(jù)公式（7）計(jì)算兩者之間的相關(guān)性。

圖4 水準(zhǔn)與時(shí)序InSAR數(shù)據(jù)之間相關(guān)性

通過(guò)式（7）計(jì)算結(jié)果可以看出兩類監(jiān)測(cè)數(shù)據(jù)的線性相關(guān)程度較高，相關(guān)系數(shù)R2為0.982，表明二者具有較強(qiáng)的一致性。為進(jìn)一步驗(yàn)證時(shí)序InSAR數(shù)據(jù)的可靠性，將2 種手段得到的監(jiān)測(cè)點(diǎn)位累積形變量進(jìn)行對(duì)比分析，結(jié)果如圖5所示。

圖5 水準(zhǔn)測(cè)量結(jié)果與時(shí)序InSAR形變結(jié)果對(duì)比

從圖5 中可以看出水準(zhǔn)測(cè)量結(jié)果與時(shí)序InSAR 形變結(jié)果在不同點(diǎn)位的最大形變差值為+9.8 mm，最小形變差值為-0.5 mm，平均差值為3.0 mm，中誤差為6.7 mm。造成二者之間差異的主要原因有：①I(mǎi)nSAR監(jiān)測(cè)結(jié)果的像元分辨率為20 m×20 m，無(wú)法在空間位置上與水準(zhǔn)點(diǎn)保持一致；②2 種監(jiān)測(cè)方式獲取的監(jiān)測(cè)結(jié)果時(shí)間上不能完全同步，有數(shù)天的差異；③時(shí)序In-SAR結(jié)果仍殘存有大氣、DEM、軌道誤差。從整體趨勢(shì)上來(lái)看，二者變化基本保持一致，時(shí)序InSAR技術(shù)監(jiān)測(cè)精度小于1 cm，具有較高的可靠性。

4 結(jié)果分析

4.1 時(shí)序結(jié)果分析

從圖3可以看出晉陽(yáng)高速、安陽(yáng)高速兩旁分布有大大小小十幾個(gè)沉陷區(qū)，均位于采空區(qū)正上方，其中面積較大、形變量級(jí)＞50 mm、距離公路2.5 km以內(nèi)的沉陷區(qū)有6個(gè)，1號(hào)沉陷區(qū)邊緣已經(jīng)侵入公路界線，須謹(jǐn)防擴(kuò)大。從2019年1月至2020年7月間最大累積形變量達(dá)到14.8 cm，位于6#采空區(qū)范圍內(nèi)。選取10個(gè)典型沉陷區(qū)的最大形變點(diǎn)做進(jìn)一步的時(shí)序分析（圖6），發(fā)現(xiàn)1、2、10號(hào)點(diǎn)形變有逐漸變緩的趨勢(shì)，其他7個(gè)沉陷區(qū)仍處于勻速或加速下沉的趨勢(shì)中。對(duì)沉陷點(diǎn)正上方的采空區(qū)開(kāi)采時(shí)間進(jìn)行統(tǒng)計(jì)分析，結(jié)果如表2所示。

圖6 沉降時(shí)序曲線

表2 采空區(qū)情況統(tǒng)計(jì)表

上述煤礦的開(kāi)采方式均采用傾斜長(zhǎng)壁采煤法，傾斜角小于10°，回采方式為后退式回采。從表2中可以看出1、2、10#煤礦開(kāi)采時(shí)間較長(zhǎng)，初采時(shí)間均在2010年前，2019年底局部區(qū)域開(kāi)采完畢，地表沉陷在半年內(nèi)趨于穩(wěn)定。3#煤礦初采時(shí)間在2008年，采深較大，監(jiān)測(cè)點(diǎn)位于2019年新掘采空區(qū)上方，形變趨勢(shì)較為平緩。4～8#采空區(qū)的初采時(shí)間為近4 a，形變速率和開(kāi)采深度正相關(guān)，采深越大，形變速率越小，且目前仍在開(kāi)采中，形變量持續(xù)增加，部分點(diǎn)還有加速下沉跡象。9#煤礦初采時(shí)間為2012 年，監(jiān)測(cè)點(diǎn)位于2020年待采區(qū)域內(nèi)，可以看出從2020年6月開(kāi)始有加速下沉跡象。

4.2 沉陷區(qū)動(dòng)態(tài)變化分析

為了深入分析采煤沉陷區(qū)的動(dòng)態(tài)變化過(guò)程，以6#、8#沉陷區(qū)作為研究對(duì)象，繪制不同方向的多個(gè)剖面進(jìn)行分析，其中剖面AA’與6#采空區(qū)多個(gè)采煤斷面掘進(jìn)方向垂直，剖面BB’與6#采空區(qū)掘進(jìn)方向平行，剖面DD’與8#采空區(qū)的掘進(jìn)方向平行。

圖7 6#采空區(qū)典型剖面累積沉降變化

6#采空區(qū)包括從東至西多組不同時(shí)期的子采空區(qū)組成，每個(gè)子采空區(qū)的掘進(jìn)方向均為由南至北，采用斜井開(kāi)拓，走向長(zhǎng)壁分層綜采，支護(hù)方式為支撐掩護(hù)式綜采支架，采深為320～470 m。為了細(xì)致地研究不同時(shí)期地面沉陷動(dòng)態(tài)變化過(guò)程，判斷開(kāi)采時(shí)間和掘進(jìn)方向與地面沉陷之間的關(guān)系，選取3 條剖面上間隔一個(gè)季度的多期累積沉降量曲線進(jìn)行分析。

從剖面AA’的沉降曲線上可以看出有3 個(gè)沉降中心，分別對(duì)應(yīng)2019 年、2018 年和2017 年3 個(gè)子采空區(qū)橫斷面，開(kāi)采時(shí)間越晚的采空區(qū)累積沉降量越大。2016 年采空斷面由于時(shí)間超過(guò)2 a，沉降趨勢(shì)基本穩(wěn)定。剖面BB’從南向北累積沉降量越來(lái)越大，對(duì)應(yīng)的開(kāi)采時(shí)間由2018 年年初至年末，初采時(shí)間越晚，被記錄到的殘余形變?cè)酱?。剖面CC’橫跨2 個(gè)子采空區(qū)，沉降曲線上相應(yīng)有2 個(gè)沉降中心，且從2019 年4—7 月第二個(gè)沉降中心開(kāi)始形成，與該區(qū)域的開(kāi)采時(shí)間基本一致，在該時(shí)間段內(nèi)第一個(gè)沉降中心也有加速趨勢(shì)，可能為第一個(gè)工作面上繼續(xù)進(jìn)行分層開(kāi)采造成的。從以上分析可以看出InSAR結(jié)果的沉降剖面上能較好地識(shí)別出不同時(shí)期采空區(qū)的沉降中心，且累積沉降量的變化趨勢(shì)能清晰反映開(kāi)采方向、開(kāi)采時(shí)間對(duì)地面沉降造成的影響。利用DD’剖面進(jìn)一步分析開(kāi)采方向上不同時(shí)間節(jié)點(diǎn)不同點(diǎn)位的沉降速率變化。8#采空區(qū)包括從南至北多組不同時(shí)期的子采空區(qū)，其中DD’剖面位于2019—2020 年采空區(qū)上，開(kāi)采方式采用斜井開(kāi)拓，長(zhǎng)壁開(kāi)采，采深350～500 m。沿掘進(jìn)方向依次選取E、F、G 3 個(gè)監(jiān)測(cè)點(diǎn)位進(jìn)行沉降速率分析，E點(diǎn)從2019年2月開(kāi)始沉降速率由小于10 mm∕a急速變大，F(xiàn) 點(diǎn)從2019 年4 月份開(kāi)始沉降速率迅速增加，G 點(diǎn)的沉降速率峰值則出現(xiàn)在2019年7 月份，最大沉降速率分別達(dá)到68 mm∕a、79 mm∕a和78 mm∕a。沉降速率超過(guò)50 mm∕a 的時(shí)間周期為2～3個(gè)月，其他時(shí)間段內(nèi)沉降速率會(huì)有小幅波動(dòng)，可能為上覆巖層穩(wěn)定性發(fā)生改變或是周邊煤礦采掘造成的影響。由于煤礦開(kāi)采造成的地表沉陷響應(yīng)較快，不同點(diǎn)位的沉降速率峰值基本能夠代表掘進(jìn)位置變化。

圖8 8#采空區(qū)不同點(diǎn)位沉降速率圖

5 結(jié) 論

本文采用時(shí)序InSAR監(jiān)測(cè)技術(shù)，基于2019 年1 月至2020 年7 月間共25 景Sentinel-1A 衛(wèi)星數(shù)據(jù)，對(duì)山西陽(yáng)城地區(qū)的采空區(qū)發(fā)育情況進(jìn)行調(diào)查研究，并利用同期水準(zhǔn)監(jiān)測(cè)結(jié)果進(jìn)行精度驗(yàn)證?？傻贸鋈缦陆Y(jié)論：

1）陽(yáng)城北部分布有10個(gè)較大規(guī)模的采空沉陷區(qū)，監(jiān)測(cè)周期內(nèi)視線向最大沉降量達(dá)到14.8 cm，現(xiàn)采空區(qū)形變沉陷趨勢(shì)較為劇烈，停采后沉陷趨勢(shì)明顯變緩。

2）InSAR監(jiān)測(cè)結(jié)果剖面能清晰識(shí)別采煤斷面，以及不同時(shí)刻采煤引起的地表沉陷變化，根據(jù)地表沉降峰值還能判斷出采煤實(shí)時(shí)掘進(jìn)位置。

3）InSAR技術(shù)可應(yīng)用于黃土地區(qū)采空區(qū)沉陷識(shí)別及監(jiān)測(cè)，為盜采、超采等隱蔽違法違規(guī)行為調(diào)查提供參考依據(jù)。