付安民,高顯連,吳發(fā)云,高金萍

(國家林業(yè)和草原局調(diào)查規(guī)劃設(shè)計院,北京 100714)

樹高是最為重要的森林定量遙感觀測參數(shù)之一,高空間分辨率、高精度的樹高遙感產(chǎn)品是進(jìn)一步開展準(zhǔn)確森林資源定量遙感調(diào)查和生態(tài)環(huán)境定量遙感評估的必要數(shù)據(jù)基礎(chǔ)[1]。目前，大范圍、高質(zhì)量的樹高遙感產(chǎn)品主要由航空激光雷達(dá)觀測數(shù)據(jù)處理得到[2],并且隨著近幾十年來激光雷達(dá)觀測和計算機(jī)技術(shù)的快速發(fā)展,航空激光雷達(dá)數(shù)據(jù)的采集和處理能力不斷得到提升[3]。激光點云數(shù)據(jù)已漸漸成為我國森林領(lǐng)域的重要數(shù)據(jù)源,服務(wù)于國家和區(qū)域性的森林資源調(diào)查、森林碳儲量估算、森林生態(tài)多樣性評估等工作[4-7]。

從激光點云到樹高的過程通常需要分別獲取林下地形(Digital Terrain Model,DTM)和林冠表面高度(Canopy Surface Model，CSM)產(chǎn)品,兩者的差值作為森林冠層高度產(chǎn)品(Canopy Height Model,CHM),因而測量樹高的能力受冠層表層和林下地形高度雙重精度限制,獲取高質(zhì)量的CHM是一項具有一定挑戰(zhàn)性的工作[8-9]。從觀測對象上看,森林觀測場景復(fù)雜多樣,不僅樹種組成多樣,而且林木高低起伏、相互遮擋、枝葉縱橫交錯,并且耦合了地形要素?！巴昝馈钡臄?shù)據(jù)采集不僅要高質(zhì)量捕捉到林下地形的起伏,同時還要準(zhǔn)確刻畫林冠大小和高低,才能保障高分辨率林下地形和林冠表層高度產(chǎn)品的準(zhǔn)確性和精度。從激光雷達(dá)載荷角度上看,不同載荷或同一載荷不同參數(shù)設(shè)置,觀測的準(zhǔn)確性和精度指標(biāo)均會受到影響,會引起觀測森林樹高的不同。由于森林是非剛體的“軟”目標(biāo),激光脈沖波束對樹冠具有一定的穿透性,直至照射至地面；激光探測器記錄透射激光脈沖的回波波形,依據(jù)波形振幅分解的點云結(jié)果進(jìn)行地面點分類濾波處理,同步獲取林木冠層和林下地形高度。對兩者的測量能力取決于激光雷達(dá)載荷波形探測器的敏感性及其波形分解算法能力、激光脈沖能量、觀測幾何條件,以及激光照射地面光斑大小等綜合因素[10]。

林冠高度產(chǎn)品CHM并非真正意義上的地面調(diào)查樹高,其物理意義需要與地面實際樹高進(jìn)行標(biāo)校釋義[8]。之后,通過進(jìn)一步的地面樹高歸一化處理,以達(dá)到聯(lián)系林木蓄積、生物量模型,實現(xiàn)森林資源調(diào)查的應(yīng)用目的。同時,受到山區(qū)林地復(fù)雜的測量條件,激光雷達(dá)觀測條件差異和處理技術(shù)方式等綜合因素的影響,CHM計量的準(zhǔn)確性及其精度指標(biāo)在測區(qū)空間分布上具有不確定性,有必要建立地面林木參照對象,對產(chǎn)品質(zhì)量進(jìn)行系統(tǒng)的檢驗分析,以促進(jìn)激光點云采集數(shù)據(jù)指標(biāo)和數(shù)據(jù)處理方法的優(yōu)化改進(jìn),為獲取高空間分辨率、高精度的樹高遙感產(chǎn)品提供保障。

