基于數(shù)學(xué)形態(tài)學(xué)的機(jī)載LiDAR采煤區(qū)沉陷信息提取

摘 要：由于采煤沉陷過程復(fù)雜和地表地形影響，機(jī)載LiDAR在采煤區(qū)沉陷監(jiān)測中不可避免地存在噪聲數(shù)據(jù)，高密度LiDAR點云中存在的噪聲容易導(dǎo)致提取的沉陷等值線出現(xiàn)鋸齒、毛邊和多邊形碎屑等問題。將數(shù)學(xué)形態(tài)學(xué)算法格網(wǎng)初始值的判定方式進(jìn)行改進(jìn)，傳統(tǒng)數(shù)學(xué)形態(tài)學(xué)算法選擇格網(wǎng)內(nèi)最低點高程作為格網(wǎng)值，改進(jìn)算法對格網(wǎng)內(nèi)所有高程值進(jìn)行平面擬合，將擬合值作為格網(wǎng)初始值。在采煤沉陷信息提取過程，增加對地面點云的改進(jìn)數(shù)學(xué)形態(tài)學(xué)算法處理，降低噪聲數(shù)據(jù)對地面DEM的影響，提高沉陷DEM精度和沉陷等值線完整度，試驗對比分析確定算法最優(yōu)參數(shù)（格網(wǎng)大小為3.5 m，結(jié)構(gòu)元素尺寸為3 m）。最后，采用該方法對研究區(qū)數(shù)據(jù)進(jìn)行處理，獲取研究區(qū)沉陷DEM，并進(jìn)行數(shù)據(jù)分析挖掘，獲取地表下沉范圍、下沉等值線、下沉面積等。結(jié)果表明：改進(jìn)算法既保證原始的地形特征和精度，又可消除沉陷等值線中出現(xiàn)的噪聲問題，為開采沉陷預(yù)計及采后環(huán)境評估提供支撐。關(guān)鍵詞：LiDAR；點云；數(shù)學(xué)形態(tài)學(xué)；沉陷監(jiān)測；沉陷DEM中圖分類號：TD 325

Subsidence information extraction of mining area by airborne LiDAR based on mathematical morphology

GAN Bin¹，ZHENG Junliang²，YAO Wanqiang²，BAI Lingxiao¹

（1.Xi’an Surveying and Mapping Institute，Xi’an 710054，China；2.College of Geomatics，Xi’an University of Science and Technology，Xi’an 710054，China）

Abstract：Due to the complexity of coal mining subsidence process and the impact of surface topography，airborne LiDAR inevitably has noise data in coal mining subsidence monitoring.The noise in LiDAR high-density point cloud was easy to trigger the problems such as serrations，burr and polygonal debris in the extracted subsidence isoline.In this paper，the method of determining the initial value of grid in mathematical morphology algorithm is improved.The traditional mathematical morphology algorithm selected the lowest point elevation in the grid as the grid value，and the improved algorithm performs plane fitting for all the elevation values in the grid，and takes the fitting value as the initial value of the grid.In the process of mining subsidence information extraction，the processing of ground point cloud by the improved mathematical morphology algorithm was added，thus reducing the impact of noise data on DEM，and promoting the precision of subsidence DEM as well as the integrity of subsidence isoline.The optimal parameters of the algorithm was determined by the experimental" analysis and comparison（the grid size is 3.5 m，and the size of structural elements is 3 m）.Finally，the data of the study area were processed by the algorithm improved in this paper，and the subsidence DEM in the study area was obtained.The subsidence DEM was examined and excavated to obtain the surface subsidence range，subsidence contour，subsidence area，etc.The results show that the improved algorithm can not only ensure the original terrain characteristics and accuracy，but also eliminate the noise in the subsidence isoline，providing support for mining subsidence prediction and post mining environment assessment.

Key words：LiDAR；point cloud；mathematical morphology；subsidence monitoring；subsidence DEM

0 引 言

