基于無人機激光雷達技術的開采沉陷監測方法與參數反演

亓立壯， 安士凱， 田超， 周大偉*， 董祥

(1.中國礦業大學環境與測繪學院， 徐州 221116； 2.平安煤炭開采工程技術研究院有限責任公司， 淮南 232033； 3.黑龍江地理信息工程院， 哈爾濱 150081； 4.浙江中測新圖地理信息技術有限公司， 湖州 313200)

煤炭資源的開采破壞了礦區的原地質結構，從而易引發地表沉陷、滑坡、泥石流、裂縫等地質災害。開采沉陷監測就是要對開采引起的地表形變狀況進行評估，進而找出沉陷變形規律，從而預測和預防沉陷災害。目前開采沉陷應用最廣泛的方法是概率積分模型，該模型需要事先以實測數據為依據獲得預計參數。因此必須在煤礦開采區域建立觀測站，收集開采引起的地表形變數據，用來反演概率積分預計參數。

目前開采沉陷監測手段主要有傳統觀測站測量、合成孔徑雷達干涉測量(interferometric synthetic aperture radar，InSAR)、三維激光掃描以及無人機攝影測量等。傳統開采沉陷監測手段包括水準儀導線測量、全站儀測量、全球導航衛星系統(global navigation satellite system，GNSS)測量等，雖然憑借其測量精度高、反映沉陷規律準確的優勢暫時無法被完全取代，但是存在工作量大、效率低，監測范圍局限等諸多問題[1]。INSAR是以雷達天線記錄的回波信號為信息源，獲取地球表面的三維地形、地表形變和地物特征變化等信息的測量技術，可以對礦區進行大面積、長時間的地表沉陷的監測，具有的全天候、大范圍覆蓋的監測優勢，但由于其原理的局限性，存在大變形區域失相干的缺點，此外高分辨率高精度數據目前比較昂貴[2]。近年來，無人機(unmanned aerial vehical,UAV)攝影測量技術發展迅速，該技術具有易操作、飛行受地形因素影響小和監測周期短等優點，利用無人機影像可以生成礦區地表數字高程模型(digital elevation model，DEM)，兩期DEM相減即可得到測區范圍內的下沉盆地，但由于受相機分辨率的影響，以及畸變、標定誤差和數據處理誤差等因素的存在，導致其高程精度還處在分米級水平，而且相片均是地物的表面信息，當植被過于密集時就不能獲取足夠多的地面信息[3]。

三維激光掃描作為變形監測領域的新手段，被廣泛應用于地形形變監測、沉陷滑坡變形監測、橋梁和土石壩變形監測以及基坑變形監測[4-6]等，數據量大且精度較高，但是傳統地面三維激光掃描單站采集范圍有限且存在地物遮擋，多站搬運步驟煩瑣、效率低下[7]。隨著無人機技術的迅速發展，特別是大疆等公司的不斷推陳出新，無人機設備性能有了極大的提升，使得采用無人機作為飛行平臺搭載三維激光掃描儀成為了現實。該方法進行開采沉陷監測具有快速靈活的特點，不需要固定測點，節約成本。由于激光可以穿透植被反射地表信息，并且獲得點云密度大，因而該方法可以獲得豐富的沉陷盆地觀測數據。已得到實際應用，湯伏全等[8]利用低空無人機LiDAR對神榆礦區某工作面進行了開采沉陷監測，對無人機三維激光點云數據建模方法進行了研究；楊凡等[9]在利用無人機激光雷達進行植被高度監測時驗證了該技術激光點云的精度可達厘米級；Ren等[10]對無人機平臺、不同傳感器及其應用領域等方面進行了探討，構建了一個無人機激光雷達變形監測應用框架。

現以實際工程為例，采用無人機激光雷達技術進行礦區開采沉陷監測及預計參數反演，探討無人機激光點云數據的處理方法。

1 研究區域概況及數據獲取

1.1 研究區域概況

