無(wú)人機(jī)LiDAR在采空區(qū)沉陷監(jiān)測(cè)中的應(yīng)用

鄭俊良 姚頑強(qiáng) 藺小虎 張聯(lián)隊(duì) 張?jiān)? 張冬

摘 要：煤炭資源高強(qiáng)度開(kāi)采對(duì)地表造成了極大破壞，嚴(yán)重威脅建（構(gòu)）筑物安全，現(xiàn)有的開(kāi)采沉陷監(jiān)測(cè)方法缺乏面域形變監(jiān)測(cè)，限制了采空區(qū)地表形變規(guī)律研究。針對(duì)傳統(tǒng)觀測(cè)方法的不足，以西部煤礦某工作面為研究對(duì)象，采用無(wú)人機(jī)LiDAR（Light Detection and Ranging）測(cè)量方法，進(jìn)行采空區(qū)面域沉陷監(jiān)測(cè)研究。針對(duì)采空區(qū)地形特征，改進(jìn)漸進(jìn)三角網(wǎng)點(diǎn)云濾波過(guò)程中種子點(diǎn)的選取方法，提高點(diǎn)云濾波精度；對(duì)研究區(qū)沉陷模型精度進(jìn)行驗(yàn)證，結(jié)果表明：沉陷模型整體誤差小于100 mm，平均誤差小于50 mm；最后，基于沉陷模型進(jìn)行數(shù)據(jù)挖掘分析，獲取下沉等值線圖、下沉面積統(tǒng)計(jì)圖、下沉分區(qū)統(tǒng)計(jì)圖等，實(shí)現(xiàn)采空區(qū)全方位沉陷監(jiān)測(cè)。為煤礦煤炭資源開(kāi)采、地表下沉監(jiān)測(cè)、土地復(fù)墾和生態(tài)環(huán)境保護(hù)等提供參考。

關(guān)鍵詞：無(wú)人機(jī)LiDAR；點(diǎn)云濾波；沉陷監(jiān)測(cè)；沉陷模型

中圖分類號(hào)：TD 171

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1672-9315（2023）04-0825-11

DOI：10.13800/j.cnki.xakjdxxb.2023.0420

Application of UAV LiDAR in monitoring subsidence of goaf

ZHENG Junliang1，YAO Wanqaing1，LIN Xiaohu1，

ZHANG Liandui2，

ZHANG Yong1，ZHANG Dong3

（1.College of Geomatics，Xian University of Science and Technology，Xian 710054，China；

2.Shaanxi Binchang Mining Group Co.，Ltd.，Xianyang 712000，China;

3.Shaanxi Yanchang Petroleum Jingbian Coal Industry Co.，Ltd.，Yulin 718500，China）Abstract：High-intensity mining of coal resources has caused great damage to the surface and seriously threatened the safety of buildings（structures）.The existing mining subsidence monitoring methods lacked overall surface deformation monitoring，which limited the study of surface deformation laws in goaf.In response to the shortcomings of traditional observation methods，a working face in the western coal mine was taken as the research object，and the unmanned aerial vehicle LiDAR（Light Detection and Ranging）survey method was used to monitor the subsidence of the goaf area.In light of the terrain characteristics of the goaf，the selection method of seed points in the process of progressive triangulation cloud filtering was improved so as to promote the precision of point cloud filtering；the accuracy of the subsidence model in the study area was verified.The results showed that the overall error of the settlement model was less than 100 mm，and the average error was less than 50 mm；through the data mining analysis of the subsidence model，the subsidence contour map，subsidence area statistical map，and subsidence zoning statistical map were obtained to achieve the all-round subsidence monitoring of the goaf.The research provides a? reference for coal mining，surface subsidence monitoring，land reclamation and ecological environment protection.

Key words：UAV LiDAR；point cloud filtering；subsidence monitoring；subsidence model

0 引 言

隨著能源生產(chǎn)和消費(fèi)革命的推進(jìn)，煤炭在能源結(jié)構(gòu)中的占比逐漸減少，但在一段時(shí)間內(nèi)煤炭仍屬于中國(guó)主體能源，具有重要戰(zhàn)略地位 ［1-2］。西部煤礦地處黃土高原中上游，地形復(fù)雜多樣，屬于典型的西部黃土溝壑區(qū)，生態(tài)環(huán)境十分脆弱。機(jī)械化綜采、垮落式頂板管理的采煤方法造成上覆巖層破裂、地表塌陷，引起地下水流失、土地沙漠化，破壞沉陷區(qū)生態(tài)環(huán)境、威脅地表建（構(gòu)）筑物安全［3-5］。因此，針對(duì)西部煤礦所在的黃土溝壑區(qū)，開(kāi)展可靠的礦區(qū)開(kāi)采沉陷監(jiān)測(cè)，準(zhǔn)確掌握開(kāi)采沉陷動(dòng)態(tài)變化過(guò)程，對(duì)開(kāi)采沉陷預(yù)計(jì)、地表建（構(gòu)）筑物損壞評(píng)價(jià)、生態(tài)修復(fù)具有重要意義。

