礦山開采損害InSAR / UAV融合監(jiān)測關(guān)鍵技術(shù)及應用

周大偉，安士凱，吳 侃，胡振琪，刁鑫鵬

(1. 中國礦業(yè)大學 環(huán)境與測繪學院，江蘇 徐州 221116；2. 平安煤炭開采工程技術(shù)研究院有限責任公司安徽省煤礦綠色低碳發(fā)展工程研究中心，安徽 淮南 232001)

0 引 言

煤炭開采在滿足國家能源需求的同時也導致嚴重的礦區(qū)地質(zhì)環(huán)境損害問題，比如礦區(qū)地表沉陷、地裂縫等地質(zhì)災害，土地損壞、植被退化、水體流失等礦區(qū)環(huán)境損害。煤炭開采引起巖層及地表變形是產(chǎn)生礦山地質(zhì)環(huán)境損害的根源。隨著國家供給側(cè)結(jié)構(gòu)性改革、我國能源戰(zhàn)略調(diào)整以及東部煤炭資源逐漸枯竭，西部地區(qū)已成為我國未來主要煤炭資源產(chǎn)地[1-3]。同時，我國采煤方法發(fā)展迅速，采煤方法已由以往的炮采、普采、綜采和綜放開采迅速向無人化、智能化開采方向變革[3-5]。這使得我國西部以大規(guī)模、高強度快速開采為主，開采引起的地表移動損害也呈現(xiàn)出變形速度快、損害程度深、波及范圍廣的特點[6-9]。由于西部生態(tài)環(huán)境脆弱，煤炭開采引起的地質(zhì)與環(huán)境損害問題成為重點關(guān)注和研究的焦點[10-13]。然而，西部礦區(qū)地表沉陷監(jiān)測數(shù)據(jù)較少，沉陷規(guī)律研究不深入，更沒有適合的沉陷參數(shù)。隨著國家對環(huán)境保護要求越來越嚴格，礦區(qū)環(huán)境損害監(jiān)測及評估已經(jīng)成為西部礦山企業(yè)面臨的巨大難題。在這種形勢下，如何快速、準確、全面地監(jiān)測西部高強度開采引起的地表沉陷與環(huán)境損害是解決問題的關(guān)鍵。而目前常規(guī)地表移動與變形監(jiān)測方法已經(jīng)無法適應我國新形勢下的礦山開采損害監(jiān)測需求[14-15]。隨著測繪技術(shù)的快速發(fā)展，如何充分利用測繪新技術(shù)，突破理論與技術(shù)瓶頸，建立新形勢下的礦山開采損害監(jiān)測理論、技術(shù)與方法，以適應西部煤炭大規(guī)模、高強度開采地質(zhì)環(huán)境損害監(jiān)測任務，是目前資源開發(fā)與生態(tài)環(huán)境保護亟須解決的重大問題。

總結(jié)現(xiàn)有礦山開采損害監(jiān)測技術(shù)的研究現(xiàn)狀，分析了各項監(jiān)測技術(shù)在礦區(qū)開采損害監(jiān)測中存在的問題、優(yōu)勢和不足?；赨AV攝影測量與InSAR技術(shù)各自的優(yōu)勢，提出了融合UAV/InSAR監(jiān)測礦山開采損害的關(guān)鍵理論與技術(shù)，并以內(nèi)蒙古鄂爾多斯王家塔煤礦為例進行了工程應用研究，取得良好的效果，研究結(jié)果可為工礦企業(yè)地表損害監(jiān)測和環(huán)境保護提供借鑒。

1 礦山開采損害監(jiān)測技術(shù)發(fā)展現(xiàn)狀

礦山開采損害監(jiān)測技術(shù)主要包括常規(guī)地表沉陷監(jiān)測技術(shù)、三維激光掃描技術(shù)、攝影測量技術(shù)(UAV攝影測量)和InSAR監(jiān)測技術(shù)等4個方面。

1.1 常規(guī)地表沉陷監(jiān)測技術(shù)

常規(guī)地表移動變形監(jiān)測方法，按照布站形式，可分為剖面線狀觀測站和網(wǎng)狀觀測站[14-15]。前者應用較廣，它是在工作面走向/傾向主斷面上布“點”成“線”組成觀測站；后者是多布設一些測點組成網(wǎng)格狀的觀測站，由于受地形、地物條件限制(如占用更多土地)，后者應用較少。2種觀測站獲取的均為有限的“點狀”數(shù)據(jù)。常規(guī)觀測站主要采用包括三角測量、導線測量(經(jīng)緯儀、全站儀)、水準測量(水準儀)以及GPS等常規(guī)測量方法；常規(guī)觀測站精度高、是目前獲取礦區(qū)地表形變最常用、最有效的手段，在礦山開采沉陷研究及實際工程中發(fā)揮了重要的作用。隨著研究的發(fā)展，出現(xiàn)一些不足：①尺度小，一般以單工作面為監(jiān)測對象，測量方法也決定了不便大范圍(大尺度)的作業(yè)；②僅以地表移動變形值為監(jiān)測內(nèi)容，無法獲得礦區(qū)生態(tài)環(huán)境損害信息(比如植被光譜、水體等)；③成本高、周期長(時效性差)、工作量大、需要埋設測點且難以長期保存等缺點；④只“點狀”觀測，數(shù)據(jù)及信息量少，不能完整反映沉陷盆地特征。這也導致了其已無法適應新形勢下西部高強度規(guī)模開采地表損害監(jiān)測任務。

