基于林內(nèi)圖像的單位面積碳儲量估計方法

王雪峰 陳珠琳 管青軍 劉嘉政 王 甜 袁 瑩

(1.中國林業(yè)科學研究院資源信息研究所 北京 100091；2.內(nèi)蒙古莫爾道嘎林業(yè)局 額爾古納 022357)

森林是陸地生態(tài)系統(tǒng)中最大的碳庫，森林植被碳儲量約占全球植被碳儲量的77%，森林土壤碳儲量約占全球土壤碳儲量的39%，森林每年固定的碳約占整個陸地生態(tài)系統(tǒng)的2/3，在全球碳循環(huán)和緩沖氣候變化過程中具有不可替代的作用(Brownetal.，1993；Johnstonetal.，2004)。據(jù)第八次全國森林資源清查結果，我國森林面積現(xiàn)已達2.08億hm2，森林覆蓋率為21.63%，其中，天然林面積1.22億hm2，占全國森林總面積的58.65%，可見，明確天然林生態(tài)系統(tǒng)碳儲量和碳庫分配特征對更好了解我國森林生態(tài)系統(tǒng)碳庫和未來森林生態(tài)系統(tǒng)碳匯均具有重要意義(魏亞偉等，2015)。

隨著溫室效應加劇，人們逐漸意識到森林在全球碳循環(huán)中的重要作用，但森林究竟是充當碳源還是碳匯角色仍存在很大爭議(Caoetal.，1998)。為了解全球碳計量，學者們對世界范圍內(nèi)不同區(qū)域、不同類型的森林碳儲量進行了估算(Bracketal.，2002；Masellietal.，2009；Rutishauseretal.，2013)，采用的方法有樣地實測法、生理生態(tài)法和遙感反演法等(Sileshi，2014)。森林碳儲量的主要載體為森林生物量，大部分研究直接或間接以生物量的現(xiàn)存量乘以植被的含碳率推算得到森林碳儲量。樣地實測法(也稱收獲法)是普遍采用的估算方法，可分為皆伐法、平均生物量法和相對生長法3類(木村允，1976；佐藤大七郎，1973)，其中，皆伐法精度較高，但要耗費大量人財物力且對環(huán)境破壞大；平均生物量法樣地選擇受人為控制，常出現(xiàn)實測結果偏高的現(xiàn)象；相對生長法通過建立林分各部分或整體的生物量回歸方程估算碳儲量，靈活多變，但需要針對不同地區(qū)、不同樹種分別考慮模型形式、參數(shù)估計和相容性等問題(Kondwanietal.，2018；Chojnackyetal.，2014)。生理生態(tài)法基于植物生理因子或生態(tài)環(huán)境因子建立數(shù)學模型，從而估算森林碳儲量，其優(yōu)點是不需要大量工作即可完成碳儲量估算，但由于到目前為止人類還不完全了解植物生理生態(tài)機制，因此該方法基本停留在研究階段。遙感反演法(Djomoetal.，2017；Main-Knornetal.，2013；Hoomanetal.，2015)在大范圍森林碳儲量估算方面具有一定優(yōu)勢，但其總體估計精度偏低，且作為數(shù)據(jù)源的衛(wèi)星影像或航片信息也非用戶隨手可得，同時不同季節(jié)、大氣效應和天氣情況等都會對估計結果產(chǎn)生影響；雖然機載激光雷達(Rossetal.，2004；Patenaudeetal.，2004)、機載高光譜技術以及融合后的多源數(shù)據(jù)能夠解決一些上述問題，但高成本等嚴重限制了其使用和推廣力度。

森林碳儲量估計是一個非常棘手的問題，人們一直在探索消耗少、精度高的碳儲量估計方法，但是研究進展緩慢。隨著圖像或視頻設備的普及，無論是通過林業(yè)現(xiàn)場的圖像傳感器還是調(diào)查員攜帶的數(shù)字相機、智能手機，人們可以非常容易獲取反映林木疏密和高矮等信息的林內(nèi)圖像，那么，如果能基于這些圖像估測森林碳儲量，則必將極大改觀現(xiàn)有碳儲量的調(diào)查方法。鑒于此，本研究從普通圖像入手，探尋適合林內(nèi)圖像的目標信息分類提取算法及與碳儲量的關系，并嘗試提出一種基于林內(nèi)圖像的單位面積森林碳儲量估計方法。

