植被季相變化對三維激光測量切溝地形的影響

王佳希, 鄧家勇, 張 巖?, 宋曉鵬, 馬 瑞

(1.北京林業(yè)大學(xué)水土保持學(xué)院，100083，北京；2.山西吉縣森林生態(tài)系統(tǒng)國家野外科學(xué)觀測研究站，100083，北京)

地形地貌是影響土壤侵蝕的重要因子。土壤侵蝕研究中，數(shù)字高程模型(digital elevation model, DEM)是重要的地形地貌基礎(chǔ)數(shù)據(jù)，可以用于提取坡度、坡向及水系特征；計算侵蝕溝形態(tài)指標(biāo)和土壤侵蝕模型中的地形因子；分析土壤侵蝕過程中的地貌變化[1-4]。DEM的制作早期只能通過實地勘測或者使用已有地形圖[5]。隨著測繪技術(shù)的不斷發(fā)展[6-7]，基于GPS[8]、衛(wèi)星和航空遙感[9]、攝影測量[10]、三維激光掃描[11]等都可制作DEM。三維激光掃描所產(chǎn)生的數(shù)據(jù)是具有空間坐標(biāo)及高程信息的點(diǎn)，稱為點(diǎn)云[11]。雖然三維激光掃描具有一定的穿透作用，但由于植被的反射作用等，點(diǎn)云中仍會存在植被等非地面點(diǎn)。所以點(diǎn)云仍需經(jīng)過濾波、插值[12-14]而最終生成高精度的DEM[13,15-19]。點(diǎn)云濾波方法有數(shù)學(xué)形態(tài)法濾波[20]、坡度濾波[21]、基于回波次數(shù)[13]、基于回波強(qiáng)度[22]等。常用的插值方法有反距離權(quán)重法(inverse distance weighted，IDW)、克里金法(Kriging)、自然鄰域法(natural neighbor，NN)[14]和不規(guī)則三角網(wǎng)(triangulated irregular network，TIN)[8]等。

點(diǎn)云濾波在平坦、植被較單一地區(qū)相對簡單，在崎嶇、植被低矮復(fù)雜地區(qū)較為困難[22]。黃土高原丘陵溝壑區(qū)土壤侵蝕嚴(yán)重，地形破碎，植被多為次生演替的灌木和草本，屬于點(diǎn)云濾波困難的地區(qū)。潘少奇等[23]將點(diǎn)云用TIN法插值生成DEM，與1∶1萬的地形圖比較，精度有所提升，但沒有植被濾波處理。馬鼎等[22]利用點(diǎn)云的激光回波強(qiáng)度衰減補(bǔ)償模型進(jìn)行了植被濾波，得到的DEM效果較好，但需要點(diǎn)云記錄有回波強(qiáng)度信息。

如上，植被對于DEM的影響研究多集中于植被濾波、插值方法的開發(fā)[24-25]，而關(guān)于植被季相變化對DEM的影響研究較少。植被非生長季時，枯萎凋落，激光腳點(diǎn)落在地面的可能性增大，DEM精度可能會提高；但是植被季相變化具體有何影響，缺乏數(shù)據(jù)支撐。筆者基于黃土高原丘陵溝壑區(qū)植物生長季 (2018年9月)、非生長季(2018年12月)2次實地測繪的點(diǎn)云，分析生長季與非生長季的點(diǎn)云以及DEM差異，討論植被季相變化是否對其具有影響，探究提高DEM精度的手段。

1 研究區(qū)概況

