礦區(qū)沉陷DEM 多重濾波方法研究

姚頑強 ，蒙延斌 ，鄭俊良 ，薛志強

（1.西安科技大學 測繪科學與技術(shù)學院，陜西 西安 710054；2.陜西彬長孟村礦業(yè)有限公司，陜西 咸陽 713602）

煤炭資源在我國的能源中占據(jù)著很重要的位置，有效地推動了國民經(jīng)濟的快速發(fā)展[1-2]。但是隨著煤炭開采強度的提高導致地表大范圍塌陷，使地表產(chǎn)生移動和變形[3-5]，給礦區(qū)環(huán)境帶來嚴重的影響[6]。在實地監(jiān)測中，現(xiàn)有的水準測量、GPSRTK 測量及InSAR 等遙感監(jiān)測方法均具有一定的局限性。前者不僅工作量大[7]，而且僅通過計算各項地表巖移參數(shù)與角量信息很難獲取到整個移動盆地的形變信息；而后者InSAR 作為1 種全新的對地觀測技術(shù)，雖然在礦區(qū)開采沉陷得到廣泛應用[8-10]，但是只能測取地表的微小形變信息，存在大梯度變形區(qū)域失相干的缺點，使得InSAR 技術(shù)在沉陷監(jiān)測中的應用受到了限制。

近些年來，無人機遙感技術(shù)發(fā)展迅速，在礦區(qū)監(jiān)測中取得了一定成效，并在多方面得到應用[11-15]。無人機LiDAR 作為1 種激光遙感技術(shù)，具有獲取數(shù)據(jù)速度快，飛行時受地形因素小的特點，可以獲取到礦區(qū)高精度、高分辨率的三維點云數(shù)據(jù)，并進行點云濾波和插值生成DEM，并將多期DEM疊加得到沉陷DEM。但是基于點云濾波[16-18]的局限性，仍存在噪聲，從而使沉陷DEM 也會包含噪聲和誤差。有學者通過不同的去噪算法來提高沉陷模型的精度，但是目前還是存在問題：①無人機所獲取的數(shù)據(jù)精度不能確定；②利用現(xiàn)在主流的點云濾波和DEM 插值算法計算結(jié)果難以驗證，影響后續(xù)DEM 的建模精度[19-20]；③沒有根據(jù)整個移動盆地的誤差分布特征進行算法的選擇，去噪效果較差，不能滿足沉陷監(jiān)測的需求。

綜上，以涼水井煤礦某工作面為研究區(qū)域，研究沉陷DEM 的去噪問題；利用無人機采集的2期點云數(shù)據(jù)進行濾波插值生成DEM，疊加后得到沉陷DEM；根據(jù)沉陷盆地的誤差分布特征選擇多重濾波組合去除噪聲，并通過實測數(shù)據(jù)對比去噪效果，得到最優(yōu)組合去噪方案，為沉陷監(jiān)測提供了好的技術(shù)支撐。

1 研究區(qū)概況

研究區(qū)位于陜西省神木市的涼水井煤礦，在陜北黃土高原北部，毛烏素沙漠的南邊，屬于丘陵區(qū)。地理坐標為：東經(jīng)111°14′22″～111°21′24″，北緯38°47′29″～38°53′24″。工作面內(nèi)沒有河流存在，地貌主要是以風積沙地貌為主，大部分區(qū)域為風積沙層所覆蓋。工作面地表的其他區(qū)域無村莊、耕地、建筑物分布。地勢由東北向東南方向逐漸降低，地表整體西部要高于東部，傾角小于1°。431303 工作面位于井田的最南端，是431 盤區(qū)的首個回采工作面，工作面為4-3煤層，4-2煤層是上覆巖層，在前期已開采完畢。煤層相對位置情況如圖1，上部為盤區(qū)的南部采空積水區(qū)，西部為煤礦輔運大巷，北部和東部主要是實煤區(qū)。礦區(qū)工作面的回采區(qū)域內(nèi)煤層屬于較穩(wěn)定煤層，煤層平均采深137.92 m。

圖1 涼水井礦雙煤層相對位置示意圖Fig.1 Diagram of the relative position of the double coal seam

實測觀測線整體布設圖如圖2。

圖2 實測觀測線布設圖Fig.2 Layout of actual observation lines

研究區(qū)工作面走向為超充分采動，走向觀測線布設在工作面的中心線上。綜合考慮工作面地質(zhì)結(jié)構(gòu)情況，工作面兩端走向線覆蓋整個工作面，并向兩端外側(cè)各延伸150 m，布設了走向觀測線、傾向觀測線以及道路觀測線，其中紅色表示地面觀測點。