1 研究區(qū)概況與數(shù)據(jù)

1.1 研究區(qū)概況

研究區(qū)位于內(nèi)蒙古大興安嶺西北麓的寒溫帶針葉原始林內(nèi)，該區(qū)森林覆蓋率94.7%，地理坐標為120°00′26″—121°19′31″E，51°03′16″—52°06′00″N，屬大陸性季風氣候區(qū)的低山丘陵地帶，坡度一般小于15°，海拔450～1 405 m，年均氣溫-5～-7 ℃，冬寒夏溫，年降水量414～528 mm，喬木以興安落葉松(Larixgmelinii)為主，另外有白樺(Betulaplatyphylla)、山楊(Populusdavidiana)、樟子松(Pinussylvestrisvar.mongolica)等；林下植物主要有興安杜鵑(Rhododendrondauricum)、杜香(Ledumpalustre)、篤斯越桔(Vacciniumuliginosum)、五味子(Schisandrachinensis)等。

1.2 研究數(shù)據(jù)

沿著防火林帶隨機設置64塊圓形樣地，樣地中心選擇在垂直距離防火帶100 m處，樣地半徑為10.3 m，中心海拔482～1 018 m。記錄樣地海拔、經(jīng)緯度、灌草、土壤等基本信息，實測樣地內(nèi)全部樹木的胸徑、樹高。使用分辨率為3 456×5 184像素的佳能單反700D相機，攝取林分縱斷面和樹冠2類圖像?？v斷面圖像采集從磁北方向開始，每隔45°縱向采集1張，每點8張，拍攝時保持過光心和成像中心的連線與攝影方向地面平行，圖1所示為實際樣地圖像。由于縱斷面圖像包含林地等信息，為避免林木外因素對預測結果的影響，采用交互式方式剔除地平線以下部分，即僅使用地平線上部分作為待用圖像。最終計算的圖像密閉度最小值為13.3%、最大值為86.8%，蓄積量分布區(qū)間為16～367 m3·hm-2，碳儲量基于張俊(2008)模型根據(jù)蓄積量計算得到。

圖1 樣地圖像

2 研究方法

林木越密集，圖像前景點數(shù)越稠密，林木越高大，成像后其所占用的縱向像素數(shù)越多，即林分縱斷面圖像隱式包含林分密度和高度2類復合信息。這種感官上的“密集”“高大”與林分碳儲量存在正相關關系。假設某空間點周圍存在虛擬環(huán)幕，外圍林木投影到環(huán)幕上，則環(huán)幕上林木點與整個環(huán)幕的比值就代表該空間的密集大小程度，區(qū)間為[0，1]，本研究將其稱為“林分密閉度”。如果在該點用攝影機環(huán)周攝影得到n幀圖像，則此n幀圖像的平均密閉度可反映林木對該點周圍空間的平均“填充”程度，如此，即將感官上林木對空間的占有程度由定性描述轉為定量刻畫，其相對客觀地反映了林木密度和大小信息，這一參數(shù)必定與單位面積林木碳儲量存在關系。

2.1 林分密閉度計算

欲得到單位面積林木碳儲量，首先要獲取林分密閉度，即需要計算某點周圍全部林木圖像的前景比例均值。

2.1.1 確定全局閾值 找到圖像中一個灰度值ts，通過如下方法將圖像變?yōu)橹挥星熬昂捅尘?個灰度級的新圖像Inew：

(1)

如果hg是灰度為g的像素數(shù)(即灰度直方圖)、a是g≤ts的灰度平均值、b是g>ts的灰度平均值，則容易滿足此約束條件。其中，

Brink等(1996)采用對稱交叉熵全局閾值，確定使下式dt達到最小的ts就能將圖像分割成只有前景和背景2個灰度級的新圖像：

(2)

對稱交叉熵法給出了一個優(yōu)化閾值，該閾值相比大津算法和最大熵算法得到的閾值更接近于最優(yōu)閾值(Li，1993)。但由于整幅圖像的等權重考慮，對于光線不均的林內(nèi)圖像常使樹體內(nèi)部出現(xiàn)空洞、非樹體部分劃歸為樹體的情況，因此必須兼顧當前像素的周邊環(huán)境信息，最終決定該像素的歸屬。