地下礦產(chǎn)資源開采會破壞區(qū)域局部受力平衡和水文條件，導(dǎo)致采空區(qū)發(fā)生塌陷，出現(xiàn)地裂縫、塌陷坑，甚至引發(fā)山體滑坡等地質(zhì)災(zāi)害^[1-3]。因此，對資源開發(fā)，特別是煤炭資源開發(fā)造成的地表塌陷進(jìn)行快速、高效監(jiān)測是綠色礦山建設(shè)的重要任務(wù)。機(jī)載LiDAR（Light Detection and Ranging）測量系統(tǒng)能夠快速獲取地表大面積高精度點云數(shù)據(jù)，相對于傳統(tǒng)水準(zhǔn)儀、全站儀、GNSS沉陷監(jiān)測，機(jī)載LiDAR具有高效、快速、面域觀測的優(yōu)勢。對井田開采引發(fā)的地表沉陷區(qū)域持續(xù)監(jiān)測，可快速取得地表多時相三維模型，進(jìn)而獲取地表沉陷分布和沉陷DEM^[4]。機(jī)載LiDAR觀測可針對工作面整體進(jìn)行動態(tài)監(jiān)測，而非傳統(tǒng)開采沉陷中走向、傾向觀測線的點式觀測。其監(jiān)測結(jié)果更為全面、直觀，并可為概率積分法預(yù)計參數(shù)求取和地表移動規(guī)律分析提供大量高精度數(shù)據(jù)^[5-7]。

機(jī)載LiDAR掃描的海量高密度點云能夠更加真實地反映地表情況，尤其是開采沉陷過后地形復(fù)雜多變，沉陷DEM能更準(zhǔn)確反應(yīng)地表沉陷狀況^[8]，也能精確提取任意點位的沉陷信息^[9-11]。但是，機(jī)載LiDAR獲取沉陷DEM需經(jīng)過點云預(yù)處理、校正、濾波、分類、插值等多個過程，數(shù)據(jù)本身和處理過程必然會產(chǎn)生一定誤差，例如機(jī)載LiDAR自身誤差、點云濾波誤差、插值計算誤差，進(jìn)而影響沉陷模型的精度。點云噪聲的存在，影響觀測面下沉信息的提取和分析，例如生成采煤下沉邊界不準(zhǔn)確，等值線存在鋸齒、毛邊、破碎線段和碎屑多邊形^[12-14]，導(dǎo)致下沉面積和體積計算不準(zhǔn)確，進(jìn)一步影響開采沉陷面域分析的其他工作。點云密度越高表達(dá)的地面信息越精細(xì)，包含的噪聲數(shù)據(jù)也會增加，在沉陷等值線提取和下沉面積統(tǒng)計時將會產(chǎn)生更大的影響。

目前，許多學(xué)者對機(jī)載LiDAR等值線提取進(jìn)行了廣泛研究。在提取LiDAR點云中的地形信息常用方法中，鄭俊良等對無人機(jī)LiDAR獲取的點云進(jìn)行分類濾波，并提取下沉等值線和下沉面積統(tǒng)計圖^[15]；牛峰明對下沉盆地剖面地形進(jìn)行擬合，該方法存在擬合過度問題，并且局限于剖面分析，即觀測線下沉分析^[16]；何睦通過特征點提取和修改分辨率的方法去除等值線中的碎屑多邊形，該方法在高密度點云的下沉盆地中存在改變地形的情況，效果較差^[17]。針對此問題，筆者提出改進(jìn)數(shù)學(xué)形態(tài)學(xué)的DEM數(shù)據(jù)處理方法，在保證精度的情況下，去除DEM噪聲，提取最優(yōu)等值線，并將改進(jìn)算法應(yīng)用于煤層開采導(dǎo)致的地表沉陷信息提取，改善噪聲導(dǎo)致的等值線破碎、不平整等問題。

1 研究區(qū)概況與數(shù)據(jù)采集

1.1 研究區(qū)概況

