鄭俊良,姚頑強,藺小虎,張聯(lián)隊,張 詠,張 冬

(1.西安科技大學(xué) 測繪科學(xué)與技術(shù)學(xué)院,陜西 西安 710054;2.陜西彬長礦業(yè)集團(tuán)有限公司,陜西 咸陽 712000;3.陜西延長石油靖邊煤業(yè)有限公司,陜西 榆林 718500)

0 引 言

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

目前,開采沉陷監(jiān)測方法包括傳統(tǒng)預(yù)埋監(jiān)測樁觀測(水準(zhǔn)測量、全站儀測量、GNSS測量等)、InSAR監(jiān)測、地面三維激光測量、傾斜攝影測量和機載LiDAR測量。傳統(tǒng)預(yù)埋監(jiān)測樁觀測方法局限于點位觀測,具有較高的點位精度(毫米級、厘米級),但不能準(zhǔn)確描述采空區(qū)地面整體的形變過程,且西部溝壑區(qū)地形結(jié)構(gòu)復(fù)雜,存在工作量大和人員無法到達(dá)的問題[6];InSAR測量能夠進(jìn)行面狀觀測,且取得較高的測量精度(毫米級),但受限于觀測周期,對于地表塌陷速度快、沉陷量大的區(qū)域,容易造成失相干現(xiàn)象[7-8];地面三維激光測量能夠取得面狀厘米級,甚至毫米級觀測精度,但由于其視野范圍、點云拼接誤差等原因,在超長超寬工作面監(jiān)測時存在工作量大、累積誤差大等問題[9];近年來,隨著影像密集特征點匹配算法的成熟,無人機傾斜攝影測量技術(shù)提高了攝影測量的觀測精度,但在地表植被覆蓋度高、地形起伏大的區(qū)域,易造成監(jiān)測結(jié)果精度低,且精度受限于像控點密度[10-11];機載LiDAR測量具有精度高、觀測范圍大的優(yōu)勢,但有人機載LiDAR測量存在設(shè)備昂貴、天氣要求高、空域申請困難等問題;隨著無人機技術(shù)發(fā)展和激光掃描設(shè)備改進(jìn),無人機LiDAR測量技術(shù)逐漸成熟,在滿足厘米級測量精度的同時能夠極大地降低成本、提高效率,彌補有人機載LiDAR測量的不足[12-13]。無人機LiDAR通過對目標(biāo)區(qū)域低空多時相掃描,可獲取高密度點云數(shù)據(jù),經(jīng)點云濾波、插值和求差計算,可獲取目標(biāo)區(qū)域高程變化。但在激光點云濾波、插值和求差計算過程中不可避免地存在分類錯誤、噪聲和誤差,一些學(xué)者通過改進(jìn)濾波算法和插值方法提高DEM精度,進(jìn)而改善沉陷模型精度[14-17],但仍存在一定問題：①現(xiàn)有成熟點云濾波算法主要針對傳統(tǒng)有人機載LiDAR點云數(shù)據(jù),相比較而言,無人機LiDAR點云具有更高的點云密度和更復(fù)雜的噪聲數(shù)據(jù),如在山區(qū)應(yīng)用效果較好的漸進(jìn)三角網(wǎng)點云濾波算法(Progressire Trianglated Irregular Network,PTIN),在無人機LiDAR點云濾波時,由于高密度和多噪聲的原因,易出現(xiàn)種子點錯誤問題,導(dǎo)致濾波錯誤;②無人機LiDAR獲取的點云、數(shù)字高程模型(Digital Elevation Model,DEM)和沉陷模型之間的精度缺乏實測數(shù)據(jù)驗證;在無人機LiDAR測量系統(tǒng)和環(huán)境條件相同情況下,沉陷模型精度是否優(yōu)于原始點云和DEM精度缺乏試驗對比;③傳統(tǒng)LiDAR監(jiān)測成果單一,對沉陷模型進(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,地勢西北高東南低。以涼水井煤礦431301工作面為研究對象,開采煤層為4-3煤層,該煤層平均埋深138 m,平均厚度1.3 m,采用長臂式綜采、頂板自然塌陷方法開采;其中上覆4-2煤層已于2013年開采完成,地表于2014年達(dá)到穩(wěn)定狀態(tài)。

圖1 研究對象位置Fig.1 Location of study object

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

圖2 研究區(qū)與監(jiān)測線布設(shè)Fig.2 Layout of study area and monitoring line

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

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

圖3 研究區(qū)與采空區(qū)位置關(guān)系Fig.3 Location relationship between study area and goaf

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

1.3 研究技術(shù)路線

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

圖4 研究技術(shù)路線Fig.4 Research technology route

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

(1)

(2)

(3)

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

2 LiDAR點云濾波與插值計算

2.1 常用濾波方法

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

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

2.2 PTIN種子點選取及濾波

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

圖5 改進(jìn)的PTIN濾波流程Fig.5 Improved PTIN filter flow chart

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

圖6 種子點選取計算示意Fig.6 Schematic diagram of seed point selection calculation

圖7 局部點云濾波效果Fig.7 Partial point cloud filtering results

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

2.3 點云插值計算

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

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

圖8 DEM誤差分布Fig.8 DEM error distribution

同時,計算不同格網(wǎng)尺寸DEM的ME,MAE,RMSE,結(jié)果見表1。

表1 DEM誤差Table 1 DEM errors

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

2.4 DEM差值計算

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

圖9 沉陷模型三維效果Fig.9 3D rendering of subsidence model

3 精度評定

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

表2 點云誤差統(tǒng)計Table 2 Point cloud error statistics

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

圖10 點云、DEM,沉陷模型誤差分布Fig.10 Errors distribution of cloud,DEM,subsidence model

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

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

圖11 下沉折線Fig.11 Subsidence lines

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

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

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

4.1 下沉量統(tǒng)計

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

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

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

表3 下沉面積統(tǒng)計Table 3 Statistics of subsidence areas

4.2 下沉剖面信息提取

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

圖13 虛擬監(jiān)測點與剖面線位置Fig.13 Location map of virtual monitoring points and profile lines

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

圖14 虛擬監(jiān)測點下沉折線Fig.14 Virtual monitoring point subsidence line

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

圖15 剖面下沉曲線Fig.15 Profile subsidence curves

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

4.3 沉陷面域分析

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

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

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

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

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

為了明晰煤層開采造成地表下沉的范圍和嚴(yán)重程度,以不同下沉量所占下沉面積為參數(shù),進(jìn)行地表損害的分等定級,將地表下沉值進(jìn)行分組,并確定其范圍和對應(yīng)面積,結(jié)果如圖18所示,可以看出,不同沉降等級區(qū)域由內(nèi)向外逐漸擴展,采空區(qū)從中心向外延伸。沉降可分為6級：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,能為工作面開采過程中地表建(構(gòu))筑物的保護(hù)、移民搬遷提供依據(jù)。

圖18 地表下沉量分級Fig.18 Classification of surface subsidence

5 結(jié) 論

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

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

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