空天地一體化監(jiān)測聯(lián)合反演開采沉陷概率積分預(yù)計參數(shù)研究

胡東升，程小凱，張雅飛，李 濤，廉旭剛

(1.華陽新材料科技集團(tuán)有限公司，山西 陽泉 045000；2.自然資源部第一大地測量隊(duì)，陜西 西安 710054；3.太原理工大學(xué) 礦業(yè)工程學(xué)院，山西 太原 030024)

煤炭開采為社會帶來經(jīng)濟(jì)效益的同時，其開采將導(dǎo)致地表沉降，給礦區(qū)以及周圍環(huán)境造成一定的負(fù)面影響。通過對開采地區(qū)長期、動態(tài)的沉陷監(jiān)測，可及時獲取礦區(qū)地表形變數(shù)據(jù)，分析移動變化規(guī)律，科學(xué)地指導(dǎo)礦區(qū)的開采活動，從而最大限度地降低因礦區(qū)開采帶來的損失。以往采用的水準(zhǔn)測量和GPS(Global Positioning System，GPS)測量等傳統(tǒng)的礦區(qū)形變監(jiān)測手段，成熟度雖高，但工作強(qiáng)度大、耗費(fèi)人力物力、易受天氣影響、工作周期長、監(jiān)測點(diǎn)易受損害、監(jiān)測范圍小，不能實(shí)現(xiàn)對礦區(qū)持續(xù)、大范圍的動態(tài)監(jiān)測[1]。

合成孔徑雷達(dá)差分干涉測量(Differential Interferometric Synthetic Aperture Radar，DInSAR)作為新興的對地觀測技術(shù)，能夠?qū)崿F(xiàn)全天時、全天候?qū)ρ芯繀^(qū)進(jìn)行觀測，且監(jiān)測精度高，覆蓋范圍廣[2]。SAR衛(wèi)星獲取的數(shù)據(jù)，不僅包括強(qiáng)度信息，還包含相位信息[3]。楊氏雙縫干涉實(shí)驗(yàn)為InSAR技術(shù)的出現(xiàn)奠定了理論基礎(chǔ)。1989年，Gabriel[4]首次使用DInSAR技術(shù)將干涉相位中的形變相位和地形相位分離，開啟了InSAR監(jiān)測地表形變的新篇章。1993年，Gabriel等DInSAR技術(shù)對地表形變進(jìn)行監(jiān)測，證實(shí)了其精度可達(dá)厘米級甚至毫米級。1996年，Carnec C[5]對ERS-1影像數(shù)據(jù)進(jìn)行處理分析，初步證明了DInSAR技術(shù)可用于礦區(qū)地表沉降監(jiān)測。國內(nèi)學(xué)者經(jīng)過多年的努力，在該領(lǐng)域也取得了重大突破。吳立新[6]等對覆蓋中國東部典型礦區(qū)的5景ERS1/2影像進(jìn)行差分處理，驗(yàn)證了DInSAR技術(shù)用于礦區(qū)沉降監(jiān)測的可行性；楊嘉威[7]等引入分布式散射體 (Distributed Scatterer，DS) 目標(biāo)，通過DS-InSAR技術(shù)對礦區(qū)地表和鐵路沿線的形變規(guī)律進(jìn)行了分析。楊澤發(fā)、朱建軍等人[8-10]分別基于單軌InSAR、基于InSAR時序形變以及聯(lián)合InSAR技術(shù)和水準(zhǔn)數(shù)據(jù)等手段對礦區(qū)開采引起的形變進(jìn)行監(jiān)測，并對其沉降規(guī)律進(jìn)行分析，為該領(lǐng)域做出了許多貢獻(xiàn)。廉旭剛[11]等結(jié)合Sentinel1A和1B數(shù)據(jù)提高DInSAR的時間分辨率，實(shí)現(xiàn)對大同礦區(qū)沉降的高精度監(jiān)測。