涼水井煤礦位于陜西省神木市以西約16 km處，井田位于陜北黃土高原北部，毛烏蘇沙漠之南。東部為黃土梁峁溝谷地貌，西部為波狀沙丘地，地勢開闊，井田南部、北部黃土沖溝發(fā)育，梁峁區(qū)及沙丘區(qū)植被覆蓋良好。地勢總體呈西高東低、中部高南北低，最高處位于西部，標(biāo)高1 326.40 m，最低處位于東南角，標(biāo)高1 100.00 m，最大高差226.40 m。以涼水井煤礦某工作面局部為研究對象，研究區(qū)位置、巷道分布、機(jī)載LiDAR航測時間、工作面回采位置如圖1所示。

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

選取大疆M600無人機(jī)搭載SZT-R250機(jī)載激光雷達(dá)對涼水井工作面局部進(jìn)行掃描。工作面從東向西推進(jìn)，根據(jù)監(jiān)測要求，共進(jìn)行3期掃描。第1期為2020年8月15日，回采1 024.6 m；第2期為2020年11月7日，回采1 384.6 m；第3期為2021年5月19日，回采1 665.5 m。每期數(shù)據(jù)均采用相同的航測參數(shù)，飛行高度為70 m，飛行速度為8 m/s，掃描頻率設(shè)定為40 kHz，獲取的點云平均密度為60點/m²左右。采集數(shù)據(jù)時天氣晴朗，采集后經(jīng)過檢驗PDOP值、觀測衛(wèi)星數(shù)以及軌跡數(shù)據(jù)緊耦合處理，認(rèn)定數(shù)據(jù)質(zhì)量合格。

2 改進(jìn)的數(shù)學(xué)形態(tài)學(xué)算法

2.1 數(shù)學(xué)形態(tài)學(xué)算法原理

數(shù)學(xué)形態(tài)學(xué)的基本思想是使用具有特定形態(tài)

的結(jié)構(gòu)元素提取圖像中的拓?fù)浣Y(jié)構(gòu)信息，達(dá)到圖像分類的目的^[18]，其基本運算包括膨脹和腐蝕2種。其中，腐蝕指的是利用結(jié)構(gòu)元素遍歷目標(biāo)圖像，得到目標(biāo)圖像中所有完全包含該結(jié)構(gòu)元素的元素集合；膨脹指的是結(jié)構(gòu)元素在遍歷目標(biāo)圖像時，得到的所有使結(jié)構(gòu)元素與目標(biāo)圖像的交集非空的元素構(gòu)成的集合^[19]。通過腐蝕和膨脹2種基本操作可以實現(xiàn)如開運算、閉運算、形態(tài)學(xué)梯度、頂帽、黑帽等多種形態(tài)學(xué)操作。

腐蝕運算定義為

f?b={x|b（x）∈f} "（1）

膨脹運算定義為

f⊕b={x|b（x）∩f≠φ} "（2）

開運算公式為

fb=f?b⊕b "（3）

式中 f為輸入的目標(biāo)圖像；b為結(jié)構(gòu)元素，開運算操作可用來消除較小物體，平滑較大物體并在狹窄處分離物體^[20]，在平滑毛刺和邊界時又不明顯改變物體總的形狀和面積。

基于形態(tài)學(xué)點云處理主要采用開運算，對點云進(jìn)行先腐蝕后膨脹操作，主要步驟如圖2所示。首先將點云劃分為row*col個小柵格，用柵格內(nèi)最低點作為每個柵格的值，如圖2（a）和（b）所示；然后對各個柵格依次進(jìn)行腐蝕操作，遍歷當(dāng)前腐蝕柵格的周圍柵格，取最小值作為當(dāng)前柵格值，如圖2（c）所示；再進(jìn)行膨脹操作，取周圍柵格最大值作為該柵格膨脹后的值，如圖2（d）所示；最后遍歷各個柵格內(nèi)的點，如圖2（e）所示，與腐蝕、膨脹后的柵格值進(jìn)行比較，大于設(shè)定值則刪除該點。

