遲臣鑫 陳偉強(qiáng) 朱鵬程 查劍鋒 殷和健 馬晨陽(yáng)

（1.中國(guó)礦業(yè)大學(xué)環(huán)境與測(cè)繪學(xué)院，江蘇 徐州 221116；2.江蘇省資源環(huán)境信息工程重點(diǎn)實(shí)驗(yàn)室，江蘇 徐州 221116；3.皖北煤電集團(tuán)有限責(zé)任公司五溝煤礦，安徽 淮北 235131）

雨水囤積與地下水匯流會(huì)在采煤塌陷區(qū)形成積水，會(huì)破壞礦區(qū)及附近生態(tài)環(huán)境［1-2］、造成耕地、道路或建筑物的破壞。據(jù)統(tǒng)計(jì)，截至2018年，我國(guó)兩淮地區(qū)的采煤塌陷區(qū)面積約300 km2，積水面積占塌陷區(qū)面積的30%～50%［3］。近年來(lái)，采煤塌陷積水區(qū)的生態(tài)修復(fù)得到了學(xué)術(shù)界的普遍關(guān)注［4-6］，精確獲取面積這一基礎(chǔ)數(shù)據(jù)，對(duì)于有效開(kāi)展采煤塌陷積水區(qū)生態(tài)修復(fù)［7］等工作具有重要意義。

采煤塌陷積水區(qū)面積獲取方法通常有RTK實(shí)測(cè)法和遙感法，后者最為常用。范忻等［8］提出了一種改進(jìn)的P-WSVM方法，用于計(jì)算塌陷積水區(qū)面積，該方法受到遙感數(shù)據(jù)分辨率和波段的影響，提取精度為87.76%；王安妮［9］利用NCIWI指數(shù)法，以Landsat影像為數(shù)據(jù)源提取了2000—2018年的沛北礦區(qū)積水區(qū)面積；胡振琪等［10］采用基于小波變換的多尺度邊緣檢測(cè)算法提取了塌陷積水區(qū)邊界。一般而言，采煤塌陷積水區(qū)具有面積現(xiàn)勢(shì)性強(qiáng)、地形復(fù)雜及測(cè)區(qū)形狀不規(guī)則的特點(diǎn)，傳統(tǒng)方法存在的不足有：①RTK實(shí)測(cè)法適用于小區(qū)域的采煤塌陷積水區(qū)面積計(jì)算，若要保證面積計(jì)算精度，則需提高RTK實(shí)地采點(diǎn)間隔，工作量會(huì)隨之增加；②考慮豐水期、枯水期時(shí)采煤塌陷積水區(qū)面積有差異，并受到遙感影像分辨率的影響［11］，遙感法適用于大范圍的采煤塌陷積水區(qū)的面積計(jì)算，且高分辨率遙感影像獲取成本較高，計(jì)算范圍較小的采煤塌陷積水區(qū)面積會(huì)產(chǎn)生數(shù)據(jù)冗余和經(jīng)費(fèi)浪費(fèi)。為此，本研究提出了一種以無(wú)人機(jī)影像為數(shù)據(jù)源，利用耦合濾波方法消除噪聲、自動(dòng)提取水域輪廓并最終獲取采煤塌陷積水區(qū)面積的方法。

1 塌陷積水區(qū)面積信息無(wú)人機(jī)采集方法

無(wú)人機(jī)采集地表信息的一般流程有航線規(guī)劃、布設(shè)像控點(diǎn)和提取信息。由于塌陷積水區(qū)隨煤炭開(kāi)采逐步變化，其面積在不同季節(jié)會(huì)有差異并且水域邊界不規(guī)則，故本研究提出了如圖1所示的面積采集思路。

