基于UVA與植被指數(shù)的排土場DEM模型構(gòu)建關(guān)鍵技術(shù)

宮傳剛，卞正富，卞和方，雷少剛,黃 赳，張周愛，郭海橋，張 浩

(1.中國礦業(yè)大學(xué) 國土資源研究所，江蘇 徐州 221116; 2.礦山生態(tài)修復(fù)教育部工程研究中心，江蘇 徐州 221116; 3.神華寶日希勒能源有限公司，內(nèi)蒙古 呼倫貝爾 021500)

隨著無人機(jī)(Unmanned Aerial Vehicle,UAV)與傳感器的不斷發(fā)展，基于無人機(jī)平臺(tái)的低空(100～1 000 m)/超低空(<100 m)攝影測量技術(shù)已顯示出獨(dú)特優(yōu)越性。無人機(jī)攝影測量技術(shù)結(jié)合了無人機(jī)駕駛飛行器技術(shù)、遙感傳感器技術(shù)、遙測遙控技術(shù)、通信技術(shù)、POS定位定姿態(tài)技術(shù)、GPS差分定位技術(shù)和遙感應(yīng)用技術(shù)，具有系統(tǒng)靈活、起降方便、成本低、質(zhì)量高、自動(dòng)化和智能化等特點(diǎn)[6-8]，在快速獲取小區(qū)域、飛行困難地區(qū)和復(fù)雜地形區(qū)域高分辨率影像方面具有明顯優(yōu)勢[9]。目前，無人機(jī)攝影測量技術(shù)已被廣泛用于地形測量[10-13]，環(huán)境監(jiān)測[14-15]，農(nóng)作物信息提取及應(yīng)用[16-17]，森林植被信息提取[18-19]，濕地、水源地信息提取[20-21]，災(zāi)害應(yīng)急救援[22-24]，電力巡視[25]，文物保護(hù)[26]，建筑物風(fēng)險(xiǎn)評(píng)估[27]等方面，其空間信息密度與精度可達(dá)機(jī)載雷達(dá)水平[28-29]。

然而，基于可見光立體像對匹配獲取的地形數(shù)據(jù)大多未考慮植被與其他地物對高程的影響，因此在高精度地形提取方面難以保證垂直精度。針對這一問題，在植被類型單一且無建筑物干擾區(qū)，可以利用可見光植被指數(shù)間接反映植被高度，利用植被指數(shù)-植被高度回歸模型剔除植被對地形的影響，從而獲得真實(shí)地形。

以寶日希勒露天煤礦外排土場北坡為研究對象，采用旋翼無人機(jī)平臺(tái)與地面GPS RTK系統(tǒng)，構(gòu)建外排土場地形測量作業(yè)平臺(tái)。建立可見光植被指數(shù)-植被高度修正模型對無人機(jī)攝影測量地形進(jìn)行修正。最后用檢查點(diǎn)對測量結(jié)果進(jìn)行精度檢驗(yàn)，同時(shí)與GPS RTK地形測量數(shù)據(jù)進(jìn)行比較分析，對該測量方法的實(shí)用性和不足進(jìn)行思考和研究，可為以后露天煤礦或其他復(fù)雜地形的測量提供參考。

1 研究區(qū)概況

寶日希勒露天礦位于內(nèi)蒙古自治區(qū)呼倫貝爾市陳巴爾虎旗寶日希勒鎮(zhèn)境內(nèi)，地理坐標(biāo)為東經(jīng):119.689°～119.761°，北緯:49.364°～49.412°。研究區(qū)位于寶日希勒露天煤礦外排土場北坡(圖1)，總面積約2.2 km2，海拔為620～720 m，坡面角度為0°～30°，區(qū)內(nèi)土地類型以草地與裸土地為主，植被以多裂葉荊芥(Schizonepeta multifida)、貝加爾針茅(Stipa baicalensis)和扁蓿豆(Pocokia ruthenica)為主。由于坡面起伏、土質(zhì)松散，在雨水沖刷下形成多條侵蝕溝，使得坡面形態(tài)復(fù)雜。

2 數(shù)據(jù)采集及預(yù)處理