目前，開(kāi)采沉陷監(jiān)測(cè)方法包括傳統(tǒng)預(yù)埋監(jiān)測(cè)樁觀測(cè)（水準(zhǔn)測(cè)量、全站儀測(cè)量、GNSS測(cè)量等）、InSAR監(jiān)測(cè)、地面三維激光測(cè)量、傾斜攝影測(cè)量和機(jī)載LiDAR測(cè)量。傳統(tǒng)預(yù)埋監(jiān)測(cè)樁觀測(cè)方法局限于點(diǎn)位觀測(cè)，具有較高的點(diǎn)位精度（毫米級(jí)、厘米級(jí)），但不能準(zhǔn)確描述采空區(qū)地面整體的形變過(guò)程，且西部溝壑區(qū)地形結(jié)構(gòu)復(fù)雜，存在工作量大和人員無(wú)法到達(dá)的問(wèn)題 ［6］；InSAR測(cè)量能夠進(jìn)行面狀觀測(cè)，且取得較高的測(cè)量精度（毫米級(jí)），但受限于觀測(cè)周期，對(duì)于地表塌陷速度快、沉陷量大的區(qū)域，容易造成失相干現(xiàn)象［7-8］；地面三維激光測(cè)量能夠取得面狀厘米級(jí)，甚至毫米級(jí)觀測(cè)精度，但由于其視野范圍、點(diǎn)云拼接誤差等原因，在超長(zhǎng)超寬工作面監(jiān)測(cè)時(shí)存在工作量大、累積誤差大等問(wèn)題［9］；近年來(lái)，隨著影像密集特征點(diǎn)匹配算法的成熟，無(wú)人機(jī)傾斜攝影測(cè)量技術(shù)提高了攝影測(cè)量的觀測(cè)精度，但在地表植被覆蓋度高、地形起伏大的區(qū)域，易造成監(jiān)測(cè)結(jié)果精度低，且精度受限于像控點(diǎn)密度［10-11］；機(jī)載LiDAR測(cè)量具有精度高、觀測(cè)范圍大的優(yōu)勢(shì)，但有人機(jī)載LiDAR測(cè)量存在設(shè)備昂貴、天氣要求高、空域申請(qǐng)困難等問(wèn)題；隨著無(wú)人機(jī)技術(shù)發(fā)展和激光掃描設(shè)備改進(jìn)，無(wú)人機(jī)LiDAR測(cè)量技術(shù)逐漸成熟，在滿足厘米級(jí)測(cè)量精度的同時(shí)能夠極大地降低成本、提高效率，彌補(bǔ)有人機(jī)載LiDAR測(cè)量的不足［12-13］。無(wú)人機(jī)LiDAR通過(guò)對(duì)目標(biāo)區(qū)域低空多時(shí)相掃描，可獲取高密度點(diǎn)云數(shù)據(jù)，經(jīng)點(diǎn)云濾波、插值和求差計(jì)算，可獲取目標(biāo)區(qū)域高程變化。

但在激光點(diǎn)云濾波、插值和求差計(jì)算過(guò)程中不可避免地存在分類錯(cuò)誤、噪聲和誤差，一些學(xué)者通過(guò)改進(jìn)濾波算法和插值方法提高DEM精度，進(jìn)而改善沉陷模型精度［14-17］，但仍存在一定問(wèn)題：①現(xiàn)有成熟點(diǎn)云濾波算法主要針對(duì)傳統(tǒng)有人機(jī)載LiDAR點(diǎn)云數(shù)據(jù)，相比較而言，無(wú)人機(jī)LiDAR點(diǎn)云具有更高的點(diǎn)云密度和更復(fù)雜的噪聲數(shù)據(jù)，如在山區(qū)應(yīng)用效果較好的漸進(jìn)三角網(wǎng)點(diǎn)云濾波算法（Progressire Trianglated Irregular Network，PTIN），在無(wú)人機(jī)LiDAR點(diǎn)云濾波時(shí)，由于高密度和多噪聲的原因，易出現(xiàn)種子點(diǎn)錯(cuò)誤問(wèn)題，導(dǎo)致濾波錯(cuò)誤；②無(wú)人機(jī)LiDAR獲取的點(diǎn)云、數(shù)字高程模型（Digital Elevation Model，DEM）和沉陷模型之間的精度缺乏實(shí)測(cè)數(shù)據(jù)驗(yàn)證；在無(wú)人機(jī)LiDAR測(cè)量系統(tǒng)和環(huán)境條件相同情況下，沉陷模型精度是否優(yōu)于原始點(diǎn)云和DEM精度缺乏試驗(yàn)對(duì)比；③傳統(tǒng)LiDAR監(jiān)測(cè)成果單一，對(duì)沉陷模型進(jìn)一步研究較少，未進(jìn)行深入的數(shù)據(jù)分析和挖掘。

1 研究區(qū)及數(shù)據(jù)獲取

1.1 研究區(qū)概況

涼水井礦位于陜西省榆林市東北部、榆神府礦區(qū)東部，處于北緯38°47′～38°53′和東經(jīng)110°14′～110°21′之間（圖1）。井田地面標(biāo)高1 100～1 326 m，地勢(shì)西北高東南低。以涼水井煤礦431301工作面為研究對(duì)象，開(kāi)采煤層為4-3煤層，該煤層平均埋深138 m，平均厚度1.3 m，采用長(zhǎng)臂式綜采、頂板自然塌陷方法開(kāi)采；其中上覆4-2煤層已于2013年開(kāi)采完成，地表于2014年達(dá)到穩(wěn)定狀態(tài)。

基于《煤礦測(cè)量規(guī)程》（2013）設(shè)置地面監(jiān)測(cè)點(diǎn)，結(jié)合臨近礦區(qū)開(kāi)采沉陷角量參數(shù)和現(xiàn)場(chǎng)實(shí)際情況，共布設(shè)3條監(jiān)測(cè)線：走向觀測(cè)線、傾向觀測(cè)線和道路觀測(cè)線，如圖2所示。地面監(jiān)測(cè)點(diǎn)通過(guò)埋設(shè)預(yù)制觀測(cè)樁，采用四等水準(zhǔn)測(cè)量，獲取監(jiān)測(cè)點(diǎn)多時(shí)相高程，實(shí)現(xiàn)地表沉陷監(jiān)測(cè)。在工作面停采線附近，選擇面積為0.4 km2的區(qū)域作為研究區(qū)，進(jìn)行無(wú)人機(jī)LiDAR試驗(yàn)數(shù)據(jù)采集和分析。431301工作面回采期間共進(jìn)行23次全面觀測(cè)，根據(jù)開(kāi)采沉陷一般規(guī)律，選取位于研究區(qū)域內(nèi)的39個(gè)地面監(jiān)測(cè)點(diǎn)作為無(wú)人機(jī)LiDAR參考數(shù)據(jù)，對(duì)無(wú)人機(jī)LiDAR監(jiān)測(cè)成果進(jìn)行精度分析，如圖2中紅色三角形所示。

1.2 數(shù)據(jù)獲取