無人機(jī)最早出現(xiàn)在1917年，是為滿足軍事上的需求。Przybilla和Wester Ebbinghaus[12]分別于1979年和1980年成功進(jìn)行固定翼和旋翼無人機(jī)的航空攝影測試。Pawe[13]等利用無人機(jī)攝影測量手段監(jiān)測到礦區(qū)的不連續(xù)形變，促進(jìn)了無人機(jī)技術(shù)應(yīng)用于礦區(qū)監(jiān)測的發(fā)展。21世紀(jì)后，國內(nèi)加強(qiáng)對無人機(jī)攝影測量技術(shù)的研究[14]。2003年，盛業(yè)華[15]等進(jìn)行地表塌陷變形試驗(yàn)并采用數(shù)字?jǐn)z影測量對其觀測。Zhou[16]等利用無人機(jī)攝影測量手段對煤礦區(qū)地表沉降進(jìn)行監(jiān)測，并反演了開采沉降參數(shù)。邱亞輝[17]等人利用無人機(jī)航測技術(shù)獲取露天采坑的DEM數(shù)據(jù)，并對多期DEM進(jìn)行差值，得到研究區(qū)相關(guān)變化量。戴嵩[18]等使用無人機(jī)技術(shù)得到尾礦壩正射影像，監(jiān)測出最大沉降范圍在0.16m之內(nèi)。廉旭剛[19]等通過無人機(jī)攝影測量技術(shù)實(shí)現(xiàn)了對陽煤礦區(qū)的地表沉陷監(jiān)測，證明了該技術(shù)在大梯度形變區(qū)域監(jiān)測精度高，監(jiān)測微小形變的精度有待改善。

目前，針對無人機(jī)與InSAR聯(lián)合監(jiān)測地災(zāi)的研究很少。大多是使用二者中的其一手段進(jìn)行觀測，借助另一種手段對監(jiān)測結(jié)果進(jìn)行空間驗(yàn)證，未能利用二者監(jiān)測精度不同的優(yōu)勢?？紤]到InSAR和無人機(jī)技術(shù)各自的優(yōu)缺點(diǎn)，本文將兩種技術(shù)的優(yōu)勢進(jìn)行互補(bǔ)，融合二者的監(jiān)測值，并使用融合后的下沉曲線對一礦81403工作面的概率積分參數(shù)進(jìn)行反演。

1 研究區(qū)概況

華陽集團(tuán)一礦81403工作面采用放頂煤綜合機(jī)械化采煤法，走向長為1345m，傾向長為226m，煤層的平均傾角是4°，平均煤厚7.24m，平均采深為446.8m。研究區(qū)位置關(guān)系如圖1所示。依據(jù)《建筑物、水體、鐵路及主要井巷煤柱留設(shè)與壓煤開采規(guī)范》的有關(guān)要求，結(jié)合本區(qū)內(nèi)觀測地表情況，采用剖面線觀測線的方式在81403工作面上方布置了三條觀測線，分別為半條走向觀測線A線和兩條傾向觀測線B線和C線。根據(jù)開采深度和設(shè)站目的，并根據(jù)實(shí)際情況，觀測線上的測點(diǎn)按30m距離埋設(shè)。設(shè)計各觀測線的總長度、分段長度以及點(diǎn)間距和點(diǎn)數(shù)見表1。地表移動觀測站的首次觀測日期為2019年7月24日，末次觀測日期為2021年11月23日，共進(jìn)行了23次地表移動觀測。本次研究主要對走向A線33個有效監(jiān)測點(diǎn)以及傾向C線17個有效監(jiān)測點(diǎn)進(jìn)行分析。

圖1 研究區(qū)相對位置

2 數(shù)據(jù)采集及處理方法

2.1 數(shù)據(jù)采集

研究選用的SAR影像為Sentinel-1A衛(wèi)星的數(shù)據(jù)，該衛(wèi)星由歐空局研發(fā)，并于2014年4月發(fā)射。衛(wèi)星上搭載C波段的傳感器，波長為5.6cm，分辨率為5m×20m，是一個全天時、全天候工作的雷達(dá)成像系統(tǒng)。通過ASF(Alaska Satellite Facility)下載平臺，覆蓋一礦81403工作面的單視復(fù)數(shù)(Single Look Complex，SLC)、升軌軌道的SAR影像數(shù)據(jù)。選取的數(shù)據(jù)時間跨度為2020年3月22日至2022年1月11日，并且為保證差分干涉結(jié)果的質(zhì)量，以衛(wèi)星的重訪周期12d為間隔進(jìn)行兩兩差分干涉，共55景影像。使用SRTM(Shuttle Radar Topography Mission)的30m分辨率DEM(Digital Elevation Model，DEM)作為差分的外部DEM。

表1 觀測線長度、點(diǎn)間距和點(diǎn)數(shù)明細(xì)