1.2 三維激光掃描技術(shù)

可分為地面激光掃描技術(shù)、機載激光雷達和航天激光雷達掃描技術(shù)。其中的機載/航天激光雷達設備較貴，使用成本較高，對于礦山企業(yè)來說使用較少。使用較多的是地面激光掃描技術(shù)。包括作者在內(nèi)的部分學者[15-18]利用地面三維激光掃描技術(shù)建立地表移動觀測站，該方法可獲得豐富的“面狀”點云數(shù)據(jù)，但是在實際觀測過程中，需要分站掃描，每一站掃描范圍有限，搬站次數(shù)多，尤其是西部山區(qū)及丘陵地區(qū)，搬站困難；因此，該方法仍存在尺度小、成本高、效率低，只獲取地表變形數(shù)據(jù)，無法獲得礦區(qū)環(huán)境損害數(shù)據(jù)等問題，限制了其在礦區(qū)的推廣應用，尤其無法滿足西部礦區(qū)開采損害監(jiān)測任務。

1.3 攝影測量技術(shù)(UAV攝影測量技術(shù))

攝影測量技術(shù)在礦山應用由來已久，航空/航天攝影測量主要用于礦區(qū)地形圖測繪[13-15]，受制于成本及精度等因素的影響，很難推廣到礦區(qū)地表沉陷監(jiān)測中；在利用地面近景攝影測量監(jiān)測礦區(qū)地表形變方面，相關(guān)學者也做了一些嘗試[16-18]，實踐證明可行，然而受相機視角的限制，其觀測區(qū)域較小，效率較低，后期數(shù)據(jù)處理復雜等因素，導致其發(fā)展緩慢。隨著近年低成本、輕型UAV技術(shù)的發(fā)展，出現(xiàn)了UAV低空攝影測量技術(shù)，成為傳統(tǒng)航空攝影測量的有力補充，可快速獲取含有豐富地物光譜信息的多(高)光譜、高分辨率光學影像，適合多尺度的監(jiān)測任務，極大推動了攝影測量技術(shù)在礦區(qū)的應用。近年UAV攝影測量在礦區(qū)應用研究得到飛速發(fā)展，主要集中于：①礦區(qū)地表DEM建立及三維建模，研究表明UAV建立DEM精度能達到厘米～分米級；如宮傳剛等[19]除去植被影響后，DEM精度達到45 mm；田帥帥等[20]利用UAV建立采煤沉陷區(qū)地表DEM，修正后的高程中誤差為12.7 cm；②礦區(qū)地表穩(wěn)態(tài)及動態(tài)沉陷盆地監(jiān)測，求取地表沉陷參數(shù)；如侯恩科等[21-22]采用輕型UAV監(jiān)測寧夏羊場灣煤礦Y120212工作面采煤沉陷量，得到最大沉陷量達6.5 m；周大偉等[23]利用UAV對內(nèi)蒙古王家塔煤礦2S201工作面地表沉陷進行監(jiān)測，得到三期地表動態(tài)下沉盆地，并求取了可靠的開采沉陷參數(shù)。③礦區(qū)環(huán)境因素監(jiān)測，包括礦區(qū)地表裂縫監(jiān)測[24-25]、土地利用分類[26]，植被分類及參數(shù)反演[27]等。如：肖武等[28-29]利用無人機多光譜影像數(shù)據(jù)，研究了采煤沉陷對玉米葉綠素和生物量的影響。

UAV低空攝影測量技術(shù)可以同時采集到地表沉陷變形與礦區(qū)生態(tài)環(huán)境信息，其高程精度為分米級，其在礦區(qū)開采沉陷大變形監(jiān)測中具有優(yōu)勢(如Y120212工作面地面6.5 m下沉量)。然而，分米級精度決定了其無法精確監(jiān)測沉陷盆地邊界區(qū)域。

1.4 InSAR技術(shù)

InSAR集合成孔徑雷達技術(shù)(SAR)與干涉測量技術(shù)于一體，用來記錄地物的后向散射強度信息及相位信息，前者反映了地表屬性(含水量、粗糙度、地物類型等)，后者則蘊含了傳感器與目標物之間的距離信息[30]。InSAR通過對同一區(qū)域兩景SAR影像進行干涉處理(即相位差)即可提取地表高程信息或分離地表形變。作為InSAR技術(shù)的延伸，差分合成孔徑雷達干涉技術(shù)(D-InSAR)則是通過引入外部DEM或三軌/四軌差分實現(xiàn)了地表變形監(jiān)測[31-32]。D-InSAR易受時空失相干、大氣延遲和軌道誤差等因素影響，尤其時空失相干是InSAR技術(shù)在礦區(qū)沉陷監(jiān)測中的最大阻礙。煤礦開采沉降量大、下沉速率快、礦區(qū)地表覆蓋和地形起伏變化快等原因?qū)е翴nSAR局部區(qū)域(如沉陷盆地中心大變形區(qū)域)相位完全失相關(guān)，無法可靠地求取沉陷盆地中心區(qū)域的大變形，然而，InSAR技術(shù)可精確探測出沉陷盆地邊緣區(qū)域變形。