無(wú)人機(jī)LiDAR型號(hào)為SZT-R250，集成了RIEGL miniVUX-1激光測(cè)距系統(tǒng)和高精度GNSS/IMU（Global Navigation Satellite System/Inertial Navigation System）系統(tǒng)，以及地面GNSS差分系統(tǒng)，掃描測(cè)程250 m。無(wú)人機(jī)有效飛行時(shí)間約20 min，試驗(yàn)飛行高度50 m，速度8 m/s，激光脈沖發(fā)射頻率為100 kHz，掃描速度為100線/s。為了保證多時(shí)相數(shù)據(jù)精度一致，在每次數(shù)據(jù)采集時(shí)均采用相同的飛行和航測(cè)參數(shù)。在2020年11月和2021年5月共進(jìn)行2次無(wú)人機(jī)LiDAR測(cè)量，數(shù)據(jù)采集范圍覆蓋該時(shí)間段地表下沉范圍，航測(cè)時(shí)間與工作面回采位置、試驗(yàn)區(qū)域如圖3所示。

圖3展示了井下巷道、工作面位置、采空區(qū)與地表航測(cè)范圍位置關(guān)系。2020年11月23日和2021年5月12日分別進(jìn)行無(wú)人機(jī)LiDAR航測(cè)，其對(duì)應(yīng)的工作面回采位置如圖3藍(lán)色豎線所示。同時(shí)航測(cè)范圍覆蓋研究期間采煤導(dǎo)致的地表下沉影響范圍，其中采空區(qū)1表示在首次數(shù)據(jù)采集時(shí)，已存在的采空區(qū)；采空區(qū)2表示第2次LiDAR數(shù)據(jù)采集時(shí)產(chǎn)生的采空區(qū)；圖中黑色三角形表示地面實(shí)測(cè)監(jiān)測(cè)點(diǎn)。

1.3 研究技術(shù)路線

無(wú)人機(jī)LiDAR開(kāi)采沉陷監(jiān)測(cè)主要步驟包括外業(yè)航測(cè)、LiDAR點(diǎn)云濾波、點(diǎn)云插值計(jì)算、DEM差值計(jì)算、沉陷模型獲取及數(shù)據(jù)挖掘分析等。同時(shí)，以地面監(jiān)測(cè)樁測(cè)量數(shù)據(jù)為依據(jù)進(jìn)行LiDAR點(diǎn)云、DEM和沉陷模型精度驗(yàn)證，技術(shù)路線如圖4所示。

采用漸進(jìn)三角網(wǎng)點(diǎn)云濾波算法進(jìn)行點(diǎn)云濾波，該算法一般由局部最低點(diǎn)作為稀疏種子點(diǎn)，創(chuàng)建初始不規(guī)則三角網(wǎng)，通過(guò)激光點(diǎn)與三角網(wǎng)的距離和角度進(jìn)行迭代分類。該算法適用性強(qiáng)、濾波效果較佳，但種子點(diǎn)的選擇是漸進(jìn)三角網(wǎng)濾波的關(guān)鍵。因此，結(jié)合礦區(qū)地表沉陷特征，對(duì)種子點(diǎn)的選取方法進(jìn)行改進(jìn)，確保種子點(diǎn)的準(zhǔn)確性。點(diǎn)云插值計(jì)算參考文獻(xiàn)［18］所述，選擇最適用于黃土溝壑區(qū)的Kriging插值算法進(jìn)行計(jì)算，獲取地表DEM。點(diǎn)云、DEM和沉陷模型精度采用平均誤差（Mean Error，ME）、平均絕對(duì)誤差（Mean Absolute Error，MAE）、中誤差（Root Mean Square Error，RMSE）驗(yàn)證，以水準(zhǔn)測(cè)量結(jié)果作為參考值，計(jì)算方法如下

式中 Rm為基于DEM提取的下沉值;Zm為地面監(jiān)測(cè)樁觀測(cè)下沉值；m為參考點(diǎn)數(shù)量。

2 LiDAR點(diǎn)云濾波與插值計(jì)算

2.1 常用濾波方法

無(wú)人機(jī)LiDAR能夠獲取海量地表點(diǎn)云數(shù)據(jù)，其中包括地面、植被、建筑物等，根據(jù)點(diǎn)云的幾何特征、色彩、高程變化等信息提取地面點(diǎn)的過(guò)程即為點(diǎn)云濾波。點(diǎn)云濾波方法分為：基于坡度濾波（Slope Based Filtering，SBF）［19］、基于漸進(jìn)三角網(wǎng)濾波［20］、形態(tài)學(xué)濾波（Morphological Filter，MF）［21］、布料模擬濾波（Cloth Simulation Filter，CSF）［22］等。SBF濾波算法采用點(diǎn)之間的坡度作為閾值參數(shù)進(jìn)行濾波，該方法對(duì)于緩坡地形效果較好，但對(duì)于地形復(fù)雜區(qū)域，由于單一閾值的限定，濾波效果不足；MF濾波算法借鑒形態(tài)學(xué)圖像處理方法進(jìn)行點(diǎn)云濾波，剔除非地面點(diǎn)時(shí)不能有效保留地形信息，易造成地形平滑；CSF濾波算法采用點(diǎn)云翻轉(zhuǎn)、模擬布料貼合倒置點(diǎn)云的方法，使用高程閾值分離地面點(diǎn)和非地面點(diǎn)，平地、緩坡和陡坡需采用不同參數(shù)，參數(shù)確定存在困難；PTIN濾波方法由Axelsson提出，根據(jù)種子點(diǎn)建立初始地形模型，對(duì)待定點(diǎn)進(jìn)行距離和角度閾值判斷，確定地面點(diǎn)或非地面點(diǎn)，逐步迭代完成點(diǎn)云濾波，該方法對(duì)于復(fù)雜的山區(qū)、森林地形具有較好的適應(yīng)性，但對(duì)種子點(diǎn)的依賴度較高。

PTIN濾波方法具有較好的適用性，特別是對(duì)于山區(qū)、溝壑區(qū)效果較好，但其種子點(diǎn)的選取是濾波的關(guān)鍵［23］。特別是對(duì)于煤層開(kāi)采導(dǎo)致地表出現(xiàn)階梯裂縫區(qū)域種子點(diǎn)的選取，如種子點(diǎn)選取在裂縫或塌陷坑底部會(huì)造成局部濾波錯(cuò)誤。因此，針對(duì)開(kāi)采沉陷區(qū)域梯度變形和偽最低噪聲點(diǎn)問(wèn)題，對(duì)種子點(diǎn)選取方法進(jìn)行改進(jìn)。

