免像控?zé)o人機(jī)攝影測量開采沉陷監(jiān)測方法研究

廉旭剛 王 站 劉曉宇 員鴻燕3

（1.太原理工大學(xué)礦業(yè)工程學(xué)院，山西 太原 030024；2.山西冶金巖土工程勘察有限公司第一分公司，山西 臨汾 041000；3.山西工程科技職業(yè)大學(xué)交通工程學(xué)院，山西 晉中 030619）

礦區(qū)開采沉陷區(qū)變形及生態(tài)環(huán)境監(jiān)測數(shù)據(jù)是礦山生態(tài)修復(fù)治理的重要依據(jù)[1]。采用水準(zhǔn)儀、全站儀、GNSS等傳統(tǒng)開采沉陷監(jiān)測手段費(fèi)時費(fèi)力，且單點(diǎn)監(jiān)測方式適用性不強(qiáng)；InSAR監(jiān)測雖然具有監(jiān)測覆蓋大、精度高等優(yōu)點(diǎn)，但其沉降監(jiān)測梯度有限[2]。隨著無人機(jī)技術(shù)及其搭載傳感器的不斷成熟，因其靈活性、效率高、高時空分辨率、低成本、操作方便等優(yōu)勢[3-4]，使得在礦區(qū)、農(nóng)業(yè)、林業(yè)監(jiān)測等領(lǐng)域得到了應(yīng)用，并取得了顯著成效[5-6]。近年來，在生態(tài)修復(fù)方面，一些學(xué)者采用無人機(jī)遙感技術(shù)開展了采煤沉陷區(qū)的耕地及積水區(qū)范圍提取[7-8]、耕地質(zhì)量評價[9]、農(nóng)作物生長狀態(tài)評估[10-11]、植被分類[12]等研究；在采動地面災(zāi)害監(jiān)測方面，無人機(jī)低空遙感手段被用于采動滑坡、崩落、地裂縫[13-14]、塌陷坑[15]及災(zāi)害預(yù)測[16]研究。

針對傳統(tǒng)開采沉陷監(jiān)測方法的不足，近年來研究人員采用無人機(jī)搭載激光雷達(dá)和高分辨率相機(jī)進(jìn)行礦區(qū)地面沉降監(jiān)測，并且取得了一定成效。無人機(jī)在礦區(qū)的應(yīng)用雖然處于起步階段，但在速度、規(guī)模和服務(wù)范圍方面發(fā)展迅速[17]，ESPOSITO等[18]、?WIAKALA 等[19]、葛林林等[20]、高冠杰等[21]、周大偉等[22]采用無人機(jī)攝影測量技術(shù)對煤礦開采地表沉陷進(jìn)行了監(jiān)測，取得了一定成效；湯伏全等[23]、張永庭等[24]基于無人機(jī)LiDAR技術(shù)對榆神礦區(qū)采煤沉陷建模方法進(jìn)行了改進(jìn)。本研究在上述成果的基礎(chǔ)上，以山西某礦區(qū)為例，采用免像控?zé)o人機(jī)攝影測量技術(shù)進(jìn)行開采沉陷監(jiān)測，進(jìn)一步提高該方法的自動化程度。

1 研究區(qū)概況

研究區(qū)位于太行山北段西側(cè)劉備山的南麓低中山區(qū)。地表地形復(fù)雜，溝谷縱橫，屬大陸性半干旱氣候，降雨主要集中在7—9 3個月，占全年總降水量的71%～91%。如圖1所示，相向開采的工作面A開采15#煤層，采深400～545 m，采厚6.5 m，長度1 200 m，寬度200 m，推進(jìn)速度3.3 m/d；工作面B開采8#煤層，采深336～506 m，采厚2.9 m，長度1 190 m，寬度220 m，推進(jìn)速度5.2 m/d。兩工作面水平方向距離32 m，垂直間距平均51 m，均采用長壁式一次采全高采煤方法，頂板全部垮落。上覆巖層以泥巖、砂質(zhì)泥巖為主，砂巖厚度為30～45 m，黃土層厚度5～10 m。本研究在圖1中布設(shè)了a、b、c3條觀測線驗(yàn)證無人機(jī)攝影測量方法的精度。

