基于偏振反射模型和隨機森林回歸的葉片氮含量反演

張子晗,晏 磊,2,劉思遠,付 瑜,姜凱文,楊 彬,劉綏華,張飛舟*

1. 北京大學地球與空間科學學院， 遙感與地理信息系統研究所空間信息集成與3S工程應用北京市重點實驗室,北京 100871 2. 桂林航空工業學院， 廣西高校無人機遙測重點實驗室,廣西 桂林 541004 3. 湖南大學電氣與信息工程學院,湖南 長沙 410082 4. 貴州師范大學地理與環境科學學院,貴州 貴陽 550001

引 言

在科研人員利用高光譜遙感數據進行植被內部重要物質含量反演之后的30多年里， 高光譜遙感技術已廣泛應用于研究與植被相關的多種重要的生物化學物質， 例如： 水分、 纖維素、 蛋白質、 葉綠素、 氮元素、 硫元素和磷元素等。 葉片氮元素含量(leaf nitrogen concentration,LNC)與植物光合呼吸作用和其他的相關生物過程密切相關， 其含量僅占葉片總質量的0.2%～6.4%， 但卻能極大程度上影響全球氣候變化過程， 并在廣域碳氮循環過程中起到主導作用， 對于解構生態環境規律有重大科學意義。 氮元素在葉片內部的主要依托為蛋白質和葉綠素， 這兩種物質在可見光近紅外(400～2 500 nm)波段有獨特的吸收特性， 因而可以利用高光譜數據進行氮元素含量反演。

利用高光譜數據進行氮含量反演的算法可以大致分為三類： 第一類方法建立在葉片反射率與葉片氮含量之間統計關系的基礎上， 以逐次線性擬合回歸(stepwise multiple linear regression,SMLR)和偏最小二乘回歸(partial least squares regression,PLSR)為代表[1]， 在高光譜波段數較多時， 這類算法的回歸計算比較耗時， 精度相比于后兩類方法更低， 但由于算法原理簡單明晰， 仍被廣泛使用； 第二類算法建立在高光譜植被指數的基礎上， 通過計算高光譜植被指數來反演葉片氮含量[2-3]， 這類算法通過提取紅邊等深層次光譜特征來構建植被指數， 對于高光譜數據的挖掘仍舊不夠深入， 但由于計算簡單， 常用于在農學相關的應用領域粗略估計葉片氮含量； 第三類算法建立在機器學習的基礎上， 以支持向量機(support vector machine,SVM)、 人工神經網絡(artificial neural network,ANN)和隨機森林(random forest,RF)為代表[4-5]， 這類算法屬于遙感領域與計算機領域學科交叉的范疇， 具有極大的精度優勢。 不同的算法選擇會產生算法固有誤差， 影響氮含量反演精度， 為了盡可能減小算法誤差， 同時盡可能利用高光譜數據蘊含的信息， 本研究中測試了多種氮含量反演方法， 最終以反演精度為指標， 選擇了第三類算法中的RF算法作為本研究中葉片氮含量反演使用的主體框架。

上述三類氮含量反演算法的輸入參數皆為高光譜反射率， 但從定量遙感分析的角度來看， 光譜測量得到的反射率包括兩個部分的貢獻： 一部分來源于入射光在葉片內部的多次散射過程， 另一部分來源于入射光在葉片表面的鏡面反射過程[6]。 在多次散射過程中， 入射光大部分透過葉片表面進入葉片內部， 與葉片內部結構中的葉綠素和蛋白質等含氮生物化學物質發生相互作用， 在葉片內部多次散射后， 攜帶氮元素信息離開葉片， 這部分來自于葉片內部多次散射的漫反射率所攜帶的光譜信息是氮含量反演密切關注的。 而在鏡面反射過程中， 入射光部分在葉表蠟質層直接發生鏡面反射， 這部分反射能量不包含與葉片內部葉綠素和蛋白質等含氮生物化學物質相關的信息， 在氮含量反演中屬于高光譜反射率自帶的數據誤差， 如果使用包含鏡面反射部分的高光譜反射率直接進行氮含量反演會產生機理誤差， 影響氮含量反演精度。 依據菲涅爾原理， 葉片反射光的鏡面反射組分是部分偏振的， 而葉片反射光的漫反射組分由于在葉片內部多次散射導致的消光作用是非偏振的。 因此， 通過偏振反射建模方法獲得雙向偏振分布函數(bidirectional polarization distribution function,BPDF)， 可模擬計算偏振反射率， 進而扣除測量獲取的葉片高光譜反射率中與氮含量反演無關的鏡面反射組分， 從源頭減少機理誤差， 提高氮含量反演精度。