2.1.2 閾值調(diào)整與前景抽出 如果g為n×m圖像中第x(x=0，…，n×m)個像素的灰度，C++通過hg++遍歷全部像素，即可得到該圖像的灰度直方圖，此時，hg是灰度為g的像素數(shù)。需要注意，在遍歷前要為從黑(0)到白(255)的每一個灰度級各分配8位內(nèi)存空間并清零。按照上述對稱交叉熵法對圖像進行全局性二值化分類，由于各像素是獨立考慮的，所以常出現(xiàn)本來是樹體像素被判斷成背景、而非樹體像素被判斷成前景的情況。客觀需求是，樹體內(nèi)部像素點無論灰度值多少均應歸為前景、而樹體外部像素點無論灰度大小均應歸于背景，一種解決策略是考慮相鄰像素屬性并采取相應措施。

在判斷焦點像素歸屬時，考慮其上下左右各r個像素，規(guī)定以當前焦點像素為中心的方形內(nèi)，小于閾值的像素數(shù)大于γ=r(3r+2)時判為前景，否則判為背景。如圖2a所示，當r=2時，a、b、c、d、o、e、f、g、h將全部被判為前景；如果r=1，則a、c被判為前景，其他被判為背景。而對于圖2b，不論r=1還是r=2，u、w、v均被判為背景。可以看出，對于物體內(nèi)部，如果考慮鄰近周邊像素屬性，適當?shù)陌霃胶团袛鄿蕜t對樹體有填充效果；而物體外部正好相反，具有抑制作用，這恰好與實際問題相符合。

圖2 考慮鄰域像素的判斷方法在樹體內(nèi)外產(chǎn)生的不同效果

在對稱交叉熵法處理基礎上，結合鄰域像素屬性進行分類通常能取得較為滿意的結果。圖3a為一幅圖像的縱剖面圖；圖3b為基于對稱交叉熵法的計算結果，此時得到的前景面積占47.09%；圖3c在獲取閾值時采用同樣方法，但在分類時不單純針對單個像素，還考慮焦點像素與鄰域像素的關系。鄰域半徑取1，當鄰域內(nèi)的前景像素數(shù)大于5歸為同類，此時前景面積占全部面積的48.74%，增加1.65%，觀察發(fā)現(xiàn)，樹體內(nèi)部應填充的部分被填充了，而樹體外部無關像素被忽略。

圖3 對稱交叉熵法處理及在此基礎上考慮鄰域像素屬性的算法比較

為了更加清晰起見，圖4給出了興安落葉松樹體部分的一個實際圖像計算結果。由于樹冠、樹葉等相互遮擋，樹體內(nèi)部光線不均，加之噪聲存在，致使樹體內(nèi)部不同區(qū)域出現(xiàn)過暗或過亮斑點(圖4a)；采用對稱交叉熵法處理，這些斑點可能區(qū)分成背景(圖4b)；加入鄰域其他像素信息后，誤判有所改善(圖4c)。這種利用與鄰域像素關系以決定當前像素歸屬的方法具有腐蝕和膨脹雙重特性，即當焦點像素處于樹體內(nèi)部時容易將該點歸為樹體，當焦點像素處于樹體外部時容易將該點歸為背景，樹體內(nèi)部被膨脹、樹體外部被腐蝕，符合客觀實際。

圖4 對稱交叉熵法及與進一步考慮鄰域像素的結果比較

需要注意的是，兼顧鄰域像素屬性的方法對鄰域半徑大小和區(qū)域內(nèi)多少像素為前景時將焦點像素判為前景的入選準則很敏感，如果這2個參數(shù)選擇不當，將出現(xiàn)較大誤判。試驗表明，鄰域半徑以不大于2為宜，入選準則選擇γ∈[3r2+2r,4r2+2r]之間的一個數(shù)較為合適。本研究中所用圖像取r=1、γ=6，即在以焦點像素為中心的3×3的鄰域內(nèi)，如果有大于6個相似像素出現(xiàn)，則將焦點像素歸為該類。

圖像被執(zhí)行前景和背景二值化歸類后，統(tǒng)計樹體像素比例，即得到該幅圖像所代表的林分密閉度，進一步計算旋轉一周各圖像林分密閉度均值，即得到攝影點的林分密閉度。很明顯，該均值算法與圖像重疊沒有關系，且采集攝影點周圍圖像越全面，越能代表該攝影點的林分密閉度。

2.2 碳儲量估計