陜西省綏德縣橋溝小流域(E 110°17′24″～110°17′43″，N 37°29′33″～37°30′13″)位于無定河下游(圖1)，屬于黃土高原丘陵溝壑區(qū)第一副區(qū)。氣候類型為溫帶半干旱大陸性季風(fēng)氣候。土壤類型均為黃綿土亞類。植被類型以草地、灌木林地為主，有少量人工林地。植物群落多為次生演替，代表性群系有白羊草(Bothriochloaischaemum)群系、白草(Pennisetumflaccidum)-鐵桿蒿(Artemisiasacrorum)群系、茭蒿(Artemisiagiraldii)群系、茵陳蒿(Artemisiacapillaris)群系和興安胡枝子(Lespedezadavurica)群系。人工林造林主要有刺槐(Robiniapsdudoacacia)、油松(Pinustabulaeformis)、側(cè)柏(Platycladusorientalis)、檸條(Caraganakorshinskii)等[26]。在橋溝小流域中選取A、B和C樣區(qū)。樣區(qū)A(1 445 m2)切溝溝頭邊緣平滑，溝壁較陡，溝道比降在1.08～1.25之間；以禾本科，菊科等草本為主，溝底有檸條，冠層高度約20～60 cm；植被覆蓋度較低。樣區(qū)B(3 177 m2)溝壁較陡，溝頭邊緣鋒利，溝道比降在0.75～1.11之間，切溝較大；除垂直溝壁裸露外，在溝底、溝間有零星刺槐分布，其高度約2.5～3.5 m；灌草覆蓋度較高，灌木冠層高度約30～50 cm。樣區(qū)C(1 386 m2)切溝密集，溝道比降較大(1.2～2.0)，但切溝寬度較??；以低矮灌木和草本為主，灌木主要分布在溝底，冠層高度20～30 cm；植被覆蓋度中等。

A、B和C為三維激光測量樣區(qū)編號。底圖來自Google 影像。Sampling sites are named A，B and C for 3D laser measurement. Background images are from Google Earth images. 圖1 研究區(qū)位置圖Fig.1 Location of the study area

2 材料與方法

2.1 數(shù)據(jù)獲取

已有研究表明黃土高原植被覆蓋在7—9月達(dá)到峰值，12月達(dá)到最低值[27]。本研究分別于植物生長季 (2018年9月)、非生長季(2018年12月)，使用拓普康IS—IMAGING STATION型三維激光掃描儀對橋溝內(nèi)切溝特征差異較大的3個樣區(qū)(A、B和C)內(nèi)11條切溝進(jìn)行了掃描。為了保證測量精度，無障礙物遮擋，儀器架設(shè)點(diǎn)選擇在樣區(qū)對面山坡視野最佳且大致水平的位置，并用手持GPS確定架設(shè)點(diǎn)地理坐標(biāo)及高程。2次測量使用了相同的架設(shè)點(diǎn)，保證2次測量之間相對誤差較小(圖1)。點(diǎn)云平面間隔為0.15 m，密度約30～50 pts/m2。所得點(diǎn)云不包括回波次數(shù)、強(qiáng)度等。另外，在2次外業(yè)時，使用激光測距儀和卷尺抽樣測量部分優(yōu)勢種群的株高，用來輔助判斷植被季相變化。

2.2 數(shù)據(jù)處理與分析

2.2.1 點(diǎn)云處理流程 第1步，通過人工判讀除去包括飛鳥、電力線路、儀器偶然誤差等造成的離群值[13]。第2步，基于Terrascan，利用坡度濾波法[21]將點(diǎn)云分類為地面點(diǎn)和非地面點(diǎn)(大部分為植被點(diǎn))，隨后去除分類中的非地面點(diǎn)，僅保留地面點(diǎn)。第3步，將濾波前后點(diǎn)云用TIN法[8]、克里金法[28]分別插值構(gòu)建DEM(柵格分辨率0.15 m)。按照上述1～3步分別處理9月、12月數(shù)據(jù)。最終分別得到2季相的4幅DEM：TIN法、濾波-TIN法、克里金法、濾波-克里法。

2.2.2 點(diǎn)云和DEM精度評價方法 對點(diǎn)云高程進(jìn)行數(shù)理統(tǒng)計，計算點(diǎn)云平均高程、標(biāo)準(zhǔn)差、最大值和最小值等，對比2次點(diǎn)云整體差異。