綜上所述， 影響葉片氮含量反演精度的誤差來源有兩個： 算法誤差和機理誤差。 選擇RF算法進行葉片氮含量反演算法可以盡可能減小算法誤差， 構建BPDF可以去除大部分機理誤差， 最終實現基于BPDF和RF的葉片氮含量精確反演。

1 實驗部分

1.1 氮含量反演數據

本研究中用于氮含量反演的數據為主數據集。 葉片氮元素反演的研究地點為德國巴伐利亞國家森林公園[7](49°3′19″N,13°12′9″E)。 根據優勢樹種的類型， 選擇使用國際地圈生物圈計劃(international geosphere-biosphere programme,IGBP)土地覆蓋分類體系， 優勢樹種與IGBP類別的對應情況見表1。

表1 巴伐利亞國家森林公園內優勢樹種對應的IGBP類別Table 1 IGBP classes of dominant vegetation inBavarian Forest National Park

由于本研究區域內的優勢樹種中闊葉樹木均為落葉闊葉品種， 針葉樹木都為常綠針葉品種， 所以在研究區域中被簡單標注為3種類型： 落葉闊葉林(IGBP04)、 常綠針葉林(IGBP01)和針闊混交林(IGBP05)， 其點位空間分布如圖1所示。

圖1 葉片氮元素反演點位空間分布Fig.1 Spatial distribution of foliar nitrogen content retrieval plot

共有26個點位作為氮元素反演研究區(8片落葉闊葉林， 8片常綠針葉林， 10片針闊混交林)。 其中， 每個采樣區在實地占地為30 m×30 m， 點位范圍用分層隨機抽樣的方法選擇8～9棵樹木采樣， 單顆樹木至少采集20片葉片， 妥善保存好后在實驗室用凱式定氮法測量氮元素含量。 高光譜數據的獲取時間為2013年7月22日， 利用HySpex傳感器系統在3 000 m高度獲取高光譜影像(HySpex傳感器系統包括2臺高光譜成像光譜儀， 一臺覆蓋400～1 000 nm的光譜范圍， 另一臺覆蓋1 000～2 500 nm的光譜范圍)， 高光譜數據共有418個窄波段， 410～2 495 nm內， 每5 nm一個通道。 氮含量反演使用的原始反射率是經過幾何校正和大氣校正處理后的26個樣區內的418個波段的雙向反射因子數據(bidirectional reflectance factor,BRF)， 除了BRF數據之外無人機系統還記錄了獲取數據時的觀測幾何(太陽天頂角、 觀測天頂角、 相對方位角)。

1.2 偏振反射建模數據

偏振反射建模數據為輔數據集。 本研究通偏振反射建模獲取BPDF， 進而估計偏振反射率。 用于偏振反射建模的BRDF-BPDF數據庫[8]由搭載在PARASOL衛星上的POLDER傳感器獲取。 本研究中為了建立合適的BPDF模型所使用的數據來自于POLDER經過校正的產品， 包含了數據獲取時的經緯度和觀測幾何、 地物的IGBP類別以及地物類型同質性、 樣區的歸一化植被指數(normalized difference vegetation index,NDVI)、 氣溶膠情況、 6個波段(490,565,670,765,865和1 020 nm)的反射率和865 nm波段偏振反射率。

該BRDF-BPDF數據庫從POLDER的連續7年的數據中挑選出質量最好的1年， 即2008年， 同時， 在POLDER數據中挑選了質量最優秀(盡可能滿足對全球范圍的覆蓋、 數據有效性高、 樣區IGBP類別勻質性高)的點位。 最終數據庫的時間跨度和點位選擇定為： 包括2008年12個月份的數據， 每個月份采取50個均勻分布在全球的點位， 在每個點位上按照不同的觀測幾何選擇高質量POLDER/PARASOL體系的數據。 本研究使用該數據庫中2008年7月份歐洲范圍內的數據。

