許昌建,劉迎春,左麗君,李建更,張 婷,韓路萌,方 宇,張 尹,王 天
(1.國家林業(yè)和草原局調(diào)查規(guī)劃設計院,北京 100714;2.北京工業(yè)大學 信息學部,北京 100124;3.中國科學院空天信息研究院,北京 100094)
森林作為陸地生態(tài)系統(tǒng)的主體,其生物量占全球陸地生態(tài)系統(tǒng)生物量的80%[1],在減緩全球氣候變化影響、維護生物多樣性和防治水土流失等方面起到至關(guān)重要的作用[2]。森林地上生物量(Above-Ground Biomass,AGB)是評價森林生產(chǎn)力和固碳速率的關(guān)鍵參數(shù),因此,快速準確地估測森林地上生物量對于量化碳儲量和了解全球碳循環(huán)是十分重要的[3]。目前,各國通用的森林地上生物量估測方法主要基于大量樣地調(diào)查數(shù)據(jù)和線性回歸模型,雖然估測精度高,但需要耗費較大的人力、物力和財力[4]。隨著遙感技術(shù)的不斷發(fā)展,激光雷達(Light Detection and Ranging,LiDAR)能夠穿透森林植被冠層,獲取森林垂直結(jié)構(gòu)信息,進而獲得較高精度森林地上生物量[5-6]。聯(lián)合星載激光雷達、機載激光雷達和地面樣地數(shù)據(jù)的北美洲寒帶針葉林地上生物量,估算結(jié)果與森林資源調(diào)查方法相差1.9%[7]。Muss等[8]利用小光斑激光雷達生成仿真大光斑激光雷達波形,對森林地上生物量進行估算,研究發(fā)現(xiàn),仿真波形法能夠構(gòu)建與生物量高度相關(guān)的模型。
在生物量估測模型的構(gòu)建方法上,大多數(shù)模型基于參數(shù)之間的線性關(guān)系構(gòu)建。如,多元線性回歸、偏最小二乘、數(shù)量回歸等方法。由于遙感因子之間往往具有較強的非線性關(guān)系,導致這些方法的擬合效果不佳[9]。人工神經(jīng)網(wǎng)絡是一種模擬人腦神經(jīng)元工作方式的建模方法,能夠較好地擬合非線性關(guān)系,并且具有一定的自組織與自學習能力[10]。近年來,人工神經(jīng)網(wǎng)絡在林業(yè)上得到了愈發(fā)廣泛的應用。例如,基于單隱含層反向傳播(Backpropagation,BP)神經(jīng)網(wǎng)絡構(gòu)建的遙感影像-森林地上生物量估測模型,在決定系數(shù)、均方根誤差以及預測精度等方面均明顯優(yōu)于多元線性回歸模型[11];以2個隱含層的多層感知器和多元線性回歸模型分別估算長白落葉松人工林地上生物量,證實了多層感知器的優(yōu)勢[12]。
盡管人工神經(jīng)網(wǎng)絡已用于森林地上生物量估測,然而目前的人工神經(jīng)網(wǎng)絡大多是淺層神經(jīng)網(wǎng)絡,即隱含層數(shù)較少。相關(guān)研究表明,與淺層神經(jīng)網(wǎng)絡相比,基于多隱含層的人工神經(jīng)網(wǎng)絡具有更強的表示能力,可挖掘出數(shù)據(jù)中更多的有用信息[13]。但更多的層數(shù)會帶來更多的模型參數(shù),需要更多的樣本訓練模型,而樣地調(diào)查數(shù)據(jù)往往不能滿足訓練樣本的數(shù)量要求。為解決訓練樣本不足的問題,本文利用機載激光雷達數(shù)據(jù)與樣地調(diào)查數(shù)據(jù)建立回歸關(guān)系,得到生物量樣本數(shù)據(jù),擴充了訓練樣本,用于多層感知器的訓練。本研究以河北省張家口市為研究區(qū)域,利用小光斑激光雷達仿真大光斑激光雷達波形,進而提取波形參數(shù),用于估算森林地上生物量,并比較了多元線性回歸模型和多層感知器模型的估測精度和擬合優(yōu)度。
研究區(qū)位于河北省西北部的張家口市,地理位置為39°30′~42°10′N,113°50′~116°30′E,屬于溫帶大陸性季風氣候區(qū),四季分明,年平均氣溫6.2℃,年均降水量400mm。張家口市南北長289.2km,東西寬216.2km,下轄6區(qū)10縣,土地面積3.68萬km2,地勢西北高、東南低,橫貫中部的陰山山脈將張家口市劃分為壩上、壩下兩部分。壩上高原區(qū)海拔1 300~1 600m,地勢平坦,草原遼闊,是典型的波狀高原景觀。壩下是華北平原和蒙古高原的過渡帶,海拔500~1 200m,地形復雜,丘陵、河谷與盆地相間分布。研究區(qū)森林覆蓋率61%,主要樹種有落葉松(Larixgmelinii)、白樺(BetulaplatyphyllaSuk.)、油松(Pinustabulaeformis)、側(cè)柏(Platycladusorientalis)、山楊(Populus
davidiana)、蒙古櫟(Quercusmongolica)等。
2.1.1機載小光斑激光雷達數(shù)據(jù)
機載小光斑激光雷達數(shù)據(jù)的采集時間為2018年9月15日—30日,天氣晴朗。采用大疆經(jīng)緯M600 Pro無人旋翼機搭載掃描鷹HS-600超輕小低空激光雷達系統(tǒng),獲取了9個1km×1km的小光斑激光雷達數(shù)據(jù)。因各個測區(qū)的面積較小且相對獨立,調(diào)查時在每個測區(qū)的一級測量控制點上架設了一臺地面基站。根據(jù)機載LiDAR航攝技術(shù)要求、測區(qū)范圍、成圖要求及HS-600航攝儀性能,制定了設備航飛參數(shù)(表1)。
表1 航飛參數(shù)
小光斑激光雷達分類。首先,采用噪聲點濾波將明顯低于地面的點或點群(低點)和明顯高于地表目標的點,以及移動地物點歸為噪點,最先分離出來。之后,使用LiDAR360軟件對點云進行自動分類,再通過人工編輯,糾正自動分類時錯分的點。最終分類后的激光雷達數(shù)據(jù)包含地面點、植被點、噪點3類(圖1)。
圖1 小光斑激光雷達分類流程
2.1.2樣地調(diào)查數(shù)據(jù)
2018年9—10月在河北省張家口市調(diào)查了16個20m×30m森林樣地。地面調(diào)查點的確定方法是:先在0.8m分辨率衛(wèi)星影像和林地一張圖上分樹種、林齡組選擇調(diào)查區(qū)域;之后,到調(diào)查區(qū)域選擇備用調(diào)查點,并在衛(wèi)星影像標記樣地起始點和矩形樣地范圍,備用調(diào)查點一般選擇有明確參考地物的地點;最后,調(diào)查人員從備選調(diào)查點選取最終調(diào)查樣地,并提供準確的樣地坐標。由于調(diào)查樣地采用標志物定位,定位誤差低于0.5m。由于無人機激光雷達數(shù)據(jù)位置精度高,且能夠清晰分辨每一株樹,調(diào)查樣地與激光雷達數(shù)據(jù)間的匹配精度高于0.5m,能夠保證二者位置的一致性。樣地位置如圖2所示。
圖2 研究區(qū)和樣地位置分布圖
這16個樣地按森林類型分為:針葉林12個和闊葉林4個;按林齡組分為:幼齡林3個、中齡林7個和成熟林6個。調(diào)查每個樣地的經(jīng)緯度、海拔、坡度,通過每木檢尺獲得樣地林木的平均胸徑,選取標準木,測量標準木樹高,采用已發(fā)表的異速生長方程(表2)計算標準木生物量,并換算成樣地生物量。用LAI-2200沿樣地對角線測定森林冠層葉面積指數(shù)(LAI)。
表2 主要樹種生物量異速生長方程
仿真波形按照以下4個步驟生成:
1) 提取直徑20m內(nèi)的小光斑激光雷達,在垂直方向上以0.15m間隔對小光斑激光雷達進行高度劃分。
2) 對小光斑激光雷達的反射強度賦予權(quán)重。一般來說,大光斑激光雷達的光斑能量分布服從二維正態(tài)分布,即光斑中心能量最高,向外逐漸衰減,光斑邊緣的能量為中心能量的1/e2[19]。點云反射強度的賦值公式見式(1):
(1)
式中:Iw,i是第i個點加權(quán)強度;Ii是該點的反射強度;xi和yi是該點的平面坐標,x0和y0是光斑中心位置的平面坐標;R是光斑半徑。
3) 對每個高度區(qū)間內(nèi)的所有點的加權(quán)強度進行累加,得到初步的仿真波形。
4) 根據(jù)仿真波形能量的最大值和最小值,對仿真波形進行歸一化處理,使其相對強度分布在(0,1)區(qū)間,得到歸一化的仿真波形(圖3),以降低點云密度對仿真波形的影響。
圖3 歸一化后的仿真大光斑激光雷達波形
提取仿真波形參數(shù)。本文提取了波形高度百分位數(shù)Qx(x=10,20,25,30,40,50,60,70,75,80,90,95)、波形內(nèi)平均高Avg、波形內(nèi)最大高Max、波形變異系數(shù)CV和波形標準差SD共16個波形參數(shù)(表3)。
表3 波形參數(shù)定義Tab.3 Definition of waveform parameters
深度學習模型所需的訓練樣本較大,樣地調(diào)查數(shù)據(jù)遠不能滿足訓練樣本的數(shù)量要求。為了擴充訓練樣本,本文以葉面積指數(shù)與冠層高度的乘積為自變量,建立與實測樣地生物量關(guān)系方程;進而以此方程推算小光斑激光雷達覆蓋區(qū)域的森林生物量,共提取1 333組生物量樣本,且兩兩樣本之間的距離大于50m。將樣本數(shù)據(jù)按照4∶1的比例,隨機地將樣本劃分為訓練集和測試集,其中訓練集1 068組數(shù)據(jù),測試集265組數(shù)據(jù)。
2.3.1葉面積指數(shù)
葉面積指數(shù)(Leaf Area Index,LAI)是單位地表面積上所有葉片表面積的一半,是表征植被冠層結(jié)構(gòu)的基本參量之一[20]。本文采用小光斑激光雷達結(jié)合樣地實測葉面積指數(shù)計算區(qū)域森林葉面積指數(shù)。
葉面積指數(shù)的計算[21]如式(2)所示:
(2)
式中:ang為平均掃描角;GF為間隙率(Gap Fraction);k為消光系數(shù),消光系數(shù)由比爾-蘭伯特方程得到。已有研究表明,大部分森林的消光系數(shù)可以為0.5[22]。
平均掃描角(ang)的計算如式(3)所示:
(3)
式中:n是小光斑激光雷達總點數(shù),anglei是第i個點的掃描角度。
間隙率(GF)的計算如式(4)所示:
(4)
式中:nground是高度值低于高度閾值的地面點數(shù)。
2.3.2冠層高度模型
冠層高度模型(Canopy Height Model,CHM)由數(shù)字表面模型(Digital Surface Model,DSM)與數(shù)字高程模型(Digital Elevation Model,DEM)做差得到。本文所使用的DSM和DEM均由分類后的點云數(shù)據(jù),在LiDAR360軟件中,通過TIN插值法生成。插值時,采用逐點插入法構(gòu)建Delaunay三角網(wǎng),從最近的臨近點組成的多個三角形共同形成的表面上提取柵格單元值。
2.4.1多元線性回歸模型
在森林地上生物量估測的相關(guān)研究中,多元線性回歸(Multiple Linear Regression,MLR)模型的應用十分廣泛[23],一般以森林地上生物量實測值作為因變量,以遙感數(shù)據(jù)參數(shù)作為自變量,通過一次多項式構(gòu)建模型,擬合自變量與因變量之間的關(guān)系(式(5)):
Y=βX+ε
(5)
式中:Y表示地上生物量;X表示自變量;β為自變量的參數(shù);ε為誤差項。
2.4.2多層感知器模型
多層感知器(Multi-Layer Perceptron,MLP),也稱作深度前饋網(wǎng)絡,是一種典型的深度學習模型[24]。多層感知器的結(jié)構(gòu)一般最左側(cè)為輸入層,中間為隱含層,最右側(cè)為輸出層(圖4)。在多層感知器中,前一層的每個節(jié)點與相鄰后一層的每個節(jié)點都是連接的,前一層接收的輸入通過矩陣運算和激活函數(shù)后輸出,作為相鄰后一層的輸入,并逐層運算至輸出層。
圖4 多層感知器結(jié)構(gòu)示意圖
(6)
式中:激活函數(shù)σ(·)為修正線性單元(Rectified Linear Unit,ReLU),其表達式為:
σ(x)=max{0,x}
(7)
本文構(gòu)建的多層感知器包含5個隱含層,每個隱含層均包含10個神經(jīng)元,輸入層神經(jīng)元個數(shù)與輸入變量個數(shù)相等,輸出層有1個神經(jīng)元。使用803組訓練樣本,學習率0.05,對模型訓練1 000輪,并使用265組驗證樣本,對模型的擬合優(yōu)度和估測精度進行評價。
2.4.3模型精度評價指標
本文從擬合優(yōu)度和估測精度兩方面對模型進行評價。評價擬合優(yōu)度的指標為決定系數(shù)(R2)和調(diào)整決定系數(shù)(Adj.R2),公式如下:
(8)
(9)
評價估測精度的指標為均方根誤差(RMSE)、相對均方根誤差(RMSEr)和平均絕對誤差(MAE)作為,公式如下:
(10)
(11)
(12)
式中:n為樣本總數(shù);ymax為樣地調(diào)查生物量最大值;ymin為樣地調(diào)查生物量最小值。
由小光斑激光雷達計算的葉面積指數(shù)(LAI),與樣地調(diào)查葉面積指數(shù)(LAI′)存在較強的線性關(guān)系:LAI′=0.893LAI+0.192,R2達到0.94(圖5)。以上述關(guān)系方程對葉面積指數(shù)進行修正,得到最終的區(qū)域葉面積指數(shù)產(chǎn)品,并與0.2m分辨率正射影像對比(圖6)。葉面積指數(shù)產(chǎn)品的空間分布與影像中植被分布趨勢一致:植被茂密區(qū)域的葉面積指數(shù)明顯高于植被稀疏區(qū)域,無植被區(qū)域的葉面積指數(shù)為0(顏色越深表示LAI值越高,0值區(qū)域為白色)。
圖5 小光斑激光雷達計算葉面積指數(shù)與樣地測量葉面積指數(shù)對比
圖6 測區(qū)葉面積指數(shù)與正射影像對比
以校正后葉面積指數(shù)(LAI')與冠層高度(H)的乘積為自變量,與樣地生物量(B)建立回歸關(guān)系(式(13)):
B=6.866×H×LAI'+2.654
(13)
式中:H從CHM中提取,R2達到0.87(圖7)。
圖7 樣地地上生物量與校正后葉面積指數(shù)(LAI')與冠層高度(H)乘積的關(guān)系
從波形參數(shù)間相關(guān)關(guān)系(圖8)可以看出,高度分位數(shù)Qx之間的相關(guān)性較高,且分位數(shù)越相近,相關(guān)系數(shù)越大;而波形平均高Avg、波形變異系數(shù)CV和波形標準差SD與其他各項參數(shù)之間的相關(guān)系數(shù)較低,相關(guān)系數(shù)均小于0.4。
由于高度分為數(shù)Qx之間的相關(guān)性較高,本文分別以Q25,Q50,Q75,Q95,Max,AVG,CV,SD單個參數(shù)進行建模,評價了單個參數(shù)對模型估測效果的影響;同時對上述參數(shù)以一定方式進行組合,觀察了不同參數(shù)組合的模型估測效果;最后用全部16個參數(shù)進行建模,從而比較多層感知器和多元線性回歸模型對參數(shù)信息的挖掘能力。
首先建立單變量線性回歸模型,觀察單個波形參數(shù)與森林地上生物量之間的關(guān)系,模型的輸入?yún)?shù)有高度百分位數(shù)Qx(x=10,25,50,75,95)、波形最大高Max、波形平均高Avg、波形變異系數(shù)CV和波形標準差SD。
注紅色表示相關(guān)系數(shù)絕對值大于0.6,黃色表示相關(guān)系數(shù)絕對值在0.4與0.6之間,綠色表示相關(guān)系數(shù)絕對值小于0.4,相關(guān)系數(shù)絕對值的數(shù)值在矩陣下半三角中列出。
表4 單變量線性回歸模型參數(shù)及擬合效果
當采用單個波形參數(shù)作為輸入時(表4),選擇波形最大高Max得到的估測效果最好。以高度百分位數(shù)Qx作為輸入變量時,x值越大,模型的估測效果越好。波形平均高Avg、波形變異系數(shù)CV和波形標準差SD與森林地上生物量的相關(guān)性較差,R2只有0.01。
在構(gòu)建多元線性回歸模型時,選取不少于3個波形參數(shù),且盡量避免關(guān)系較弱的參數(shù),建立多元線性回歸模型(表5)。
各多元線性回歸模型的評估指標見表6豎線左側(cè)??梢钥闯?當參與建模的變量增多時,模型的擬合效果得到提升,選擇全部16個變量進行建模時,模型的決定系數(shù)與調(diào)整決定系數(shù)均為0.60。
表6中豎線右側(cè)列出了各個多層感知器模型的擬合優(yōu)度與估測精度,為了使比較結(jié)果更直觀,估測結(jié)果更優(yōu)的指標用黑體表示。
表5 13種波形參數(shù)組合得到的多元線性回歸模型
表6 基于13種波形參數(shù)組合的多元線性回歸模型和多層感知器模型估算森林地上生物量效果比較
從訓練集與驗證集的均方誤差(圖9)可以看出,在前80輪訓練中,模型在訓練集和驗證集上的均方誤差迅速下降。在第80到200輪訓練中,模型在訓練集上的均方誤差緩慢下降,但在驗證集上的均方誤差并不穩(wěn)定。在200輪訓練之后,模型在訓練集上的均方誤差緩慢下降,在驗證集上的均方誤差逐漸穩(wěn)定,不再隨訓練輪數(shù)的增加而下降。訓練完成后,模型在訓練集與驗證集上的均方誤差相差不大,沒有出現(xiàn)明顯的過擬合。
從表6可以看出,使用相同的變量組合進行建模時,除參數(shù)組合2中多元線性回歸模型的RMSE和RMSEr優(yōu)于多層感知器模型外,其余各項評價指標多層感知器模型均優(yōu)于多元線性回歸模型。
圖9 多層感知器(參數(shù)組合12)估算地上生物量的均方根誤差隨訓練輪數(shù)變化趨勢
當參與建模的參數(shù)增加時,多元線性回歸模型和多層感知器模型的估測效果均得到提升。但隨著輸入?yún)?shù)的增多,多元線性回歸模型的估測效果得到的提升越來越有限,參數(shù)組合12比參數(shù)組合13的建模參數(shù)多7個,但多元線性回歸模型估測結(jié)果的決定系數(shù)相同,其他評價指標也相近。從參數(shù)組合9,10,11,12和13看出,隨著建模參數(shù)的增加,多層感知器模型的估測效果能夠得到持續(xù)的提升,說明多層感知器模型的擬合能力更強,能夠從新增的參數(shù)中挖掘出更多的信息。
以參數(shù)組合12為例,對多層感知器模型和多元線性回歸模型估測16個樣地地上生物量的結(jié)果進行對比。對多元線性回歸模型而言,估測的生物量均值為48.46t/hm2,比實測均值51.52t/hm2低5.93%(圖10(a))。6個樣地的預測值比調(diào)查值高,10個樣地的預測值比調(diào)查低,偏差幅度為-34.96~23.28t/hm2(圖10(b))。對多層感知器模型而言,估測的生物量均值為49.43t/hm2,比實測均值低4.06%(圖10(c))。6個樣地的預測生物量比調(diào)查值高,10個樣地的預測生物量比調(diào)查值低,偏差幅度為-19.09~20.19t/hm2(圖10(d))。因此,多層感知器的預測結(jié)果比多元線性回歸模型更接近實測值。
本文利用樣地數(shù)據(jù)和小光斑激光雷達數(shù)據(jù)擴充生物量樣本,進而以小光斑激光雷達仿真大光斑激光雷達。雖然這會引入一些誤差,但由于很多區(qū)域無法同時獲取樣地數(shù)據(jù)和大光斑激光雷達數(shù)據(jù),可以將機載小光斑激光雷達數(shù)據(jù)作為樣地數(shù)據(jù)與大光斑激光雷達數(shù)據(jù)的中間尺度數(shù)據(jù)。本文以仿真大光斑激光雷達波形參數(shù)作為模型輸入?yún)?shù),對比了不同參數(shù)組合的多元線性回歸模型和多層感知器模型的森林地上生物量估測效果。實驗發(fā)現(xiàn),以單個參數(shù)構(gòu)建模型時,波形最大高Max的估測效果最好。以高度百分位數(shù)參數(shù)Qx作為模型輸入?yún)?shù)時,x的值越大,模型的估測效果越好。波形平均高Avg、變異系數(shù)CV和標準差SD作為參數(shù)時,模型的估測效果不好。
通過對比多元線性回歸模型和多層感知器模型發(fā)現(xiàn),當建模參數(shù)較少時,通過增加建模參數(shù),兩種模型的生物量估測效果均可得到提升。但參與建模的參數(shù)較多時,多元線性回歸模型已不能夠通過增加建模參數(shù)來明顯提升估測效果,而多層感知器模型仍然能夠從增加的參數(shù)中挖掘信息,提升模型的估測效果。當兩種模型的建模參數(shù)相同時,多層感知器模型的擬合優(yōu)度與估測精度顯著高于多元線性回歸模型。
本文工作也存在有待完善之處。樣本數(shù)據(jù)由樣地實測數(shù)據(jù)建立回歸關(guān)系而得到,這樣做雖然大大擴充了訓練模型的樣本,使訓練多層感知器模型成為可能,但由樣地實測數(shù)據(jù)到樣本數(shù)據(jù)的轉(zhuǎn)化中,可能引入新的誤差,并對模型的估測結(jié)果產(chǎn)生影響,這種影響有待進一步的研究和分析。
目前,一些學者聯(lián)合Lidar、可見光、高光譜、SAR等數(shù)據(jù),對森林地上生物量進行估測,發(fā)現(xiàn)能夠獲得更好的生物量估測效果,這些數(shù)據(jù)特征提取和多源數(shù)據(jù)融合方式,對本研究的建模參數(shù)提取具有一定指導意義。未來應嘗試更加多樣的參數(shù)組合,將不同數(shù)據(jù)源的遙感數(shù)據(jù)進行融合,進一步提高模型的生物量估測能力。