采用中誤差和判定系數(shù)表征DEM和點(diǎn)云的符合程度。中誤差，即均方根誤差(root mean square error，RMSE)。它是最常用的體現(xiàn)DEM偏移情況的指標(biāo)，計算公式為：

(1)

式中：σ為中誤差，m；R為點(diǎn)云高程，m；Z為插值高程，m；n為測得的總點(diǎn)數(shù)；m為第m個點(diǎn)。中誤差越小，DEM越符合點(diǎn)云。

判定系數(shù)(R2)用以表征插值高程與實測高程總體吻合程度。計算公式為：

(2)

圖2 A、B和C樣區(qū)點(diǎn)云高程頻數(shù)分布圖Fig.2 Frequency distribution of point cloud elevation in the site A，B and C

本研究共選取7條切溝縱剖面，29條橫剖面來提取溝長、溝寬和溝深，用來代表不同方法生成的DEM地貌形態(tài)差異。另外，通過三維可視化、坡度分析和切溝縱剖面圖對比不同季相，不同濾波插值方法對生成DEM的影響。三維可視化可以直觀判斷DEM質(zhì)量，地形地貌還原程度，植被覆蓋特征等。坡度圖可以直觀地顯示出DEM坡度變化，植被在其中多表現(xiàn)為“鼓包”，另外從坡度圖中也很容易識別出溝沿線等地貌特征。本研究區(qū)地面粗糙度不大，而三維可視化圖、坡度圖中粗糙度較大區(qū)域，也一定程度上反映了植被表面的影響。

3 結(jié)果與分析

3.1 植被季相變化

3.1.1 點(diǎn)云高程統(tǒng)計特征 點(diǎn)云高程統(tǒng)計分析表明，樣區(qū)A、B和C生長季平均高程比非生長季分別高0.22、0.47和0.07 m。實地抽樣方面，共測量優(yōu)勢種灌木、多年生草本3種，9和12月分別測184和189株。其平均株高分別為0.73和0.68 m。測量的9株喬木高度無明顯變化。非生長季點(diǎn)云標(biāo)準(zhǔn)差普遍比生長季略低，且中位數(shù)、眾數(shù)也有不同程度減小，說明非生長季高程分布更趨于集中，且整體偏低。從圖2可看出3個樣區(qū)類似的偏態(tài)分布規(guī)律，且非生長季比生長季的高值部分略少，其他部分明顯增多。樣區(qū)B點(diǎn)云頻率分布特征差異最大。植被覆蓋度方面，樣區(qū)A約0.60，樣區(qū)B約0.85，樣區(qū)C約0.75，說明點(diǎn)云高程變化和植被覆蓋度有一定關(guān)系，植被覆蓋度越高，點(diǎn)云高程季相變化越大。

3.1.2 植被季相變化與濾波效果 3個樣區(qū)12月所測點(diǎn)云高程平均標(biāo)準(zhǔn)差比9月有所降低，說明非生長季的點(diǎn)云高程離散趨勢減弱。12月所測數(shù)據(jù)濾波后保留的點(diǎn)云更多(表1)。從點(diǎn)云剖面圖可看出，9月所測點(diǎn)云中高大植被在12月幾乎消失，這說明經(jīng)過坡度濾波后，高大喬灌被去除同時地面點(diǎn)盡可能地被保留(圖3)。樣區(qū)A、B和C點(diǎn)云平均高程分別降低了0.25、0.56和0.44 m(9月)。

表1 4種濾波插值方法所得DEM精度統(tǒng)計

深色點(diǎn)云為濾波前，淺色點(diǎn)云為濾波。 The dark point cloud is the one before filtering, and the light point cloud is the one after filtering圖3 濾波前后切溝橫剖面點(diǎn)云分布圖Fig.3 Distribution of point cloud of gully profiles before and after filtering

3.2 不同方法生成DEM對比