2.2 PTIN種子點(diǎn)選取及濾波

針對(duì)礦區(qū)，特別是淺埋煤層、大采高開(kāi)采造成的階梯塌陷地表，無(wú)人機(jī)LiDAR漸進(jìn)三角網(wǎng)點(diǎn)云濾波過(guò)程中種子點(diǎn)的選取至關(guān)重要。因此，針對(duì)礦區(qū)階梯狀地表，在漸進(jìn)三角網(wǎng)濾波算法的基礎(chǔ)上，進(jìn)行種子點(diǎn)選取方法的改進(jìn)，去除錯(cuò)誤的種子點(diǎn)，提高種子點(diǎn)質(zhì)量，計(jì)算過(guò)程如圖5所示。

格網(wǎng)劃分主要是用規(guī)則的虛擬格網(wǎng)覆蓋全部點(diǎn)云數(shù)據(jù)，將點(diǎn)云按照平面坐標(biāo)分配到相應(yīng)的虛擬網(wǎng)格內(nèi)。種子點(diǎn)的選取是濾波的關(guān)鍵，采用改進(jìn)的擴(kuò)展局部最小值方法進(jìn)行種子點(diǎn)選取，原理如圖6所示。

傳統(tǒng)種子點(diǎn)選取主要采用格網(wǎng)內(nèi)最低點(diǎn)作為種子點(diǎn)，圖6中采用傳統(tǒng)方法選取的種子點(diǎn)應(yīng)為P1。但P1相較于格網(wǎng)內(nèi)其他點(diǎn)，高程明顯過(guò)于小，屬于最低點(diǎn)噪聲數(shù)據(jù)。文中方法將格網(wǎng)內(nèi)激光點(diǎn)按高程值由小到大進(jìn)行排序，并計(jì)算各自與后續(xù)點(diǎn)之間的高程差值（圖6）。將高差平均值作為閾值。若當(dāng)前點(diǎn)與后點(diǎn)的高差小于高差閾值，即

di＜（i=1，2，…）

，則當(dāng)前點(diǎn)為種子點(diǎn)，若大于則計(jì)算第2最低點(diǎn)與第3最低點(diǎn)高差，逐步計(jì)算，直到獲取滿足高差閾值要求的數(shù)據(jù)點(diǎn)，并將該點(diǎn)作為種子點(diǎn)。該方法能夠確保選擇的最低點(diǎn)為地面點(diǎn)，避免噪聲點(diǎn)、異常點(diǎn)的影響，如圖6中噪聲點(diǎn)P1。采用文中算法選取種子點(diǎn)，然后采用PTIN算法進(jìn)行地面點(diǎn)提取，結(jié)果如圖7所示。

對(duì)2020年11月和2021年5月點(diǎn)云數(shù)據(jù)分別進(jìn)行點(diǎn)云濾波，并對(duì)濾波結(jié)果進(jìn)行統(tǒng)計(jì)，地面點(diǎn)云占比分別為93%和83%，由于不同時(shí)期地表植被覆蓋度不同，導(dǎo)致了地面點(diǎn)云占比不同。

2.3 點(diǎn)云插值計(jì)算

插值算法優(yōu)劣直接影響地面DEM精度，目前常用的插值算法包括反距離權(quán)重（Inverse Distance Weighted，IDW）、克里金（Kriging）、自然鄰域（Natural Neighbor，NN）、樣條函數(shù)（Spline）、不規(guī)則三角網(wǎng)（Triangulated Irregular Network，TIN）等［18］。IDW插值以采樣點(diǎn)與插值點(diǎn)間的距離為權(quán)重計(jì)算加權(quán)平均值，離插值點(diǎn)越近的采樣點(diǎn)權(quán)重越大，超過(guò)一定距離時(shí)權(quán)重忽略不計(jì)，適用于采樣點(diǎn)均勻且密集的區(qū)域［24］；Kriging插值與IDW類似，該方法采用半變異函數(shù)作為權(quán)重進(jìn)行無(wú)偏最優(yōu)估計(jì)，不僅考慮距離權(quán)重，還兼顧空間分布關(guān)系，適用于空間相關(guān)性較好的采樣點(diǎn)，但容易出現(xiàn)邊緣效應(yīng)［25］；自然鄰域插值是基于Voronoi結(jié)構(gòu)的插值計(jì)算，該方法對(duì)局部特征的繼承性較好，但局部點(diǎn)缺失時(shí)容易造成失真［26］；Spline以最小曲率面逼近各采樣點(diǎn)，以獲取最佳平滑曲面，顧及大范圍的采樣點(diǎn)，適用于復(fù)雜高曲率曲面，但平滑易造成地形失真［26］；不規(guī)則三角網(wǎng)（TIN）利用采樣點(diǎn)構(gòu)建不規(guī)則三角形，能夠在地形復(fù)雜區(qū)域保留地形特征細(xì)節(jié)，但其邊緣梯度可能存在突然變化［25］。

文獻(xiàn)［18］的研究成果表明，在黃土高原溝壑區(qū)，Kriging插值在無(wú)人機(jī)LiDAR點(diǎn)云插值計(jì)算中具有較好效果。因此，采用Kriging插值算法對(duì)地面點(diǎn)云進(jìn)行插值計(jì)算，并將插值格網(wǎng)大小對(duì)結(jié)果精度的影響進(jìn)行試驗(yàn)。采用DEM誤差分布和誤差值大小對(duì)比的方法進(jìn)行結(jié)果分析，不同格網(wǎng)大小引起的DEM誤差分布如圖8所示。

同時(shí)，計(jì)算不同格網(wǎng)尺寸DEM的ME，MAE，RMSE，結(jié)果見(jiàn)表1。