1.3 數據一致性說明

使用兩組相對獨立的數據集， 主數據集用于氮含量反演， 輔數據集用于建立偏振反射模型， 使用主數據集和輔數據集的流程圖如圖2所示。 從數據獲取的時間來看， 考慮到偏振反射建模具有較強的季節周期性特征， 為了與巴伐利亞國家森林公園的氮含量反演高光譜數據的獲取時間契合， 在進行偏振反射建模時僅使用BRDF-BPDF數據庫數據中與機載數據同月份的數據， 以保證時間一致性。 從數據獲取的地理位置來看， 由于氮含量反演數據的研究位置在德國， 同時考慮到在BRDF-BPDF數據庫數據中獲取的點位數量不能太少， 最終根據BRDF-BPDF數據庫數據的經緯度， 篩選出歐洲范圍內的數據進行偏振反射建模， 以保證地理位置一致性。 從兩組數據集的空間分辨率來看， POLDER數據的分辨率為6.5 km， 無人機高光譜數據的分辨率為1.65 m(400～1 000 nm)和3.3 m(1 000～2 500 nm)， 但由于BRDF-BPDF數據庫中收錄的POLDER數據勻質性大于75%， 可認為在數據匹配時受空間尺度的影響較小。

圖2 數據使用流程圖Fig.2 Flowchart of data utilization

1.3 方法

1.3.1 偏振反射建模

研究中獲取的高光譜數據為418個窄波段的BRF。 根據菲涅爾定律， 地表反射是部分偏振的， 其中的偏振部分主要為線偏振， 也是由觀測幾何(太陽天頂角、 觀測天頂角、 相對方位角)決定的， 可用BPDF表征其多角度偏振反射特性。 有文獻表明， 在可見光近紅外波段， BPDF可視為觀測幾何穩定時的光譜不變量， 在不同波段的偏振反射率為定值[9]。 BPDF的決定因素是地物類型和觀測幾何， 為了去除地表反射率的偏振部分， 需要利用BRDF-BPDF數據庫中的7月份歐洲范圍的數據， 建立3個分別適應IGBP01， IGBP04和IGBP05這3種不同地表IGBP類別的BPDF模型以模擬估計偏振反射率(polarization bidirectional reflectance factor,PBRF)， 進而從BRF中去扣除具有部分偏振特性的鏡面反射組分——偏振反射率PBRF， 最終獲得漫反射反射率(diffuse bidirectional reflectance factor,DBRF)， 以消除機理誤差的影響。

目前精度較高的BPDF模型幾乎全部為半經驗模型， 純粹的經驗模型雖然能夠在局部較好擬合偏振反射率， 但由于缺乏物理機理， 模型參數嚴重依賴訓練數據， 且過擬合明顯。 半經驗模型通過較少的經驗參數對包含物理意義的模型進行校正調整， 其結果是更具有普適性。 本研究選擇了5個應用最廣、 精度最高的BPDF模型用于模擬IGBP01， IGBP04和IGBP05這3種不同地表IGBP類別下的偏振反射率PBRF， 分別為： Nadal模型[10]、 Waquet模型[11]、 Maignan模型[12]、 Litvinov模型[13]和Diner模型[14]。