試驗區位于內蒙古唐家會礦區，該礦區處在內蒙古自治區準格爾煤田東孔兌普查區的西南部，行政區劃隸屬準格爾旗大路鎮管轄，具體位于準格爾旗人民政府所在地薛家灣鎮西北約4 km處，其地理坐標為：111°10′27″E～111°14′34″E，39°52′45″N～ 39°57′22″N。測區地處鄂爾多斯高原，地形多種多樣，屬典型的侵蝕性丘陵地貌，區域內溝壑縱橫，呈樹枝狀分布，地形起伏較大，植被較為稀疏，主要以低矮植物為主。唐家會煤礦該礦在衛星影像上的位置如圖1所示。選擇61202工作面上方地表作為觀測地點，該工作面尺寸為240 m×1 960 m，平均傾角2°，平均采深538 m，平均采厚19 m，采用長壁后退式全垮綜放采煤法開采，設計采放比1∶3。

圖1 唐家會礦測區位置Fig.1 Location of survey area in Tangjiahui mine

1.2 UAV激光雷達地表沉陷監測的方法

1.2.1 UAV激光雷達地表沉陷監測原理

UAV激光雷達是以多旋翼無人機平臺為載體，一體化集成高精度激光掃描儀、全球定位系統(global positioning system,GPS)、慣性測量單元(inertial measurement unit,IMU)等傳感器，通過配備的數據處理和應用軟件，快速生成數字地表模型(digital surface model,DSM)、DEM，具有重量輕、攜帶方便、成果處理效率高的優點。

在工作面開采上方地表劃定適當范圍的測區，在工作面推進的兩個時刻分別對該測區進行觀測，并利用采集到的激光點云數據建立兩期數字高程模型DEM，這兩期DEM相減即為這兩個觀測時間點之間由于工作面開采所產生的地表下沉盆地。基于概率積分法動態求參原理，利用全盆地的地表下沉數據，再結合工作面開采相關信息，即可獲得該地區開采沉陷預計參數。

1.2.2 UAV激光雷達地表沉陷監測技術流程

UAV激光雷達地表沉陷監測流程主要為外業采集、求取下沉盆地、預計參數動態反演。其中外業采集的關鍵步驟為測區選擇、航線規劃、實地觀測；求取下沉盆地需要先進行數據預處理，經過POS定位定姿系統(position and orientation system,POS)組合解算、數據融合解算之后生成具有較高精度三維坐標的地面激光點云，再經過分類濾波得到地面點云，最后構建兩期DEM，DEM相減得到下沉盆地；預計參數動態反演主要是在全盆地均勻選點，根據概率積分原理進行預計參數求取。主要技術流程如圖2所示。

圖2 無人機激光雷達監測地表沉陷技術流程Fig.2 Technical process of UAV lidar monitoring surface subsidence

1.2.3 數據采集

采用ARS-450i激光測量系統，并利用大疆M600Pro無人機作為飛行平臺。該設備集成了激光掃描儀、GPS、慣導等傳感器，可同步獲取三維激光點云和定位定姿數據，儀器參數如表1、表2所示。

表1 智喙ARS-450i激光測量系統主要參數Table 1 Main parameters of ARS-450i laser measurement system

表2 大疆M600Pro無人機飛行平臺主要參數Table 2 Main parameters of DJI M600Pro UAV flight platform

首先根據測區現場實際情況以及無人機續航能力等因素，設計飛行方案。本實驗測區面積2.36 km2，分6架次全覆蓋掃描。該設備推薦相鄰航線重疊率為30%，按照儀器參數設計飛行高度為70 m。分別于2020年8月和2021年3月對上述測區進行了兩期數據采集。圖3為操作無人機進行外業采集的現場照片，圖4為外業采集航線設計示意圖。

圖3 無人機激光雷達數據采集Fig.3 UAV lidar data acquisition

圖4 唐家會測區無人激光雷達航線示意圖Fig.4 Route map of unmanned lidar in Tangjiahui survey area

現場實測前需要在測區附近的已知點位(或者現場設置木樁，并測量樁位三維坐標)架設GPS基站，該步驟的目的是采集GPS靜態數據與無人機激光雷達采集的POS數據進行組合解算，得到有較高精度的定位定姿數據。確認基站接收數據后，先開啟激光雷達靜態對齊5 min，目的是對慣性導航系統進行初始位置對準。然后按照飛行方案進行飛行和數據采集，數據采集完畢飛機安全降落后，再次靜態對齊5 min，確認數據無誤后方可關閉基站。