由表1可知，ME，MAE和RMSE在0.05～1 m區(qū)間內(nèi)變化不大，當(dāng)格網(wǎng)尺寸大于1 m時(shí)，誤差急劇增大。同時(shí)，結(jié)合圖8綜合考慮，最終選擇0.1 m的格網(wǎng)尺寸進(jìn)行插值計(jì)算，獲取較高精度的研究區(qū)DEM。

2.4 DEM差值計(jì)算

沉陷模型是由多時(shí)相地面DEM差值計(jì)算獲得，與DEM結(jié)構(gòu)和功能類似，能夠反映一段時(shí)間內(nèi)地表高程變化，即地表沉陷量，是一種高精度地表沉陷表示模型［13］。通過(guò)不同時(shí)相DEM之間的差值計(jì)算，可獲取該段時(shí)間內(nèi)地表的下沉變化。以首次觀測(cè)時(shí)獲取的DEM為基礎(chǔ)，與開(kāi)采結(jié)束后獲取的DEM做減法操作，可得到該段時(shí)間內(nèi)地表的沉陷模型，其三維效果如圖9所示，由于地表人為破壞、建筑施工等，導(dǎo)致沉陷量部分異常值過(guò)大或過(guò)小。

3 精度評(píng)定

基于地面監(jiān)測(cè)樁實(shí)測(cè)數(shù)據(jù)對(duì)地面點(diǎn)云、DEM和沉陷模型進(jìn)行精度評(píng)定。由于點(diǎn)云由一系列離散點(diǎn)構(gòu)成，無(wú)法進(jìn)行點(diǎn)對(duì)點(diǎn)差值直接計(jì)算，因此采用局部平均值的方法進(jìn)行差值計(jì)算，計(jì)算監(jiān)測(cè)樁附近局部點(diǎn)云的高程平均值，然后與實(shí)測(cè)數(shù)據(jù)對(duì)比，結(jié)果見(jiàn)表2，2020年11月LiDAR測(cè)量點(diǎn)云中誤差為60.6 mm，2021年5月為59.9 mm。

其次，基于監(jiān)測(cè)樁實(shí)測(cè)數(shù)據(jù)對(duì)DEM和沉陷模型誤差進(jìn)行計(jì)算，誤差分布如圖10所示。

從圖10可以看出，原始地面點(diǎn)云、DEM和沉陷模型三者中DEM的精度較差，沉陷模型精度最優(yōu)。在無(wú)人機(jī)LiDAR航測(cè)過(guò)程中，多時(shí)相數(shù)據(jù)采集時(shí)使用相同航測(cè)參數(shù)，能夠消除一定的系統(tǒng)誤差，如激光掃描設(shè)備系統(tǒng)誤差、飛行高度、地形變化等因素的影響，提高沉陷模型的精度。

基于沉陷模型提取地面監(jiān)測(cè)樁的下沉量，與地表監(jiān)測(cè)點(diǎn)監(jiān)測(cè)值進(jìn)行對(duì)比，繪制下沉折線圖，如圖11所示。

在2020年11月至2021年5月期間，地面水準(zhǔn)測(cè)量獲取的最大下沉值為1 872 mm，最小下沉值為0 mm，如圖11中折線所示；基于沉陷模型獲取的最大下沉值為1 860 mm；除A203監(jiān)測(cè)點(diǎn)實(shí)測(cè)數(shù)據(jù)下沉異常，實(shí)測(cè)下沉折線圖與基于沉陷模型獲取的下沉折線圖走勢(shì)和拐點(diǎn)基本一致。相同位置實(shí)測(cè)下沉量與沉陷模型下沉量的誤差絕對(duì)值最大為82 mm，最小為3 mm，平均為47 mm，中位數(shù)為42.5 mm，其中小于50 mm的占比60%，誤差絕對(duì)值整體小于100 mm，誤差分布相對(duì)集中。

4 基于沉陷模型的數(shù)據(jù)挖掘與分析

沉陷模型包含豐富的地表下沉信息，包括地表平面位置、任意位置下沉量、下沉變化等，基于沉陷模型進(jìn)行數(shù)據(jù)挖掘分析，能夠獲取開(kāi)采沉陷引起的下沉范圍、邊界角、下沉面積、不同程度損害區(qū)域面積等。

4.1 下沉量統(tǒng)計(jì)

為了準(zhǔn)確獲取地表不同位置的下沉量和不同下沉量的分布情況，對(duì)沉陷模型進(jìn)行下沉量分級(jí)分類統(tǒng)計(jì)，結(jié)果如圖12所示。

圖12表示不同下沉量區(qū)域分布情況，下沉區(qū)域呈橢圓形分布，最大下沉區(qū)位于采空區(qū)中心，由中心逐漸向外擴(kuò)散，下沉值逐漸減小，相應(yīng)的下沉面積逐漸增大。圖中圓形標(biāo)注區(qū)域出現(xiàn)異常下沉，經(jīng)實(shí)地調(diào)查發(fā)現(xiàn)是由于煤層開(kāi)采期間，人為施工造成。

對(duì)下沉區(qū)域面積進(jìn)行分級(jí)計(jì)算，以大于100 mm的下沉值作為開(kāi)采沉陷引起的地表下沉，計(jì)算結(jié)果見(jiàn)表3。下沉值大于100 mm的下沉總面積為180 075 m2，采空區(qū)面積為49 789 m2。為準(zhǔn)確確定開(kāi)采沉陷引起的不同下沉量所占面積，暫定以下沉值大于100 mm的沉陷量作為煤層開(kāi)采引起的地表下沉。結(jié)果表明，下沉值小于300 mm的面積占比超過(guò)一半以上，隨著下沉量的增大，對(duì)應(yīng)的下沉面積逐漸減小。

4.2 下沉剖面信息提取

傳統(tǒng)開(kāi)采沉陷監(jiān)測(cè)采用預(yù)埋監(jiān)測(cè)樁的方式，監(jiān)測(cè)煤層開(kāi)采過(guò)程中地表變化，該方法獲取的監(jiān)測(cè)數(shù)據(jù)少、范圍小。試驗(yàn)基于沉陷模型采用虛擬監(jiān)測(cè)點(diǎn)和剖面線的形式，分別對(duì)沉陷區(qū)域進(jìn)行走向和傾向方向下沉信息提取，虛擬監(jiān)測(cè)點(diǎn)與剖面線位置如圖13所示。

