王 瑞，薛玉玲，萬(wàn)士楨，鄭玉琳 ，張 婷

(1.贛南科技學(xué)院，江西 贛州 341000；2.中國(guó)礦業(yè)大學(xué) 環(huán)境與測(cè)繪工程學(xué)院，江蘇 徐州 221116;3.江西理工大學(xué) 應(yīng)用科學(xué)學(xué)院，江西 贛州 341000)

0 引言

我國(guó)礦產(chǎn)資源比較豐富，礦產(chǎn)資源開(kāi)采在帶來(lái)可觀的經(jīng)濟(jì)效益的同時(shí)，也給生態(tài)環(huán)境以及人類(lèi)的生產(chǎn)生活帶來(lái)了一定危害。采空區(qū)上部巖層在重力作用下會(huì)發(fā)生一系列變形，造成礦區(qū)地表下沉，甚至引發(fā)一系列環(huán)境地質(zhì)問(wèn)題。

隨著衛(wèi)星遙感技術(shù)的發(fā)展，合成孔徑雷達(dá)干涉測(cè)量(InSAR)技術(shù)逐漸成為我國(guó)形變監(jiān)測(cè)的主要手段。該技術(shù)是一種全新的空間大地測(cè)量技術(shù)，其具有全天觀測(cè)、實(shí)時(shí)觀測(cè)且觀測(cè)范圍廣的優(yōu)點(diǎn)，能夠獲取高精度、高分辨率、大范圍的觀測(cè)信息，且無(wú)需在地面布置觀測(cè)點(diǎn)，因此該技術(shù)被廣泛應(yīng)用于地表及各類(lèi)地物形變監(jiān)測(cè)。在監(jiān)測(cè)地震同震變形場(chǎng)時(shí)，差分合成孔徑雷達(dá)干涉測(cè)量(DInSAR)技術(shù)能夠獲得雷達(dá)視線方向上的形變量[1-2]。

地下煤炭開(kāi)采常會(huì)引起地表沉陷，從小規(guī)模小尺度的沉陷逐漸發(fā)展到區(qū)域性大尺度沉陷，沉陷量及范圍主要取決于地質(zhì)條件、巖石力學(xué)性質(zhì)以及開(kāi)采技術(shù)方法[3]。地表沉陷并不是在短時(shí)間內(nèi)發(fā)生的，而是隨著開(kāi)采面的增大而逐漸顯現(xiàn)的，因此持續(xù)監(jiān)測(cè)是必要的。地下煤炭開(kāi)采引起的大尺度形變是一個(gè)復(fù)雜的時(shí)空過(guò)程，通過(guò)地表監(jiān)測(cè)對(duì)礦區(qū)地質(zhì)災(zāi)害的防治具有重要意義。針對(duì)采煤引起的地表沉陷，推廣程度最高、應(yīng)用范圍最廣的地表沉陷信息獲取技術(shù)為常規(guī)觀測(cè)站法，即在擬開(kāi)采工作面上方沿工作面走向和傾向每隔 5～25 m布置一系列觀測(cè)點(diǎn)，在工作面開(kāi)采的不同階段采用高精度水準(zhǔn)儀、全站儀或GPS進(jìn)行實(shí)地監(jiān)測(cè)，經(jīng)數(shù)據(jù)處理獲得沿煤層走向和傾向的下沉量及變形情況。該方法是逐點(diǎn)監(jiān)測(cè)，空間分辨率低、效率不高；在長(zhǎng)期監(jiān)測(cè)過(guò)程中，觀測(cè)站易受礦區(qū)地表復(fù)雜條件影響。

隨著測(cè)量技術(shù)的發(fā)展，InSAR技術(shù)為解決區(qū)域大范圍沉陷監(jiān)測(cè)問(wèn)題提供了一種有效方法，其可以提取存檔數(shù)據(jù)分析沉陷前后的規(guī)律[4-7]。本文以西部某礦區(qū)的工作面為研究對(duì)象，利用Radarsat-2圖像作為實(shí)驗(yàn)數(shù)據(jù)，通過(guò)圖像解譯獲取沉陷盆地邊界。

1 InSAR技術(shù)原理及DEM修正

1.1 InSAR技術(shù)

InSAR技術(shù)是一種圖像處理技術(shù)，可以從SAR傳感器的兩個(gè)天線干涉通道生成數(shù)字高程模型(DEM)和相對(duì)相干性。隨著InSAR技術(shù)的發(fā)展，其應(yīng)用領(lǐng)域不斷擴(kuò)大，如通過(guò)同一區(qū)域兩顆衛(wèi)星通道測(cè)量衛(wèi)星和地球表面之間的相位差，此相位差的產(chǎn)生原因主要有：①兩次或多次衛(wèi)星軌跡之間的位置差異，但是差異遠(yuǎn)小于衛(wèi)星到地球表面的距離，因此可以應(yīng)用于地形地貌形變監(jiān)測(cè)；②兩次或多次衛(wèi)星采集觀測(cè)到的區(qū)域位移，此位移可能由地下資源的開(kāi)采、地震及火山運(yùn)動(dòng)等因素造成。因此，當(dāng)SAR系統(tǒng)對(duì)同一區(qū)域進(jìn)行兩次或多次觀測(cè)時(shí)，如該區(qū)域的幾何位置相對(duì)于傳感器發(fā)生了變化，則視為地表產(chǎn)生了形變。