為了驗證無人機激光雷達數據精度，兩期實驗分別還在測區內均勻采集了若干檢查點的三維坐標。檢查點主要選在道路等交通便利、無植被的裸露地面，考慮到測區內地形復雜交通不便以及工作效率，實驗利用唐家會礦區的連續運行衛星定位服務參考站(continuously operating reference stations,CORS)系統對檢查點進行了實時動態(real-time kinematic,RTK)測量。

2 數據處理及沉陷盆地構建

2.1 激光點云數據處理

2.1.1 數據預處理

無人機激光雷達數據采集整理過后，并不能直接用于DEM構建，而是需要經過POS組合解算、數據融合之后，才能獲得包含較高精度三維坐標信息的激光點云數據。具體步驟如下。

步驟1激光雷達外業采集的位姿數據包括機載GPS天線數據、慣性導航數據以及里程計數據，POS組合結算是指將上述三者與地面基站收集的同時段GPS靜態數據進行聯合結算，獲得具有較高精度的定位定姿數據。解算方法本文采用了緊組合方式進行組合解算[11]，該方法是運用卡爾曼濾波原理解算信號殘差，對GPS與IMU坐標速度之差進行校正，因其耦合程度較高目前得到了最廣泛的應用。上述步驟主要在Inertial Explorer軟件里實現。

步驟2數據融合解算是指利用步驟1得到的高精度的POS數據，對激光雷達采集的原始三維激光點云數據進行融合處理，得到實驗所需坐標系下的具有較高精度三維坐標信息的激光點云數據。因為每個架次的飛行參數一致，所以采集到的數據都是在同一個坐標系下，可以自動拼接為一整個測區的完整點云數據。本實驗利用檢查點的RTK數據進行四參數求取和坐標轉換，最終得到54坐標系下的測區激光點云。

2.1.2 地面點提取

監測區域位于荒漠化草原和草原之間的過渡地帶，地面存在人工種植的苔草、沙棘、灌木等植被，此外還存在零星建筑物、車輛、電線桿等，這些會對DEM的構建造成影響，圖5為兩期原始點云數據。要想獲得高精度的地面DEM，就必須對初始點云數據進行濾波處理僅留取地面點云。目前常用的濾波算法有基于布料模擬濾波算法、形態學濾波算法、基于不規則三角網濾波算法等[12-13]。Zhang等[14]提出的布料模擬激光點云濾波算法(cloth simulation filtering，CSF)是將原始點云數據看作一個整體，將其進行倒置，再通過模擬虛擬布料下降的物理過程實現地面點云的過濾。如圖6所示，筆者在測區內選取尺寸為100 m的正方形區域作為實驗數據，用手動濾波的方式制作標準地面點云，再用多種濾波算法分別進行濾波與之進行對比，具體比較結果如表3所示，在比較多種濾波算法結果后發現布料模擬濾波算法在測區內最為適用，故本實驗激光點云數據處理均采用布料模擬算法進行濾波。最終得到兩期地面激光點云(圖7)。

圖5 唐家會激光雷達兩期原始激光雷達數據Fig.5 Original lidar data of two phases of Tangjiahui lidar

2.2 DEM構建及精度分析

利用兩期地面激光點云數據，以1 m的間隔進行克里金插值，分別構建兩個時間點的測區地表DEM。再利用測區內的檢查點數據對兩期DEM分別進行精度評定。兩期DEM如圖8所示，DEM精度評定結果如圖9所示。通過計算得出第一期DEM高程中誤差為0.034 m，第二期DEM高程中誤差為0.037 m。

2.3 下沉盆地建立及分析

用第二期DEM減去第一期DEM，得到2020年8月15日—2021年3月13日唐家會測區61202工作面開采引起的地表下沉盆地，如圖10所示。測區內最大下沉值約為6.9 m，最大下沉點位于工作面推進約475 m處，下沉盆地東北方向矸石山因堆積升高最大26 m，測區西北部挖土坑深13 m。