首先，9月所測點(diǎn)云經(jīng)過濾波之后，點(diǎn)云數(shù)目大幅度下降。12月所測點(diǎn)云經(jīng)過濾波之后，點(diǎn)云數(shù)目也有所下降，但幅度不及9月。尤其樣區(qū)A，可能是由于其以草本為主，且植被蓋度較低?？死锝鹋cTIN法所得點(diǎn)云高程均值幾乎相同，但濾波之后高程有所降低：9月所測點(diǎn)云大概0.4 m，12月所測點(diǎn)云在0.2 m左右。TIN法與克里金法交叉驗證的R2都在0.9以上，且克里金和TIN法中誤差在0.30～0.50 m之間，濾波后插值中誤差降低至0.20～0.10 m(表1)。TIN法中誤差普遍比克里金法要小，說明TIN法所得高程更接近于實測點(diǎn)云，地形地貌刻畫精度較高。

由圖4可見， b比a插值結(jié)果更平滑，粗糙度降低，但是c，d粗糙度差異不大。4種方法中溝沿線都十分清晰，a保留的地貌和植被細(xì)節(jié)更多。a、b中的高大植被在c、d中已被去除。a3中2次測得剖面相比，中上部基本吻合，下部12月比9月明顯更低，并且12月下部與中上部更為平滑，更符合實地特征。這是由于溝底的喬灌在12月落葉枯萎造成的。b3比a3相比，其不能去除9月植被且平滑了溝頭坡度。經(jīng)過濾波之后，c3，d3分別比a3、b3更為平滑，并且2次測得的結(jié)果基本一致，但是d3溝頭部分仍有明顯差異。

a為TIN法，b為克里金法，c為濾波-TIN法，d為濾波-克里金法；1為三維可視化疊加山體陰影效果圖，2為坡度圖，3為對應(yīng)三維可視化疊加山體陰影效果圖中的白色剖線。a refers to TIN, b refers to Kriging, c refers to filtering-TIN, and d refers to filtering-Kriging. 1 refers to 3D model with hill shade effect, 2 refers to slope map, and 3 refers to profiles of white lines in corresponding 3D model.圖4 4種方法生成的DEM比較(以樣區(qū)B為例)Fig.4 Comparison of DEMs by four processes (Taking site B as an sample)

在2次三維激光掃描之間，研究區(qū)內(nèi)沒有發(fā)生強(qiáng)烈的侵蝕性降雨，不存在劇烈侵蝕的情況。DEM中切溝形態(tài)參數(shù)變化主要是由植被變化引起的。切溝形態(tài)參數(shù)統(tǒng)計表明(表2)，克里金法所得的12月平均溝長、溝寬和溝深都更大，相差0.05～0.35 m。而TIN法所得12月平均溝長和溝寬比9月分別減小0.72和0.26 m，但溝深更大。濾波-TIN法所得12月平均溝長、溝寬和溝深比9月分別增大0.03、0.25和0.10 m，說明具有去除植被的效果。由于經(jīng)過濾波后，12月保留更多地貌特征點(diǎn)，即非生長季比生長季更接近于真實地貌。濾波-克里金法所得12月份的溝長和溝深分別比9月份小0.36和0.34 m，說明克里金法的平滑作用使得插值結(jié)果稍偏離實際地貌形態(tài)。

表2 4種處理方法提取的切溝形態(tài)參數(shù)比較

4 討論

4.1 植被季相變化對濾波效果的影響

點(diǎn)云平均高程和標(biāo)準(zhǔn)差等都表明12月點(diǎn)云高程總體上低于9月，高程分布更為集中。點(diǎn)云高程高值部分在9月數(shù)量比12月多，而中低值9月比12月數(shù)量要少(圖2)。這是由于非生長季植被枯萎凋零，覆蓋度降低，植被平均高度降低，使激光腳點(diǎn)落在地面的概率增大。由于三維激光掃描儀無法測量到植株的最高點(diǎn)，所以點(diǎn)云中的植被點(diǎn)總小于最大植株高度[29]，即不能精確反映植被高度，加之抽樣株高不能完全反映植被總體情況，所以9月與12月點(diǎn)云平均高程差(0.10 m)，略高于抽樣平均株高差(0.05 m)。總體上看，植被變化與點(diǎn)云變化規(guī)律相符合，表明植被變化對于點(diǎn)云具有一定影響。

