基于UVA與植被指數的排土場DEM模型構建關鍵技術

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

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

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

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

以寶日希勒露天煤礦外排土場北坡為研究對象，采用旋翼無人機平臺與地面GPS RTK系統，構建外排土場地形測量作業平臺。建立可見光植被指數-植被高度修正模型對無人機攝影測量地形進行修正。最后用檢查點對測量結果進行精度檢驗，同時與GPS RTK地形測量數據進行比較分析，對該測量方法的實用性和不足進行思考和研究，可為以后露天煤礦或其他復雜地形的測量提供參考。

1 研究區概況

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

2 數據采集及預處理

數據采集空中飛行單元采用旋翼無人機飛行平臺;測量作業單元由機載數碼相機、機載數據采集存儲子系統組成;地面控制單元由地面站和數傳電臺組成;像控點測量單元、植被實測量單元由GPS RTK和鋼尺組成。

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

2.1 像控點、檢查點和植被實測點量測

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

干預前兩組心理狀態評分比較，差異無統計學意義（P＞0.05）。干預后，研究組SAS，SDS評分明顯低于對照組，差異有統計學意義（P＜0.05）。見表1。

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

2.2 影像獲取

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

表1 無人機及飛行參數
Table 1 UAV and ground survey details

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

2.3 數據處理

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

3 可見光植被指數與地形修正模型

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

3.1 可見光植被指數-植被高度相關性分析

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

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

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

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

注:ρgreen,ρred,ρblue分別為無人機影像在綠、紅、藍波段的反射率。

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

表3 植被指數與植被高度的相關系數
Table 3 Correlation coefficient between vegetation index and vegetation height

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

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

(1)

(2)

式中，E為模型預測值與真實值的差值;n為驗證樣本數。

表4 植被指數-植被高度回歸模型與精度檢驗
Table 4 Vegetation index-vegetation height regression model and precision test

植被指數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可知，不同植被指數—植被高度回歸模型的R2,MAD,RMSE表現出明顯差異。其中,R2最高的為EXG模型(R2=0.952)，其次為RGBVI模型(R2=0.950)，最低的為NGRDI模型(R2=0.217)。驗證樣本中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種評價指標及3.2節中植被指數與植被高度相關性分析結果可知，EXG模型與RGBVI模型能有效的反映研究區植被高度。下面以RGBVI模型為例估算研究區植被高度，進而得到研究區修正后數字高程模型(Digital Elevation Model，DEM)。

3.2 地形修正模型

利用RGBVI模型指數估算得到研究區內植被高度分部如圖3所示。

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

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

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

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

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

4 結果與分析

4.1 影像精度評估

為確定修正后地形數據的準確性，采用像控點誤差統計和檢查點驗證2種方式進行精度評估。像控點誤差統計見表6，檢查點高程值與影像解算對應高程比較見表7。

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

表6 像控點誤差統計
Table 6 Photo-control-points error statistics

點號高程誤差/m投影誤差/像素校準影像數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 檢查點與解算點高程比較
Table 7 Comparison between checkpointsand solution points elevation

檢查點號修正前影像高程/m修正后影像高程/m檢查點高程/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))對檢查點的可靠性進行評估[40]，檢查點均方根誤差由修正前的0.111 m變為0.045 m，精度顯著提高，根據《數字航空攝影測量空中三角測量規范》(GB/T 23236—2009)的要求，測量結果符合1∶500地形數據制作要求。

4.2 基于GPS RTK的DEM數據精度驗證

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

圖5 研究區DEM數據比較Fig.5 DEM data comparison in the study area

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

表8 隨機點地面植被高度統計
Table 8 Vegetation height statistics of groundrandom points

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

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

圖6 隨機點高程差正態分布曲線Fig.6 Difference of elevation normal curve of random points

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

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

5 結論與討論

(1)利用無人機攝影測量技術，采用正射與傾斜相結合的影像獲取方式保證影像質量，結合地面像控點保證影像精度，可獲取研究區高分辨率DSM數據。

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

(3)通過檢查點檢驗，修正后DEM數據精度顯著提高，滿足1∶500地形數據制作要求。通過隨機點和剖面線對比已有GPS RTK測量地形數據，修正后DEM數據在地形復雜多變區可更好的展示真實地形細節。

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