具體實(shí)施步驟為：①圍繞采煤塌陷積水區(qū)布設(shè)像控點(diǎn)并用無(wú)人機(jī)拍攝正射影像［12］，保證航向、旁向重疊率在70%以上；②對(duì)獲得的正射影像進(jìn)行預(yù)處理，包括灰度化、去噪和二值化處理等；③提取并處理正射影像上積水區(qū)的輪廓，將輪廓點(diǎn)位坐標(biāo)轉(zhuǎn)換為實(shí)地坐標(biāo)，采用向量叉乘法計(jì)算積水區(qū)面積；④對(duì)采煤塌陷積水區(qū)面積的計(jì)算精度進(jìn)行評(píng)估。

2 塌陷積水區(qū)面積信息提取的關(guān)鍵技術(shù)

2.1 影像數(shù)據(jù)處理

2.1.1 灰度化

灰度化是影像二值化處理的前提，為保留采煤塌陷積水區(qū)影像灰度化后的對(duì)比度，采用文獻(xiàn)［13］提出的基于最大加權(quán)投影求解的灰度化算法進(jìn)行影像灰度化處理。

2.1.2 噪聲處理

常見(jiàn)的圖像噪聲有高斯噪聲和椒鹽噪聲等。高斯噪聲是指概率密度函數(shù)服從高斯分布的噪聲，椒鹽噪聲是指圖像傳感器、傳輸信道及解碼處理在影像獲取、傳輸?shù)冗^(guò)程中產(chǎn)生的黑白相間的亮暗點(diǎn)噪聲［14］。圖像去噪的經(jīng)典方法有中值濾波、高斯濾波、均值濾波、高通濾波等。高斯濾波適用于處理高斯噪聲［15］；均值濾波同樣對(duì)高斯噪聲的去除效果較好，但會(huì)損失圖像的邊緣信息［16］；高通濾波可以抑制低頻噪聲，對(duì)圖像進(jìn)行銳化；中值濾波較適合去除椒鹽噪聲。目前，已有較多學(xué)者對(duì)上述濾波方法進(jìn)行了改進(jìn)［17-20］，無(wú)人機(jī)采集的采煤塌陷積水區(qū)影像中存在多種噪聲，單一濾波方法降噪效果有限，采用濾波耦合去噪具有較好的效果。本研究耦合去噪采用中值濾波、高通濾波算法實(shí)現(xiàn)。其中，中值濾波選擇5×5的像素模板，是因?yàn)椴捎幂^小尺寸模板（3×3）雖然能大量保留圖像細(xì)節(jié)信息，但圖像平滑效果較差；采用較大尺寸模板（7×7）雖然平滑效果好，但會(huì)丟失圖像細(xì)節(jié)信息；高通濾波算法通過(guò)sobel算子實(shí)現(xiàn)，可達(dá)到突出輪廓的目的。圖2為采用不同濾波方法得到的采煤塌陷積水區(qū)輪廓。

分析圖2可知：采用本研究耦合濾波去噪思路處理后，圖像中噪聲得到了有效抑制，水體輪廓信息提取較為精準(zhǔn)，效果優(yōu)于單純采用高通濾波算法、中值濾波算法以及不進(jìn)行任何濾波處理情況下的效果。圖2中，除了圖2（b）中噪聲點(diǎn)較多，無(wú)法計(jì)算輪廓面積外，其余3類算法處理后的輪廓面積分別為674.34 km2、779.10 km2、726.39 km2，只有濾波耦合去噪算法得到的輪廓面積與依據(jù)礦區(qū)實(shí)測(cè)數(shù)據(jù)的Pix4Dmapper軟件計(jì)算值（716.82 km2）最為接近，可見(jiàn)，本研究提出的濾波算法有助于減少圖像噪聲，突出積水區(qū)的輪廓信息。

2.1.3 二值化處理

二值化處理可以達(dá)到減少數(shù)據(jù)量并突出采煤塌陷積水區(qū)輪廓信息的目的，常用的閾值選取方法有［21］：三角法求閾值［22］、大津法［23］、平均值閾值法以及固定閾值法。采用該類方法生成的二值化圖如圖3所示。