數(shù)學(xué)形態(tài)學(xué)算法能較好去除噪聲數(shù)據(jù)，但是在去除點云冗余信息過程中仍存在一定問題，例如格網(wǎng)及結(jié)構(gòu)元素參數(shù)的確定，格網(wǎng)內(nèi)點云分類的判定方式等。傳統(tǒng)數(shù)學(xué)形態(tài)學(xué)算法在對點云進(jìn)行分類時，計算點膨脹后的高程值和點原始高程值之差的絕對值，大于閾值時則刪除該點，其主要用于點云濾波。而試驗點云為地面點，將格網(wǎng)內(nèi)最低點作為柵格值，且與膨脹后柵格的最大高程值做差會對原始地形產(chǎn)生改變，除此之外其參數(shù)選取的規(guī)則和范圍也不符合工作需求，因此筆者從這兩方面對數(shù)學(xué)形態(tài)學(xué)算法進(jìn)行改進(jìn)。實現(xiàn)在地形形態(tài)完整情況下，去除小范圍內(nèi)凸起或凹陷的噪聲。

2.2 數(shù)學(xué)形態(tài)學(xué)算法的改進(jìn)

2.2.1 改進(jìn)原理

針對數(shù)學(xué)形態(tài)學(xué)算法參數(shù)選取，程序設(shè)計為可任意調(diào)整如格網(wǎng)大小、殘差范圍、結(jié)構(gòu)元素尺寸等參數(shù)。根據(jù)經(jīng)驗任意調(diào)整數(shù)據(jù)處理參數(shù)可以避免傳統(tǒng)形態(tài)學(xué)參數(shù)閾值不適用于去除研究區(qū)地形冗余信息的工作。

將格網(wǎng)高程值的判定方式進(jìn)行改進(jìn)。傳統(tǒng)數(shù)學(xué)形態(tài)學(xué)算法是選擇格網(wǎng)內(nèi)最低點高程作為格網(wǎng)的值來進(jìn)行開運算，難以去除凹陷處的冗余地形。改進(jìn)算法將格網(wǎng)內(nèi)的點進(jìn)行平面擬合，將擬合值作為格網(wǎng)的初始值，再通過開運算得到膨脹后的高程。

2.2.2 改進(jìn)算法流程

改進(jìn)的數(shù)學(xué)形態(tài)學(xué)提取地表沉陷信息方法，基于VS 2017平臺，應(yīng)用QT框架開發(fā)，并優(yōu)化數(shù)據(jù)結(jié)構(gòu)，提高算法效率。算法流程如圖3所示。

1）將地面點云按照一定格網(wǎng)大小存儲并建立索引，同時對柵格中的空洞進(jìn)行修補(bǔ)。

2）對柵格內(nèi)的點進(jìn)行平面擬合，將擬合后平面的高程作為柵格的值。

3）根據(jù)地形特征選擇結(jié)構(gòu)元素的尺寸和形狀，用結(jié)構(gòu)元素遍歷整個柵格來進(jìn)行先腐蝕后膨脹的開運算。

4）對經(jīng)過開運算后的柵格值與初始柵格值做差，設(shè)定高差范圍，刪除超出該范圍的點。

5）對多期點云數(shù)據(jù)分別進(jìn)行步驟1）～4），實現(xiàn)地面點云的開運算，進(jìn)而構(gòu)建地面格網(wǎng)DEM。

6）按高程平均值法計算格網(wǎng)高程，將兩期DEM疊加做差，獲取沉陷DEM。

7）基于初始沉陷DEM，提取下沉等值線并進(jìn)行等值線平滑。

8）通過調(diào)整殘差大小、格網(wǎng)大小和結(jié)構(gòu)元素的形態(tài)，對比各控制點的真實下沉值，并將沉陷等值線進(jìn)行可視化分析，確定最優(yōu)參數(shù)。

3 改進(jìn)算法試驗

3.1 參數(shù)對比及算法評價

改進(jìn)的形態(tài)學(xué)算法中4個參數(shù)需要調(diào)整，包括格網(wǎng)大小、殘差范圍、結(jié)構(gòu)元素的尺寸和結(jié)構(gòu)元素的形態(tài)。合理的參數(shù)選擇可消除DEM上因隨機(jī)誤差出現(xiàn)的破碎多邊形，平滑大多數(shù)多邊形的邊界，使沉陷等值線更加符合實際地形和沉陷特征。該試驗通過分析各參數(shù)之間的聯(lián)系，對參數(shù)進(jìn)行分組試驗，對比分析。