本文以航空激光雷達(dá)森林資源調(diào)查數(shù)據(jù)工程質(zhì)量控制為目的,設(shè)計了航空激光雷達(dá)CHM高度產(chǎn)品評估規(guī)程方法,旨在標(biāo)稱CHM的“真實”物理意義,分析CHM高度產(chǎn)品的準(zhǔn)確性和精度指標(biāo)水平,檢驗樣地數(shù)據(jù)與航空激光點云數(shù)據(jù)的一致性。

1 研究區(qū)和數(shù)據(jù)

1.1 研究區(qū)

2018年9月1日—9月27日,在吉林省東部東北虎豹國家公園部分地區(qū)(寧安、汪清和穆棱境內(nèi)),地理坐標(biāo)(43°19′～44°10′N,129°4′～130°21′E),東西長108km,南北寬105km范圍內(nèi),開展了總面積 5 600km2的航空飛行遙感森林調(diào)查綜合實驗,如圖1所示。實驗以復(fù)雜山區(qū)地形條件下航空森林資源調(diào)查為主要目的,同步開展地面森林資源數(shù)據(jù)調(diào)查。

研究區(qū)屬于溫帶大陸性季風(fēng)氣候,氣候特征表現(xiàn)為春季多風(fēng)少雨干旱,夏季炎熱短促,秋季冷涼降溫迅速,冬季寒冷漫長。年降水量變化在450～750mm,降水集中在5—9月。研究區(qū)以森林生態(tài)系統(tǒng)為主體,主要是溫帶針闊葉混交林,其中,次生林分布廣泛,原生性紅松闊葉混交林僅呈零星分布[11]。

注:方框為312幅激光雷達(dá)文件數(shù)據(jù),紅點表示地面森林調(diào)查樣地。

1.2 數(shù)據(jù)

1.2.1地面林木調(diào)查數(shù)據(jù)

地面森林調(diào)查結(jié)合歷史森林資源規(guī)劃調(diào)查數(shù)據(jù),進(jìn)行了10個類型:7個純林(云杉(PiceaasperataMast.)、冷杉(Abiesfabri(Mast.)Craib)、落葉松(Larixgmelinii(Rupr.)Kuzen.)、樺木(Betula)、楊樹(PopulusL.)、椴樹(TiliatuanSzyszyl.)、櫟類(QuercusL.))和3個混交林(針葉混、針闊混、闊葉混)的抽樣設(shè)計,共計完成439個有效樣地調(diào)查工作,分布情況如圖1所示。樣地為圓形,直徑大小為30m。樣地內(nèi)胸徑大于5cm的所有樣木均劃為檢尺范圍,所有檢尺樣木記錄樹種、胸徑、樹高,部分樣木記錄冠幅、枝下高參數(shù)。采用微波、激光測距儀進(jìn)行樹高和冠幅的量測記錄,采用胸徑尺進(jìn)行胸徑(樹干高度1.3m處)測量；相關(guān)林木指標(biāo)如表1所示。

樣地中心和單樣木水平位置,首先協(xié)同衛(wèi)星定位服務(wù)連續(xù)運行參照站(Continuously Operating Reference Station,CORS)利用差分技術(shù)實現(xiàn)雙點絕對定位；然后再架設(shè)全站儀利用后方交會測量方法實現(xiàn)單木精密水平定位,名義上樣地中心和單樣木水平相對精度可優(yōu)于0.5m,絕對精度優(yōu)于2m；相關(guān)林木定位指標(biāo)如表1所示。

表1 樣地樣木野外測量指標(biāo)

1.2.2航空激光雷達(dá)數(shù)據(jù)

航空激光雷達(dá)數(shù)據(jù)采集利用機(jī)載激光雷達(dá)掃描系統(tǒng)聯(lián)合IMU/DGPS輔助航空測量技術(shù)開展,具體選用塞斯納208B飛機(jī)搭載RIEGL-VQ-1560i激光雷達(dá)載荷,集成慣性導(dǎo)航(Inertial Measurement Unit,IMU)和全球定位導(dǎo)航系統(tǒng)(Global Navigation Satellite System,GNSS),并聯(lián)合CORS站點同步觀測,具體設(shè)備及其飛行指標(biāo)如表2所示。