1.5 研究現(xiàn)狀分析

綜上，常規(guī)觀測站存在尺度小、僅能觀測地表形變(不能獲得環(huán)境因素)，成本高效率低且數(shù)據(jù)量少等不足；地面三維激光掃描技術(shù)雖然可獲得大量“面狀”點云數(shù)據(jù)，但是也存在尺度小、成本高效率低的問題；在我國煤炭開采向西部轉(zhuǎn)移，大規(guī)模、高強度開采成為主流的新形勢下，這些監(jiān)測方法無法適應多尺度開采損害監(jiān)測任務。

雖然UAV與InSAR技術(shù)在礦山開采變形監(jiān)測應用中也各有不足，單獨采用其中一種技術(shù)目前均無法完整地監(jiān)測整個地表沉陷盆地信息。然而，二者優(yōu)勢互補：UAV彌補了InSAR無法提取沉陷盆地大變形的問題，InSAR彌補了UAV在監(jiān)測沉陷盆地邊緣誤差大的問題；UAV獲取的多(高)光譜影像數(shù)據(jù)，在礦區(qū)環(huán)境因素監(jiān)測方面具有很大優(yōu)勢，彌補了SAR影像中缺少地物光譜特性的問題；兩種影像融合，可以獲得同時包涵地物后向散射強度、相位信息及地物光譜特性的影像，即可完整提取沉陷盆地整體形態(tài)，也可提高礦區(qū)環(huán)境損害監(jiān)測能力?；诖耍岢鋈诤蟄AV/InSAR技術(shù)進行礦區(qū)多尺度開采損害監(jiān)測，實現(xiàn)礦區(qū)關(guān)鍵地質(zhì)與環(huán)境損害信息精準提取，建立新形勢下集地表沉陷(地表形變)與環(huán)境損害為一體的綜合監(jiān)測理論、技術(shù)與方法體系，以適應西部大規(guī)模、高強度開采地質(zhì)環(huán)境損害監(jiān)測需求。

2 UAV/InSAR融合監(jiān)測礦區(qū)地表損害關(guān)鍵理論與技術(shù)

UAV/InSAR融合監(jiān)測礦區(qū)地表損害關(guān)鍵技術(shù)體系主要包括：UAV/InSAR多尺度礦山開采損害協(xié)同監(jiān)測技術(shù)；UAV/InSAR跨尺度融合理論與方法；基于融合數(shù)據(jù)的礦區(qū)開采沉陷與關(guān)鍵環(huán)境因素精細提取。首先，UAV/InSAR進行協(xié)同監(jiān)測(時間、空間和參數(shù)協(xié)同[33])，建立多尺度新型開采損害觀測站，包括不同尺度觀測站的布設模式(如工作面尺度、采區(qū)尺度和礦區(qū)尺度等)、觀測點合理布設和觀測時間間隔等；其次，UAV與InSAR獲取影像是跨尺度的異質(zhì)數(shù)據(jù)，顧及異質(zhì)數(shù)據(jù)的空間差異的情況下，充分利用各自的優(yōu)勢，建立跨尺度異質(zhì)遙感數(shù)據(jù)融合理論與方法；最后，基于融合數(shù)據(jù)精準提取地表沉陷變形和關(guān)鍵環(huán)境因素，建立基于融合數(shù)據(jù)的礦區(qū)開采沉陷與關(guān)鍵環(huán)境因素精細提取的理論與方法。

2.1 UAV/InSAR多尺度協(xié)同監(jiān)測技術(shù)

UAV/InSAR多尺度協(xié)同監(jiān)測獲取遙感影像是研究的技術(shù)支撐和數(shù)據(jù)基礎。主要包括2個方面：建立基于UAV/InSAR的多尺度的“空(InSAR)-天(UAV)-地(傳統(tǒng)觀測點)”一體化新型觀測站模式；InSAR時-空間分辨率約束下的UAV/InSAR協(xié)同觀測。