為了使公式簡潔明了， 在BPDF模型的公式表達中使用一些簡記。 其中， 太陽光方向與天頂方向(其夾角稱為太陽天頂角θS， 觀測方向與天頂方向的夾角稱為觀測天頂角θV， 太陽光方向在地表的投影與觀測方向在地表的投影之間的夾角稱為相對方位角φ， 鏡面反射過程中的法線與天頂方向的夾角稱為半天頂角θh， 入射方向與反射方向之間的夾角的一半為入射角αI， 折射角為αT， 入射方向與反射方向之間方向的變化為方向散射角γ， 折射率為N。 部分角度的余弦值簡記為

μS=cosθS

(1)

μV=cosθV

(2)

μI=cosαI

(3)

μT=cosαT

(4)

由菲涅爾公式Fp(γ,N)得到的變量之間的相關關系見式(5)

5種模型的偏振反射建模公式表達見表2。

表2 5種偏振反射模型的公式表達Table 2 Formulas of 5 BPDF models

利用輔數據BRDF-BPDF數據庫中的數據分別建立IGBP01，IGBP04和IGBP05這3種地表覆蓋類型下的5種BPDF模型， 分別擬合得到自由參數。 通過十折交叉驗證的數據分組策略， 以均方根誤差(rootmeansquareerror,RMSE)和可決系數(coefficientofdetermination,RSQ)為指標， 挑選3種地表覆蓋類型下的最優BPDF模型來模擬偏振反射率PBRF。

1.3.2 氮含量反演

氮含量反演研究基于德國巴伐利亞國家森林公園的實測氮含量主數據、 無人機高光譜主數據和偏振反射建模獲得模擬仿真結果， 研究主要包括三個步驟： 首先， 通過BPDF模擬計算偏振反射率PBRF， 進而計算漫反射率DBRF， 獲取氮含量反演的算法輸入對照組(原始反射率BRF和去偏漫反射率DBRF)； 然后， 設置RF算法參數， 利用整體散射系數和去偏后的漫反射散射系數作為算法輸入對照組進行葉片氮含量反演； 最后， 橫向對比RF算法和其他算法的表現， 并對偏振反射率在氮含量反演中的影響進行定量評估。

1.3.3 精度評估

偏振反射建模使用POLDER/PARASOL體系的數據庫中數據量較大， 采用十折交叉驗證的思路進行分組后驗證精度， 重復10次十折交叉驗證， 取均方根誤差RMSE和可決系數RSQ的均值作為BPDF模型精度的度量。 氮含量反演使用26個點位的高光譜數據受到數據量的限制， 采用留一法交叉驗證的思路進行數據分組后驗證精度， 同樣以均方根誤差RMSE和可決系數RSQ作為葉片氮含量反演精度的度量， 見式(6)和式(7)。

2 結果與討論

2.1 BPDF模型選擇結果與分析

首先， 按照氮含量反演研究時間、 區域進行輔數據集BRDF-BPDF數據庫的篩選， 選擇了POLDER數據庫中2008年7月采集的歐洲范圍內的數據； 然后， 利用十折交叉驗證的數據劃分思路將數據庫中的數據隨機劃分為訓練集和驗證集， 利用訓練集數據擬合獲得5個BPDF模型的自由參數； 獲得參數組后利用中位數方法求得單樣區最優模型參數組， 根據該模型參數組預測偏振反射率PBRF； 將預測得到的偏振反射率與包含于數據庫中的實測865 nm偏振反射率數據進行對比， 最后利用均方根誤差RMSE和可決系數RSQ來評價BPDF模型的效果， 5種模型的評價結果見表3。

表3 IGBP01/04/05地表覆蓋類型下5種BPDF模型精度評價Table 3 Accuracy assessment of 5 BPDF models on IGBP01/04/05

最終， 以可決系數RSQ為主要評判指標， 均方根誤差RMSE為輔助評判指標， 在RSQ盡可能大的情況下選擇RMSE盡量小的BPDF模型。 從整體來看， 在同一地表覆蓋類型下不同BPDF模型的模擬精度差異不大(保留3位小數后， 均方根誤差RMSE差異極小， 可決系數RSQ差異相對明顯)， 而在不同地表覆蓋類型下， IGBP04地表覆蓋類型的BPDF模型表現顯著優于IGBP01和IGBP05， 說明闊葉樹木的偏振效應相較于針葉樹木而言更易模擬， 推測可能與起偏介質的面積大小有關。 模型選擇結果為： 在IGBP01地表覆蓋類型下最優BPDF模型為Nadal模型， 在IGBP04地表覆蓋類型下最優BPDF模型為Litvinov模型， 在IGBP05地表覆蓋類型下最優BPDF模型為Maignan模型。 由于不同地表覆蓋類型的起偏特性不同， 最優BPDF模型并不統一， 分別利用選出的3種BPDF模型模擬對應的3種地表覆蓋類型的偏振反射率。

2.2 RF模型參數設置與葉片氮反演精度橫向對比分析

在偏振反射建模的基礎上利用原始反射率BRF和漫反射率DBRF作為RF算法輸入進行葉片氮含量反演。 RF的核心參數為決策樹數量， 本研究分析比對了去偏前后決策樹數量為5， 10， 50和100的RF算法的表現， 見表4。 研究結果說明， 在決策樹數量一定時， 使用去偏數據進行葉片氮含量反演普遍能獲得精度提升且在RF決策樹數量設定為10時獲得最優氮含量預測模型。 利用RF算法進行氮含量反演能獲取的最高精度為： 可決系數達到0.803， 均方根誤差達到0.252。

表4 不同決策樹數量的RF算法精度評價Table 4 RF algorithm accuracy assessment ofdifferent estimator numbers

2.3 偏振反射率占比統計分析

在26個采樣點位， 偏振反射率PBRF在整體反射率BRF中的平均占比為3.580%～11.649%， 各個點位內的偏振反射率占比箱形圖見圖3， 箱形圖形象地展示了偏振反射率PBRF占比的算數均值、 最小值、 下四分位數、 中位數、 上四分位數和最大值。 該分布表明偏振反射率PBRF在整體反射率BRF中平均占比為3.350%， 最高占比16.519%， 說明偏振反射率的影響是不可忽略的， 機理誤差足以影響反演精度， 在高光譜葉片氮含量反演中有去除偏振反射率PBRF的必要性。

圖3 偏振反射率在整體反射率中占比的分布圖中箱體內的×表示算數均值， 箱體內的黑色橫線表示中位數Fig.3 Ratio of PBRF to BRFIn this figure,the × in the box represents the mean value,and the transverse line in the box represents the median value

2.4 葉片氮反演精度橫向對比分析

橫向對比了PLSR算法、 主成分回歸(principal component regression,PCR)算法、 支持向量回歸(support vector regression,SVR)算法、 K-近鄰(K-nearest neighbor,KNN)回歸算法和RF回歸算法在氮含量反演中的表現， 利用留一法進行交叉驗證， 精度評價結果見表5。 其中， 由于偏最小二乘算法應用最為廣泛[1]， 將其作為基線算法進行整體精度評估。 研究結果說明， 在利用上述方法進行氮含量反演的過程中， 去除偏振后葉片氮含量反演的精度普遍都有提升(KNN算法在選擇的最鄰近元數目K足夠大時， 最鄰近元的選擇在去除偏振前后相同， 因而模型預測結果相同， 無精度提升)。 以均方根誤差RMSE為核心指標評估提升幅度， 平均提升幅度為4.244%， 其中RF算法去除偏振后葉片氮含量反演精度提升幅度最大， 為13.103%， 且精度最優(達到了上述5種反演算法中的最高可決系數0.803和最低均方根誤差0.252)， 充分說明了RF算法利用高光譜數據進行葉片氮含量反演的潛力。 相較于作為基線算法的PLSR(未進行光譜去偏)， 通過偏振反射建模去除大部分光譜機理誤差并通過選擇隨機森林模型減小算法誤差后， 葉片氮含量反演精度整體提高32.440%， 證明了去除機理誤差和減少算法誤差的必要性， 實現了高精度葉片氮含量反演。

表5 5種葉片氮含量反演算法精度評價Table 5 Accuracy assessment of 5 LNC retrieval algorithms

3 結 論

研究旨在提高葉片氮含量反演精度， 從機理誤差和算法誤差兩個誤差源入手， 通過偏振反射率建模去除了大部分機理誤差， 通過選擇精度最優的隨機森林算法并調整參數盡可能降低了算法誤差， 最終實現了葉片氮含量精確反演， 相比于基線方法整體精度提升為32.440%(可決系數RSQ達到0.803， 均方根誤差RMSE達到0.252)。 同時， 在討論中定量評估了偏振反射率PBRF在整體反射率BRF中的占比分布并對比了去除偏振前后多種氮含量反演算法的精度變化， 體現了隨機森林算法在高光譜信息挖掘中的優勢， 也證明了去除不攜帶氮含量信息的偏振反射率的必要性。

致謝：本研究使用的德國巴伐利亞國家森林公園葉片氮含量實測數據和無人機高光譜數據由威斯康星大學麥迪遜分校的王智慧博士提供， 在此表示衷心的感謝！