數(shù)據(jù)采集空中飛行單元采用旋翼無人機(jī)飛行平臺(tái);測量作業(yè)單元由機(jī)載數(shù)碼相機(jī)、機(jī)載數(shù)據(jù)采集存儲(chǔ)子系統(tǒng)組成;地面控制單元由地面站和數(shù)傳電臺(tái)組成;像控點(diǎn)測量單元、植被實(shí)測量單元由GPS RTK和鋼尺組成。

圖1 研究區(qū)地理位置示意Fig.1 Location map of study area

2.1 像控點(diǎn)、檢查點(diǎn)和植被實(shí)測點(diǎn)量測

無人機(jī)上搭載的GPS精度較低，因此需對航拍添加地面像控點(diǎn)以及POS差分?jǐn)?shù)據(jù)以保證精度[30]。此外，在研究區(qū)布置更多的控制點(diǎn)(檢查點(diǎn))，可以對解算的數(shù)據(jù)進(jìn)行精度檢驗(yàn)。像控點(diǎn)按照GB/T 7931—2008《1∶500，1∶1 000，1∶2 000 地形圖航空攝影測量外業(yè)規(guī)范》采用內(nèi)業(yè)人員設(shè)計(jì)、外業(yè)人員布設(shè)的方式進(jìn)行。

干預(yù)前兩組心理狀態(tài)評(píng)分比較，差異無統(tǒng)計(jì)學(xué)意義（P＞0.05）。干預(yù)后，研究組SAS，SDS評(píng)分明顯低于對照組，差異有統(tǒng)計(jì)學(xué)意義（P＜0.05）。見表1。

研究區(qū)共布設(shè)像控點(diǎn)19個(gè)，檢查點(diǎn)8個(gè)，植被高度實(shí)測點(diǎn)43個(gè)(圖1)。像控點(diǎn)均布設(shè)在無植被覆蓋區(qū)且保證航向及旁向6片重疊范圍內(nèi)可見;檢查點(diǎn)與植被高度實(shí)測點(diǎn)布置在坡頂、坡腳、平臺(tái)及坡面等區(qū)域，其中1個(gè)檢查點(diǎn)位于裸土地上、7個(gè)檢查點(diǎn)位于植被區(qū);植被高度實(shí)測點(diǎn)采用5點(diǎn)取樣法在10 cm×10 cm樣方內(nèi)測得。所有點(diǎn)位置信息均采用GPS RTK進(jìn)行測量，測量模式為“地面控制點(diǎn)模式”，單點(diǎn)測量180次取平均值。

2.2 影像獲取

研究區(qū)東西長約2 400 m，南北寬約900 m，在設(shè)計(jì)航線時(shí)需根據(jù)航向重疊、旁向重疊、地面分辨率的要求來設(shè)定航高以及拍攝間隔，再由航高和相機(jī)參數(shù)來確定航線地面最低點(diǎn)分辨率[31]。無人機(jī)及相機(jī)等參數(shù)見表1，飛行中由GPS飛行控制系統(tǒng)控制相機(jī)快門進(jìn)行定點(diǎn)曝光，相機(jī)設(shè)置為固定無窮遠(yuǎn)對焦、固定光圈以保證統(tǒng)一物鏡畸變參數(shù)。除此之外，由于研究區(qū)地勢起伏大，為減小影像畸變，重點(diǎn)區(qū)域采用傾斜拍攝與垂直拍攝相結(jié)合的影像獲取方式，充分利用側(cè)視影像來獲取地物的側(cè)面紋理信息，保證整體精度。實(shí)驗(yàn)數(shù)據(jù)采集于2017-08-30，天氣晴朗無云，微風(fēng)，共采集航片1 282張。

表1 無人機(jī)及飛行參數(shù)
Table 1 UAV and ground survey details

飛行參數(shù)型號(hào)及數(shù)值傳感器FC300X_3.6_4000x3000(RGB)航線相對高度/m110航向重疊度/%80旁向重疊度/%75采集照片數(shù)量1 282地面像控點(diǎn)數(shù)量19地面檢查點(diǎn)數(shù)量8影像分辨率/(cm·pix-1)5