無人機(jī)數(shù)據(jù)使用飛馬智能航測系統(tǒng)D2000采集，搭載D-LIDAR2000激光雷達(dá)模塊進(jìn)行LiDAR點(diǎn)云數(shù)據(jù)的采集，并且保證與SAR數(shù)據(jù)時間上的統(tǒng)一。該采集系統(tǒng)聯(lián)合飛馬系統(tǒng)高精度的組合導(dǎo)航產(chǎn)品，測量精度可達(dá)5cm(50m航高)，當(dāng)?shù)V區(qū)地表沉降量為分米級時，該手段可滿足沉陷監(jiān)測的需求。實(shí)地采集數(shù)據(jù)時，依據(jù)研究區(qū)的地形起伏等情況進(jìn)行航線規(guī)劃，重疊度不低于50%。使用變高航線，當(dāng)相對地面的飛行高度為177m時，可以采集到密度為64點(diǎn)/m2的點(diǎn)云。數(shù)據(jù)采集前的航線規(guī)劃如圖2所示。

圖2 機(jī)載LiDAR航線規(guī)劃

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

2.2.1 DInSAR數(shù)據(jù)處理

使用GAMMA軟件的Diff&Geo (Differential Interferometry and Geocoding) 差分干涉與地理編碼處理模塊。對原始SAR數(shù)據(jù)進(jìn)行多視處理，多視比為4∶1，以此抑制影像的斑點(diǎn)噪聲。將主影像與副影像兩景數(shù)據(jù)進(jìn)行配準(zhǔn)后，使用外部DEM進(jìn)行二軌差分干涉。通過自適應(yīng)濾波和去平地效應(yīng)等步驟，將形變相位外的其余相位去除，再進(jìn)行相位解纏以及地理編碼，經(jīng)雷達(dá)視線向(Line of Sight，LOS)形變轉(zhuǎn)豎直向形變得到最終的沉降結(jié)果。二軌法DInSAR處理流程如圖3所示。

圖3 二軌法DInSAR處理流程

處理數(shù)據(jù)時將數(shù)據(jù)按時間序列組成54個干涉像對，以保證每個干涉像對的時間間隔最小，最后利用ArcGIS軟件將差分結(jié)果進(jìn)行疊加得到2020年3月22日至2022年1月11日的差分干涉圖，如圖4所示，各時間為疊加的最終時間。

對差分干涉結(jié)果圖進(jìn)行分析時，會發(fā)現(xiàn)DInSAR手段在監(jiān)測大梯度變形時會產(chǎn)生失相干。DInSAR技術(shù)雖然可以監(jiān)測到整體下沉盆地的范圍，在沉降盆地邊緣處表達(dá)能力強(qiáng)，但是無法監(jiān)測到盆地中間的大梯度變形，由于沉降盆地中心沉降量級大、突變等原因，造成干涉相位的不連續(xù)，出現(xiàn)失相干及監(jiān)測結(jié)果的錯誤，使得盆地中心的監(jiān)測值與實(shí)際下沉情況不符合。

2.2.2 無人機(jī)機(jī)載 LiDAR處理

通過機(jī)載 LiDAR 獲取原始數(shù)據(jù)的預(yù)處理，即對采集結(jié)果進(jìn)行數(shù)據(jù)解算與轉(zhuǎn)換，生成標(biāo)準(zhǔn)las點(diǎn)云數(shù)據(jù)文件。采用ATIN算法對點(diǎn)云數(shù)據(jù)進(jìn)行點(diǎn)云濾波、反距離加權(quán)作插值處理，得到2020年3月26日與2022年1月16兩期相應(yīng)的DEM數(shù)據(jù)。使用ArcGIS將兩期DEM數(shù)據(jù)進(jìn)行相減，得到該時間段內(nèi)的無人機(jī)監(jiān)測的地表沉降，結(jié)果如圖5所示。

圖4 DInSAR差分結(jié)果

采用無人機(jī)監(jiān)測開采沉陷能夠較為全面地反映出礦區(qū)沉陷影響范圍，在對大梯度形變區(qū)域進(jìn)行監(jiān)測時，測量精度高，可以監(jiān)測出最大下沉。但監(jiān)測精度達(dá)不到傳統(tǒng)監(jiān)測毫米級精度，難以對礦區(qū)邊緣進(jìn)行高精度監(jiān)測。

圖5 無人機(jī)首尾兩期DEM相減

3 概率積分法參數(shù)反演

3.1 融合DInSAR與無人機(jī)技術(shù)監(jiān)測值

DInSAR技術(shù)在大量級、大梯度區(qū)域的監(jiān)測精度差，更適用于沉降盆地邊緣沉降值的監(jiān)測。無人機(jī)技術(shù)監(jiān)測沉降盆地中心大變形時精度較高，在沉降盆地邊緣監(jiān)測精度低。因此把研究區(qū)的無人機(jī)監(jiān)測結(jié)果與DInSAR差分干涉得到的沉降值進(jìn)行融合，得到更加完整、準(zhǔn)確的監(jiān)測結(jié)果。