2 無人機(jī)攝影測量數(shù)據(jù)采集及處理

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

數(shù)據(jù)采集使用某型號無人機(jī)航測系統(tǒng)，搭載索尼A6000相機(jī)，有效像素2 430萬，航向和旁向重疊分別為80%和60%，采用仿地飛行，飛行速度14 m/s，地面像素分辨率為4 cm。在2019年9月、2019年11月、2020年1月、2020年5月、2020年9月進(jìn)行了5次數(shù)據(jù)采集，共采集全彩色照片3 400幅。地表觀測線a、b、c的測點(diǎn)間距為30 m，采用全站儀進(jìn)行了4次（對應(yīng)于無人機(jī)前4次采集時間）沉降觀測。

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

數(shù)據(jù)處理流程見圖2，采用實(shí)時動態(tài)載波相位差分技術(shù)（RTK）融合GNSS動態(tài)后處理差分技術(shù)（PPK）輔助空三解算，提高空三加密的精度，以實(shí)現(xiàn)免像控的目的。盡管RTK可以實(shí)現(xiàn)高精度的單點(diǎn)定位，但對于快速飛行的無人機(jī)而言易造成衛(wèi)星失鎖現(xiàn)象，融合PPK技術(shù)可以通過反向卡爾曼濾波，解決失鎖時間段內(nèi)的歷元數(shù)據(jù)，從而提高固定速率和定位精度。

采用無人機(jī)攝影測量技術(shù)進(jìn)行地表沉降監(jiān)測的關(guān)鍵是獲取高精度數(shù)字高程模型（DEM），而實(shí)際的地表存在建筑物及植被等非地面點(diǎn)，即為數(shù)字表面模型（DSM）。從DSM到DEM需要進(jìn)行點(diǎn)云濾波，點(diǎn)云濾波的原始數(shù)據(jù)即為該區(qū)域的多視影像密集匹配點(diǎn)云。使用Terrascan軟件中的自動濾波工具獲得初始濾波結(jié)果，并通過手動分類工具對初始結(jié)果進(jìn)行改進(jìn)，以保證DEM質(zhì)量。通過可視化點(diǎn)云的橫截面手工檢測相應(yīng)的點(diǎn)云，并糾正錯誤分類的點(diǎn)。研究區(qū)域局部的點(diǎn)云濾波前后效果對比如圖3所示。

3 開采沉陷結(jié)果分析

3.1 動態(tài)沉陷盆地及其精度分析

選取一期點(diǎn)云數(shù)據(jù)濾波后進(jìn)行插值算法對比試驗(yàn)，分別使用克里金插值、反距離權(quán)重插值、自然鄰點(diǎn)插值和線性插值4種方法繪制了等值線圖，如圖4所示。可見，自然鄰點(diǎn)插值和線性插值在圖像邊界處對于地形表達(dá)的真實(shí)性不及克里金插值法和反距離權(quán)重插值法。

為了進(jìn)一步驗(yàn)證各種插值方法的精度，在實(shí)地使用RTK采集檢查點(diǎn)的三維坐標(biāo)，以此高程數(shù)據(jù)作為測量的真值與利用平面坐標(biāo)查找的等值線插值高程進(jìn)行比較，進(jìn)而計(jì)算精度[25-27]。通過插值比較得出，利用不同插值方法得出的高程不同，其中克里金插值和反距離權(quán)重插值的高程較接近實(shí)際觀測值，可以確定這兩種插值方法的數(shù)據(jù)比較可靠。但在復(fù)雜地區(qū)且無人機(jī)影像匹配密集點(diǎn)云數(shù)據(jù)量較大的情況下，克里金插值計(jì)算量較大，較為耗時。綜上分析，本研究采用反距離權(quán)重插值法。