2.3 數(shù)據(jù)處理

將外業(yè)采集影像與像控點(diǎn)導(dǎo)入Pix4D Mapper，在軟件中對相機(jī)參數(shù)進(jìn)行檢校，并將每個(gè)像控點(diǎn)轉(zhuǎn)刺到5～8張影像上。利用軟件進(jìn)行空中三角測量與密集匹配[32]、SfM攝影測量處理、多視圖三維立體重建和幾何模型建立等處理，進(jìn)而生成高精度三維點(diǎn)云數(shù)據(jù)和數(shù)字地表模型(Digital Surface Model，DSM)，再對影像進(jìn)行數(shù)字微分糾正、拼接、鑲嵌等處理，可獲得整個(gè)研究區(qū)數(shù)字正射影像(Digital Orthophoto Map，DOM)。其中，DSM數(shù)據(jù)為包含地面植被高度信息的地面高程模型，DOM數(shù)據(jù)為含有RGB波段的正射影像(圖2)，其空間分辨率均為5 cm。

3 可見光植被指數(shù)與地形修正模型

由于Pix4D Mapper解算獲得的DSM數(shù)據(jù)中高程值包含了地表植被高度，為獲取研究區(qū)真實(shí)地形(DEM數(shù)據(jù))，需建立可見光植被指數(shù)-植被高度模型修正DSM數(shù)據(jù)。

3.1 可見光植被指數(shù)-植被高度相關(guān)性分析

植被指數(shù)是地物反射波段間的不同組合方式，目前已公開發(fā)表的植被指數(shù)模型超過150種[33]，其中基于紅光、近紅外波段反射率差異構(gòu)建的植被指數(shù)NDVI得到廣泛應(yīng)用。同樣，綠色植被在可見光波段的光譜差異性也有大量研究，其部分常見可見光植被指數(shù)及其變形式見表2。

圖2 研究區(qū)數(shù)字正射影像與數(shù)字地表模型示意Fig.2 Orthophotomap and DSM data of study area

表2 主要可見光植被指數(shù)
Table 2 Major vegetation indices of visible bands

植被指數(shù)名稱計(jì)算公式綠紅比值指數(shù)[34](Green-Red Ratio Index,GRRI)ρgreenρred綠藍(lán)比值指數(shù)[35](Green-Blue Ratio Index,GBRI)ρgreenρblue歸一化綠紅差異指數(shù)[36](Normalized Green-Red Difference Index,NGRDI)ρgreen-ρredρgreen+ρred歸一化綠藍(lán)差異指數(shù)[37](Normalized Green-Blue Difference Index,NGBDI)ρgreen-ρblueρgreen+ρblue超綠指數(shù)[38](Excess Green,EXG)2ρgreen-ρred-ρblue紅綠藍(lán)植被指數(shù)[39](Red Green Blue Vegetation Index,RGBVI)ρ2green-(ρredρblue)ρ2green+(ρredρblue)

注:ρgreen,ρred,ρblue分別為無人機(jī)影像在綠、紅、藍(lán)波段的反射率。

表2中，為了使所有植被指數(shù)與植被高度呈正相關(guān)關(guān)系，取原文獻(xiàn)中GRRI與GBRI植被指數(shù)的倒數(shù)。利用表2所示公式分別計(jì)算出研究區(qū)地面實(shí)測樣本的6種植被指數(shù)，利用SPSS建立植被指數(shù)-植被高度相關(guān)性模型，其相關(guān)性見表3。由表3可知，研究區(qū)植被高度與植被指數(shù)均呈正相關(guān)關(guān)系，其中植被高度與RGBVI的相關(guān)性系數(shù)最高(R2=0.942)，其次為EXG(R2=0.908)，最低的為GRRI(R2=0.531)。在所列植被指數(shù)之間，GRRI與NGRDI,GBRI 與NGBDI相關(guān)性較高，絕對值均大于0.99。對比發(fā)現(xiàn)，與植被高度相關(guān)性較高的植被指數(shù)均由紅、綠、藍(lán)3個(gè)波段通過不同組合方式計(jì)算所得，如RGBVI和EXG。