(1)從儀器自身角度來看：掃描儀測角、測距本身就存在一定的誤差，GPS定位存在軌道誤差、鐘差、大氣延遲、多路徑效應等誤差，IMU的陀螺儀和加速計也存在自身產生的偏差；機載GPS與IMU之間的相對位置關系也會影響測量系統的定姿定位，雖然廠商會提供相關的參數，但在實際使用中難免存在磕碰、松動等不利因素，導致傳感器安置誤差的產生。

(2)從數據采集的角度來看：內蒙古地區多風沙，雖然是在儀器允許條件下起飛，但是數據的采集質量依然受到天氣條件的影響，這主要是受飛行平臺的性能(如穩定性、抗干擾能力等)限制；航高也是影響數據精度的一個重要因素，且在一定范圍內激光點云誤差和航高呈線性關系，點云密度也會隨著航高增加而降低。

圖6 不同濾波方式誤差分布Fig.6 Error distribution of different filtering methods

表3 不同濾波方式比較Table 3 Comparison of different filtering methods

圖7 地面點云Fig.7 Ground point cloud

圖8 唐家會激光雷達DEMFig.8 Tangjiahui lidar DEM

圖9 唐家會激光雷達DEM精度統計Fig.9 DEM accuracy statistics of Tangjiahui lidar

圖10 無人機激光雷達監測結果Fig.10 UAV lidar monitoring results

(3)從數據處理的角度來看：無人機激光雷達數據處理過程比較煩瑣，尤其是濾波很難做到完美提取地面點，總是會存在噪聲誤差；利用地面點構建DEM進行插值時也難免會產生誤差。

3 概率積分參數反演及分析

在煤礦開采監測領域中，為了研究復雜的沉降過程以及地表破壞程度，獲取沉降值和沉降參數至關重要。預測開采沉陷也是核心研究內容之一[15-17]，對煤田生產和地面防護具有重要意義，可用于確定地層結構是否受采動影響及其相應程度。此外，預測開采沉陷是維修、加固、重建或采取地下措施的基礎。預計參數是決定沉陷預測是否準確的關鍵因素，因此，求取準確的預計參數有重要意義。傳統的常規求參方法需要設置觀測線，但存在點位較少，需要獲取最終下沉，且容易丟失破壞測點、野外勞動量大、監測周期長等諸多問題。無人機激光雷達無需獲取完整采動過程的下沉值，只需對監測區域某一時間段或多個時間段內的兩期或多期數據進行采集，得到高密度全盆地動態沉陷數據，然后通過建立動態求參模型進行參數求取，可以獲取更加精確的沉陷參數。

3.1 概率積分法動態預計原理

在開采沉陷預計方法中，基于隨機介質理論的概率積分法在中國的應用最為廣泛[15]，該方法是同時考慮了地表沉降的空間位置和時間的相關性，即用概率積分法穩態預測模型乘以時間效應函數，以獲取在煤礦開采下不同時間的地表動態下沉值；而概率積分法穩態預測模型是指工作面開采結束后所引起地表沉降的最終狀態，可以用來說明影響范圍內任意點的預測下沉，其單位開采所造成地表沉陷盆地的計算公式為

(1)

式(1)中：We(x)為x位置的地表下沉值；r為主要影響半徑。

由于工作面開采和地表變形是在三維空間內的，則開采一個長為t寬為S的范圍所造成的任一點A(x,y)的下沉可表示為

(2)

當工作繼續推進，直至開采完畢時，由全部開采導致的A點下沉值的預測公式為

(3)

式(3)中：l為走向長度；L為傾向寬度；W0為達到充分開采的最大下沉值，W0=mqcosα，其中，q為下沉系數；m為采煤厚度；α為煤層傾角；主要影響半徑r=H/tanβ，其中，tanβ為主要影響正切，H為開采深度。

因此，若要根據式(1)～式(3)推導出地表其他位置的下沉值，需要確定預測參數。

開采所造成的地表沉陷是一個連續的時空過程，與空間位置、時間、開采速度、覆巖性質等密切相關；其移動和變形可以看作連續的時空函數。動態預測模型也是建立沉陷與空間位置和時間的相關性。用Wt表示t時刻地表點的下沉值，則Wt=W0f(t)，其中f(t)為反映開采沉陷動態變化的時間效應函數。

Knothe時間函數[16]反映了任意時刻的下沉至與最終穩態下沉值之差和該時刻的下沉速度的比例關系，其最終形式為