由林分密閉度求算單位面積碳儲量，最直接有效的方法就是使用模型。林業(yè)統(tǒng)計模型眾多，從中篩選出效果較好且使用較多的8個模型列于表1。模型中，x為林分密閉度，y為林分單位面積碳儲量。式(3)～(6)為2個參數(shù)，式(4)～(6)可首先由線性形式給出初始參數(shù)，式(7)～(10)為3個參數(shù)。需要說明的是，盡管增加參數(shù)可能帶來更高的模擬精度，但未必是最佳模型；不過，由于生物統(tǒng)計模型中的廣泛使用，本研究仍選擇使用頻率較高的4個3參數(shù)模型，以便于比較。

表1 碳儲量預測模型選取

首先,根據(jù)建模數(shù)據(jù)計算決定系數(shù)、平均殘差和均方根誤差；然后，利用其他獨立數(shù)據(jù)對模型進行驗證，計算的參數(shù)包括平均誤差、平均絕對誤差、平均百分比誤差和平均絕對百分比誤差等，計算公式參見文獻(李鳳日，2004)；最后，綜合考慮這些數(shù)據(jù)，評價模型擬合效果。

3 結果與分析

從64塊樣地中隨機抽取50塊樣地數(shù)據(jù)用于建模，剩余14塊樣地數(shù)據(jù)作為獨立數(shù)據(jù)用于模型檢驗。模型擬合使用IBM SPSS Statistics 21軟件，見表2、3。

表2 模型擬合參數(shù)

從表3可以看出，2參數(shù)模型中線性模型(3)的決定系數(shù)最大，3參數(shù)模型中(8)～(10)結果接近，以邏輯斯蒂模型(9)擬合效果最好。從平均誤差結果看，線性模型的絕對值最小，表明林分密閉度與林分單位面積碳儲量之間可能存在線性關系，但綜合考慮其他指標，邏輯斯蒂模型擬合效果相對較優(yōu)。為了更清晰了解擬合效果，從2參數(shù)和3參數(shù)模型中分別選出擬合精度較高的前2個模型，繪制實測值和預測值對應的散點圖(圖5)，粗略地看，幾個模型的預測結果接近，但仔細觀察發(fā)現(xiàn)，邏輯斯蒂模型擬合效果在較大值和較小值的兩端表現(xiàn)更加優(yōu)良。

表3 模型擬合精度

雖然圖5的總體趨勢較為理想，但各點間較分散，說明還存在沒有考慮到的其他因素影響擬合效果。由于大多情況下海拔越高樹高越低，從而影響碳儲量，為此引入一個虛擬變量k代表海拔因素。

圖5 預測值與實測值的對應效果

考慮到大興安嶺北部地區(qū)海拔差異不是特別大，同時不至于使模型難于應用，以本次樣地中心海拔平均值750 m為界，如果海拔高小于該數(shù)值定義k=0，否則k=1。由于不分類時線性模型平均誤差最小，而邏輯斯蒂模型在均方根誤差和決定系數(shù)等指標上表現(xiàn)最優(yōu)，因此本研究僅測試這2個模型。

嘗試將虛擬變量加到模型不同參數(shù)位置上，最終顯示，表4中模型設置方法得到的模擬精度最高。觀察表4結果，考慮海拔因素后，2個模型的模擬精度均得到很大改觀，如邏輯斯蒂模型，決定系數(shù)從0.859升至0.949，平均絕對誤差從9.738降至4.229。

表4 加入海拔后線性和邏輯斯蒂模型的形式、參數(shù)和精度

圖6 加入虛擬變量前后直線與邏輯斯蒂模型的殘差比較

4 討論

4.1 碳儲量模型形式

從本次數(shù)據(jù)模擬結果看，基于縱向圖像預估單位面積碳儲量，邏輯斯蒂曲線精度最高，表現(xiàn)出較好的適應性；但是僅有2個參數(shù)的線性方程表現(xiàn)出與3參數(shù)型曲線很接近的結果，由此推測，圖像密閉度與單位面積碳儲量之間可能存在線性關系。如果二者確實是線性關系，那么不僅給模型參數(shù)估計帶來方便，同時也更容易理解，因為固定空間內(nèi)隨著樹木不斷生長即碳儲量不斷增加，林分密閉度也不斷加大，二者之間呈線性關系最容易理解；當達到土地的最大承載力時，單位面積碳儲量達到極限值，此時林分密閉度為1。以本研究模型比如海拔高度大于750 m的樣地(k=1)為例，當密閉度達到1時，林地上單位面積碳儲量達到133.436 2 t，與邏輯斯蒂曲線的最大值131.952 2 t接近。所以有必要繼續(xù)測試圖像密閉度與單位面積碳儲量間的關系，如果確是線性，沒必要引用更復雜的多參數(shù)方程。