虛擬監(jiān)測(cè)點(diǎn)以15 m間隔在走向和傾向方向均勻分布，如圖13中黑色三角形所示；剖面線分為走向剖面線和傾向剖面線，位置如圖13中灰色實(shí)線所示?；诔料菽Ｐ驮谔摂M監(jiān)測(cè)點(diǎn)位提取地表下沉值，并進(jìn)行折線圖繪制，結(jié)果如圖14所示。虛擬監(jiān)測(cè)點(diǎn)最大下沉值為1 761 mm，走向下沉折線呈半盆地分布，傾向下沉折線圖呈全盆地分布，符合開(kāi)采沉陷下沉規(guī)律。

根據(jù)沉陷模型分辨率大小，下沉剖面線能夠提取大量監(jiān)測(cè)點(diǎn)數(shù)據(jù)，相比離散監(jiān)測(cè)點(diǎn)具有更好的魯棒性。因此，在沉陷區(qū)域設(shè)置一條長(zhǎng)740 m的走向剖面線和長(zhǎng)413 m的傾向剖面線，基于剖面線以沉陷模型分辨率0.1 m為間隔，提取下沉值。對(duì)提取結(jié)果進(jìn)行散點(diǎn)繪制，并提取中心線，結(jié)果如圖15所示。

基于沉陷模型提取的走向下沉曲線呈盆地形狀，說(shuō)明在走向方向達(dá)到超充分采動(dòng)，符合開(kāi)采沉陷一般規(guī)律；傾向方向下沉曲線呈山谷對(duì)稱分布，其中存在部分異常值，經(jīng)實(shí)地探查，造成該異常的主要原因是由于地表人為開(kāi)挖施工，去除異常值，整體符合開(kāi)采沉陷地表下沉特征。

4.3 沉陷面域分析

傳統(tǒng)開(kāi)采沉陷監(jiān)測(cè)主要對(duì)固定點(diǎn)進(jìn)行監(jiān)測(cè)，基于沉陷模型不僅能夠進(jìn)行傳統(tǒng)點(diǎn)監(jiān)測(cè)、剖面線監(jiān)測(cè)，獲取下沉值和下沉折線圖。同時(shí)也能夠進(jìn)行面域監(jiān)測(cè)分析，進(jìn)一步提取下沉等值線、下沉范圍、地表?yè)p害區(qū)域圖等。

等值線是一種由相同值連接構(gòu)成的曲線，如等高線、等溫線、等深線等。文中采用等值線表示地表下沉值，能夠清晰地表示地表下沉情況。等值線的精度一般取決于等值線的間隔值，間隔值越小，圖形越精細(xì)、等值線越密、表達(dá)的效果越好。但是等值線間隔值越小，越容易受異常值和離散數(shù)的影響，圖形越復(fù)雜、異常數(shù)據(jù)影響越明顯。因此，合適的等值線間隔能夠完美表達(dá)地表下沉情況，同時(shí)能夠降低異常值的影響?；诠ぷ髅嫦鲁林荡笮『筒煌鲁林档姆植紖^(qū)間，采用多組、分梯度的方法設(shè)置等值線間隔，分別選用30，50，100 mm和200 mm間隔值進(jìn)行下沉等值線繪制，結(jié)果如圖16所示。

從圖16可以看出，不同等值距的下沉等值線所表達(dá)的特征詳細(xì)度不同。較小的等值線區(qū)間值對(duì)應(yīng)較詳細(xì)的沉陷特征，同時(shí)也易受噪聲和異常值的影響。通過(guò)綜合比較發(fā)現(xiàn)，沉降間隔為100 mm的等值線既可以保留采煤沉陷特征細(xì)節(jié)，也能夠降低噪聲和異常值的影響。

傳統(tǒng)開(kāi)采沉陷以邊界角確定煤層開(kāi)采影響范圍，但由于西部煤礦地形復(fù)雜，不同地形下沉特征有所區(qū)別，走向和傾向觀測(cè)線獲取的邊界角難以準(zhǔn)確描述所有地形條件下的影響邊界。因此，基于沉陷模型進(jìn)行開(kāi)采邊界的提取能夠彌補(bǔ)傳統(tǒng)測(cè)線的不足。但由于地表自然變化（如雨水沖刷、風(fēng)沙侵蝕等）和無(wú)人機(jī)LiDAR測(cè)量精度的影響，獲取10 mm下沉邊界存在一定困難。因此，考慮自然因素影響，以100 mm下沉值作為煤層開(kāi)采引起的下沉邊界，提取地表下沉邊界，如圖17所示。

從圖17可以看出，采空區(qū)A開(kāi)采引起的沉陷區(qū)域面積140 271 m2，小于表4中計(jì)算的下沉面積180 075 m2，造成該結(jié)果的原因是由于研究區(qū)內(nèi)存在人工開(kāi)挖區(qū)域和遠(yuǎn)離采空區(qū)的下沉面積，該面積并不是由于煤層開(kāi)采造成的地表下沉。下沉區(qū)域整體偏向采空區(qū)B，造成該結(jié)果的原因是在2020年11月無(wú)人機(jī)LiDAR航測(cè)前，采空區(qū)B已存在，其對(duì)地表造成的影響并未結(jié)束，因此后續(xù)地表下沉范圍偏向采空區(qū)B。無(wú)人機(jī)LiDAR監(jiān)測(cè)期間，地表下沉面積與采空區(qū)面積比例為2.8∶1，該比例能夠?yàn)楹罄m(xù)煤層開(kāi)采地表影響范圍的確定提供參考。