由圖3可知：采用三角法確定閾值時(shí)，水體輪廓清晰，邊緣平滑且無(wú)明顯噪聲點(diǎn)，與實(shí)際水域的輪廓最為接近；采用大津法確定閾值時(shí)，輪廓邊緣噪聲點(diǎn)較多，無(wú)法準(zhǔn)確區(qū)分水域邊緣的雜草、濕地等干擾信息，提取的水域輪廓面積較實(shí)際偏大；采用平均值法選取閾值時(shí)，輪廓邊緣不平滑，有較多的噪聲點(diǎn)被提取出；采用固定閾值法提取的輪廓邊緣較為平滑，輪廓提取速度較快，噪聲點(diǎn)也較少，提取的輪廓面積與實(shí)際面積相接近，但對(duì)于不同的圖像，較難選取一個(gè)通用的固定閾值。故本研究采用三角法進(jìn)行閾值選取。

2.2 塌陷積水區(qū)面積計(jì)算方法

2.2.1 提取水域輪廓及坐標(biāo)

OpenCV中的Findcontours函數(shù)具有提取影像輪廓的功能，將處理后的二值化正射影像輸入，則會(huì)生成完整的嵌套輪廓層次結(jié)構(gòu)并返回Contours列表。

2.2.2 坐標(biāo)轉(zhuǎn)換

輪廓坐標(biāo)為圖上坐標(biāo)，要計(jì)算輪廓面積，需將圖上坐標(biāo)轉(zhuǎn)換為實(shí)地坐標(biāo)。圖上坐標(biāo)轉(zhuǎn)換為實(shí)地坐標(biāo)的公式為

式中，trans(0)、trans(3)代表圖像左上角點(diǎn)的實(shí)地坐標(biāo)；trans(1)代表像元寬度；trans(5)代表像元高度；若圖像指北，則參數(shù)trans(2)、trans(4)默認(rèn)為0；m、n為水域輪廓上各個(gè)點(diǎn)的像素坐標(biāo)；x、y為水域輪廓上各個(gè)點(diǎn)的實(shí)地坐標(biāo)。

2.2.3 面積計(jì)算

三角形的面積Area(i)等于其兩邊向量積（叉乘）的模的1/2，本研究將任意n邊形分解為n-2個(gè)三角形并求和，來(lái)計(jì)算不規(guī)則多邊形的面積Y，計(jì)算公式為

式中，Area(i)為第i個(gè)三角形的面積；n為輪廓的頂點(diǎn)數(shù)（邊數(shù)）。

讀取水域輪廓上各個(gè)點(diǎn)的x、y坐標(biāo)，并將其轉(zhuǎn)換為頂點(diǎn)Zp(p=1,2,3,…,n)，即水體的輪廓是頂點(diǎn)為Zp的任意多邊形，頂點(diǎn)坐標(biāo)按逆時(shí)針?lè)较蚺帕?，分別為(x1,y1),(x2,y2),…,(xn,yn)。由于提取的坐標(biāo)點(diǎn)數(shù)量較多，故采煤塌陷積水區(qū)可近似為不規(guī)則多邊形，面積計(jì)算公式為

式中，x,y為輪廓點(diǎn)的坐標(biāo)。

2.3 塌陷積水區(qū)面積計(jì)算精度評(píng)估方法

塌陷積水區(qū)面積計(jì)算精度受到無(wú)人機(jī)正射影像獲取精度、正射影像處理效果、水域輪廓坐標(biāo)點(diǎn)提取精度和面積計(jì)算精度等因素影響。因水域輪廓坐標(biāo)點(diǎn)較多，故其精度對(duì)面積計(jì)算的影響最大。本研究依據(jù)輪廓坐標(biāo)點(diǎn)的精度來(lái)評(píng)估面積計(jì)算精度。實(shí)現(xiàn)步驟有：①生成N(N>100)組輪廓坐標(biāo)點(diǎn)的高斯隨機(jī)數(shù)，賦予輪廓坐標(biāo)點(diǎn)不同的方差；②計(jì)算N組帶有隨機(jī)方差的輪廓面積；③比較每組面積的期望值與依據(jù)相關(guān)算法計(jì)算值的差異。不同坐標(biāo)點(diǎn)位測(cè)量誤差對(duì)應(yīng)的面積誤差如圖4所示。