表3 植被指數(shù)與植被高度的相關(guān)系數(shù)
Table 3 Correlation coefficient between vegetation index and vegetation height

注:**表示在置信度(雙測)為0.01時(shí)顯著相關(guān)。

利用植被高度實(shí)測樣本(隨機(jī)布置在草地和裸土地上，共43個(gè))分為兩部分，隨機(jī)選擇70%(30個(gè))用于建立植被指數(shù)-植被高度回歸模型，30%(13個(gè))用于檢驗(yàn)。利用SPSS回歸分析模塊擬合得到植被高度與各種植被指數(shù)之間最佳擬合曲線為線性擬合。統(tǒng)計(jì)回歸模型的相關(guān)性系數(shù)R2和檢驗(yàn)樣本的平均絕對誤差(Mean Absolute Deviation，MAD)及均方根誤差(Root Mean Squared Error，RMSE)，公式如下，統(tǒng)計(jì)結(jié)果見表4。

(1)

(2)

式中，E為模型預(yù)測值與真實(shí)值的差值;n為驗(yàn)證樣本數(shù)。

表4 植被指數(shù)-植被高度回歸模型與精度檢驗(yàn)
Table 4 Vegetation index-vegetation height regression model and precision test

植被指數(shù)VI模型R2MAD/mRMSE/mGRRIy=0.757x-0.6650.2740.075 0.097GBRIy=0.511x-0.4620.7670.046 0.056NGRDIy=1.796x+0.0730.2170.080 0.102NGBDIy=1.424x+0.0260.7890.040 0.053EXGy=0.007x+0.0230.9520.020 0.025RGBVIy=1.280x+0.0220.9500.022 0.026

由表4可知，不同植被指數(shù)—植被高度回歸模型的R2,MAD,RMSE表現(xiàn)出明顯差異。其中,R2最高的為EXG模型(R2=0.952)，其次為RGBVI模型(R2=0.950)，最低的為NGRDI模型(R2=0.217)。驗(yàn)證樣本中MAD和RMSE最小值為EXG模型(MAD=0.02 m，RMSE=0.025 m)，其次為RGBVI模型(MAD=0.02 m，RMSE=0.025 m)，最大值為NGRDI模型(MAD=0.08 m，RMSE=0.102 m)。

綜合分析3種評(píng)價(jià)指標(biāo)及3.2節(jié)中植被指數(shù)與植被高度相關(guān)性分析結(jié)果可知，EXG模型與RGBVI模型能有效的反映研究區(qū)植被高度。下面以RGBVI模型為例估算研究區(qū)植被高度，進(jìn)而得到研究區(qū)修正后數(shù)字高程模型(Digital Elevation Model，DEM)。

3.2 地形修正模型

利用RGBVI模型指數(shù)估算得到研究區(qū)內(nèi)植被高度分部如圖3所示。

由圖3可知，研究區(qū)內(nèi)裸土地(植被高度為0)面積占總面積32.7%，植被高度低于15 cm區(qū)域面積占總面積60.9%，植被高度高于15 cm區(qū)域面積占總面積6.4%。整體來看研究區(qū)東部比西部植被高度高，坡頂和平臺(tái)比坡面高度高。圖3中A,B，C,D為研究區(qū)4塊典型區(qū)域，其植被高度估算結(jié)果細(xì)節(jié)如圖4所示。

由圖4可知，在研究區(qū)不同區(qū)域，植被高度對地形的影響程度與實(shí)際情況基本一致。區(qū)域A位于排土場坡頂，植被高度差異性明顯，部分區(qū)域植被高度較高，對地形測量影響明顯。區(qū)域B接近坡腳，內(nèi)有沖蝕溝發(fā)育，地形起伏較大。區(qū)域C位于坡頂與邊坡過渡區(qū)，內(nèi)有防護(hù)林與少量人工建筑，植被高度差異性大。區(qū)域D位于坡腳，區(qū)內(nèi)有公路、花壇等地物，其植被高度為0且輪廓非常清晰。