f(t)=1-e-ct

(4)

式(4)中：c為時間常數。

3.2 基于動態預計的反演方法

運用概率積分法進行開采沉陷預計需要事先獲取預計參數，這些參數可以根據實測數據運用數學方法進行反演估計[17]。目前參數反演方法主要有模矢法、線性最小二乘法、遺傳算法、粒子群算法等。其中模矢法應用最為廣泛，該方法的基本步驟如下。

步驟1構建目標函數(預計值與實際值之差的平方和)，給定一組參數初始值作為基點，代入計算函數值。

步驟2按照一定步長在參數初始值兩側取值，分別代入目標函數得到兩組新函數值。

步驟3比較3個函數值大小，選擇其中最小的那組參數作為矢點，然后沿坐標系搜索最佳點作為新的基點。

步驟4繼續迭代直到求得滿足要求的解，該步驟可以在開采沉陷預測預報系統中實現，該軟件是由吳侃開發的一款基于概率積分原理的預計模擬軟件，可以進行煤礦開采沉陷預計與參數反演。

3.3 求參結果及分析

以2.3節所述整個下沉盆地作為原始數據，運用本章節介紹的反演方法，進行唐家會測區工作面開采沉陷預計參數求取，主要是求取下沉系數q和主要影響正切tanβ。此外，為了對反演結果進行驗證，收集唐家會礦區61101工作面的傳統觀測線數據，該數據是在工作面正上方沿走向和傾向布設兩條觀測線，采用GNSS技術，觀測了該工作面開采前至穩沉的地表移動數據。以相同方法進行參數求取并進行對比，參數求取結果如表4所示。

兩組結果對比，運用無人機激光雷達觀測的全盆地數據求參結果和傳統觀測線求取結果差距較小，其中下沉系數相差13%，主要影響正切基本一致。對于下沉系數，傳統觀測站由于其線狀的布站形式，極有可能捕捉不到最大下沉點，而無人機激光雷達則可以獲取整個盆地的下沉數據，且2.3小結分析無人機激光雷達求下沉盆地精度約為0.05 m，僅占測區開采沉陷最大值6.9 m的0.7%，因此可以認為該方法反演下沉系數是可靠的；對于主要影響正切，雖然與傳統觀測線數據求參結果基本一致，但由于中國一般以下沉10 mm為沉陷盆地邊界，無人機激光雷達受其精度限制無法捕捉盆地邊界，故認為運用該數據求取主要影響正切理論上可能存在一定誤差。文獻[3]對全盆地求取預計參數的可靠性做了詳細分析，當下沉數據中誤差占最大下沉值的比例小于7%時，計算的下沉系數、主要影響正切是可靠的。

表4 預計參數Table 4 Estimated parameters

4 討論

無人機激光雷達應用于開采沉陷監測其最大的優勢就在于可以快速高效的礦區地表三維信息的獲取和建模。具體表現在：①搭載的激光雷達可以獲取較高精度的、高密度的礦區地物激光點云，由于激光具有穿透性，經過濾波處理即可獲得大量的地面信息；②以無人機作為飛行平臺可以節省大量的人力物力，相當于把測站設在了空中，克服了地面三維激光掃描儀視線遮擋的問題；無人機與三維激光掃描儀的結合可以發揮各自優勢，掃描面積更大、獲得的信息更全面，且數據處理的流程已經較為成熟，可以形成一整套無人機激光雷達數據采集、處理技術流程。中國煤炭開采正不斷向大規模高強度發展，開采活動對地表的影響更加明顯，與礦區居民和企業的安全發展密切相關，傳統的監測手段存在獲取信息量少、工作量大的問題，無法適應煤炭開采的發展趨勢，在這種形勢下，無人機激光雷達作就成為一種較為合理的技術選擇。

