無人機LiDAR在采空區沉陷監測中的應用

鄭俊良,姚頑強,藺小虎,張聯隊,張 詠,張 冬

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

0 引 言

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

目前,開采沉陷監測方法包括傳統預埋監測樁觀測(水準測量、全站儀測量、GNSS測量等)、InSAR監測、地面三維激光測量、傾斜攝影測量和機載LiDAR測量。傳統預埋監測樁觀測方法局限于點位觀測,具有較高的點位精度(毫米級、厘米級),但不能準確描述采空區地面整體的形變過程,且西部溝壑區地形結構復雜,存在工作量大和人員無法到達的問題[6];InSAR測量能夠進行面狀觀測,且取得較高的測量精度(毫米級),但受限于觀測周期,對于地表塌陷速度快、沉陷量大的區域,容易造成失相干現象[7-8];地面三維激光測量能夠取得面狀厘米級,甚至毫米級觀測精度,但由于其視野范圍、點云拼接誤差等原因,在超長超寬工作面監測時存在工作量大、累積誤差大等問題[9];近年來,隨著影像密集特征點匹配算法的成熟,無人機傾斜攝影測量技術提高了攝影測量的觀測精度,但在地表植被覆蓋度高、地形起伏大的區域,易造成監測結果精度低,且精度受限于像控點密度[10-11];機載LiDAR測量具有精度高、觀測范圍大的優勢,但有人機載LiDAR測量存在設備昂貴、天氣要求高、空域申請困難等問題;隨著無人機技術發展和激光掃描設備改進,無人機LiDAR測量技術逐漸成熟,在滿足厘米級測量精度的同時能夠極大地降低成本、提高效率,彌補有人機載LiDAR測量的不足[12-13]。無人機LiDAR通過對目標區域低空多時相掃描,可獲取高密度點云數據,經點云濾波、插值和求差計算,可獲取目標區域高程變化。但在激光點云濾波、插值和求差計算過程中不可避免地存在分類錯誤、噪聲和誤差,一些學者通過改進濾波算法和插值方法提高DEM精度,進而改善沉陷模型精度[14-17],但仍存在一定問題：①現有成熟點云濾波算法主要針對傳統有人機載LiDAR點云數據,相比較而言,無人機LiDAR點云具有更高的點云密度和更復雜的噪聲數據,如在山區應用效果較好的漸進三角網點云濾波算法(Progressire Trianglated Irregular Network,PTIN),在無人機LiDAR點云濾波時,由于高密度和多噪聲的原因,易出現種子點錯誤問題,導致濾波錯誤;②無人機LiDAR獲取的點云、數字高程模型(Digital Elevation Model,DEM)和沉陷模型之間的精度缺乏實測數據驗證;在無人機LiDAR測量系統和環境條件相同情況下,沉陷模型精度是否優于原始點云和DEM精度缺乏試驗對比;③傳統LiDAR監測成果單一,對沉陷模型進一步研究較少,未進行深入的數據分析和挖掘。

1 研究區及數據獲取

1.1 研究區概況

涼水井礦位于陜西省榆林市東北部、榆神府礦區東部,處于北緯38°47′～38°53′和東經110°14′～110°21′之間(圖1)。井田地面標高1 100～1 326 m,地勢西北高東南低。以涼水井煤礦431301工作面為研究對象,開采煤層為4-3煤層,該煤層平均埋深138 m,平均厚度1.3 m,采用長臂式綜采、頂板自然塌陷方法開采;其中上覆4-2煤層已于2013年開采完成,地表于2014年達到穩定狀態。

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

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

圖2 研究區與監測線布設Fig.2 Layout of study area and monitoring line

1.2 數據獲取

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

圖3 研究區與采空區位置關系Fig.3 Location relationship between study area and goaf

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

1.3 研究技術路線

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

圖4 研究技術路線Fig.4 Research technology route

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

(1)

(2)

(3)

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

2 LiDAR點云濾波與插值計算

2.1 常用濾波方法

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

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

2.2 PTIN種子點選取及濾波

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

圖5 改進的PTIN濾波流程Fig.5 Improved PTIN filter flow chart

格網劃分主要是用規則的虛擬格網覆蓋全部點云數據,將點云按照平面坐標分配到相應的虛擬網格內。種子點的選取是濾波的關鍵,采用改進的擴展局部最小值方法進行種子點選取,原理如圖6所示。

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

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

對2020年11月和2021年5月點云數據分別進行點云濾波,并對濾波結果進行統計,地面點云占比分別為93%和83%,由于不同時期地表植被覆蓋度不同,導致了地面點云占比不同。

2.3 點云插值計算

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

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

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

同時,計算不同格網尺寸DEM的ME,MAE,RMSE,結果見表1。

表1 DEM誤差Table 1 DEM errors

由表1可知,ME,MAE和RMSE在0.05～1 m區間內變化不大,當格網尺寸大于1 m時,誤差急劇增大。同時,結合圖8綜合考慮,最終選擇0.1 m的格網尺寸進行插值計算,獲取較高精度的研究區DEM。

2.4 DEM差值計算

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

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

3 精度評定

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

表2 點云誤差統計Table 2 Point cloud error statistics

其次,基于監測樁實測數據對DEM和沉陷模型誤差進行計算,誤差分布如圖10所示。

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

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

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

圖11 下沉折線Fig.11 Subsidence lines

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

4 基于沉陷模型的數據挖掘與分析

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

4.1 下沉量統計

為了準確獲取地表不同位置的下沉量和不同下沉量的分布情況,對沉陷模型進行下沉量分級分類統計,結果如圖12所示。

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

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

表3 下沉面積統計Table 3 Statistics of subsidence areas

4.2 下沉剖面信息提取

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

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

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

圖14 虛擬監測點下沉折線Fig.14 Virtual monitoring point subsidence line

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

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

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

4.3 沉陷面域分析

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

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

從圖16可以看出,不同等值距的下沉等值線所表達的特征詳細度不同。較小的等值線區間值對應較詳細的沉陷特征,同時也易受噪聲和異常值的影響。通過綜合比較發現,沉降間隔為100 mm的等值線既可以保留采煤沉陷特征細節,也能夠降低噪聲和異常值的影響。

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

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

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

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

5 結 論

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

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

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