1)基于UAV/InSAR協(xié)同的多尺度觀測站的建立模式。收集目標礦區(qū)的地質(zhì)資料(巖層性質(zhì)、鉆孔柱狀圖、井上下對照圖等)及采礦資料(開采時間、開采速度、采深、采厚等)，結(jié)合開采時間和開采速度，優(yōu)選出合理的InSAR數(shù)據(jù)，包括確定極化方式、時間分辨率和空間分辨率。根據(jù)不同目的與服務對象，建立如下多尺度的開采損害觀測站形式：工作面尺度：采用“點-線-面”結(jié)合布設形式，在沉陷盆地邊界點、拐點和最大下沉點等關(guān)鍵位置設置少量觀測點(9～20個觀測點不等)，這些點也作為UAV低空攝影測量時的像控點，該觀測站的目的是求取開采沉陷參數(shù)、研究地表移動規(guī)律服務的。采區(qū)尺度和礦區(qū)尺度的“面狀”觀測站，采用UAV攝影測量進行“面狀”測量，規(guī)劃航攝范圍，獲得影像數(shù)據(jù)，前者為采區(qū)建(構(gòu))筑物損害、沉陷規(guī)律及環(huán)境損害評估服務，后者為礦區(qū)地質(zhì)環(huán)境災害監(jiān)測、環(huán)境評估及預警服務，如圖1所示。

圖1 多尺度觀測站建立模式Fig. 1 Multi-scale observation station establishment model

2)InSAR時-空分辨率約束下的UAV/InSAR協(xié)同觀測：在優(yōu)選出InSAR數(shù)據(jù)后，InSAR的時間和空間分辨率已定；UAV監(jiān)測周期(時間分辨率)及空間分辨率均可靈活調(diào)整。在UAV航攝測量時，根據(jù)InSAR的時空分辨率，對其航攝時間、地面分辨率、飛行高度、像片重疊度等航攝方案中所涉及的參數(shù)進行優(yōu)化設計，達到UAV/InSAR協(xié)同觀測的目的。分析UAV與InSAR時間分辨率的最佳關(guān)系；研究UAV/InSAR達到最佳融合效果對應的二者空間分辨率之間的關(guān)系；分析UAV地面空間分辨率的影響因素，通過UAV光學影像時間-空間分辨率的調(diào)整，從而實現(xiàn)二者時間協(xié)同、空間協(xié)同和參數(shù)協(xié)同的目的，為給后面數(shù)據(jù)融合和信息提取提供最佳的遙感影像數(shù)據(jù)。

2.2 UAV/InSAR融合理論與方法

UAV獲得光學影像分辨率為5～10 cm，InSAR獲得的衛(wèi)星雷達影像分辨率為3～100 m，2種不同平臺傳感器獲取的數(shù)據(jù)，在點-線-面觀測尺度上差異較大，二者融合是屬于跨尺度異質(zhì)遙感數(shù)據(jù)融合，目前的研究成果相對較少。如果要二者融合，需要對以下兩方面做研究：①顧及異質(zhì)數(shù)據(jù)的時空差異，挖掘它們之間的相關(guān)性與映射關(guān)系，建立UAV/InSAR影像的精確配準策略；②配準后兩類影像的像素級和特征級融合，分析融合結(jié)果的評價指標體系，建立融合精度和質(zhì)量定量評價方法，形成UAV/InSAR跨尺度融合理論與方法體系。

UAV/InSAR融合方法分為像素級融合和特征級融合；根據(jù)不同的目的采用不同的融合方法。

1)提取礦區(qū)地表裂縫及生態(tài)環(huán)境因素，采用像素級融合。像素級融合方法又具體分為IHS變換法、主成分變換法(PCA)和加權(quán)平均法等方法，各方法實施具體如下：①IHS變換法影像融合。IHS變換是將圖像RGB彩色空間的紅、綠和藍變換為IHS彩色空間的明度、色調(diào)和飽和度。變換方法為：將UAV影像進行IHS變換，再利用InSAR圖像的灰度代替I分量，進行IHS逆變換得到融合圖像，如圖2所示。②主成分分析法(PCA)影像融合。PCA是在統(tǒng)計特征基礎上的一種常用降維變換融合方法，其基本原理是將原始存在較高相關(guān)性的N維特征的高維數(shù)據(jù)，通過PCA變換投影到M維正交子空間(M≤N)，用新子空間中少數(shù)幾個不相關(guān)的、方差大的、信息量大的特征來表示原始數(shù)據(jù)信息。主要步驟：首先對UAV光學影像進行PCA變換，得到n個獨立主分量；再對高分辨率InSAR圖像進行拉伸處理，使其灰度的均值與方差和UAV影像的第i(i=1, 2, …,n)主分量一致，替換該主分量；最后進行PCA逆變換，得到融合影像，如圖3所示。③加權(quán)平均法影像融合。該方法的基本原理是對經(jīng)配準的兩幅影像中各對應像素做加權(quán)平均處理，融合成一幅新圖像。設A和B分別為UAV像片與SAR圖像，圖像大小為M×N,F為融合后的圖像，則加權(quán)平均過程可以表示為：

F(m,n)=w1A(m,n)+w2B(m,n)

(1)

式中：m為圖像中像素的行號，m=1,2,…,M；n為圖像中像素的列號，m=1,2,…,N；w1、w2為權(quán)重，通常w1+w2=1。

圖2 IHS變換融合方法流程Fig. 2 Flow of IHS transformation and fusion method

圖3 PCA變換融合方法流程Fig. 3 Flow of PCA transform fusion method