為了明晰煤層開(kāi)采造成地表下沉的范圍和嚴(yán)重程度，以不同下沉量所占下沉面積為參數(shù)，進(jìn)行地表?yè)p害的分等定級(jí)，將地表下沉值進(jìn)行分組，并確定其范圍和對(duì)應(yīng)面積，結(jié)果如圖18所示，可以看出，不同沉降等級(jí)區(qū)域由內(nèi)向外逐漸擴(kuò)展，采空區(qū)從中心向外延伸。沉降可分為6級(jí)：100～300，300～600，600～900，900～1 200，1 200～1 500 mm和1 500 mm以上，相應(yīng)的沉降區(qū)面積為65 101，29 083，16 354，13 178，11 123 m2和5 432 m2。沉降面積之比為46∶21∶12∶9∶8∶4，能為工作面開(kāi)采過(guò)程中地表建（構(gòu)）筑物的保護(hù)、移民搬遷提供依據(jù)。

5 結(jié) 論

1）在煤礦地形復(fù)雜區(qū)域，采用無(wú)人機(jī)LiDAR進(jìn)行煤層開(kāi)采沉陷監(jiān)測(cè)，針對(duì)沉陷區(qū)地形特征，對(duì)PTIN濾波算法種子點(diǎn)的選取方法進(jìn)行改進(jìn)，剔除了最低點(diǎn)噪聲，實(shí)現(xiàn)無(wú)人機(jī)LiDAR地面點(diǎn)云的穩(wěn)健提取。采用Kriging插值算法對(duì)地面點(diǎn)云進(jìn)行內(nèi)插生成DEM，對(duì)格網(wǎng)尺寸進(jìn)行試驗(yàn)，確定最優(yōu)的格網(wǎng)尺寸大?。?.1 m），保證DEM誤差分布和精度均處于最優(yōu)，獲取較高精度的地面DEM和沉陷模型。

2）無(wú)人機(jī)LiDAR測(cè)量的高程精度可達(dá)25～75 mm，受地面噪聲、飛行條件、地形地貌的影響，DEM的精度可達(dá)25～100 mm，由于相似誤差的抵償，沉陷模型精度可達(dá)3～82 mm。因此，研究區(qū)基于無(wú)人機(jī)LiDAR獲取的沉陷模型整體精度優(yōu)于100 mm，平均誤差小于47 mm。

3）根據(jù)沉陷模型下沉特征進(jìn)行多方位分析，對(duì)沉陷模型進(jìn)行挖掘和分析，獲取地表下沉折線圖、下沉等值線圖、地表?yè)p害分布圖、地表下沉面積與采空區(qū)比例關(guān)系等數(shù)據(jù)，能夠?yàn)殚_(kāi)采沉陷預(yù)計(jì)和地表建構(gòu)筑物保護(hù)提供依據(jù)。

參考文獻(xiàn)（References）：

［1］ 王雙明，段中會(huì)，馬麗，等.西部煤炭綠色開(kāi)發(fā)地質(zhì)保障技術(shù)研究現(xiàn)狀與發(fā)展趨勢(shì)［J］.煤炭科學(xué)技術(shù)，2019，47（2）：1-6.WANG Shuangming，DUAN Zhonghui，MA Li，et al.Ｒesearch status and future trends of geological assurance technology for coal green development in Western China［J］.Coal Science and Technology，2019，47（2）：1-6.

［2］王冠，劉靜.煤炭經(jīng)濟(jì)發(fā)展的綠色向度與綠色道路［J］.人民論壇，2020（29）：70-71.WANG guan，LIU Jing.Green dimension and green road of coal economic development［J］.Peoples Tribune，2020（29）：70-71.

［3］蘆家欣，湯伏全，趙軍儀，等.黃土礦區(qū)開(kāi)采沉陷與地表?yè)p害研究述評(píng)［J］.西安科技大學(xué)學(xué)報(bào)，2019，39（5）：859-866.LU Jiaxin，TANG Fuquan，ZHAO Junyi，et al．Review of study on mining subsidence and ground surface damage in loess mining area［J］．Journal of Xian University of Science and Technology，2019，39（5）：859-866．

［4］余學(xué)義，王昭舜，楊云.大采深綜放開(kāi)采地表移動(dòng)變形規(guī)律［J］.西安科技大學(xué)學(xué)報(bào)，2019，39（4）：555-563.YU Xueyi，WANG Zhaoshun，YANG Yun．Surface movement and deformation law of fully-mechanized caving mining with large mining depth［J］．Journal of Xian University of Science and Technology，2019，39（4）：555-563．

［5］許家林，錢(qián)鳴高.綠色開(kāi)采的理念與技術(shù)框架［J］.科技導(dǎo)報(bào)，2007，25（7）：61-65.XU Jialin，QIAN Minggao.Concept of green mining and its technical framework［J］.Science and Technology Herald，2007，25（7）：61-65.

［6］張海君，吳奕樞，王飛，等.準(zhǔn)格爾煤田龍王溝煤礦特厚煤層綜放開(kāi)采地表沉陷規(guī)律［J］.西安科技大學(xué)學(xué)報(bào)，2022，42（5）：874-883.ZHANG Haijun，WU Yishu，WANG Fei，et al.Surface subsidence law in fully-mechanized caving mining of extra-thick seam in Longwanggou coal mine of Zhungeer coalfield［J］.Journal of Xian University of Science and Technology，2022，42（5）：874-883.

［7］湯伏全，董龍凱，王宗良，等.基于沉陷對(duì)稱特征的近水平煤層開(kāi)采InSAR三維位移反演模型［J］.煤炭學(xué)報(bào)，2019，44（1）：210-220.TANG Fuquan，DONG Longkai，WANG Zongliang，et al．A 3D inversion model for InSAＲ detected displacements based on ground subsidence symmetry induced by horizontal coal mining［J］．Journal of China Coal Society，2019，44（1）：210-220．

［8］李海濱，唐國(guó)茜，趙桂娟，等.三維激光掃描技術(shù)在路堤沉降觀測(cè)中的應(yīng)用［J］.西安科技大學(xué)學(xué)報(bào)，2021，41（1）：87-93.LI Haibin，TANG Guoxi，ZHAO Guijuan，et al.Application of 3D laser scanning technology in the observation of embankment settlements［J］．Journal of Xian University of Science and Technology，2021，41（1）：87-93．

［9］MATWIJ W，W GRUSZCZYSKI，PUNIACH E，et al.Determination of underground mining-induced displacement field using multi-temporal TLS point cloud registration［J］.Measurement，2021（4）：109482.