對點(diǎn)云濾波之后得到的地面點(diǎn)云統(tǒng)一進(jìn)行反距離加權(quán)插值處理得到各期DEM，將2019年9月首次無人機(jī)觀測獲得的DEM分別減去2019年11月、2020年1月、2020年5月和2020年9月觀測獲得的DEM，得到如圖5所示的工作面開采引起的地表下沉盆地的動態(tài)發(fā)展情況，結(jié)果分別對應(yīng)于圖5（a）、圖5（b）、圖5（c）和圖5（d）。在監(jiān)測時段內(nèi)，兩個工作面及附近地表均有不同程度的沉降，隨著工作面的推進(jìn)，A-15#和B-8#工作面地表最大下沉值逐漸增大，下沉盆地范圍不斷擴(kuò)大。B-8#和A-15#兩個工作面沉陷區(qū)域分別向南和北方向延伸，與工作面開采方向一致。由圖 5（c）到圖 5（d）的過程中，在A-15#工作面中心最大下沉值和整體沉陷盆地范圍變化不明顯。

本研究通過無人機(jī)相應(yīng)時間段DEM的高程減去兩個時段相應(yīng)點(diǎn)的全站儀高程，計(jì)算無人機(jī)沉降盆地?cái)?shù)據(jù)的均方誤差（M），公式為

式中，n為用于對比的點(diǎn)數(shù)，個；Δ為DEM與全站儀沉陷值的差值，m。

計(jì)算結(jié)果表明：沉降值的最大均方誤差為0.17 m，最小均方誤差為0.07 m，平均均方誤差為0.15 m（表1）。

不同時間段內(nèi)全站儀與無人機(jī)獲取的下沉曲線對比見圖6，其中不連續(xù)部分是由于部分監(jiān)測點(diǎn)遭受破壞而未獲得全站儀觀測數(shù)據(jù)所致。

由圖6可知：圖6（b2）、圖6（b3）、圖6（c3）中左側(cè)的矩形下沉量較小的區(qū)域，無人機(jī)數(shù)據(jù)跳點(diǎn)較多，與全站儀實(shí)測數(shù)據(jù)偏差較大，在圖6（b2）、圖6（b3）、圖6（c2）、圖6（c3）中右側(cè)虛線矩形區(qū)域中最大下沉值附近偏差減小，兩條曲線趨勢也較為接近。傾向b線和c線下沉值與下沉范圍隨著工作面的推進(jìn)下沉幅度逐步增大，呈現(xiàn)兩側(cè)平緩、邊緣較陡、底部較狹窄的“V”型曲線，走向a線由于對稱布設(shè)半條，未能全部顯示工作面的下沉曲線發(fā)展情況。

3.2 誤差來源及分析

統(tǒng)一點(diǎn)云濾波算法和插值算法所構(gòu)建的無人機(jī)沉陷DEM仍然存在明顯誤差，綜合圖5和圖6可知，無人機(jī)沉陷監(jiān)測中DEM誤差主要來源于以下幾方面：

（1）無人機(jī)容易受到一些內(nèi)部或外部因素的影響。儀器傳感器量化過程中帶來的系統(tǒng)誤差以及天氣狀況對飛行器姿態(tài)和成像質(zhì)量的影響產(chǎn)生的誤差會影響無人機(jī)數(shù)據(jù)質(zhì)量。因此每次無人機(jī)采集數(shù)據(jù)前都需要考慮天氣情況，選擇空氣潔凈、光照充足的時間段，選用成像質(zhì)量好的相機(jī)進(jìn)行航攝，數(shù)據(jù)后處理解算差分POS數(shù)據(jù)提高定位精度。

（2）濾波算法未去除的非地面點(diǎn)引起的模型誤差。點(diǎn)云濾波算法的局限性導(dǎo)致非地面點(diǎn)不能完全去除，以及數(shù)據(jù)采集覆蓋不同的季節(jié)，研究區(qū)植物的生長狀況會有差異，植被茂密時會降低所獲取DEM的精度?？赏ㄟ^多項(xiàng)式擬合或鄰域平均等方法減少噪聲來降低該類誤差，也可以選取冬季或春季采集數(shù)據(jù)，盡量減小植被對無人機(jī)數(shù)據(jù)精度的影響。