實際飛行數(shù)據(jù)達(dá)到激光脈沖密度優(yōu)于10pulse/m2。飛行任務(wù)規(guī)劃為南、北兩個測區(qū),飛行高度分別為2 100m和2 300m,飛行速度均為260km/h,飛行航帶重疊度優(yōu)于25%,激光測距儀載荷的參數(shù)設(shè)置保持一致。南、北兩測區(qū)分別設(shè)置了構(gòu)架航帶,構(gòu)架線航測區(qū)域布設(shè)校準(zhǔn)面,對飛行航線進(jìn)行整體校準(zhǔn)。

表2 航空激光雷達(dá)數(shù)據(jù)采集設(shè)備及其指標(biāo)

2 原理和技術(shù)方法

2.1 樣地樣木數(shù)據(jù)處理

樣地樣木數(shù)據(jù)處理過程主要分為4個步驟:1)GNSS數(shù)據(jù)進(jìn)行差分定位處理；2)協(xié)同全站儀后方交會測量數(shù)據(jù)和GNSS差分結(jié)果的單木定位和樣地中心定位；3)樣地林木統(tǒng)計；4)建立森林資源樣地和樣木數(shù)據(jù)庫。其中GNSS絕對定位精度優(yōu)于2m；全站儀單木相對定位精度優(yōu)于0.5m。平面坐標(biāo)為CGCS 2000坐標(biāo)系,高程采用1985國家高程基準(zhǔn),投影采用高斯克里格,3度分帶,中央子午線129°。

其中樣木的水平定位需要特別注意,開展林區(qū)定位測量,GNSS接收機(jī)接收衛(wèi)星信號易受到地形、林木遮擋,經(jīng)常發(fā)生衛(wèi)星信號多路徑和無法直接接收情況,定位精度受到影響,如圖2所示。通常需要延長GNSS接收機(jī)的觀測時間達(dá)到15min以上,才能達(dá)到絕對定位精度優(yōu)于2m的需求。此外,事后GNSS差分軟件報告精度與實際位置精度也沒有相關(guān)性,無法通過軟件報告的精度水平,推斷實際的定位水平[12]。因而激光點云冠層高度模型數(shù)據(jù)CHM與地面樣木進(jìn)行高度對比時,需要進(jìn)行水平一致性檢驗,并對偏移明顯的樣地進(jìn)行修正,使得兩者的位置匹配精度優(yōu)于2m。

圖2 GNSS衛(wèi)星接收信號多路徑效應(yīng)和無法直接

2.2 激光雷達(dá)數(shù)據(jù)及DTM CHM產(chǎn)品

激光雷達(dá)數(shù)據(jù)處理主要分為4個步驟:1)利用POSPac軟件將GNSS和IMU慣導(dǎo)數(shù)據(jù),聯(lián)合地面參照站數(shù)據(jù)進(jìn)行差分定位處理,逐航帶進(jìn)行定位定向精結(jié)算；2)利用RiPROCESS軟件將精解算姿態(tài)和定位數(shù)據(jù)賦給每個激光脈沖數(shù)據(jù),并根據(jù)架構(gòu)航線對各航帶激光雷達(dá)數(shù)據(jù)進(jìn)行航帶間拼接和整體平差；3)剔除明顯的噪聲點；4)區(qū)域分幅,對激光點云數(shù)據(jù)進(jìn)行濾波和分類,區(qū)分地面點和非地面點。東北虎豹國家公園機(jī)載激光雷達(dá)數(shù)據(jù)2,3級產(chǎn)品融合一體,共計312幅文件,總數(shù)據(jù)量2.7TB,如圖1和表3所示。

在激光雷達(dá)2,3級數(shù)據(jù)的基礎(chǔ)上,首先,利用地面分類點,構(gòu)建TIN網(wǎng)空間插值生成1m空間分辨率的DTM數(shù)據(jù)；其次,利用首次回波點云(First Return)數(shù)據(jù),根據(jù)1m網(wǎng)格高程最大值統(tǒng)計方法,以非空間插值方式,獲取CSM數(shù)據(jù)；最后,兩者差值獲取1m空間分辨率的CHM數(shù)據(jù)。

2.3 CHM質(zhì)量分析方法