在實(shí)際的采煤沉陷研究中，一般取沉陷值為10mm的點(diǎn)作為邊界點(diǎn)。然而，DInSAR技術(shù)的觀測精度可以達(dá)到毫米級。因此，為了獲取礦區(qū)更完整的高精度監(jiān)測，使用DInSAR手段監(jiān)測沉降盆地邊緣的小變形，使用無人機(jī)技術(shù)監(jiān)測沉降盆地中心的大梯度變形。

DInSAR干涉的相干性決定其相位解纏結(jié)果的準(zhǔn)確性，因此使用DInSAR平均相干系數(shù)對DInSAR監(jiān)測值進(jìn)行篩選。Baran等通過經(jīng)驗(yàn)統(tǒng)計等方法對DInSAR可檢測的形變梯度與相干性的關(guān)系進(jìn)行了深入研究，得到了以下數(shù)學(xué)模型[20]。

Dmax=dmax+0.002(γ-1)

(1)

式中，Dmax為最大可監(jiān)測到的形變梯度；γ為差分干涉圖的相干系數(shù)，該值介于0到1之間；dmax為DInSAR理論上可監(jiān)測到的最大形變梯度，等于衛(wèi)星傳感器波長與2倍影像分辨率之比。

由式(1)可看出相干性增大，可監(jiān)測的最大形變梯度也隨之增大。實(shí)際中受許多種去相關(guān)因素的影響，DInSAR可監(jiān)測的形變梯度一般遠(yuǎn)小于理論值。對于采用 C 波段、分辨率為20m的Sentinel-1A的影像，通過式(1)可知，若DInSAR相干系數(shù)小于0.3，則其最大可監(jiān)測的形變梯度等于0。

基于以上分析，選取相干系數(shù)低于0.3并最靠近研究區(qū)工作面沉降盆地邊緣的監(jiān)測點(diǎn)，以該點(diǎn)作為DInSAR監(jiān)測結(jié)果的篩選閾值。通過提取各監(jiān)測點(diǎn)的相干系數(shù)，經(jīng)對比發(fā)現(xiàn)，從A9點(diǎn)之后的A12監(jiān)測點(diǎn)開始出現(xiàn)相干系數(shù)小于0.3的情況，且越靠近沉陷盆地中心，相干性系數(shù)越低甚至出現(xiàn)失相干現(xiàn)象。因此確定對于A9號監(jiān)測點(diǎn)以及其之前的點(diǎn)采用DInSAR監(jiān)測值，其余監(jiān)測點(diǎn)采用無人機(jī)監(jiān)測值。對于C9號與C24號之間的監(jiān)測點(diǎn)，其監(jiān)測值采用無人機(jī)監(jiān)測值，其余采用DInSAR監(jiān)測值。

3.2 概率積分預(yù)計參數(shù)反演

當(dāng)前開采沉陷預(yù)計使用較為廣泛的一種方法為概率積分法(Probability Integral Method，PIM)，由我國學(xué)者劉寶琛、廖國華等在隨機(jī)介質(zhì)理論的基礎(chǔ)上發(fā)展而得[21]。概率積分法共有5個與沉陷相關(guān)的參數(shù)，各自為下沉系數(shù)q、主要影響角正切tanβ，水平移動系數(shù)b，拐點(diǎn)偏移距Si(i=1，2，3，4)和開采影響傳播角θ0[22]。一般通過在工作面上方布置觀測站來獲取實(shí)測數(shù)據(jù)，以此數(shù)據(jù)為基礎(chǔ)確定概率積分法最佳參數(shù)[23]。

本次研究使用InSAR技術(shù)與無人機(jī)攝影測量技術(shù)的融合監(jiān)測值來反演概率積分參數(shù)。準(zhǔn)備好預(yù)計所需的工作面及觀測線數(shù)據(jù)，根據(jù)工作面實(shí)際情況設(shè)置初始概率積分參數(shù)，采用《山區(qū)煤礦開采地表移動變形預(yù)計系統(tǒng)(MMSPS)》計算軟件進(jìn)行計算。將計算結(jié)果與融合后的曲線進(jìn)行對比，分析結(jié)果的差異以調(diào)整預(yù)計參數(shù)，再次進(jìn)行模擬計算，重復(fù)該步驟直到取得最佳的模擬效果，該結(jié)果下的預(yù)計參數(shù)就是最優(yōu)模擬參數(shù)。經(jīng)多次模擬計算得到的概率積分參數(shù)結(jié)果見表2。

表2 模擬求參計算結(jié)果