對濾波前后的點(diǎn)云和DEM進(jìn)行了對比分析，結(jié)果表明坡度濾波法對于去除植被具有一定效果。在植被濾波完全去除植被的理想狀態(tài)下，不同季相所得濾波結(jié)果應(yīng)完全一致。但本文9月與12月植被濾波后的點(diǎn)云并不一致，樣區(qū)A 2次點(diǎn)云平均高程相差在0.11 m左右，樣區(qū)B和C相差了0.3 m左右。所以即使濾波之后，植被季相變化對點(diǎn)云高程的影響仍然存在，進(jìn)而影響DEM精度。

4.2 插值方法對DEM精度的影響

前人比較克里金法與最鄰近法、反距離權(quán)重法等將點(diǎn)云插值生成DEM的精度，結(jié)果表明克里金法效果較好[28]。本研究將克里金法和TIN法進(jìn)行比較，發(fā)現(xiàn)TIN法比克里金法R2更高，σ更小，且σ均符合國家測繪局1∶500 DEM一級精度要求(坡度25°以上，柵格分辨率為0.5 m的，一級精度σ<0.70，坡度6°～25°，一級精度σ<0.50)。TIN法生成的DEM與點(diǎn)云擬合效果更好，精度更高，對濾波前后的地形變化更為敏感。而克里金法有一定的平滑效果，一方面有使平滑高大植被的效果，使其高程一定程度地降低(9月，σ=0.46)，但另一方面，也會出現(xiàn)“牛眼”現(xiàn)象[28]，造成與實際地貌的偏離。另外，未經(jīng)濾波處理直接生成DEM時，克里金法比TIN法得到的平均溝長溝寬更大，但平均溝深更淺，而經(jīng)過濾波后，克里金法比TIN法生成的溝長溝寬和溝深更小?？赡苁怯捎诳死锝鸱ú逯禃r，溝沿線附近、溝底的部分點(diǎn)云插值權(quán)重過低，導(dǎo)致溝沿線附近、溝底略緩，出現(xiàn)DEM平滑現(xiàn)象。

總之，TIN法生成的DEM溝長溝寬等參數(shù)更接近于真實的切溝地貌形態(tài)。經(jīng)過濾波之后比未濾波生成DEM的切溝形態(tài)參數(shù)有所增大，說明具有一定地植被去除效果。

5 結(jié)論

1)植被季相變化對于三維激光掃描具有一定影響，生長季測量的點(diǎn)云高程均值高于非生長季0.25 m，從而進(jìn)一步影響到DEM的精度和地形地貌刻畫。

2)基于濾波-TIN法生成的DEM平均溝長、溝寬和溝深，2次相比，分別增大0.03、0.25和0.10 m，結(jié)合2次實地抽樣的植物株高，說明坡度濾波能一定程度上去除植被。但濾波后9月份(生長季)點(diǎn)云高程均值仍然比12月份(非生長季)高出0.15 m，兩者并不完全一致，表明植被季相變化對DEM所造成的影響無法完全去除。

3)與克里金法相比，TIN法生成的DEM精度較高。前者存在平滑現(xiàn)象，植被生長季時能粗略地去除高大植被，而在非生長季時平滑溝道邊緣。TIN法比克里金法的插值結(jié)果更接近實際地貌。

黃土高原丘陵溝壑區(qū)植被季相變化對DEM精度具有影響，且僅通過植被濾波算法很難完全去除。選擇植被非生長季時進(jìn)行測繪，并使用適合的濾波插值方法，才能有效地提高DEM的精度。