CHM質(zhì)量分析方法主要有6個步驟。1)確定樣木空間尺度大小,假定單木冠幅直徑不小于5m。在此假設(shè)條件下,以樣木水平定位坐標(biāo)為中心,建立2m半徑的圓形緩沖區(qū),用來分析樣木彼此的空間臨近關(guān)系,執(zhí)行空間疊加運算分析。2)優(yōu)化選擇地面參照對比樣木:為避免林木高低起伏,相互遮擋影響,運用空間拓?fù)浞治龊途植繕涓咦罡咴瓌t,逐個樣地選擇具備較好觀測天窗的上層林木,且確保周邊沒有其他林木遮擋。3)空間疊加分析:將地面樣木數(shù)據(jù)與CHM產(chǎn)品空間疊加,逐個樣木獲取緩沖區(qū)內(nèi)高度最大值CHMmax,賦予每株樣木；建立CHM樹高、地面調(diào)查樹高的雙樣本集。4)從樣木總體看CHM樹高:根據(jù)雙樣本散點圖和差值直方圖統(tǒng)計量,反復(fù)迭代去除粗差(均值μ及2倍標(biāo)準(zhǔn)差2σ為參變量)；利用優(yōu)化樣本總體確立通用型CHM樹高歸一化模型,對標(biāo)CHM高度“真實”物理意義,評估其準(zhǔn)確性。5)從分樹種組樣木總體看CHM樹高:根據(jù)雙樣本散點圖和差值直方圖統(tǒng)計量,查看各樹種組CHM高度“真實”物理意義和準(zhǔn)確性分異情況。6)評估CHM產(chǎn)品質(zhì)量:制定精度水平分級指標(biāo),統(tǒng)計抽樣參照對比樣木的CHM測量精度分級情況,指示CHM質(zhì)量狀況。

表3 激光點云數(shù)據(jù)產(chǎn)品分級表

圖3 CHM產(chǎn)品生產(chǎn)與精度評價流程圖

3 結(jié)果與分析

3.1 CHM樹高準(zhǔn)確性和精度指標(biāo)分析

從地面林木調(diào)查數(shù)據(jù)中選取上層林優(yōu)勢木作為“真實”數(shù)據(jù),用于CHM樹高準(zhǔn)確性和精度的對標(biāo)分析。首先,從地面每株林木定位精度優(yōu)于2m條件出發(fā),構(gòu)建每一株樣木的2m半徑圓形緩沖區(qū)；其次,利用空間分析工具確定局部優(yōu)勢的高大喬木,即挑選的每株林木在與其緩沖區(qū)交疊的所有鄰近林木中具有優(yōu)勢高度；剔除周邊低矮林木。圖4顯示了初步優(yōu)化后在中密度林和高密度林中的典型抽樣結(jié)果。最終,從439個樣地的全體25 467株林木中,挑選了參照樣木7 873株。

將確立的參照樣木按照緩沖區(qū)分析范圍提取最高值CHMmax作為機(jī)載激光雷達(dá)樹高,建立雙樣本對比數(shù)據(jù),圖5(a)和(b)顯示了兩者的散點圖和差值直方圖(參照樣木樹高減去CHMmax)。從圖5(a)散點圖上看,兩者存在明顯的線性關(guān)系,R2=0.68,同時中低林木觀測存在顯著的CHM樹高大于地面參照樣木的現(xiàn)象,大多是被臨近高大林木遮擋導(dǎo)致的,如圖4所示。從圖5(b)樹高差值直方圖上看,在-2.5～4m的區(qū)間范圍呈現(xiàn)明顯的高斯峰分布；同時還存有700余株的CHM樹高顯著大于地面參照樣木現(xiàn)象(5m以上)。綜上,CHM樹高和地面參照樣木雙樣本呈現(xiàn)顯著的線性相關(guān)性,且CHM樹高略低于地面參照樹高,均值約為1m左右；但是也存在一定數(shù)量的粗差樣本對。

注:黑色數(shù)字表示地面調(diào)查樹高,紅色數(shù)字表示地面樹高與CHMmax的差值,高亮圓圈表示受鄰近高大林木影響導(dǎo)致的粗差樣本。

圖5 地面參照樣木樹高與相應(yīng)2m半徑圓形緩沖區(qū)CHMmax的散點圖和差值直方圖(迭代去粗差前后對比)