圖3 基于RGBVI的研究區(qū)植被高度估算結(jié)果Fig.3 Estimation results of vegetation height in the study area based on RGBVI model

根據(jù)以上分析，可建立地形修正模型以除去研究區(qū)內(nèi)植被高度對地形的影響。假設(shè)無人機(jī)攝影測量所獲DSM數(shù)據(jù)任意點(diǎn)O(x,y)高程值為H0，利用植被指數(shù)-植被高度修正模型估算該點(diǎn)植被高度為h0，則修正后該點(diǎn)真實(shí)坐標(biāo)為

A′(x,y,z)=A(x,y,H0-h0)

4 結(jié)果與分析

4.1 影像精度評(píng)估

為確定修正后地形數(shù)據(jù)的準(zhǔn)確性，采用像控點(diǎn)誤差統(tǒng)計(jì)和檢查點(diǎn)驗(yàn)證2種方式進(jìn)行精度評(píng)估。像控點(diǎn)誤差統(tǒng)計(jì)見表6，檢查點(diǎn)高程值與影像解算對應(yīng)高程比較見表7。

由表6，7可知，所有像控點(diǎn)和檢查點(diǎn)均有效，未出現(xiàn)異常值。其中像控點(diǎn)投影誤差均控制在0.5個(gè)像素之內(nèi)，精度滿足設(shè)計(jì)要求;修正后影像高程大多小于修正前，點(diǎn)CheckP-5處為裸土地，植被指數(shù)為0，修正前后影像高程值不變。

表6 像控點(diǎn)誤差統(tǒng)計(jì)
Table 6 Photo-control-points error statistics

點(diǎn)號(hào)高程誤差/m投影誤差/像素校準(zhǔn)影像數(shù)GCP-01-0.0230.4817GCP-020.0210.4275GCP-030.0390.4435GCP-040.0190.3695GCP-05-0.0240.1455GCP-060.0260.2635GCP-070.0180.3067GCP-080.0220.2655GCP-090.0140.1945GCP-10-0.0060.4745GCP-110.0190.4385GCP-12-0.0180.4575GCP-130.0190.3155GCP-14-0.0570.3625GCP-15-0.0610.2166GCP-16-0.0770.4795GCP-17-0.0250.4778GCP-180.0210.4708GCP-190.0210.3498

表7 檢查點(diǎn)與解算點(diǎn)高程比較
Table 7 Comparison between checkpointsand solution points elevation

檢查點(diǎn)號(hào)修正前影像高程/m修正后影像高程/m檢查點(diǎn)高程/mCheckP-1651.110651.014651.038CheckP-2625.712625.614625.591CheckP-3687.992687.771687.861CheckP-4631.266631.254631.221CheckP-5716.697716.697716.667CheckP-6697.706697.614697.631CheckP-7720.503720.411720.352CheckP-8698.755698.547698.579

使用均方根誤差式(式(2))對檢查點(diǎn)的可靠性進(jìn)行評(píng)估[40]，檢查點(diǎn)均方根誤差由修正前的0.111 m變?yōu)?.045 m，精度顯著提高，根據(jù)《數(shù)字航空攝影測量空中三角測量規(guī)范》(GB/T 23236—2009)的要求，測量結(jié)果符合1∶500地形數(shù)據(jù)制作要求。

4.2 基于GPS RTK的DEM數(shù)據(jù)精度驗(yàn)證

為進(jìn)一步評(píng)估修正后無人機(jī)攝影測量地形數(shù)據(jù)的可靠性，將修正后地形數(shù)據(jù)與同期使用GPS RTK所測地形數(shù)據(jù)進(jìn)行比較和分析(圖5)。由于無人機(jī)攝影測量數(shù)據(jù)由點(diǎn)云生成，而GPS RTK地形測量數(shù)據(jù)由多點(diǎn)插值所得，因此無人機(jī)攝影測量(圖5(a))解算數(shù)據(jù)比已有GPS RTK地形測量數(shù)據(jù)(圖5(b))更加精細(xì)，在測量復(fù)雜地形時(shí)無人機(jī)攝影測量體現(xiàn)出明顯優(yōu)勢。同時(shí)，將兩組DEM數(shù)據(jù)相重疊(圖5(c))可知，兩組數(shù)據(jù)在不同區(qū)域呈現(xiàn)出不同高程值，下面進(jìn)行具體論述。