3.1.1 結(jié)構(gòu)元素的形狀

數(shù)學(xué)形態(tài)學(xué)算法的本質(zhì)就是使用結(jié)構(gòu)元素對圖像進(jìn)行先腐蝕后膨脹的開運算處理，刪除影響地形特征的隨機(jī)誤差點。結(jié)構(gòu)元素可分為非對稱元素和對稱元素，經(jīng)試驗測試，非對稱元素會導(dǎo)致DEM發(fā)生一定的偏移情況，因此選擇對稱元素中的正方形作為結(jié)構(gòu)元素。

3.1.2 格網(wǎng)大小、結(jié)構(gòu)元素尺寸、殘差范圍

試驗?zāi)康氖峭ㄟ^沉陷DEM提取出準(zhǔn)確的沉陷等值線，并確定開采沉陷下沉影響范圍。因此，需去除導(dǎo)致等值線中出現(xiàn)的破碎獨立多邊形和鋸齒狀等值線的噪聲點。格網(wǎng)大小應(yīng)根據(jù)地形起伏和沉陷DEM下沉量大小設(shè)定。格網(wǎng)尺寸過大會造成地形特征點丟失；格網(wǎng)過小地形更加復(fù)雜精細(xì)，同時噪聲導(dǎo)致的DEM不確定誤差也會對沉陷等值線產(chǎn)生更大影響，噪聲點去除效果差。結(jié)構(gòu)元素尺寸與格網(wǎng)尺寸效果類似。殘差范圍決定地形的起伏狀況，殘差范圍越大保留的點越多，地形信息就越詳細(xì)；殘差范圍越小，點云刪除率越大，地形越平坦。試驗對研究區(qū)內(nèi)兩期濾波后獲取的地面點云進(jìn)行改進(jìn)數(shù)學(xué)形態(tài)學(xué)處理，獲取地面信息。參考相關(guān)文獻(xiàn)，經(jīng)多次嘗試，最終選取6組效果較好的參數(shù)進(jìn)行數(shù)據(jù)分析，對比所生成沉陷等值線結(jié)構(gòu)合理性和實測參考點高程的平均絕對誤差（MAE）和均方根誤差（RMSE）。統(tǒng)計結(jié)果見表1和圖4所示。

采用半盆地觀測的形式，選取下沉盆地右側(cè)作為試驗區(qū)，兩期各有約400萬點云。通過試驗發(fā)現(xiàn)，格網(wǎng)大小對沉陷等值線影響最大，結(jié)構(gòu)元素尺寸的影響最小。圖4（a）為未經(jīng)過改進(jìn)算法處理的沉陷等值線，等值線過于零碎，存在的噪聲數(shù)據(jù)較大，等值線無法滿足沉陷監(jiān)測分析需求。圖4（b）至圖4（g）分別表示表1中序號Ⅰ-Ⅵ組參數(shù)所生成的沉陷等值線。綜合對比，經(jīng)過改進(jìn)形態(tài)學(xué)算法提取出地表信息后獲取的等值線均較初始情況有所改善。其中第Ⅰ、Ⅱ、Ⅲ組參數(shù)所生成的沉陷等值線展現(xiàn)的地表信息更為詳細(xì)。但是，其中存在大量破碎多邊形，相比較而言，第Ⅱ組中效果較好。于是在Ⅱ的基礎(chǔ)上嘗試?yán)^續(xù)增大格網(wǎng)或減小殘差范圍。Ⅵ在Ⅱ的基礎(chǔ)上增大格網(wǎng)，發(fā)現(xiàn)生成的等值線走勢和整體形態(tài)有較大改變，由于去除過多地面點，而對地表形態(tài)造成改變，如圖4（g）所示。Ⅴ中減小殘差范圍，降低噪聲點的去除比率，造成等值線更加尖銳，等值線精度過低，如圖4（f）所示。綜合Ⅱ和Ⅵ的試驗結(jié)果，調(diào)整結(jié)構(gòu)尺寸和格網(wǎng)大小，采用Ⅳ參數(shù)進(jìn)行試驗，結(jié)果如圖4（e）所示，在保證精度的同時能夠取得較好的等值線效果。