像素級融合能夠充分應用原始數(shù)據(jù)中包含的數(shù)據(jù)和信息量，獲得更高質(zhì)量的含SAR強度、相位信息和UAV光譜信息的數(shù)據(jù)，為礦區(qū)地表植被、土地利用和地表裂縫提供服務。

2)提取礦區(qū)地表沉陷變形值，采用UAV/InSAR的特征級融合方法。特征級融合是指從不同數(shù)據(jù)中分別提取相關(guān)特征，通過對提取的特征進行融合處理，生成新的特征或特征矢量，以便于后續(xù)地物解譯。采用2種方法進行UAV與InSAR特征級融合。①第1種融合方法采用相位融合：由于礦區(qū)開采沉陷快速大變形的特點，造成InSAR解纏過程中相位“整周數(shù)”的丟失，使得監(jiān)測結(jié)果出現(xiàn)較大偏差(最大下沉量嚴重偏小)；利用UAV獲得礦區(qū)大變形沉陷信息修正SAR影像處理過程中丟失的“整周數(shù)”相位，以提高InSAR技術(shù)監(jiān)測礦區(qū)沉陷量級。具體步驟如下：首先，將獲取的N景SAR影像進行干涉和解纏處理，獲取N-1幅解纏相位，通過“疊加”和重纏繞處理獲得復數(shù)纏繞相位；其次，UAV航攝及數(shù)據(jù)處理獲得礦區(qū)DEM，兩期DEM相減獲得礦區(qū)沉陷大變形；再次，將UAV獲得“大形變”根據(jù)雷達波長轉(zhuǎn)換為絕對干涉相位，并轉(zhuǎn)換為纏繞相位。最后，通過InSAR獲得纏繞相位與UAV獲得的含有大變形的纏繞相位融合，修正SAR相位中的“整周數(shù)”，得到含有大變形的融合相位，通過解算處理，最終得到完整、準確的礦區(qū)沉陷盆地變形。②第2種融合方法的基本思想：首先利用UAV攝影測量和InSAR技術(shù)分別提取同一個區(qū)域的地表沉陷盆地信息，通過邊界提取算法提取UAV沉陷盆地精確的外邊緣和InSAR監(jiān)測的沉陷盆地內(nèi)邊緣區(qū)域，對提取的邊界進行融合處理，得到完整、精確的沉陷盆地，如圖4所示。

圖4 UAV與InSAR沉陷盆地融合方法思路Fig.4 Idea of UAV and InSAR subsidence basin fusion method

2.3 礦山開采損害信息精準提取

礦山開采損害可分為礦山地質(zhì)損害與礦山生態(tài)環(huán)境損害，以地表沉陷變形(包括下沉和水平移動)與地表裂縫為代表的礦山地質(zhì)損害和以土地、植被與地表水體等為關(guān)鍵因素代表的礦山環(huán)境損害。從3個方面開展研究：①地表沉陷變形的獲取，研究UAV/InSAR融合精確獲取礦區(qū)動態(tài)及穩(wěn)態(tài)沉陷盆地的方法以及如何提取礦區(qū)水平移動值；②礦區(qū)地表裂縫的提取，研究分別從融合影像中提取地表裂縫的方法；③研究礦區(qū)地表土地利用、植被覆蓋及水體等關(guān)鍵環(huán)境因素的變化檢測方法，對比分析性能及優(yōu)缺點，優(yōu)選最佳提取方法。從UAV/InSAR融合數(shù)據(jù)中精準提取礦山地表形變與環(huán)境損害信息，建立面向礦山損害監(jiān)測應用的信息提取理論與方法。包括3個方面：UAV/InSAR融合地表沉陷變形精確獲取(下沉和水平移動)；UAV/InSAR融合礦區(qū)地表裂縫的提取方法；UAV/InSAR融合礦區(qū)關(guān)鍵環(huán)境因素變化檢測方法。

1)UAV/InSAR融合地表沉陷變形精確獲取，包括地表下沉盆地的提取和地表水平移動的提取。①下沉盆地的提取方法。地表下沉盆地提取有2種方法：方法1的思路是利用UAV獲得礦區(qū)大變形修正InSAR影像處理過程中相位的“整周數(shù)”，從而獲得完整的下沉盆地；方法2的思路是利用UAV和InSAR技術(shù)分別監(jiān)測礦區(qū)沉陷盆地，將UAV大變形區(qū)域與InSAR的邊界區(qū)域融合，從而得到完整沉陷盆地數(shù)據(jù)，如圖4所示。求取的下沉盆地，根據(jù)開采時間不同，可以分為動態(tài)下沉盆地和穩(wěn)態(tài)下沉盆地。②水平移動提取方法。方法1：可以采用UAV航攝的像控點求取水平移動值。方法2：由于UAV攝影測量和InSAR測量時無固定的目標點，結(jié)合UAV/InSAR像素級融合后的影像，通過融合影像精確提取礦區(qū)地表穩(wěn)定的特征點、線(如道路)或面(建筑物頂面)，求取礦區(qū)水平移動。以提取線的中點和面的重心為基準，在精確配準后，以第一期數(shù)據(jù)為準，計算第i期提取的點、線(中點)或面(重心)與第一期的距離，多期距離相減求取地表水平移動。