圖5 研究區(qū)DEM數(shù)據(jù)比較Fig.5 DEM data comparison in the study area

利用ArcGIS在研究區(qū)內(nèi)取1 000個(gè)隨機(jī)點(diǎn)(圖5(c))，隨機(jī)點(diǎn)處地面植被高度見表8。分別將修正后無人機(jī)解算高程數(shù)據(jù)和GPS RTK地形測量高程數(shù)據(jù)提取至隨機(jī)點(diǎn)，建立同名點(diǎn)高程屬性對照表，得出兩組數(shù)據(jù)同名點(diǎn)高程差。除去異常值后，統(tǒng)計(jì)隨機(jī)點(diǎn)高程差正態(tài)分布直方圖如圖6所示。

表8 隨機(jī)點(diǎn)地面植被高度統(tǒng)計(jì)
Table 8 Vegetation height statistics of groundrandom points

植被高度/cm0(裸地)0～55～1010～20>20所占比例/%35.9024.7019.9015.603.90

由圖6可知，無人機(jī)攝影測量解算高程數(shù)據(jù)與GPS RTK地形測量高程數(shù)據(jù)高程差主要集中在0.1 m內(nèi)，但也有部分點(diǎn)的差值超過0.1 m。

圖6 隨機(jī)點(diǎn)高程差正態(tài)分布曲線Fig.6 Difference of elevation normal curve of random points

為進(jìn)一步研究誤差來源，在研究區(qū)中心區(qū)隨機(jī)取3條剖面線(圖7)進(jìn)行地形對比。根據(jù)實(shí)地調(diào)查，研究區(qū)排土場以平臺(tái)加斜坡的形式堆放，其中斜坡呈20°～30°，在從平臺(tái)過度到斜坡時(shí)地形變化劇烈、棱角分明。由圖7可知，修正后無人機(jī)解算地形與現(xiàn)實(shí)相符，而GPS RTK測量地形數(shù)據(jù)只能在整體上體現(xiàn)地形起伏。在地形平坦處兩組數(shù)據(jù)高度重合，在地形突變處差異較大。

圖7 無人機(jī)攝影測量解算數(shù)據(jù)與已有GPS RTK測量數(shù)據(jù)剖面對比示意Fig.7 Cross-sections through modelled surface with perpendicular distance between the surface obtained from UVA photogrammetric data and the reference surface

5 結(jié)論與討論

(1)利用無人機(jī)攝影測量技術(shù)，采用正射與傾斜相結(jié)合的影像獲取方式保證影像質(zhì)量，結(jié)合地面像控點(diǎn)保證影像精度，可獲取研究區(qū)高分辨率DSM數(shù)據(jù)。

(2)基于紅綠藍(lán)波段的植被指數(shù)(RGBVI和EXG)回歸模型對植被高度具有較好的模擬和預(yù)測精度，可用于研究區(qū)內(nèi)植被高度的估算。利用回歸模型可將DSM數(shù)據(jù)修正為DEM數(shù)據(jù)。

(3)通過檢查點(diǎn)檢驗(yàn)，修正后DEM數(shù)據(jù)精度顯著提高，滿足1∶500地形數(shù)據(jù)制作要求。通過隨機(jī)點(diǎn)和剖面線對比已有GPS RTK測量地形數(shù)據(jù)，修正后DEM數(shù)據(jù)在地形復(fù)雜多變區(qū)可更好的展示真實(shí)地形細(xì)節(jié)。

(4)相對于傳統(tǒng)地形測量手段，利用無人機(jī)平臺(tái)可大幅提高工作效率和數(shù)據(jù)精度。但同時(shí)，由于可見光波段信息較少，本文中建立的可見光植被指數(shù)-植被高度回歸模型不具有普適性。在植被覆蓋復(fù)雜區(qū)，可先對研究區(qū)進(jìn)行植被類型分類，對不同類型植被區(qū)分別建模以提高擬合精度。