［10］LIAN X，LIU X，GE L，et al.Time-series unmanned aerial vehicle photogrammetry monitoring method without ground control points to measure mining subsidence［J］.Journal of Applied Remote Sensing，2021，15（2）：24505.

［11］詹容若，段亮，羅曉容，等.無(wú)人機(jī)多點(diǎn)位航拍高分辨率三維數(shù)字露頭建模［J］.西安科技大學(xué)學(xué)報(bào)，2021，41（6）：1050-1058.ZHAN Rongruo，DUAN Liang，LUO Xiaorong，et al.3D digital outcrop modeling with high resolution using drone-based multi-point photography［J］．Journal of Xian University of Science and Technology，2021，41（6）：1050-1058．

［12］何柯璐，湯伏全，韋書(shū)平，等.機(jī)載激光點(diǎn)云構(gòu)建西部礦區(qū)開(kāi)采沉陷模型研究［J］.測(cè)繪科學(xué)，2021，46（2）：130-138.HE Kelu，TANG Fuquan，WEI Shuping，et al.Model research on mining subsidence of western mining area based on airborne laser point cloud［J］.Science of Surveying and Mapping，2021，46（2）：130-138.

［13］ZHENG J，YAO W，LIN X，et al.An accurate digital subsidence model for deformation detection of coal mining areas using a UAV-based LiDAR［J］.Remote Sen-sing，2022，14（2）：421.

［14］ALI M，TAHA E D，MOHAMED M，et al.Generation of digital terrain model from multispectral LiDar using different ground filtering techniques［J］.Egyptian Journal of Remote Sensing and Space Science，2021，24（2）：181-189.

［15］王競(jìng)雪，張雪洋，洪紹軒，等.結(jié)合形態(tài)學(xué)和TIN三角網(wǎng)的機(jī)載LiDAR點(diǎn)云濾波算法［J］.測(cè)繪科學(xué)，2019，44（5）：151-156，183.WANG Jingxue，ZHANG Xueyang，HONG Shaoxuan，et al.Aerial LiDAR point cloud filtering algorithm combining mathematical morphology and TIN［J］.Science of Surveying and Mapping，2019，44（5）：151-156，183.

［16］甘斌，鄭俊良，姚頑強(qiáng)，等.基于數(shù)學(xué)形態(tài)學(xué)的機(jī)載LiDAR采煤區(qū)沉陷信息提取［J］.西安科技大學(xué)學(xué)報(bào)，2023，43（1）：175-182.GAN Bin，ZHENG Junliang，YAO Wanqiang，et al.Subsidence information extraction of mining area by airborne LiDAR based on mathematical morphology［J］.Journal of Xian University of Science and Technology，2023，43（1）：175-182.

［17］HUI Z，JIN S，XIA Y，et al.A mean shift segmentation morphological filter for airborne LiDAR DTM extraction under forest canopy［J］.Optics & Laser Technology，2021，136（1）：106728.

［18］李朋飛，張曉晨，嚴(yán)露，等.復(fù)雜地形中機(jī)載LiDAR點(diǎn)云構(gòu)建DEM的插值算法對(duì)比［J］.農(nóng)業(yè)工程學(xué)報(bào)，2021，37（15）：146-153.LI Pengfei，ZHANG Xiaochen，YAN Lu，et al.Comparison of interpolation algorithms for DEMs in topographically complex areas using airborne LiDAR point clouds［J］.Transactions of the Chinese Society of Agricultural Engineering，2021，37（15）：146-153.

［19］VOSSELMAN G.Slope based filtering of laser altimetry data［J］.International Archives of Photogrammetry & Remote Sensing，2000，33（B3/2；PART3）：935-942.

［20］隋立春，張熠斌，張碩，等.基于漸進(jìn)三角網(wǎng)的機(jī)載LiDAR點(diǎn)云數(shù)據(jù)濾波［J］.武漢大學(xué)學(xué)報(bào)（信息科學(xué)版），2011，36（10）：1159-1163.SUI Lichun，ZHANG Yibin，ZHANG Shuo，et al.Filtering of airborne LiDAR point cloud data based on progressive TIN［J］.Geomatics and Information Science of Wuhan University，2011，36（10）：1159-1163.

［21］MORGAN M，TEMPFLI K.Automatic building extraction from airborne laser scanning data［J］.International Archives of Photogrammetry & Remote Sensing，2000

，33（B3/2；PART3）：616-623.

［22］ZHANG W，QI J，WAN P，et al.An easy-to-use Aairborne LiDAR data filtering method based on cloth simulation［J］.Remote Sensing，2016，8（6）：501.

［23］李卉，李德仁，黃先鋒，等.一種漸進(jìn)加密三角網(wǎng)LIDAR點(diǎn)云濾波的改進(jìn)算法［J］.測(cè)繪科學(xué)，2009，34（3）：39-40，216.LI Hui，LI Deren，HUANG Xianfeng，et al.Advanced adaptive TIN filter for LIDAR point clouds data［J］.Science of Surveying and Mapping，2009，34（3）：39-40，216.

［24］張孟丹，余鐘波，谷黃河，等.無(wú)定河流域降水量空間插值方法比較研究［J］.人民黃河，2021，43（4）：30-37，99.ZHANG Mengdan，YU Zhongbo，

GU Huanghe，et al.Comparative study on spatial interpolation methods of precipitation in Wuding River basin［J］.Peoples Yellow River，2021，43（4）：30-37，99.

［25］LI J，HEAP A D.Spatial interpolation methods applied in the environmental sciences：A review［J］.Environmental Modelling & Software，2014，53（Mar）：173-189.

［26］馮波，陳明濤，岳冬冬，等.基于兩種插值算法的三維地質(zhì)建模對(duì)比［J］.吉林大學(xué)學(xué)報(bào)（地球科學(xué)版），2019，49（4）：1200-1208.FENG Bo，CHEN Mingtao，YUE Dongdong，et al.Comparison of 3D geological modeling based on two interpolation algorithms［J］.Journal of Jilin University（Earth Science Edition），2019，49（4）：1200-1208.

（責(zé)任編輯：高佳）