（3）采用點(diǎn)云內(nèi)插算法生成格網(wǎng)DEM引起的模型誤差。點(diǎn)云內(nèi)插會造成精度損失，尤其是在地形變化復(fù)雜和點(diǎn)云密度較低的地區(qū)，可通過改進(jìn)點(diǎn)云插值算法或幾種插值算法優(yōu)選減小這類誤差。

（4）影像分辨率的限制。影像分辨率越高，影像密集匹配和空中三角測量過程中精度就會越高。但是要求無人機(jī)降低飛行高度、減小飛行速度及設(shè)置更小的地面采用距離，會導(dǎo)致總體航攝效率會下降，增加無人機(jī)飛行的安全風(fēng)險，可通過適度提高影像地面分辨率來解決。

3.3 主斷面沉降特征分析

為了更直觀地反映本研究沉陷DEM監(jiān)測情況，基于2019-11—2019-09、2020-01—2019.09和2020-05—2019-09 3期沉陷DEM提取A-15#和B-8#兩個工作面走向主斷面下沉剖面數(shù)據(jù)，對3期數(shù)據(jù)進(jìn)行Savitzky-Golay平滑去噪[28]。該濾波器是一個基于多項(xiàng)式擬合方法的低通濾波器，相對于低頻信號變化 緩慢，波形平滑，相對于高頻信號變化非?？?，波形有突變。無人機(jī)數(shù)據(jù)中存在因飛行過程中姿態(tài)不穩(wěn)引起的高程跳點(diǎn)較多的情況，因此Savitzky-Golay平滑適合于剔除無人機(jī)數(shù)據(jù)中的高程突變噪聲。

?

該方法對于提取的無人機(jī)數(shù)據(jù)有較好的平滑效果，盡可能保留了沉陷盆地底部有效的沉陷信息。Savitzky-Golay平滑濾波的效果隨著選取窗口寬度的不同而不同，可以滿足多種不同場合的需求[25]。

無人機(jī)沉陷DEM兩個工作面走向下沉主斷面的下沉曲線如圖7所示。剖面路徑起點(diǎn)位置和目標(biāo)位置確定，方向均由北到南，圖中最大下沉值分別向靠近原點(diǎn)和遠(yuǎn)離原點(diǎn)方向偏移，這與兩個工作面的開采方向一致。圖7（a）中2020-05—2019-09曲線最大下沉值為2 912 mm，實(shí)測最大值為2 892 mm，誤差為20 mm，占實(shí)測走向最大下沉值的0.69%。總體來看，采用無人機(jī)監(jiān)測礦區(qū)沉陷能夠較為全面地反映出礦區(qū)沉陷影響范圍，最大下沉值的監(jiān)測精度較高。

4 結(jié)論

（1）本研究提出了一種基于免像控?zé)o人機(jī)攝影測量的礦區(qū)沉陷監(jiān)測方法，闡述了無人機(jī)攝影測量礦區(qū)監(jiān)測作業(yè)的實(shí)際工作流程，構(gòu)建了基于差分DEM的動態(tài)沉陷盆地模型，并驗(yàn)證了無人機(jī)沉陷監(jiān)測的精度，討論了隨著工作面推進(jìn)沉陷盆地發(fā)展情況和盆地主斷面下沉規(guī)律。

（2）免像控?zé)o人機(jī)攝影測量技術(shù)為在礦區(qū)受到地形限制或人員難以進(jìn)入的區(qū)域提供了一種安全、可靠、靈活的地表沉陷監(jiān)測方法；無人機(jī)攝影測量產(chǎn)品精度受到數(shù)據(jù)集質(zhì)量、研究區(qū)域特征、飛行參數(shù)設(shè)置、濾波插值過程引入的模型誤差以及使用的處理工具等因素影響。

（3）未來需對無人機(jī)多時相數(shù)據(jù)的可重復(fù)性進(jìn)行更全面地評估。根據(jù)不同的表面類型進(jìn)行多重誤差評估以產(chǎn)生更真實(shí)的地形變化的空間表示。