SAR系統(tǒng)在采集數(shù)據(jù)過(guò)程中，反射信號(hào)會(huì)受噪聲、地形、大氣及地表運(yùn)動(dòng)等因素的影響。DInSAR原理及地表形變幾何關(guān)系示意圖見(jiàn)圖 1。

圖1 DInSAR原理及地表形變幾何關(guān)系圖

地表形變干涉相位的計(jì)算步驟為：干涉相位φ表示同一個(gè)目標(biāo)P由衛(wèi)星S1和S2在兩次不同采集的時(shí)間內(nèi)，P點(diǎn)發(fā)生了形變而移至P′點(diǎn)。干涉相位φ受多個(gè)因素的影響，其計(jì)算式如下[8]：

(1)

(2)

φInt=φTopo+φFlat+φDefo+φAtm+φNoise，

(3)

(4)

h=H-R1cosθ，

(5)

h′=H-R2cosθ，

(6)

式中：λ為雷達(dá)衛(wèi)星的波長(zhǎng)，R1、R2分別為衛(wèi)星兩次通過(guò)目標(biāo)P、P′時(shí)的距離，φS1、φS2分別為主衛(wèi)星和副衛(wèi)星位置的干涉相位，φTopo為受地形影響的相位，φFlat為參考面因素引起的相位，φDefo為沿視線向(LOS)地表形變因素引起的相位，φAtm為大氣延遲因素引起的相位，φNoise為噪聲因素引起的相位。式(4)、式(5)和式(6)揭示了干涉相位φInt與高程h、h′之間的函數(shù)關(guān)系。已知天線位置的高H、基線長(zhǎng)度B、基線水平角α、往返雙程相位差N，可得到雷達(dá)成像視角θ，再計(jì)算形變前后點(diǎn)P的高程h和h′，最后獲得形變量。

采用GAMMA的傳統(tǒng)連續(xù)干涉測(cè)量方法解算有效的影像。ZEBKER等[9]的研究成果表明，為能夠準(zhǔn)確有效解算出雷達(dá)視線的形變量，單位像素所能解算出的形變量為波長(zhǎng)的一半。RadarSat-2的波長(zhǎng)為0.056 m，故干涉圖所能解算出的最大單位像素理論變形值為0.028 m。當(dāng)任意兩幅影像之間相鄰單位像素的形變大于波長(zhǎng)的二分之一時(shí)，則可能無(wú)法在時(shí)域內(nèi)正確地解纏。

1.2 DEM修正

何秀風(fēng)等[10]的研究結(jié)果表明，在采用DInSAR技術(shù)進(jìn)行沉陷解算時(shí)，因垂直基線不會(huì)為0，所以DEM精度對(duì)兩種通法的解算具有關(guān)鍵性作用。因此本文在12.5m DEM的基礎(chǔ)上，融合無(wú)人機(jī)(Unmanned Aerial Vehicle,UAV)獲取的開(kāi)采前地表DEM數(shù)據(jù)，分別對(duì)由原始DEM和UAV獲取的DEM數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)差計(jì)算，再據(jù)此計(jì)算兩種融合DEM的權(quán)重[11-12]。融合模型如式(7)和式(8)所示，融合結(jié)果如圖2所示。

圖2 融合DEM

(7)

hr=Jhy1+(1-J)hy2，

(8)

式中：J為數(shù)據(jù)融合權(quán)重，取值范圍在0～1；s1和s2為融合DEM標(biāo)準(zhǔn)差；l為融合影像誤差的相關(guān)系數(shù)；hr為融合后的DEM數(shù)據(jù)；hy1和hy2為融合前兩種原始DEM數(shù)據(jù)。

2 研究區(qū)域概況及數(shù)據(jù)

2.1 研究區(qū)域概況

本文以西部某礦區(qū)為研究對(duì)象，其位于內(nèi)蒙古自治區(qū)鄂爾多斯市。礦區(qū)為不規(guī)則多邊形，南北長(zhǎng)14.40 km，東西最寬6.77 km，標(biāo)高1 150～1 420 m，地表起伏較大。礦區(qū)煤田廣闊、資源豐富、煤質(zhì)優(yōu)良；煤層埋深較小，約為150～456 m。

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