分析無人機激光雷達的DEM建模平均高程中誤差為3.5 cm，結合現有文獻[4-9]可以確定運用該方法進行地表建模高程精度在5 cm之內。未作固定觀測站來檢核下沉盆地精度，利用誤差傳播定律得到下沉盆地精度為5.0 cm，事實上兩期DEM相減可以消除部分系統誤差，下沉盆地精度理論上應該更高。中國一般以下沉10 mm等值線作為開采沉陷盆地邊界，因此無人機激光雷達的精度達不到監測沉陷盆地邊界及其他小變形的要求，但是與該目標已經十分接近，并且由于其全盆地監測的優勢，沉陷數據全面，在求取概率積分預計參數方面較為可靠。綜合分析，無人機激光雷達可以在以下方面提高檢測精度：一是進一步提高激光雷達的精密度，儀器本身為多傳感器集成，激光掃描儀、慣導、GPS接收機均存在一定誤差，三者的位置關系是否固定也影響系統定姿定位；二是進一步提高無人機性能，無人機的抗風能力和穩定性影響掃描姿態進而影響激光點云數據的精度；三是進一步改進數據處理方法，改進濾波算法以及DEM建模方法以獲取更為精確的地面信息。

無人機除搭載激光雷達之外，還能搭載相機進行攝影測量。無人機攝影測量與無人機激光雷達相比各有優勢與不足[18-21]，分析如下。

(1)精度方面，無人機攝影測量由于其相機性能、相片幾何畸變、重疊度、飛行平臺穩定性等因素的制約，目前還處于分米級的階段，與無人機激光雷達還存在一定差距，但二者都能反映全盆地的沉陷信息，在進行預計參數反演尤其是下沉系數的求取時都具有相當高的可靠性。

(2)數據方面，無人機攝影測量獲取的實時圖像數據，首先由于相片本身是地物表面信息，因此當監測區域植被過于密集時，難以捕捉足夠的地物信息，濾波過后存在較多空白區域；其次當地形起伏較大時，會造成連續的影像之間重疊度損失，實際工程常常被迫設置高重疊度從而增加工作任務，再者當大氣質量較差或光線較暗時，相機難以獲取高清晰度的影像信息，而激光雷達是主動發射激光接受反射信息，以其強大的穿透性可以穿過大氣、植被直接反射地表信息，亦不受環境亮度的影響。

(3)工作效率方面，受點云密度的制約，搭載激光雷達的無人機往往比搭載相機的無人機飛行速度慢；激光雷達設備重量大，對飛機性能尤其是電池續航能力要求更高，實際應用中往往需要頻繁起降更換電池；激光點云反射距離有限，使得激光雷達無人機飛行高度一般要低于攝影測量無人機；攝影測量需要在測區內布設若干像控點，并且隨測區增大所需點數亦增大，激光掃描則不需要像控點的布設，因而節省了一部分人力物力。

綜上所述，當需要快速大面積監測作業，且對精度要求較低時，選擇無人機攝影測量更為合適，而測量范圍相對較小，且要求厘米級精度時，無人機機光雷達則更具優勢。隨著無人機技術的繼續發展和激光測量系統的進一步改進，無人機激光雷達在數據采集效率和精度上還有極大的進步空間，長遠來看，該技術在開采沉陷監測領域更具發展優勢和潛力。

5 結論

(1)目前主流的開采沉陷監測技術均存在著精度與效率兩者不可兼得的問題，無人機激光雷達充分集合了三維激光掃描數據量大精度高與無人機操作方便工作迅速的優點，既能快速獲取大范圍的地表激光點云數據，又能有效節約時間與勞動成本，尤其適合地形復雜多樣的西部礦區進行開采沉陷監測。

(2)對原始激光數據分別進行組合解算、數據融合、分類濾波之后，建立兩期測區數字高程模型DEM，通過與檢查點對比得到中誤差分別為0.034 m和0.037 m；兩期DEM相減得到觀測時段內的開采沉陷盆地，最大下沉值約為6.9 m，根據誤差傳播定律得到下沉盆地中誤差約為0.050 m。

(3)對下沉盆地進行均勻取點作為原始觀測數據，采取基于動態預計的反演方法，求取該礦區概率積分預計參數，其中下沉系數主要影響正切比較可靠。收集同一礦區內臨近工作面的傳統觀測線數據進行參數反演，得到的預計參數結果與無人機激光雷達全盆地求參基本一致。綜上所述，利用無人機雷達進行開采沉降監測是可行的，利用無人機激光雷達監測所得完整沉陷盆地進行預計參數反演是可靠的。