去除粗差樣本,優(yōu)化上層林木確立“真實”數(shù)據(jù)樣本集,用于樹高物理意義的對標(biāo)分析。以μ±2σ為限值(樣本均值μ和標(biāo)準(zhǔn)差σ),經(jīng)4次去粗差迭代,去除粗差樣本。重新獲取的優(yōu)化上層林木樣本總計5 856株,如,圖5(c)和圖5(d)顯示了總體散點圖和樹高差值統(tǒng)計直方圖。兩者線性相關(guān)系數(shù)R2=0.97,樹高差均值μ為0.7m,標(biāo)準(zhǔn)差σ為±1.0m。表4顯示了各個樹種(組)的專題統(tǒng)計結(jié)果(1)吉林省林業(yè)廳.吉林省立木材積、出材率表.2016.。從各個樹種(組)的樹高差均值μ看,均呈現(xiàn)CHM樹高比上層林木樹高要低的現(xiàn)象,偏差范圍0.3～1.3m不等,其中偏差較大的為人工針葉林Ⅰ的1.0m、人工針葉林Ⅲ的1.3m；從各樹種(組)樹高差變量標(biāo)準(zhǔn)差上看,大小范圍±0.7～±1.2m不等。以上分析表明,CHM樹高與地面真實樹高存在系統(tǒng)偏差,各個樹種略有差異。

表4 CHM樹高和精度標(biāo)校分析表

3.2 CHM樹高產(chǎn)品精度評價

CHM樹高產(chǎn)品精度評價就是篩選適合的上層林木作為“真實”樹高,開展對比分析,定量化反映產(chǎn)品的精度水平及其空間分異情況。其中,地面調(diào)查林木的單株水平定位精度優(yōu)于2m；同時,地面樹高測量精度受林地復(fù)雜觀測條件限制,通常誤差范圍允許最大達(dá)到“真實”樹高的10%。這樣的條件表明,CHM數(shù)據(jù)總體精度指標(biāo)不能由個別樹高決定,而是應(yīng)從樣地中抽取一定數(shù)量上層林參照樣木,通過數(shù)據(jù)統(tǒng)計,更為可靠地指示CHM數(shù)據(jù)準(zhǔn)確性和精度誤差；在局部樣地尺度,也應(yīng)形成一定數(shù)量團(tuán)組抽樣模式,達(dá)到指示樣地單元精度水平的目的。

從圖4可以看出,2m半徑的圓形緩沖區(qū)內(nèi)有時會涵蓋鄰近高大林木的樹冠,因而CHMmax的算法設(shè)計導(dǎo)致了有些林木樹高提取不正確。顯然,需進(jìn)一步優(yōu)化上層林樣木,主要采取以下兩個步驟:一是去除樣地外圍高大林木的遮擋影響；二是提高上層林木緩沖區(qū)窗口尺寸,盡量剔除周邊低矮林木。具體設(shè)置參數(shù)為:依次采用2,3,4,5m半徑的圓形緩沖區(qū)窗口；設(shè)置原樣地范圍邊界向內(nèi)縮減5m。隨著緩沖區(qū)大小/樣地半徑變化,從表5看出,地面參照林木的株數(shù)依次為7 873株(5),3 366株(4),1 601株(3),1 017株(2),714株(1),逐步快速遞減。這與保障每個樣地保持一定數(shù)量抽檢樣木如5株以上的目標(biāo)矛盾。為平衡樣地尺度數(shù)量問題,在方案(5)基礎(chǔ)上,每個樣地抽取5株最高林木,作為補(bǔ)充抽樣方案(6)。