基于改進(jìn)的數(shù)學(xué)形態(tài)學(xué)算法提取沉陷DEM地形信息，既可以很大程度地保留地形特征，保證下沉盆地的準(zhǔn)確性，又可以盡可能地消除碎屑多邊形。在保證精度的同時能夠取得較好的等值線效果，如圖4（e）所示。格網(wǎng)尺寸為3.5 m，結(jié)構(gòu)元素尺寸為3 m時，沉陷DEM的殘差為0.1 m，平均絕對誤差和均方根誤差分別達(dá)到0.134，0.152 m。用此方法提取出的沉陷等值線在保證精度情況下，符合工程測量與繪圖要求。

3.2 試驗分析

結(jié)合涼水井煤礦工作面局部進(jìn)行試驗分析，采用3.1節(jié)中獲取的最優(yōu)參數(shù)對整個研究區(qū)地面點云數(shù)據(jù)進(jìn)行沉陷等值線提取。未采用改進(jìn)算法的等值線提取結(jié)果如圖5所示，采用文中改進(jìn)算法處理后的等值線效果如圖6所示。

通過對比發(fā)現(xiàn)，改進(jìn)的數(shù)學(xué)形態(tài)學(xué)算法處理后獲取的沉陷等值線，能夠在保留地形特征的同時，盡可能地去除噪聲數(shù)據(jù)的影響，但仍會存在部分等值線碎屑等噪聲數(shù)據(jù)。使用改進(jìn)的數(shù)學(xué)形態(tài)學(xué)算法對復(fù)雜地形提取沉陷地形信息的效果明顯優(yōu)于原始沉陷等值線。

參考《煤礦測量規(guī)程》（2013），地表沉陷盆地以下沉10 mm為下沉邊界。因此，需要對沉陷DEM提取下沉10 mm的邊界線，由于沉陷DEM邊緣區(qū)域下沉誤差具有隨機(jī)性和正、負(fù)抵償特性^[21-23]，搜索出下沉值接近10 mm的柵格，認(rèn)定其中心為地表下沉的邊界，圖7中黃色邊界所示。研究區(qū)2020年11月和2021年5月下沉區(qū)域及等值線變化如圖7所示，采用0.25 m為沉陷等值線間隔，對下沉盆地內(nèi)的下沉區(qū)域進(jìn)行分級，由外向內(nèi)沉陷值分別為0.01，0.30，0.55，0.80，1.05，1.30，1.55，1.80 m。

根據(jù)下沉邊界以及沉陷等值線，繪制沉陷區(qū)域面積統(tǒng)計圖，如圖8所示，圖中統(tǒng)計各沉陷區(qū)間對應(yīng)的沉陷面積。下沉量最大區(qū)域為開采區(qū)域中部，并依次向工作面邊界部分遞減，形成下沉盆地。綜合分析，沉陷盆地形態(tài)，保持著從內(nèi)向外延伸的特征，與工作面回采進(jìn)度相吻合。

除沉陷等值線圖和沉陷面積統(tǒng)計圖外，還可以對沉陷面積進(jìn)行分級統(tǒng)計，如圖9所示，對下沉面積進(jìn)行分級統(tǒng)計后可以看出各個下沉量對應(yīng)下沉面積在兩次航飛期間的變化情況，為相似地形下的開采活動和采后環(huán)境狀況評估提供依據(jù)。

圖10為使用改進(jìn)算法處理與未使用算法生成的三維下沉盆地對比情況。可以發(fā)現(xiàn)在使用改進(jìn)算法后，有效去除噪聲數(shù)據(jù)造成的沉陷盆地局部毛刺，使得下沉盆地更加準(zhǔn)確。

4 結(jié) 論