實驗所采用的機載LiDAR 系統(tǒng)型號為SZT-250，掃描的測程可以達到250 m，具有高精度的導航定位系統(tǒng)和激光掃描系統(tǒng)。數(shù)據(jù)采集時通過移動測量系統(tǒng)操控軟件來獲取點云的姿態(tài)和導航定位信息。激光掃描儀的測距精度為15 mm，視場角為360°，掃描速度在10～100 線/s，在進行多目標探測時，每束激光可以達到5 次回波，每期無人機航飛都采集出2 組數(shù)據(jù)，這樣可以保證所采集的點云具有足夠的密度，點云的質(zhì)量也可以得到很好的保證。所掃描出的每2 期點云數(shù)據(jù)的平均密度為60～700 脈沖/m2，經(jīng)過檢驗后，點云的質(zhì)量達到標準，可進行點云數(shù)據(jù)的處理工作。

2 數(shù)據(jù)獲取與誤差分析

2.1 沉陷DEM 獲取

在本次數(shù)據(jù)采集中,無人機LiDAR 對涼水井礦區(qū)工作面分別在2020 年8 月15 日（回采位置1 024.6 m）、2020 年11 月7 日（回采位置1 384.6 m）和2021 年5 月19 日（回采完畢）一共進行了3 次觀測,無人機飛行時工作面開采位置如圖3。

圖3 無人機飛行時工作面開采位置Fig.3 Working face mining position at the time of the drone flight

點云處理及成果分析流程如圖4。通過對期不同的點云數(shù)據(jù)進行濾波、插值,然后利用第1 期和第3 期的DEM 數(shù)據(jù)進行相減，得到累積的沉陷DEM。并視為這個DEM 為原始的沉陷DEM，由于利用現(xiàn)有的濾波算法不能準確的分離出地面點和非地面點，并且對點云進行插值的時候,還會存在插值誤差，所以使得初始沉陷模型也會包含誤差，對后期的實驗處理結(jié)果會有很大的影響。將沉陷DEM 利用ArcScene 進行場景屬性設置，基于范圍進行計算，將沉陷DEM 高程方向拉伸40 倍，以突出沉陷模型的誤差特征。沉陷DEM三維模型如圖5。

圖5 去噪前沉陷DEM 三維視圖Fig.5 3D view of DEM before noise removal

2.2 誤差分析

通過對原始點云的預處理，利用點云插值和濾波生成的DEM 仍然存在誤差，主要來源有3 個方面：

1）濾波算法時未除去的非地面點誤差。對點云進行濾波時，不同的濾波算法都會有一定的限制性，不能完全剔除一些誤差對于模型的影響。在不同的季節(jié)，植被的生長，從而使一些植被點不重疊，導致沉陷DEM 會出現(xiàn)噪聲。

2）不同點云插值生成DEM 引起的模型誤差。每一期所采集的點云數(shù)據(jù)都會有所差異，例如點云的分布位置，點云的厚度和寬度，在內(nèi)插生成DEM 時,會對精度造成一定的損失。采集點云數(shù)據(jù)的區(qū)域情況，當?shù)乜赡軙霈F(xiàn)地形起伏,坡度和坡向變化的情況。這樣會使得沉陷DEM 有很大的模型誤差,要想減小這種誤差,可以利用最優(yōu)的點云插值算法或者去噪方法,同時也可以適當增加點云的寬度和厚度，讓點云的分布更均勻,方便數(shù)據(jù)的處理。

3）無人機LiDAR 外業(yè)掃描引起的誤差。無人機在外業(yè)掃描時，會受到外界環(huán)境和人為因素的影響。例如無人機的飛行姿態(tài)時刻都在變化，以及無人機在掃描地面時會和地面有1 個傾角范圍，角度過大或過小會陷入盲區(qū)，導致有些地物或者植被掃描不到，點云也會有缺失和分布稀疏的情況。不同的航線之間采集數(shù)據(jù)具有重疊度，在沒有誤差的情況下,2 個航帶所采集的數(shù)據(jù)是相互交融到一起的。這樣在進行兩期DEM 影像相減的時候也會出現(xiàn)誤差。

然而減小這種誤差的途徑有很多，操作人員可以掌控無人機和地面的飛行高度，通過增加掃描范圍獲取更多點云的數(shù)量，利用TerraSolid 軟件可以將范圍線導入點云中，看范圍線是否被包含在內(nèi)，檢查點云的質(zhì)量。可以利用PCL 點云庫對點云進行配準，根據(jù)源點云點數(shù)，目標點云點數(shù)進行重采樣，最后輸出變換矩陣，得到最終的配準點云數(shù)量，減小系統(tǒng)性誤差和產(chǎn)生的一些隨機性誤差。

根據(jù)所分析的幾種誤差可知,不同的誤差類型都可以影響沉陷DEM 的精度,所以為了后續(xù)不再對沉陷模型產(chǎn)生誤差，直接對沉陷DEM 進行去噪。