3 實(shí)例分析

3.1 塌陷積水區(qū)概況

采煤塌陷積水區(qū)位于某礦開(kāi)采區(qū)，如圖5所示，采煤引起地表塌陷，雨水囤積與地下水匯流導(dǎo)致該區(qū)域形成了3處采煤塌陷積水區(qū)。3處塌陷積水區(qū)內(nèi)部均無(wú)水草覆蓋，水體邊緣有少量水草和樹(shù)木，水域1、2輪廓形狀不規(guī)則，用作農(nóng)田灌溉蓄水池，水域3為連通復(fù)雜輪廓，用作魚(yú)塘養(yǎng)殖。為滿足該礦對(duì)塌陷積水區(qū)生態(tài)修復(fù)及安全防護(hù)的要求，應(yīng)用所提方法對(duì)3處積水區(qū)的面積進(jìn)行了計(jì)算分析。

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

在以3處積水區(qū)為中心的矩形區(qū)域上分別布設(shè)了8個(gè)像控點(diǎn)，采用大疆Phantom4 RTK無(wú)人機(jī)獲取塌陷積水區(qū)的正射影像，用Pix4Dmapper及Global Mapper軟件提取并處理正射影像信息。

3.3 塌陷積水區(qū)面積提取

根據(jù)本研究方法實(shí)現(xiàn)思路，使用python語(yǔ)言編寫(xiě)程序，依據(jù)無(wú)人機(jī)影像自動(dòng)計(jì)算采煤塌陷積水區(qū)的面積。以水域2為例，該處水域輪廓如圖6所示，水域2輪廓點(diǎn)位提取結(jié)果如圖7所示。

3.4 塌陷積水區(qū)面積測(cè)量精度分析

本研究分別采用Pix4Dmapper軟件、RTK實(shí)測(cè)法以及本研究方法計(jì)算了塌陷積水區(qū)的面積，結(jié)果如表1所示。

由表1可知：本研究所提方法與其他兩種方法計(jì)算的水域面積相對(duì)誤差在6%以內(nèi)；所提方法計(jì)算的水域1面積大于Pix4Dmapper軟件，該方法計(jì)算的水域2、水域3的面積均小于RTK實(shí)測(cè)法和Pix4Dmap?per軟件，原因在于水域2、水域3的水體邊緣有較多樹(shù)木、雜草等干擾，水域1的水體邊緣樹(shù)木、雜草等干擾較少。所提方法在提取水體輪廓時(shí)會(huì)將水體邊緣的樹(shù)木、雜草等干擾信息排除，故計(jì)算的面積偏?。划?dāng)水體邊緣無(wú)樹(shù)木、雜草等干擾信息時(shí)，輪廓提取較為精細(xì)，故計(jì)算的面積偏大。

4 結(jié)論

（1）為滿足采煤塌陷積水區(qū)生態(tài)修復(fù)需求，提出了一種基于無(wú)人機(jī)影像的采煤塌陷積水區(qū)面積采集方法。結(jié)合某礦工程實(shí)例，對(duì)該方法的適用性進(jìn)行了分析。結(jié)果表明：該方法與商用化軟件Pix4Dmapper及RTK實(shí)測(cè)法相比，采煤塌陷積水區(qū)面積計(jì)算誤差在6%以內(nèi)，表明該方法對(duì)于采煤塌陷積水區(qū)面積計(jì)算有一定的適用性，可為礦山生態(tài)修復(fù)等工作提供可靠依據(jù)。

（2）由于無(wú)人機(jī)航攝存在續(xù)航時(shí)間較短、像控點(diǎn)布設(shè)數(shù)量與測(cè)區(qū)面積成正比的不足，故該方法適用于中小面積的采煤塌陷積水區(qū)面積計(jì)算。相對(duì)于RTK實(shí)測(cè)法和遙感法，該方法工作量較少，且時(shí)效性較強(qiáng)。