3.1章節(jié)中初步明確了正常情況下大部分地面參照林木與CHM樹高的總體精度水平,即CHM樹高與地面參照林木樹高存在系統(tǒng)誤差μ=0.7m,標(biāo)準(zhǔn)方差σ=±1m。據(jù)此采用標(biāo)準(zhǔn)差分級方法,以算術(shù)平均值μ作為中間0點,以±σ,±2σ,±3σ,±3σ以上/下為分界點,進(jìn)行分組計數(shù)統(tǒng)計,從而指示CHM的總體精度水平和顯著異常數(shù)據(jù)。從表5看,方案(6)最為有效的抑制了林木遮擋影響,-3σ以下的分組占比最低,僅為0.7%,可作為評價CHM樹高產(chǎn)品精度水平的最優(yōu)方案。由此得到示范區(qū)CHM產(chǎn)品的精度水平狀況:其中0～±σ的百分比為57.6%,0～±2σ的百分比為79.4%,-2σ以下占比2.0%,+2σ以上占比18.6%。圖6具體展示了表5中抽樣方案(5)和(6)的前后樣地尺度抽檢樣木的空間分布和高度誤差的對比情況。方案(6)避免了過多粗差樣木的干擾,集中通過5株上層林木準(zhǔn)確反映了樣地尺度的樹高誤差情況。

表5 地面參照林木抽樣方案及其標(biāo)準(zhǔn)差分級統(tǒng)計表

(a)樣地參照樣木全部合格情況

(b)樣地參照樣木全部不合格情況

(c)樣地參照樣木僅部分合格情況

注:黑色數(shù)字表示地面調(diào)查樹高,紅色數(shù)字表示地面樹高與CHMmax的差值。假定CHM樹高與真實樹高差值介于μ～±2σ為合格抽檢樣木,

依據(jù)表5抽樣方案(6)精度指標(biāo)分級統(tǒng)計結(jié)果,+2σ以上占比為18.6%。說明CHM產(chǎn)品主要問題是與實測樹高比較CHM樹高偏低。雖然激光雷達(dá)測量樹高的理論十分簡單,但是樹高的準(zhǔn)確性和精度受到儀器性能、采集規(guī)范和數(shù)據(jù)處理諸多因素的影響,十分復(fù)雜,具有一定的不確定性。CHM樹高誤差Etree的3個方面來源[8]分別是:1)冠頂高度偏低誤差Etop,其程度大下是林木自身形態(tài)如樹冠大小和形狀、枝葉疏密與激光脈沖能量、脈寬、傳感器靈敏度等相互共同作用的結(jié)果；2)林下地形DTM數(shù)據(jù)高估誤差EDTM,其程度與采集過程中的地面點質(zhì)量和地面點濾波處理的有效性密切相關(guān)；3)林木自身傾斜導(dǎo)致的誤差ELean,其程度與林木傾斜角度和局地地形坡度大小相關(guān)。除了以上誤差因素外,地面調(diào)查數(shù)據(jù)的錯誤也參雜其中,如,樹高記錄錯誤、林木定位錯誤等。上述誤差和錯誤可依據(jù)指標(biāo)分級的指示,逐個樣地進(jìn)行人工目視判讀,識別誤差來源和錯誤情況如圖6所示。

4 結(jié)論

航空激光雷達(dá)數(shù)據(jù)及其產(chǎn)品已成為森林管理和森林動態(tài)監(jiān)測的重要技術(shù)手段,但是其精度水平目前主要以平坦區(qū)域或明顯的人工地物控制點進(jìn)行評價；山區(qū)林地由于復(fù)雜地形和森林測量條件成為精度評估的“盲區(qū)”。隨著航空激光雷達(dá)在森林調(diào)查和管理工作中的廣泛應(yīng)用,大面積林區(qū)CHM數(shù)據(jù)精度指標(biāo)缺失的問題日漸突出。

本研究從現(xiàn)有地面林木水平定位精度優(yōu)于2m、測高精度為林木“真實”樹高10%的實際工作條件出發(fā),利用空間分析方法和統(tǒng)計分析策略,從工程應(yīng)用角度制定了CHM樹高產(chǎn)品的檢驗規(guī)程。在我國東北虎豹國家公園航空激光雷達(dá)森林資源調(diào)查示范區(qū)的應(yīng)用結(jié)果表明,該方法準(zhǔn)確、可靠地反映了地面調(diào)查樣地、樣木與CHM高度數(shù)據(jù)的空間一致性情況:一是驗證了山地林區(qū)CHM產(chǎn)品精度水平,指示出偏差異常情況,從而為查找數(shù)據(jù)采集和處理過程中存在問題提供線索；二是能夠檢驗地面調(diào)查樣地的樣木定位和高度測量的準(zhǔn)確性。