4.2 郁閉度與碳儲量的關系

由密閉度一個參數(shù)計算單位面積碳儲量，應用時僅攝取縱向圖像即可，不需要額外參數(shù)就可以估測林分碳儲量，從而給應用帶來方便，這也是本研究的出發(fā)點。郁閉度是森林調(diào)查時一個常用的調(diào)查因子，在遙感應用中，基于郁閉度估測林分蓄積量或林內(nèi)因子也是常用做法。研究表明，碳儲量隨郁閉度增加呈上升趨勢(劉茂秀等，2011)，因此，如果將郁閉度作為碳儲量估計的一個自變量或分類變量能夠進一步提高碳儲量估計精度，則可以考慮繼續(xù)建立基于密閉度和郁閉度的“雙度”碳儲量預測模型，類似于一元和二元材積表，這將給用戶提供更加豐富的選擇權力。由于在采集林木縱向圖像時，將鏡頭轉向冠層就可以獲得冠層圖像，由冠層圖像計算郁閉度簡單且精度高，而這不需要增加太多的內(nèi)外業(yè)工作量，因此，繼續(xù)研究模型中是否包含郁閉度情況下，對林分碳儲量估計帶來多大影響有重要應用價值。

5 結論

森林碳儲量估計工作量巨大。本研究基于實際感知經(jīng)驗，提出了一種基于林內(nèi)圖像的單位面積碳儲量估計方法，并以興安落葉松為例，給出了具體估計方法、步驟和精度分析。

1)基于林內(nèi)圖像能夠?qū)崿F(xiàn)對單位面積森林碳儲量的較高精度估計。表面上看，本研究屬于圖像理解探索內(nèi)容，但實際上是基于林內(nèi)圖像求算林分密度和高度進而估計碳儲量的一種方法。很明顯，單位面積森林碳儲量與其林木密度和高度直接相關，而林分縱向圖像正好隱式包含這2類信息，同時圖像方法本身又能相對準確地獲取這些信息，因此，估計碳儲量也僅僅是一個度量轉換問題。

2)采用對稱交叉熵法獲得初始值，進一步考慮鄰域像素屬性以決定當前像素所屬，能夠減弱林內(nèi)光線不均對圖像灰度造成的影響，提高林分分類的準確性。盡管對稱交叉熵法是一種較好的圖像二值化算法，但是由于算法的整體統(tǒng)一處理，致使計算局部光照不均具有普遍性的林內(nèi)圖像時出現(xiàn)較大偏差。本研究提出利用與鄰近像素關系以決定當前像素歸屬的方法，既保留了全局算法的優(yōu)點，又增加了對具體像素的處理信息，從而大大增強了適應性，提高了分類精度。

3)將圖像作為輸入變量，計算出林分密閉度后根據(jù)邏輯斯蒂模型就可以估計林地上單位面積碳儲量，為碳儲量估計提供了一種解決之路。隨著林業(yè)物聯(lián)網(wǎng)技術的不斷應用，林內(nèi)圖像等各類傳感數(shù)據(jù)不斷傳回，如何抽取此大量數(shù)據(jù)中所隱含的實用價值，特別是如何理解其中的海量圖像，已成為林業(yè)大數(shù)據(jù)應用的難點之一，本研究方法可以從實時圖像中獲取監(jiān)測地碳儲量信息，對于林業(yè)大數(shù)據(jù)應用是一個借鑒。

一直以來，林業(yè)上使用材積表估計林地蓄積量，對應起見，本研究基于林內(nèi)圖像求算碳儲量(或蓄積量)的方法稱為圖像表法。能夠看到，通過建立的圖像表也可以實現(xiàn)對單位面積森林碳儲量(或蓄積量)的估計，因為該方法通過手機或照相機掃描一周圖像即可得到該點為中心的單位面積碳儲量(或蓄積量)，效率要遠高于目前的測徑方法。因此，通過建立某一樹種廣泛樣地的圖像表，以此為基礎估計該類型單位面積林地碳儲量(或蓄積量)可能成為未來一種全新的快速測算方法。