2)UAV/InSAR融合礦區(qū)地表裂縫提取方法。煤礦開采引起的地表裂縫的發(fā)育規(guī)律為(圖5)：裂縫均位于工作面外側(cè)拉伸變形區(qū)，均由拉伸變形和正曲率變形引起，其基本形態(tài)是以工作面為中心的圓弧形或橢圓形。因此，煤礦地表裂縫的發(fā)育區(qū)域和形態(tài)特征均可作為先驗知識，基于此，建立先驗知識模型，以含有光譜特性的UAV光學影像為主，含有地物強度信息的InSAR影像為輔助修正，通過二者融合提高裂縫提取精度；利用UAV/InSAR融合影像數(shù)據(jù)，通過掩膜技術(shù)去除干擾，采用特征提取算子提取礦山地表裂縫，形成基于先驗知識模型的礦區(qū)地表裂縫提取方法。具體步驟：首先，根據(jù)地裂縫與其他地物在光譜、幾何特征和區(qū)域特征上的差異，采用掩膜技術(shù)去除部分干擾；然后，基于掩膜結(jié)果，利用邊緣檢測算子(如Canny算子和SIFT算子)提取圓弧形或橢圓形線性特征，去除不具有該線性特征的地表形態(tài)；最后，根據(jù)裂縫提取結(jié)果，對于散亂小裂縫等進行人工修正，得到符合發(fā)育規(guī)律的礦區(qū)地表裂縫信息。

圖5 煤礦開采地表裂縫動態(tài)發(fā)展過程和基本形態(tài)Fig.5 Dynamic development process and basic form of surface cracks in coal mining

3)UAV/InSAR融合礦區(qū)關(guān)鍵環(huán)境因素變化提取。以UAV/InSAR像素級融合影像為基礎，通過像元級和對象級變化監(jiān)測。①像元級變化檢測方法。提取多個時段的影像的亮度、NDVI特征。采用差值法得到差異圖，該差異圖集中了絕大部分的變化信息, 而濾除了影像中相同的背景部分；在差異圖中進行地物分類，提取礦區(qū)關(guān)鍵環(huán)境因素的變化信息。②面向?qū)ο蟮淖兓瘷z測方法。UAV獲得光學影像具有豐富的地物光譜特性；由于不同地物散射特性的差異, 在SAR圖像上會表現(xiàn)出不同的亮度(灰度)和紋理；因此，灰度特征和紋理特征是SAR影像用于分類的主要特征。面向?qū)ο蟮淖兓瘷z測方法是根據(jù)同一對象同質(zhì)性進行分類；利用UAV/SAR融合影像中綜合利用地物的光譜特征、灰度特征和紋理特征進行同質(zhì)點選擇，利用面向?qū)ο蟮姆诸惼鬟M行地物分類，分類器采用支持向量機(SVM)分類算法以及深度學習算法實現(xiàn)，從分類結(jié)果中統(tǒng)計礦區(qū)關(guān)鍵環(huán)境因素(比如土地、植被和水體等)的面積、形狀和數(shù)量等基礎屬性信息。通過多期數(shù)據(jù)進行對比差值運算獲得礦區(qū)關(guān)鍵環(huán)境因素的變化信息，為礦區(qū)環(huán)境損害研究提供數(shù)據(jù)支持。需要說明，雖然均是利用融合數(shù)據(jù)提取地物，不同地物提取方法的側(cè)重點略有不同，比如植被的提取以含有地物光譜特征的UAV光學影像為主，InSAR影像的強度信息作為輔助修正，以提高提取精度。而對于建筑土地或者水體，由于其反射強度信息占主，因此需以InSAR影像為主，UAV光學影像為輔助修正。

3 工程應用

針對上述UAV/InSAR融合監(jiān)測礦區(qū)損害理論體系，以內(nèi)蒙古鄂爾多斯王家塔煤礦地表沉陷變形監(jiān)測為例，主要以地表沉陷為監(jiān)測內(nèi)容，說明上述方法在西部礦區(qū)的應用。

3.1 研究區(qū)域概況

以位于內(nèi)蒙古自治區(qū)鄂爾多斯市的王家塔煤礦2S201工作面為研究區(qū)域，工作面如圖6所示。該工作面傾向長度260 m，走向長度1 253 m，煤層平均傾角為2°，平均采深約為200 m，平均采厚為3.26 m。該工作面于2018年7月11日正式回采，2018年10月25日收尾。

3.2 工作面尺度“點-線-面”結(jié)合模式觀測站