由于采用InSAR與無人機(jī)技術(shù)監(jiān)測開始于2020年3月，而現(xiàn)場實(shí)測開始于2019年7月24日，為了使時間統(tǒng)一，對實(shí)測數(shù)據(jù)進(jìn)行了處理，將2022年1月的累計實(shí)測下沉減去2020年3月的實(shí)測累計沉降，得到2020年3月至2022年1月的實(shí)測下沉值。做差后的實(shí)測下沉曲線與融合曲線整體趨勢相同。

A線和C線融合下沉曲線結(jié)果與模擬下沉曲線結(jié)果對比如圖6所示。

圖6 81403工作面A線和C線融合下沉曲線與模擬下沉曲線對比

使用融合DInSAR監(jiān)測值與無人機(jī)監(jiān)測值后的下沉曲線，對一礦81403工作面概率積分法參數(shù)進(jìn)行反演，反演結(jié)果為：下沉率q=0.8，水平移動系數(shù)b=0.3，主要影響角正切tanβ=2.5，開采影響傳播角θ0=85°，拐點(diǎn)偏移距S1=S2=S3=S4=-20m。通過對比，反演的參數(shù)結(jié)果與實(shí)測數(shù)據(jù)所擬合的參數(shù)結(jié)果基本接近，證明了獲得的概率積分法參數(shù)的合理性和可靠性。經(jīng)過概率積分法求參，只得到了非充分采動情況下的預(yù)計參數(shù)。依據(jù)《建筑物、水體、鐵路及主要井巷煤柱留設(shè)與壓煤開采規(guī)程》所述，要獲得充分采動時的預(yù)計參數(shù)，需要對非充分采動條件下的預(yù)計參數(shù)進(jìn)行修正，其中水平移動系數(shù)和開采影響傳播角與回采尺寸之間的關(guān)系不明顯，可以不予修正。按照規(guī)程中下沉系數(shù)、主要影響角正切、拐點(diǎn)偏移距等參數(shù)與回采尺寸間的關(guān)系，再結(jié)合研究區(qū)覆巖性質(zhì)，對概率積分預(yù)計參數(shù)進(jìn)行修正。陽泉礦區(qū)覆巖性質(zhì)為中硬，修正后的參數(shù)：下沉系數(shù)q=0.82，水平移動系數(shù)b=0.3，主要影響角正切tanβ=2.6，開采影響傳播角θ0=85°，拐點(diǎn)偏移距S1=S2=S3=S4=-33.3m。該研究獲取的DInSAR與無人機(jī)監(jiān)測值均為一維下沉值，反演出水平移動系數(shù)這個參數(shù)的參考意義不大。需要說明的是，因?yàn)橐坏V81403工作面的觀測站于2019年7月24開始觀測，而無人機(jī)與InSAR技術(shù)監(jiān)測沉降是于2020年3月開始，此時該工作面已經(jīng)開采了300多米，所以造成A線預(yù)計下沉曲線比無人機(jī)與DInSAR融合后的下沉曲線靠右。

4 結(jié) 論

1)通過對研究區(qū)的Sentinel-1A數(shù)據(jù)進(jìn)行差分處理，結(jié)果表明：由于沉降量、突變等原因使得影像出現(xiàn)失相干的現(xiàn)象，對DInSAR監(jiān)測結(jié)果影響較大，沉降邊緣沉降量小的點(diǎn)，誤差??；而在沉降中心沉降量大的點(diǎn)，DInSAR在該區(qū)域獲取的沉降值仍十分微小，未能有效監(jiān)測出正確的下沉值。符合DInSAR監(jiān)測微小形變的特點(diǎn)。

2)由無人機(jī)機(jī)載LiDAR得到的監(jiān)測值結(jié)果表明：采用無人機(jī)監(jiān)測可以監(jiān)測出最大下沉，較全面的反映出礦區(qū)沉陷影響范圍。但是難以對礦區(qū)邊緣進(jìn)行高精度監(jiān)測，邊緣表達(dá)能力差，無法監(jiān)測微小沉陷變形。

3)融合DInSAR和無人機(jī)兩種非等精度的監(jiān)測值，一定程度上可以彌補(bǔ)DInSAR手段在大量級、大梯度形變區(qū)失相干的問題以及無人機(jī)技術(shù)在監(jiān)測邊緣微小形變時精度低的不足。并利用融合后的下沉曲線對一礦81403工作面的概率積分參數(shù)進(jìn)行反演。結(jié)果表明，反演的參數(shù)結(jié)果與實(shí)測參數(shù)結(jié)果基本接近。該方法對礦區(qū)地質(zhì)災(zāi)害評估、預(yù)防礦山地質(zhì)災(zāi)害具有一定的參考價值。