姚頑強(qiáng) ，蒙延斌 ，鄭俊良 ，薛志強(qiáng)

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

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

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

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

1 研究區(qū)概況

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

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

實(shí)測(cè)觀測(cè)線整體布設(shè)圖如圖2。

圖2 實(shí)測(cè)觀測(cè)線布設(shè)圖Fig.2 Layout of actual observation lines

研究區(qū)工作面走向?yàn)槌浞植蓜?dòng)，走向觀測(cè)線布設(shè)在工作面的中心線上。綜合考慮工作面地質(zhì)結(jié)構(gòu)情況，工作面兩端走向線覆蓋整個(gè)工作面，并向兩端外側(cè)各延伸150 m，布設(shè)了走向觀測(cè)線、傾向觀測(cè)線以及道路觀測(cè)線，其中紅色表示地面觀測(cè)點(diǎn)。

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

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

2.1 沉陷DEM 獲取

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

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

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

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

2.2 誤差分析

通過對(duì)原始點(diǎn)云的預(yù)處理，利用點(diǎn)云插值和濾波生成的DEM 仍然存在誤差，主要來源有3 個(gè)方面：

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

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

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

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

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

3 沉陷DEM 濾波及實(shí)驗(yàn)分析

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

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

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

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

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

3.2 濾波方法對(duì)比

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

PSNR 值是1 種圖像客觀的評(píng)價(jià)指標(biāo)，基于誤差敏感的圖像質(zhì)量評(píng)價(jià)，其PSNR 的計(jì)算公式為：

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

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

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

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

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

3.3 與實(shí)測(cè)下沉量對(duì)比

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

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

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

3.4 去噪后誤差分析和精度驗(yàn)證

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

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

因?yàn)楦咚乖肼暤姆秶鷷?huì)遍布到所有的灰度級(jí)，以及其他噪聲分布也遵循一定的規(guī)律，所以為了達(dá)到更好地去除噪聲效果，針對(duì)不同噪聲采用不用的濾波器。相比而言，維納濾波+中值濾波的對(duì)于高斯噪聲的效果最好，維納濾波對(duì)于服從高斯分布的突變誤差更好，而其他幾種方法不但消除了突變誤差，還消除了邊緣誤差，非地面點(diǎn)誤差，以及非開采沉陷因素所引起的誤差。使得沉陷DEM的邊緣信息更好，沉陷模型的精度相對(duì)于4 種獨(dú)立去噪方法的精度提高了很多，去噪效果的好壞的順序?yàn)榫S納濾波+中值濾波＞維納濾波+中值濾波+高斯濾波＞均值濾波+維納濾波＞中值濾波+均值濾波。

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

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

表2 多重濾波去噪誤差指標(biāo)Table 2 Multi-filter denoising error indicators

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

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

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

4 結(jié)語

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