采用工作面尺度的“點-線-面”結(jié)合模式建立觀測站。在研究區(qū)域范圍內(nèi)，在沉陷盆地邊界點、拐點和最大下沉點等關(guān)鍵位置布設了15個測點(也作為UAV攝影測量的像控點和檢查點)；根據(jù)SAR衛(wèi)星經(jīng)過研究區(qū)域的時間(表1)，選定InSAR數(shù)據(jù)為2018-06-11—2018-09-04之間8景Sentinel-1A的IW SLC影像數(shù)據(jù)，影像參數(shù)見表1。2018-06-09和2018-09-04兩個時間采用UAV攝影測量對研究區(qū)域進行了兩次觀測，UAV測量參數(shù)及數(shù)據(jù)基本信息見表2。同時，采用全站儀/水準測量的方法觀測地表測點。因此，同時獲得該區(qū)域的UAV光學影像、InSAR影像和主斷面關(guān)鍵點水準數(shù)據(jù)，建立了工作面尺度的 “點-線-面”結(jié)合模式觀測站即為傳統(tǒng)觀測站與UAV/InSAR測量新技術(shù)相結(jié)合的新型觀測模式，具體如圖6所示。

表1 Sentinel-1A影像相關(guān)參數(shù)及S1A干涉對信息

表2 UAV攝影測量航攝信息

圖6 研究區(qū)域、全站儀監(jiān)測點及航帶分布Fig.6 Excavation view of study area and layout of image con-trol points

3.3 UAV/InSAR數(shù)據(jù)融合方法

本次主要提取礦區(qū)地表沉陷變形值，采用了UAV/InSAR的特征級融合方法中的第2種融合思路(圖4)，具體步驟如圖7所示。

圖7 UAV/InSAR沉陷盆地特征級融合具體步驟Fig. 7 Specific steps of UAV/InSAR subsidence basin feature-level fusion

1)利用UAV攝影測量和InSAR技術(shù)分別提取該工作面地表沉陷盆地。

2)確定沉陷盆地邊界區(qū)域和盆地中心區(qū)域：由于已知InSAR在外邊緣區(qū)域精度較高，UAV可以監(jiān)測沉陷盆地中心大變形區(qū)域。因此，根據(jù)煤礦開采沉陷學主要影響半徑的定義[7]，選取0.16Wmax和0.84Wmax(Wmax為最大下沉值)的點為邊界，小于0.16Wmax的沉陷區(qū)域采用InSAR監(jiān)測數(shù)據(jù)，大于0.84Wmax的沉陷區(qū)域采用無人機監(jiān)測數(shù)據(jù)。

3)確定0.16Wmax和0.84Wmax之間的區(qū)域：通過概率積分法擬合迭代選取，即首先固定盆地邊界和中心實測數(shù)據(jù)，然后通過概率積分模型進行擬合循環(huán)迭代，每次迭代一次，提取與擬合下沉盆地(即為擬合曲面)差值最小的格網(wǎng)點，以擬合誤差中誤差最小作為目標函數(shù)，通過多次擬合迭代可以獲取0.16Wmax和0.84Wmax之間區(qū)域最佳監(jiān)測數(shù)據(jù)。再聯(lián)合提取的邊緣區(qū)和盆地中心區(qū)可以得到完整的融合沉陷盆地。

3.4 UAV/InSAR融合提取地表沉陷盆地

1)UAV解算下沉盆地。DEM的構(gòu)建是生成地表沉陷盆地的關(guān)鍵步驟，UAV數(shù)據(jù)處理得到該區(qū)域不同期DEM，多期DEM相減獲得監(jiān)測區(qū)域的地表沉陷。通過對UAV原始影像進行畸變校正、相對定向、絕對定向、光束法區(qū)域網(wǎng)平差以及影像密集匹配等處理步驟得到監(jiān)測區(qū)域的點云，對點云進行濾波、去噪和壓縮處理，得到去除非地面點之后的點云。根據(jù)數(shù)字高程模型的生成算法，建立20 m分辨率方格網(wǎng)，每個格網(wǎng)內(nèi)僅保留一個最低點作為地面點，通過內(nèi)插格網(wǎng)節(jié)點高程構(gòu)建DEM[12-13]，圖8所示為建立的第1期DEM和第2期DEM。兩期DEM相減獲得該區(qū)域由2018-06-09—2018-09-04的動態(tài)地表移動盆地，如圖9a所示，圖9b和圖9c為走向(A—A′)和傾向(B—B′)剖面下沉曲線。UAV監(jiān)測到該期間的最大下沉值為2 487 mm，通過與水準測量結(jié)果對比，評估UAV下沉盆地的下沉值中誤差為81 mm；UAV沉陷盆地誤差對于沉陷盆地邊界影響較大，對盆地中間大變形區(qū)域影響相對較??；剖面下沉曲線也可看出，相比于沉陷盆地中間區(qū)域(區(qū)域2)，沉陷盆地邊界區(qū)域(區(qū)域1)波動相對更大，在后面求參數(shù)時主要影響角正切的誤差偏大也能說明該問題。

圖8 UAV攝影測量解算的兩期DEM和DOMFig.8 Two phases of DEM and DOM solved by UAV photogrammetry method

圖9 UAV攝影測量解算結(jié)果Fig.9 Result of UAV photogrammetvy