利用SAR影像作為實(shí)驗(yàn)解算數(shù)據(jù)，雷達(dá)影像采用加拿大太空總署發(fā)射的RadarSat-2(搭載C波段傳感器，重訪周期24 d)采集差分干涉測(cè)量的影像數(shù)據(jù)。干涉影像的采集日期及相關(guān)參數(shù)如表1所示。

表1 SAR影像參數(shù)

3 解算結(jié)果及精度分析

3.1 解算結(jié)果分析

該礦區(qū)某工作面自2018年7月12日開(kāi)始回采，至2018年10月25日回采完畢，共耗時(shí)105 d，工作面總長(zhǎng)度約1 263.7 m。工作面地表環(huán)境復(fù)雜，地勢(shì)起伏較大，地表土壤屬于沙地，表面有許多被雨水沖刷而成的溝壑，傳統(tǒng)測(cè)量方法難以監(jiān)測(cè)地表變形。因此，本文采用DInSAR技術(shù)構(gòu)建完整沉陷盆地模型，用以分析沉陷規(guī)律。

圖3為相鄰兩期圖像之間的時(shí)間序列形變圖。由圖3可知，受2S201工作面東側(cè)的2S202工作面開(kāi)采的殘余形變影響，2S202工作面自2018年6月9日至2019年1月11日地表出現(xiàn)形變，形變量隨著時(shí)間的推移不斷減小，此后逐漸消失；2S201工作面于2018年7月12日開(kāi)始回采，隨著工作面的推進(jìn)，地表開(kāi)始出現(xiàn)形變，并且形變沿著工作面回采的方向移動(dòng)；2018年7月12日至2018年10月25日為工作面強(qiáng)回采期，地表形變SAR影像失相干較為嚴(yán)重，僅能識(shí)別出較小形變；隨著回采的結(jié)束，地表形變逐漸減小，DInSAR技術(shù)可識(shí)別出沉陷盆地邊緣，最大識(shí)別沉陷量出現(xiàn)在2018年11月24日至2019年1月11日期間，為0.144 m。

圖4為2018年7月12日至2019年4月16日期間的時(shí)間序列累計(jì)形變圖。由圖4可知，在此時(shí)期內(nèi)視線方向所能識(shí)別的累計(jì)最大沉陷量為0.342 m，隨著工作面的推進(jìn)和時(shí)間的推移，受回采的影響沉陷盆地逐漸增大。圖4a-圖4c為回采時(shí)期的形變圖，該時(shí)期內(nèi)地表沉陷量較大，SAR影像失相干嚴(yán)重，識(shí)別沉陷量較小。

圖3 相鄰兩期圖像之間的時(shí)間序列形變圖圖4 DInSAR累計(jì)沉陷監(jiān)測(cè)

3.2 精度分析

在2S201工作面上部布置監(jiān)測(cè)點(diǎn)進(jìn)行精度驗(yàn)證。采用三角高程中間觀測(cè)方法進(jìn)行數(shù)據(jù)采集，根據(jù)坐標(biāo)位置提取DInSAR時(shí)序累計(jì)沉陷量。從監(jiān)測(cè)點(diǎn)中選取具有代表性的10個(gè)監(jiān)測(cè)點(diǎn)進(jìn)行誤差統(tǒng)計(jì)，結(jié)果如表2所示。

表2 DInSAR誤差統(tǒng)計(jì) 單位：m

DInSAR技術(shù)無(wú)法真實(shí)解算大量級(jí)形變，不適合大沉陷監(jiān)測(cè)，沉陷越小的區(qū)域誤差越小。由表2可知：對(duì)于小尺度形變，DInSAR技術(shù)的精度可在厘米級(jí)；當(dāng)沉陷量增大時(shí)，誤差為分米級(jí)和米級(jí)。2S201工作面自南向北回采，地表由南向北沉陷量遞減，因此在累計(jì)沉陷量監(jiān)測(cè)中北部大于南部。DInSAR技術(shù)對(duì)于沉陷盆地周邊小區(qū)域的沉陷識(shí)別較為敏感。

4 結(jié)論

本文采用DInSAR技術(shù)提取因井下煤炭開(kāi)采引起的地表變形數(shù)據(jù)，在數(shù)據(jù)解算過(guò)程中引入多源DEM數(shù)據(jù)，計(jì)算每種DEM數(shù)據(jù)標(biāo)準(zhǔn)差并建立不同權(quán)值融合模型。提取西部某礦區(qū)某工作面2018年7月12日至2019年4月16日的完整地表形變場(chǎng)進(jìn)行分析，結(jié)果表明，在地表?xiàng)l件特殊且地面觀測(cè)站較少的情況下，DInSAR技術(shù)可以有效分析地表形變規(guī)律，其雖受相干影響而無(wú)法識(shí)別大尺度真實(shí)形變，但可以識(shí)別沉陷邊界。