3 沉陷DEM 濾波及實驗分析

3.1 經(jīng)典濾波和多重濾波原理

在圖像中，由于成像系統(tǒng)、傳輸介質(zhì)的不完善，圖像在形成過程中往往會受到多種噪聲的污染，擾亂圖像的可觀測的信息。沉陷DEM 與圖像類似，采用圖像濾波方法可實現(xiàn)沉陷DEM 中噪聲的去除，細節(jié)特征更加明顯。沉陷DEM 去噪流程圖如圖6。

圖6 沉陷DEM 去噪流程Fig.6 Subsidence DEM denoising process

經(jīng)典濾波目前常用的方法有高斯濾波、維納濾波、中值濾波、均值濾波。4 種濾波都有不同的特性。均值濾波對高斯噪聲有良好的去噪能力，但是在消除噪聲的同時也會對圖像的高頻細節(jié)成分造成破壞和損失。高斯濾波根據(jù)圖像的邊緣特征，既含有低頻分量又含有高頻分量，平滑時不會被不需要的高頻信號污染，保留了大部分信號。維納濾波對高斯噪聲有明顯抑制作用，但是容易失去邊緣信息。中值濾波讓圖像周圍的像素值接近真實值，消耗孤立的噪聲點，能保持良好的邊緣特性，但是會失真過多。

多重濾波是根據(jù)誤差的分布情況和算法特點選擇不同的組合算法，以及根據(jù)沉陷DEM 的誤差分布特性，并結(jié)合地表沉陷盆地本身的特征分布進行算法的選擇疊加，本研究主要是選擇了中值濾波+維納濾波，維納濾波+中值濾波+高斯濾波，均值濾波+維納濾波，中值濾波+均值濾波這4 種組合方式，并對算法進行疊加，對地表沉陷未除去的噪聲誤差繼續(xù)實驗驗證。

3.2 濾波方法對比

在MATLAB 的工作環(huán)境下,為了驗證幾種濾波去除噪聲效果的好壞，從主觀和客觀的角度綜合進行分析評價。分別使用中值濾波、維納濾波、均值濾波、高斯濾波進行濾波處理。對比分析4種濾波的去噪效果；得到的中值濾波算法、維納濾波算法、均值濾波算法、高斯濾波算法的PSNR 值分別為27.526 8、28.439 1、28.614 2、26.992 7。

PSNR 值是1 種圖像客觀的評價指標，基于誤差敏感的圖像質(zhì)量評價，其PSNR 的計算公式為：

式中：PSNR 為圖像的峰值信噪比；MSE 為原圖像與處理圖像之間均方誤差。

式中：I、K分別為原始影像和處理后影像；i、j為像素值的行和列；m、n為表示大小的尺寸值。

在幾種濾波方法中，對于含有噪聲的沉陷DEM，PSNR 值有所不同，對于突變噪聲，高斯和維納濾波可以取得更好地去除效果，但是各種濾波算法對噪聲傷害都很大，在盡量保留圖像細節(jié)特征的條件下對噪聲進行去除。而每種濾波特性有所不同，處理不好會使圖像質(zhì)量下降，失真，PSNR 值越大，表示圖像的失真越少。

通過4 組沉陷DEM 去噪實驗，得到的去噪后沉陷DEM 三維視圖如圖7。由圖7 可以看出：去噪后的沉陷DEM 相對于初始的沉陷DEM 去噪效果更好,并且在沉陷DEM 中的非沉陷區(qū)域中服從高斯分布的突變噪聲信號被去除了,保留了下沉盆地的一些真實的細節(jié)特征。也在一定程度上反映出了礦區(qū)沉陷模型的真實情況。

圖7 去噪后沉陷DEM 三維視圖Fig.7 3D view of DEM subsidence after noise removal

3.3 與實測下沉量對比

在開采沉陷期間,研究區(qū)可采煤層包括4-2煤層和4-3煤層，工作面實際煤層采高與觀測下沉結(jié)果之間不符合開采沉陷的客觀規(guī)律，是由于研究區(qū)多煤層重復采動異常造成的。4-2煤層雖然在開采完畢后經(jīng)歷長時間塌陷，但仍可能存在一定的裂隙和空洞。上覆巖層在短期內(nèi)很難達到穩(wěn)定狀態(tài)。4-3煤層的開采對4-2煤層上覆巖層造成進一步擾動，破壞進一步加劇，下沉量更大。為了更直觀地反映本研究不同方法進行沉陷DEM 去噪的效果，基于處理前、后的沉陷DEM 沿下沉盆地走向線提取下沉數(shù)據(jù)，以橫坐標表示控制點點號,縱坐標表示下沉量，繪制下沉量曲線圖。走向觀測線不同濾波去噪下沉情況如圖8。

圖8 走向觀測線不同濾波去噪下沉情況Fig.8 De-noising subsidence of strike observation line with different filters