1）改進(jìn)數(shù)學(xué)形態(tài)學(xué)算法可以消除機(jī)載LiDAR地面點云中噪聲數(shù)據(jù)引起的沉陷等值線鋸齒、毛邊、破碎線段及碎屑多邊形情況，其平均絕對誤差和均方根誤差達(dá)到0.134，0.152 m，在保證其精度的同時，極大改善所提取等值線、下沉面和三維盆地效果。

2）基于改進(jìn)數(shù)學(xué)形態(tài)學(xué)算法提取的沉陷等值線，可進(jìn)行邊界提取、沉陷分布圖繪制、沉陷量和沉陷面積統(tǒng)計等工作，突出機(jī)載LiDAR在面域監(jiān)測上的優(yōu)勢，為開采沉陷巖移參數(shù)獲取和沉陷預(yù)計等工作奠定基礎(chǔ)。

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

[1] 李愛國．基于時序差分干涉測量的黃土溝壑礦區(qū)地表沉降監(jiān)測研究[D].西安：長安大學(xué)，2014．

LI Aiguo．Monitoring of surface subsidence in Loess Gully mining area based on time-series differential interferometry[D]．Xi’an：Chang’an University，2014．

[2]戴華陽．巖層與地表移動變形量的時空關(guān)系及描述方法[J]．煤炭學(xué)報，2018，43（S2）：450-459．

DAI Huayang．Mining subsidence variables and their time-space relationship description［J］．Journal of China Coal Society，2018，43（S2）：450-459．

[3]范立民，馬雄德，李永紅，等．西部高強(qiáng)度采煤區(qū)礦山地質(zhì)災(zāi)害現(xiàn)狀與防控技術(shù)[J]．煤炭學(xué)報，2017，42（2）：276-285．

FAN Limin，MA Xiongde，LI Yonghong，et al．Geological disasters and control technology in high intensity mining area of western China［J］．Journal of China Coal Society，2017，42（2）：276-285．

[4]于海洋，楊禮，張春芳，等．基于LiDAR DEM不確定性分析的礦區(qū)沉陷信息提取[J]．金屬礦山，2017（10）：1-7．

YU Haiyang，YANG Li，ZHANG Chunfang，et al．Mining subsidence information extraction based on LiDAR DEM uncertainty analysis[J]．Metal Mine，2017（10）：1-7．

[5]杜斌，司馬俊杰，王琦，等．基于DEM的礦山開采沉陷分析方法研究[J]．北京測繪，2015（1）：14-17．

DU Bin，SIMA Junjie，WANG Qi，et al．Research on mining subsidence analysis method based on DEM[J]．Beijing Surveying and Mapping，2015（1）：14-17．

[6]于彩霞，董箭，李改肖，等．基于LiDAR點云的高精度DEM構(gòu)建方法[J]．地理空間信息，2015，13（3）：22-24，8．

YU Caixia，DONG Jian，LI Gaixiao，et al．High precision DEM construction method based on LiDAR point cloud[J]．Geospatial Information，2015，13（3）：22-24，8．

[7]曹建軍，那嘉明，湯國安，等．DEM隨機(jī)噪聲誤差去除方法研究——以黃土溝沿線地形剖面為例[J]．地理與地理信息科學(xué)，2017，33（2）：12-18，127．CAO Jianjun，NA Jiaming，TANG Guoan，et al．Study on DEM random noise error Removal method：A case study of topographic profile along Loess Gully[J]．Geography and Geo-Information Science，2017，33（2）：12-18，127．

[8]HUI Z，JIN S，XIA Y，et al．A mean shift segmentation morphological filter for airborne LiDAR DTM extraction under forest canopy[J]．Optics amp; Laser Technology，2021，136（1）：106728．

[9]袁養(yǎng)林．基于機(jī)載LiDAR點云數(shù)據(jù)生成DEM和等高線研究[D]．西安：長安大學(xué)，2020．YUAN Yanglin．Study on DEM and contour generation based on airborne LiDAR point cloud data[D]．Xi’an：Chang’an University，2020．

[10]馮梅，鐘斌．基于LiDAR點云自動生成等高線的方法研究[J]．測繪與空間地理信息，2012，35（6）：87-90，93．

