宮傳剛,卞正富,卞和方,雷少剛,黃 赳,張周愛,郭海橋,張 浩
(1.中國(guó)礦業(yè)大學(xué) 國(guó)土資源研究所,江蘇 徐州 221116; 2.礦山生態(tài)修復(fù)教育部工程研究中心,江蘇 徐州 221116; 3.神華寶日希勒能源有限公司,內(nèi)蒙古 呼倫貝爾 021500)
隨著無(wú)人機(jī)(Unmanned Aerial Vehicle,UAV)與傳感器的不斷發(fā)展,基于無(wú)人機(jī)平臺(tái)的低空(100~1 000 m)/超低空(<100 m)攝影測(cè)量技術(shù)已顯示出獨(dú)特優(yōu)越性。無(wú)人機(jī)攝影測(cè)量技術(shù)結(jié)合了無(wú)人機(jī)駕駛飛行器技術(shù)、遙感傳感器技術(shù)、遙測(cè)遙控技術(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)勢(shì)[9]。目前,無(wú)人機(jī)攝影測(cè)量技術(shù)已被廣泛用于地形測(cè)量[10-13],環(huán)境監(jiān)測(cè)[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]。
然而,基于可見光立體像對(duì)匹配獲取的地形數(shù)據(jù)大多未考慮植被與其他地物對(duì)高程的影響,因此在高精度地形提取方面難以保證垂直精度。針對(duì)這一問題,在植被類型單一且無(wú)建筑物干擾區(qū),可以利用可見光植被指數(shù)間接反映植被高度,利用植被指數(shù)-植被高度回歸模型剔除植被對(duì)地形的影響,從而獲得真實(shí)地形。
以寶日希勒露天煤礦外排土場(chǎng)北坡為研究對(duì)象,采用旋翼無(wú)人機(jī)平臺(tái)與地面GPS RTK系統(tǒng),構(gòu)建外排土場(chǎng)地形測(cè)量作業(yè)平臺(tái)。建立可見光植被指數(shù)-植被高度修正模型對(duì)無(wú)人機(jī)攝影測(cè)量地形進(jìn)行修正。最后用檢查點(diǎn)對(duì)測(cè)量結(jié)果進(jìn)行精度檢驗(yàn),同時(shí)與GPS RTK地形測(cè)量數(shù)據(jù)進(jìn)行比較分析,對(duì)該測(cè)量方法的實(shí)用性和不足進(jìn)行思考和研究,可為以后露天煤礦或其他復(fù)雜地形的測(cè)量提供參考。
寶日希勒露天礦位于內(nèi)蒙古自治區(qū)呼倫貝爾市陳巴爾虎旗寶日希勒鎮(zhèn)境內(nèi),地理坐標(biāo)為東經(jīng):119.689°~119.761°,北緯:49.364°~49.412°。研究區(qū)位于寶日希勒露天煤礦外排土場(chǎng)北坡(圖1),總面積約2.2 km2,海拔為620~720 m,坡面角度為0°~30°,區(qū)內(nèi)土地類型以草地與裸土地為主,植被以多裂葉荊芥(Schizonepeta multifida)、貝加爾針茅(Stipa baicalensis)和扁蓿豆(Pocokia ruthenica)為主。由于坡面起伏、土質(zhì)松散,在雨水沖刷下形成多條侵蝕溝,使得坡面形態(tài)復(fù)雜。
數(shù)據(jù)采集空中飛行單元采用旋翼無(wú)人機(jī)飛行平臺(tái);測(cè)量作業(yè)單元由機(jī)載數(shù)碼相機(jī)、機(jī)載數(shù)據(jù)采集存儲(chǔ)子系統(tǒng)組成;地面控制單元由地面站和數(shù)傳電臺(tái)組成;像控點(diǎn)測(cè)量單元、植被實(shí)測(cè)量單元由GPS RTK和鋼尺組成。
圖1 研究區(qū)地理位置示意Fig.1 Location map of study area
無(wú)人機(jī)上搭載的GPS精度較低,因此需對(duì)航拍添加地面像控點(diǎn)以及POS差分?jǐn)?shù)據(jù)以保證精度[30]。此外,在研究區(qū)布置更多的控制點(diǎn)(檢查點(diǎn)),可以對(duì)解算的數(shù)據(jù)進(jìn)行精度檢驗(yàn)。像控點(diǎn)按照GB/T 7931—2008《1∶500,1∶1 000,1∶2 000 地形圖航空攝影測(cè)量外業(yè)規(guī)范》采用內(nèi)業(yè)人員設(shè)計(jì)、外業(yè)人員布設(shè)的方式進(jìn)行。
干預(yù)前兩組心理狀態(tài)評(píng)分比較,差異無(wú)統(tǒng)計(jì)學(xué)意義(P>0.05)。干預(yù)后,研究組SAS,SDS評(píng)分明顯低于對(duì)照組,差異有統(tǒng)計(jì)學(xué)意義(P<0.05)。見表1。
研究區(qū)共布設(shè)像控點(diǎn)19個(gè),檢查點(diǎn)8個(gè),植被高度實(shí)測(cè)點(diǎn)43個(gè)(圖1)。像控點(diǎn)均布設(shè)在無(wú)植被覆蓋區(qū)且保證航向及旁向6片重疊范圍內(nèi)可見;檢查點(diǎn)與植被高度實(shí)測(cè)點(diǎn)布置在坡頂、坡腳、平臺(tái)及坡面等區(qū)域,其中1個(gè)檢查點(diǎn)位于裸土地上、7個(gè)檢查點(diǎn)位于植被區(qū);植被高度實(shí)測(cè)點(diǎn)采用5點(diǎn)取樣法在10 cm×10 cm樣方內(nèi)測(cè)得。所有點(diǎn)位置信息均采用GPS RTK進(jìn)行測(cè)量,測(cè)量模式為“地面控制點(diǎn)模式”,單點(diǎn)測(cè)量180次取平均值。
研究區(qū)東西長(zhǎng)約2 400 m,南北寬約900 m,在設(shè)計(jì)航線時(shí)需根據(jù)航向重疊、旁向重疊、地面分辨率的要求來(lái)設(shè)定航高以及拍攝間隔,再由航高和相機(jī)參數(shù)來(lái)確定航線地面最低點(diǎn)分辨率[31]。無(wú)人機(jī)及相機(jī)等參數(shù)見表1,飛行中由GPS飛行控制系統(tǒng)控制相機(jī)快門進(jìn)行定點(diǎn)曝光,相機(jī)設(shè)置為固定無(wú)窮遠(yuǎn)對(duì)焦、固定光圈以保證統(tǒng)一物鏡畸變參數(shù)。除此之外,由于研究區(qū)地勢(shì)起伏大,為減小影像畸變,重點(diǎn)區(qū)域采用傾斜拍攝與垂直拍攝相結(jié)合的影像獲取方式,充分利用側(cè)視影像來(lái)獲取地物的側(cè)面紋理信息,保證整體精度。實(shí)驗(yàn)數(shù)據(jù)采集于2017-08-30,天氣晴朗無(wú)云,微風(fēng),共采集航片1 282張。
表1 無(wú)人機(jī)及飛行參數(shù)
Table 1 UAV and ground survey details
飛行參數(shù)型號(hào)及數(shù)值傳感器FC300X_3.6_4000x3000(RGB)航線相對(duì)高度/m110航向重疊度/%80旁向重疊度/%75采集照片數(shù)量1 282地面像控點(diǎn)數(shù)量19地面檢查點(diǎn)數(shù)量8影像分辨率/(cm·pix-1)5
將外業(yè)采集影像與像控點(diǎn)導(dǎo)入Pix4D Mapper,在軟件中對(duì)相機(jī)參數(shù)進(jìn)行檢校,并將每個(gè)像控點(diǎn)轉(zhuǎn)刺到5~8張影像上。利用軟件進(jìn)行空中三角測(cè)量與密集匹配[32]、SfM攝影測(cè)量處理、多視圖三維立體重建和幾何模型建立等處理,進(jìn)而生成高精度三維點(diǎn)云數(shù)據(jù)和數(shù)字地表模型(Digital Surface Model,DSM),再對(duì)影像進(jìn)行數(shù)字微分糾正、拼接、鑲嵌等處理,可獲得整個(gè)研究區(qū)數(shù)字正射影像(Digital Orthophoto Map,DOM)。其中,DSM數(shù)據(jù)為包含地面植被高度信息的地面高程模型,DOM數(shù)據(jù)為含有RGB波段的正射影像(圖2),其空間分辨率均為5 cm。
由于Pix4D Mapper解算獲得的DSM數(shù)據(jù)中高程值包含了地表植被高度,為獲取研究區(qū)真實(shí)地形(DEM數(shù)據(jù)),需建立可見光植被指數(shù)-植被高度模型修正DSM數(shù)據(jù)。
植被指數(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分別為無(wú)人機(jī)影像在綠、紅、藍(lán)波段的反射率。
表2中,為了使所有植被指數(shù)與植被高度呈正相關(guān)關(guān)系,取原文獻(xiàn)中GRRI與GBRI植被指數(shù)的倒數(shù)。利用表2所示公式分別計(jì)算出研究區(qū)地面實(shí)測(cè)樣本的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)性較高,絕對(duì)值均大于0.99。對(duì)比發(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
注:**表示在置信度(雙測(cè))為0.01時(shí)顯著相關(guān)。
利用植被高度實(shí)測(cè)樣本(隨機(jī)布置在草地和裸土地上,共43個(gè))分為兩部分,隨機(jī)選擇70%(30個(gè))用于建立植被指數(shù)-植被高度回歸模型,30%(13個(gè))用于檢驗(yàn)。利用SPSS回歸分析模塊擬合得到植被高度與各種植被指數(shù)之間最佳擬合曲線為線性擬合。統(tǒng)計(jì)回歸模型的相關(guān)性系數(shù)R2和檢驗(yàn)樣本的平均絕對(duì)誤差(Mean Absolute Deviation,MAD)及均方根誤差(Root Mean Squared Error,RMSE),公式如下,統(tǒng)計(jì)結(jié)果見表4。
(1)
(2)
式中,E為模型預(yù)測(cè)值與真實(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)。
利用RGBVI模型指數(shù)估算得到研究區(qū)內(nèi)植被高度分部如圖3所示。
由圖3可知,研究區(qū)內(nèi)裸土地(植被高度為0)面積占總面積32.7%,植被高度低于15 cm區(qū)域面積占總面積60.9%,植被高度高于15 cm區(qū)域面積占總面積6.4%。整體來(lái)看研究區(qū)東部比西部植被高度高,坡頂和平臺(tái)比坡面高度高。圖3中A,B,C,D為研究區(qū)4塊典型區(qū)域,其植被高度估算結(jié)果細(xì)節(jié)如圖4所示。
由圖4可知,在研究區(qū)不同區(qū)域,植被高度對(duì)地形的影響程度與實(shí)際情況基本一致。區(qū)域A位于排土場(chǎng)坡頂,植被高度差異性明顯,部分區(qū)域植被高度較高,對(duì)地形測(cè)量影響明顯。區(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)植被高度對(duì)地形的影響。假設(shè)無(wú)人機(jī)攝影測(cè)量所獲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)
為確定修正后地形數(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)高程值與影像解算對(duì)應(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))對(duì)檢查點(diǎn)的可靠性進(jìn)行評(píng)估[40],檢查點(diǎn)均方根誤差由修正前的0.111 m變?yōu)?.045 m,精度顯著提高,根據(jù)《數(shù)字航空攝影測(cè)量空中三角測(cè)量規(guī)范》(GB/T 23236—2009)的要求,測(cè)量結(jié)果符合1∶500地形數(shù)據(jù)制作要求。
為進(jìn)一步評(píng)估修正后無(wú)人機(jī)攝影測(cè)量地形數(shù)據(jù)的可靠性,將修正后地形數(shù)據(jù)與同期使用GPS RTK所測(cè)地形數(shù)據(jù)進(jìn)行比較和分析(圖5)。由于無(wú)人機(jī)攝影測(cè)量數(shù)據(jù)由點(diǎn)云生成,而GPS RTK地形測(cè)量數(shù)據(jù)由多點(diǎn)插值所得,因此無(wú)人機(jī)攝影測(cè)量(圖5(a))解算數(shù)據(jù)比已有GPS RTK地形測(cè)量數(shù)據(jù)(圖5(b))更加精細(xì),在測(cè)量復(fù)雜地形時(shí)無(wú)人機(jī)攝影測(cè)量體現(xiàn)出明顯優(yōu)勢(shì)。同時(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。分別將修正后無(wú)人機(jī)解算高程數(shù)據(jù)和GPS RTK地形測(cè)量高程數(shù)據(jù)提取至隨機(jī)點(diǎn),建立同名點(diǎn)高程屬性對(duì)照表,得出兩組數(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可知,無(wú)人機(jī)攝影測(cè)量解算高程數(shù)據(jù)與GPS RTK地形測(cè)量高程數(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)一步研究誤差來(lái)源,在研究區(qū)中心區(qū)隨機(jī)取3條剖面線(圖7)進(jìn)行地形對(duì)比。根據(jù)實(shí)地調(diào)查,研究區(qū)排土場(chǎng)以平臺(tái)加斜坡的形式堆放,其中斜坡呈20°~30°,在從平臺(tái)過度到斜坡時(shí)地形變化劇烈、棱角分明。由圖7可知,修正后無(wú)人機(jī)解算地形與現(xiàn)實(shí)相符,而GPS RTK測(cè)量地形數(shù)據(jù)只能在整體上體現(xiàn)地形起伏。在地形平坦處兩組數(shù)據(jù)高度重合,在地形突變處差異較大。
圖7 無(wú)人機(jī)攝影測(cè)量解算數(shù)據(jù)與已有GPS RTK測(cè)量數(shù)據(jù)剖面對(duì)比示意Fig.7 Cross-sections through modelled surface with perpendicular distance between the surface obtained from UVA photogrammetric data and the reference surface
(1)利用無(wú)人機(jī)攝影測(cè)量技術(shù),采用正射與傾斜相結(jié)合的影像獲取方式保證影像質(zhì)量,結(jié)合地面像控點(diǎn)保證影像精度,可獲取研究區(qū)高分辨率DSM數(shù)據(jù)。
(2)基于紅綠藍(lán)波段的植被指數(shù)(RGBVI和EXG)回歸模型對(duì)植被高度具有較好的模擬和預(yù)測(cè)精度,可用于研究區(qū)內(nèi)植被高度的估算。利用回歸模型可將DSM數(shù)據(jù)修正為DEM數(shù)據(jù)。
(3)通過檢查點(diǎn)檢驗(yàn),修正后DEM數(shù)據(jù)精度顯著提高,滿足1∶500地形數(shù)據(jù)制作要求。通過隨機(jī)點(diǎn)和剖面線對(duì)比已有GPS RTK測(cè)量地形數(shù)據(jù),修正后DEM數(shù)據(jù)在地形復(fù)雜多變區(qū)可更好的展示真實(shí)地形細(xì)節(jié)。
(4)相對(duì)于傳統(tǒng)地形測(cè)量手段,利用無(wú)人機(jī)平臺(tái)可大幅提高工作效率和數(shù)據(jù)精度。但同時(shí),由于可見光波段信息較少,本文中建立的可見光植被指數(shù)-植被高度回歸模型不具有普適性。在植被覆蓋復(fù)雜區(qū),可先對(duì)研究區(qū)進(jìn)行植被類型分類,對(duì)不同類型植被區(qū)分別建模以提高擬合精度。