2)InSAR解算沉陷盆地。通過對8景InSAR影像數(shù)據(jù)處理得到圖10a中的7個干涉像對的形變圖，干涉對信息見表1。將圖中由D-InSAR解算得到的各干涉像對的形變量進行疊加解算，得到2S201工作面觀測時間從2018-06-11—2018-09-03的動態(tài)地表沉降量，動態(tài)沉陷盆地如圖10b所示，圖10c為走向剖面(A—A′)下沉曲線。根據(jù)水準測量結(jié)果，該時段內(nèi)地表最大下沉值為2 415 mm。而InSAR監(jiān)測的最大下沉值僅為110 mm，遠小于實際下沉量。因此，InSAR在監(jiān)測礦區(qū)形變梯度大時，受時空失相干因素的影響，無法準確獲取大變形區(qū)域沉降。

圖10 D-InSAR解算結(jié)果Fig. 10 Result of D-InSAR

3)UAV/InSAR融合獲得精確下沉盆地。根據(jù)圖7中的融合步驟對InSAR和UAV監(jiān)測的下沉盆地進行融合，得到融合下沉盆地如圖11a所示，做走向(A-A’)和傾向(B-B’)剖面線，得到融合盆地主斷面下沉曲線，如圖11b和圖11c所示。從圖11b和圖11c中可看出，邊界區(qū)域的相對較平滑，通過與水準測量結(jié)果對比，邊界區(qū)域精度達到6 mm，融合下沉盆地整體中誤差為56.3 mm。與UAV沉陷盆地相比，融合下沉盆地整體精度和邊界區(qū)域精度均有提高。

圖11 UAV/InSAR融合結(jié)果Fig. 11 Result of UAV/InSAR fusion method

3.5 基于融合數(shù)據(jù)的沉陷參數(shù)反演

首先，利用UAV沉陷盆地數(shù)據(jù)進行參數(shù)反演，采用概率積分法結(jié)合時間函數(shù)的動態(tài)求參方法[14,39]，擬合中誤差為181 mm，占最大下沉值的7.3%，實測與預計擬合效果如圖12a所示；其次，采用融合下沉盆地數(shù)據(jù)進行求參，擬合中誤差為46.5 mm，占最大下沉值的1.9%，擬合效果如圖12b所示；兩者求參結(jié)果見表3。通過與水準測量求參結(jié)果相比：

1)UAV下沉盆地求參結(jié)果中，下沉系數(shù)q與水準測量數(shù)據(jù)求參結(jié)果相符，相對誤差僅為1.4%，主要影響角正切tanβ與之相比偏差較大，偏離約20%。tanβ是開采深度H與主要影響半徑r的比值(tanβ=H/r)，在充分開采的條件下，r的范圍處在沉陷盆地邊界區(qū)域，由于UAV下沉盆地邊界區(qū)域的誤差較大，導致求出的主要影響角正切值誤差偏大。

2)融合數(shù)據(jù)求取的下沉系數(shù)與水準求參結(jié)果一致，主要影響角正切與水準求參結(jié)果的相對誤差僅為5%。融合數(shù)據(jù)參數(shù)反演的結(jié)果能夠更好地預計礦區(qū)地表沉陷。

3)與單獨UAV下沉盆地求參結(jié)果相比，融合下沉盆地求出的參數(shù)更準確。由于InSAR解算的下沉盆地的邊界區(qū)域精度較高，所以用融合求參很好地解決了單一UAV求參時主要影響角正切偏差較大的問題。

表3 無人機數(shù)據(jù)求取的測量參數(shù)

圖12 利用UAV下沉盆地和融合下沉盆地求參擬合效果Fig.12 Fitting effect of measured subsidence and predicted subsidence by parameters inversion with UAV subsidence basin and the fused subsidence basin

4 結(jié) 論

1)提出了融合UAV攝影測量/InSAR監(jiān)測礦山開采損害的思路。論述了基于UAV/InSAR融合的新型多尺度地表損害觀測站的建立模式，分析了InSAR時-空間分辨率約束下的UAV/InSAR協(xié)同觀測方法；探討了UAV/InSAR跨尺度異質(zhì)遙感數(shù)據(jù)融合及礦區(qū)地質(zhì)環(huán)境損害信息精細提取方法。

2)以內(nèi)蒙古王家塔煤礦為例進行了應用，結(jié)果表明：UAV測量誤差對于沉陷盆地邊界影響較大；與UAV沉陷盆地相比，InSAR/UAV融合下沉盆地整體精度和邊界區(qū)域精度均更高。與單獨UAV求參相比，融合下沉盆地數(shù)據(jù)求出的參數(shù)更準確，主要影響角正切的相對誤差僅為5%，融合數(shù)據(jù)很好解決了單一UAV求參時tanβ誤差偏大的問題。

3)實際工程應用表明了UAV/InSAR融合技術(shù)在礦區(qū)地表沉陷變形監(jiān)測中具有顯著的優(yōu)勢；UAV彌補了InSAR無法提取沉陷盆地大變形的問題，InSAR彌補了UAV在監(jiān)測沉陷盆地邊緣誤差大的問題，且適合礦區(qū)多尺度的監(jiān)測任務；可為我國西部煤礦開采地表損害監(jiān)測提供技術(shù)支撐。