FENG Mei，ZHONG Bin．Automatic contour generation method based on LiDAR point cloud[J]．Geomatics and Spatial Information Technology，2012，35（6）：87-90，93．

[11]盧遙，余濤，盧小平，等．基于高差分析的點云數(shù)據(jù)提取礦區(qū)地表沉陷信息方法[J]．測繪通報，2013（3）：22-25．

LU Yao，YU Tao，LU Xiaoping，et al．Method of extracting surface subsidence information from point cloud data based on elevation difference analysis[J].Bulletin of Surveying and Mapping，2013（3）：22-25．

[12]STOCKDONF H F，SALLENGER J A H，LIST J H，et al．Estimation of shoreline position and change using airborne topographic lidar data[J]．Journal of Coastal Research，2002，18（3）：502-513．

[13]ROBERTSON W，WHITMAN D，ZHANG K，et al．Mapping shoreline position using airborne laser altimetry[J]．Journal of Coastal Research，2004，20（3）：884-892．

[14]于彩霞，黃文騫，鄭義東，等．一種基于ArcGIS的LiDAR點云數(shù)據(jù)生成等值線的方法[J]．海洋測繪，2015，35（2）：52-54．YU Caixia，HUANG Wenqian，ZHENG Yidong，et al．A method for generating contour lines from LiDAR point cloud data based on ArcGIS[J]．Marine Surveying and Mapping，2015，35（2）：52-54．

[15]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 Sensing，2022，14（2）：421．

[16]牛峰明．機(jī)載LiDAR礦區(qū)沉陷信息提取方法研究[D]．鄭州：河南理工大學(xué)，2016．NIU Fengming．Research on subsidence information extraction method of airborne LiDAR Mining area[D]．Zhengzhou：Henan Polytechnic University，2016．

[17]何睦，李夢珍．基于LiDAR點云數(shù)據(jù)的等值線繪制方法[J]．計算機(jī)系統(tǒng)應(yīng)用，2018，27（6）：158-164．

HE Mu，LI Mengzhen．Contour rendering method based on LiDAR point cloud data[J]．Application of Computer Systems，2018，27（6）：158-164．

[18]ZHANG K，CHEN S C，WHITMAN D，et al．A progressive morphological filter for removing nonground measurements from airborne LIDAR data[J]．IEEE Transactions on Geoscience and Remote Sensing，2003，41（4）：872-882．

[19]文福林，張凱，蒲鋒，等．一種基于噪聲點鄰域的形態(tài)學(xué)濾波算法研究[J]．電子技術(shù)應(yīng)用，2021，47（3）：71-74．

WEN Fulin，ZHANG Kai，PU Feng，et al．Research on morphological filtering algorithm based on noise point neighborhood[J]．Application of Electronic Technique，2021，47（3）：71-74．

[20]周欽坤，岳建平，李子寬，等．一種改進(jìn)的形態(tài)學(xué)LiDAR點云濾波算法[J]．甘肅科學(xué)學(xué)報，2020，32（5）：1-5．ZHOU Qinkun，YUE Jianping，LI Zikuan，et al．An improved morphological LiDAR point cloud filtering algorithm[J]．Journal of Gansu Science，2020，32（5）：1-5．

[21]湯伏全，蘆家欣，韋書平，等．基于無人機(jī)LiDAR的榆神礦區(qū)采煤沉陷建模方法改進(jìn)[J]．煤炭學(xué)報，2020，45（7）：2655-2666．

TANG Fuquan，LU Jiaxin，WEI Shuping，et al．Improvement of mining subsidence modeling method based on UAV LiDAR in Yushen Mining area[J]．Journal of China Coal Society，2020，45（7）：2655-2666．

[22]何柯璐，湯伏全，韋書平，等．機(jī)載激光點云構(gòu)建西部礦區(qū)開采沉陷模型研究[J]．測繪科學(xué)，2021，46（2）：130-138．

HE Kelu，TANG Fuquan，WEI Shuping，et al．Research on mining subsidence model of western mining area constructed by airborne laser point cloud[J]．Science of Surveying and Mapping，2021，46（2）：130-138．

[23]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．

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