在研究區(qū)選取54 個觀測點，通過與實測數(shù)據(jù)對比可以看出：從控制點A173 到A182 的曲線走勢情況來看,在這段區(qū)間地表有明顯的下沉,形成1 個下沉盆地。綜上，沉陷模型在去噪前的誤差基礎上，地表下沉盆地在邊緣的部分相對于來說較平緩，在礦區(qū)工作面的中心區(qū)域下沉具有明顯的突變特征?？傮w來預判，高斯濾波和維納濾波的去噪效果相對較好。

3.4 去噪后誤差分析和精度驗證

通過上述下沉量曲線可以看出，高斯濾波和維納濾波有著很好地去噪效果。這說明沉陷DEM中的一些誤差服從高斯分布，去除了沉陷邊緣的突變噪聲。而均值濾波和中值濾波的效果卻一般。分析得出，根據(jù)沉陷DEM 的誤差特性，后兩者的去噪范圍并沒有將此誤差包含在內(nèi)，不符合地表沉陷的分布特征，從而使得去除噪聲效果較差。所以需要對4 種濾波算法組合去除噪聲，將沉陷模型精度進一步提高。濾波后的直方圖如圖9。

圖9 多重濾波直方圖Fig.9 Multi-filter histograms

因為高斯噪聲的范圍會遍布到所有的灰度級，以及其他噪聲分布也遵循一定的規(guī)律，所以為了達到更好地去除噪聲效果，針對不同噪聲采用不用的濾波器。相比而言，維納濾波+中值濾波的對于高斯噪聲的效果最好，維納濾波對于服從高斯分布的突變誤差更好，而其他幾種方法不但消除了突變誤差，還消除了邊緣誤差，非地面點誤差，以及非開采沉陷因素所引起的誤差。使得沉陷DEM的邊緣信息更好，沉陷模型的精度相對于4 種獨立去噪方法的精度提高了很多，去噪效果的好壞的順序為維納濾波+中值濾波＞維納濾波+中值濾波+高斯濾波＞均值濾波+維納濾波＞中值濾波+均值濾波。

對多重濾波的去噪結(jié)果進行驗證與分析，走向觀測線多重濾波下沉情況如圖10，經(jīng)典濾波去噪誤差指標見表1，多重濾波去噪誤差指標見表2，不同濾波算法誤差分布如圖11。

表1 經(jīng)典濾波去噪誤差指標Table 1 Classical filter denoising error indexes

表2 多重濾波去噪誤差指標Table 2 Multi-filter denoising error indicators

圖10 走向觀測線多重濾波下沉情況Fig.10 Multiple filter sinking of strike observation line

圖11 不同濾波算法誤差分布圖Fig.11 Error distribution of different filtering algorithms

實驗表明，根據(jù)沉陷DEM 的誤差分布特征，所組合起來的幾種濾波去噪方法有很大的成效，相比于上述的其他4 種單獨濾波方法不僅去除了沉陷DEM 服從高斯分布的突變誤差，還有非地面點誤差，植被覆蓋度所引起的重疊誤差，以及邊緣誤差，使邊緣輪廓更加清晰。沉陷區(qū)域的下沉量也明顯減少了很多。沉陷DEM 的精度得到了改善，表1 和表2 中的誤差指標也符合要求，均方根誤差和平均絕對誤差都有所減小，也體現(xiàn)出了去噪前后誤差的分布情況，與預想是一致的，并且在上述的去噪實驗中，維納+中值濾波去噪效果是最好的，也證實了本研究對于沉陷DEM 去噪的效果具有較高的可靠性。

4 結(jié)語

利用無人機LiDAR 技術(shù)進行礦區(qū)地表沉陷監(jiān)測時，可以獲取到高精度的點云數(shù)據(jù)，但是目前主流的點云插值與濾波算法生成沉陷DEM 時會有精度不足的現(xiàn)象。其中，非地面點噪聲、點云插值誤差以及外業(yè)掃描所引起的點云密度不足等誤差都會使沉陷DEM 有噪聲出現(xiàn)?；诘V區(qū)沉陷DEM 的建模精度和誤差分布情況，根據(jù)幾種單獨濾波和多重濾波結(jié)果表明，多重濾波不僅去除了服從高斯分布的突變誤差，還去除了非沉陷區(qū)域中的邊緣誤差，非地面點誤差等誤差。其中維納+中值濾波去噪效果最佳，下沉盆地的細節(jié)突出更加明顯，說明基于沉陷DEM 誤差所采取的去噪方法具有一定的合理性和可靠性。但是在本研究中，所得到的實驗結(jié)果是在榆神礦區(qū)得到的，地理環(huán)境相對較復雜，還要其他未知誤差的來源，所以利用無人機LiDAR 對開采沉陷還需要進一步研究。