地表反照率機器學習估算方法

吳錦超，吳永靜，林超，竇寶成，劉銳

(1．廣東省國土資源技術中心，廣州 510000；2．北京吉威數源信息技術有限公司，北京 100043)

0 引言

地表反照率是廣泛應用于地表能量平衡、中長期天氣預報和全球變化研究中的重要參數之一[1-3]，其定義為：地表向半球空間反射的所有輻射能量與所有入射能量之比[4]。地表反照率作為能量平衡方程中的重要參量之一，反映了地球表面反射太陽輻射的能力，其時空變化受到自然過程(如降雪)以及人類活動(如森林破壞)等的影響，是全球環境變化的指示因子之一[5-10]。地球系統科學和氣候變化研究對多種尺度、長時間序列、空間連續且高質量的地表反照率數據有著迫切的應用需求，因此，對衛星遙感估算地表反照率開展研究具有重要科學意義。隨著對地觀測技術的快速發展，傳感器種類日趨豐富，衛星遙感方法為全球尺度、長時間序列地表反照率的估算提供了重要途徑[11-14]。

目前根據不同衛星遙感數據特點結合多種數學物理模型發展起來的地表反照率反演方法眾多，根據反演算法中使用模型的方式以及對遙感數據的利用方式，國際上主流的反照率生成算法主要包括基于反演二向性反射模型的反照率估算方法和直接估算反照率方法。為了支持全球長時間序列環境和氣候變化的研究，2012年GLASS團隊發布了第一版多種地表參數的長時間序列產品，2014年又更新到了第三版，其中的全球地表反照率產品具有時間序列長、時空連續一致的特點。GLASS全球反照率產品基于直接估算方法綜合利用AVHRR、MODIS等遙感數據獲取的目前全球最長時間序列的地表反照率產品[15-16]，已被應用在多個區域和全球尺度的研究中，包括：理解快速城鎮化過程中的區域輻射強迫，中國東北地區的森林擾動，半干旱內陸河流域的蒸散發估算，辨識氣候模式像元的空間變化及其影響[17]，驗證、標定和改進氣候模式中的模擬和參數化[18]，揭示全球森林的時空變冷和變暖效應[19]，估算日間凈輻射[20]等。

業務化運行的遙感產品生成中往往采用半經驗或經驗線性模型來兼顧反演的精度和效率，GLASS反照率產品算法便采用多元線性回歸的方法。近年來，機器學習方法在廣泛的研究和工業領域備受青睞，很多優秀的機器學習方法和框架涌現出來。這些機器學習方法依賴現在強大的計算能力可以建立大數據集上或簡單或復雜的各種回歸模型，強大的模型為傳統上復雜的遙感參數建模提供了一種新的可能。此外，對于傳統上較簡單的遙感參數建模，機器學習方法可以建立復雜度可控的模型來對反演精度調優。梯度提升決策樹(gradient boost decision tree，GBDT)是集成學習中boosting算法中基于梯度下降迭代的回歸樹算法，適用于回歸和分類問題。梯度提升決策樹是使用大量的簡單決策樹對訓練數據的不同特征建立模型，然后共同決策預測值。在回歸中使用平方誤差作為損失函數，每一棵決策樹迭代上一輪所有決策樹的預測結果并計算殘差，通過迭代加性訓練方法，在最小梯度上逐漸減小殘差值，從而較快地獲取較高的預測精度。梯度提升決策樹方法目前在遙感上的應用還較少，Fan等[21]使用梯度提升決策樹方法進行了空氣污染的時空預測。近年來，基于梯度提升決策樹方法的機器學習框架XGBoost[22]和LightGBM[23]以其快速高精度的特點在各大機器學習競賽和工業應用中備受關注，它們也給大數據量遙感參數快速生成提供了一種新的思路。本研究以基于MODIS反射率數據估算反照率模型為例，采用梯度提升決策樹算法，基于POLDER多角度觀測數據集實現該模型。

1 數據

POLDER-3是可以在全球進行BRDF觀測的星載傳感器，其利用視場重疊獲取多角度數據，單次過境時每個像元最多可有16個不同角度的觀測，單個像元一個月累積的觀測角度最多可達到300多個，基本可以實現全方位的角度觀測。POLDER-3數據地表類型覆蓋廣泛，本研究采用基于POLDER-3 BRDF數據集構建的MODIS波段BRDF數據集進行反照率反演建模。POLDER-3 BRDF數據集經過篩選和插值生成各角度網格的地表方向反射率，然后進行波段轉換得到對應于MODIS波段的植被、裸土和部分冰雪覆蓋地表的訓練數據集。模型模擬的冰雪BRDF數據集采用物理模型模擬不同參數下的純冰雪像元BRDF數據。地表寬波段反照率數據是先對BRDF進行半球積分獲得窄波段反照率，再通過窄波段反照率向寬波段反照率的轉換公式得到地表寬波段反照率。角度格網模擬數據按太陽天頂角、觀測天頂角和相對方位角生成：太陽天頂角范圍是0～80°，每2°間隔進行劃分，共分為41個間隔，格網中心分別為0°、2°、4°等。觀測天頂角范圍是0～64°，每2°間隔進行劃分，共分為33個間隔，格網中心分別為0°、2°、4°等。相對方位角范圍是0～180°，每5°間隔進行劃分，共分為37個間隔，格網中心分別為0°、5°、10°等。因此，角度格網在太陽/觀測角度空間共分成41×36×37=50 061個格網。

為了評估算法效果，利用地面站點觀測的反照率同衛星反演反照率進行分析評價，地面參考反照率數據來自北美AmeriFlux、SURFRAD和ARM等觀測網絡的28個站點，站點信息如表1所示。

表1 反照率地面驗證站點信息

2 方法

2.1 基于多元線性回歸的反照率估算方法

GLASS反照率估算采用基于MODIS地表反射率數據直接估算反照率的算法(AB1)，可以利用MODIS傳感器每天的地表方向反射率數據(已經過大氣校正)直接反演地表寬波段反照率，生成日地表反照率中間產品。AB1算法的核心思路是建立MODIS地表方向反射率與地表寬波段反照率之間分格網的多元線性回歸關系，首先采用POLDER-3/PARASOL BRDF數據集和模型模擬數據生成各種地表類型的MODIS地表方向反射率和地表寬波段反照率，然后進行分格網的回歸。GLASS反照率估算是在反射率角度格網上進行多元線性回歸，建立每個格網上MODIS地表方向反射率與地表寬波段反照率的回歸關系，即轉換系數，如式(1)、式(2)所示。

(1)

(2)

式中：αws是寬波段白空反照率；αbs(θs(k))是寬波段黑空反照率；θs(k)是太陽天頂角，取值0°～80°間隔5°，即k取值1，2，3，…，17；i取值對于MODIS數據而言為1，…，7，對于AVHRR數據而言為1，2，分別代表MODIS和AVHRR的短波窄波段；m0和n0(k)是回歸表達式的常數項；mi和ni(k)是表達式的回歸系數，回歸系數按照角度格網和地表類型建立；ρi(θs，θv，φ)是方向反射率。

2.2 基于梯度提升決策樹的反照率估算方法

不同于直接估算方法中采用的多元線性回歸方法，本研究采用擬合能力更強的梯度下降決策樹算法建立MODIS地表方向反射率與地表寬波段反照率的回歸關系。在角度格網的反照率直接估算方法中，算法基于大量樣本BRDF訓練數據集建立角度網格上的線性回歸關系。但實際上數據集可能存在非線性特性，此時使用非線性回歸算法代替線性回歸可以得到更加精確的結果，甚至可以建立角度連續的估算模型，從而緩解角度格網化離散降低精度的問題。這里以GLASS中基于MODIS的反照率估算為例，基于GLASS的多地表類型BRDF數據集采用機器學習方法建立MODIS地表方向反射率到寬波段反照率的估算模型。如果在所有角度上，對所有樣本進行訓練得到估算模型，相當于模型不依賴于角度格網查找表，可以直接估算任意角度的反照率。本文機器學習方法選擇的是LightGBM框架的梯度提升決策樹算法，模型的自我評價采用的是留余法，訓練和測試數據集都是采用“標記數據+輸入數據”的格式，標記數據為訓練數據中的寬波段反照率，輸入數據為窄波段反射率和角度等。

1)模型的構建策略。根據模型建立的粒度(即訓練數據集的粒度)，模型的構建包括以下三種尺度。

(1)在單一角度格網上，對單一地表類型數據訓練得到估算模型。這相當于還是依賴原有的格網查找表，但替換了原有的多元線性回歸建立的估算模型。

(2)在部分或所有角度上，對單一地表類型數據進行訓練得到估算模型。這相當于不依賴查找表，并且角度連續。

(3)在所有角度上，對所有地表類型數據進行訓練得到估算模型。這相當于不依賴于查找表，角度連續且隱含類型，但此時訓練數據集會異常龐大，一般機器難以運算。

2)訓練和預測數據集的構建。GLASS BRDF數據集的反射率數據是角度格網化記錄的，其中反照率數據首先存儲白空反照率，然后存儲黑空反照率數據，其中黑空反照率按照太陽天頂角遞增的順序排列，范圍是0°～80°，間隔為4°，短波波段、可見光波段和近紅外波段反照率依次存儲。在進行學習訓練時，每種反照率需要單獨構建訓練和預測數據集。

為了進行模型精度的自評價，采用了留余法。具體來說，首先將樣本按等間距抽樣劃分為10組，具體操作為將1，11，21，…，n抽樣為第1組，然后依次抽樣剩下9組。訓練和檢驗分10次進行，每次選擇其中9組作為訓練數據，剩余1組作為檢驗數據，10次精度評價結果的平均作為總體訓練的精度。

3)模型的訓練和預測應用。模型訓練選擇的是LightGBM工具和GBDT方法。為了能夠客觀評價模型的精度和泛化能力，模型訓練和預測是在相互獨立的訓練和預測數據集上進行的。考慮到模型訓練數據集大小的可操作性，模型的訓練主要在兩個層次上進行：一是分角度格網進行單獨訓練，二是按照角度格網的三個維度(太陽天頂角、觀測天頂角和觀測方位角)分別進行訓練。模型的預測應用首先讀取MODIS反射率數據集構建預測數據集，然后進行預測獲取反照率。

3 結果與討論

分別基于訓練數據集采用MLR方法和GBDT方法建立模型，使用預測數據集的輸入估算反照率，并以預測數據集的標記反照率為參考計算RMSE和R2，評價兩種方法的精度。

3.1 分格網訓練結果評價

以太陽主平面觀測(RAA=0°)時的33個觀測天頂角乘以41個太陽天頂角為例，分別使用多元線性回歸(簡稱MLR)和GBDT方法，在每個格網尺度上建立模型并統計模型精度情況在角度格網上的分布，如圖1所示。其中圖1(a)和圖1(b)分別是MLR和GBDT方法估算反照率的RMSE，圖1(c)和圖1(d)分別是MLR和GBDT方法估算反照率的R2，圖中橫軸是平面的太陽天頂角維，縱軸是平面的觀測天頂角維。可以看出，太陽天頂角和觀測天頂角均較大的格網相比于其他角度格網，兩種方法估算的誤差均較大、相關性較差，其中MLR方法估算的RMSE高達0.05，可決系數低于0.5，但GBDT方法在大角度時的精度衰減相比于其他角度不太顯著。同時在太陽天頂角等于觀測天頂角時的格網，也存在類似的精度衰減情況。在大部分的角度格網上，GBDT方法都比MLR方法具有更低的RMSE和更高的R2。

圖1 太陽主平面下角度格網的精度分布比較

圖2 太陽主平面下所有角度格網的RMSE和R2頻率分布

進一步，對RMSE按照0.01間隔，對R2按照0.1間隔統計MLR方法和GBDT方法的RMSE和R2分布。圖2給出了太陽主平面下所有角度格網的RMSE和R2的頻率分布。在RMSE小于0.02的格網統計中，GBDT方法格網占比94.5%，而MLR方法格網僅占69.6%；在RMSE小于0.01的格網統計上，GBDT方法格網占比近半數47.8%，而MLR方法格網僅占24%，GBDT方法相比MLR方法有比較明顯的精度提升。在R2大于0.9的格網統計中，GBDT方法格網占比97.2%，MLR方法格網占比84.9%，GBDT方法有更好的相關性。

3.2 單類不同粒度格網組合訓練結果評價

不同于原來的MLR方法，GBDT方法可以在不同的數據集粒度上建立統一模型，因此在單一角度格網、單一平面單一觀測天頂角、單一平面和多平面基礎上分別建立GBDT模型并進行精度對比。圖3給出了不同粒度下建立的GBDT模型與格網MLR方法精度比較的結果，圖中估算反照率為兩種方法估算的反照率，參考反照率為預測數據集給出的參考反照率。四種粒度數據集包括：(a)RAA=0，VZA=0，SZA=0，也即單一角度格網；(b)RAA=0，VZA=0，也即單一平面單一觀測天頂角，太陽主平面天底觀測；(c)RAA=0，也即單一平面，太陽主平面；(d)RAA=0～50，也即多平面。在較多數據集的統計結果上，MLR方法的RMSE為0.017到0.02，GBDT方法的RMSE為0.009到0.013；MLR方法的R2為0.934到0.956，GBDT方法的R2為0.972到0.987。可以看出，在四種粒度的數據集上建立的GBDT方法無論在RMSE還是R2上均優于MLR方法。

學習算法的精度依賴于模型對于數據集的表達能力。從圖3可以看出，單一平面訓練數據集GBDT的結果(圖3(c2)，樣本數量超過4 500萬)和多平面訓練數據集GBDT的結果(圖3(d2)，樣本數量超過5億)相比單一平面單一觀測天頂角數據集機器學習的結果(圖3(b2))，模型的精度有所衰減，這主要是由于模型對于數據集的表達能力不足導致的，因此應該增加學習模型的復雜度。多決策樹方法可以通過增加模型中樹的數量進而增加模型的復雜度，對單一和多平面數據集時，分別設置5 000棵樹代替原來設置的1 000棵樹重新訓練評價，兩種設置的估算結果如圖4所示。圖4(a)和圖4(b)分別是單一平面下1 000棵樹和5 000棵樹的估算結果，圖4(c)和圖4(d)分別是多平面下1 000棵樹和5 000棵樹的估算結果。可以看出，5 000棵樹模型相比于1 000棵樹模型，在單一平面數據集上將R2由0.972提升到0.979，將RMSE由0.013降到0.011；在多平面數據集上將R2由0.975提升到0.982，將RMSE由0.012降到0.01。可見，增加決策樹數量可以在大數據量訓練時提升估算模型的精度。

圖3 MLR和GBDT方法在不同粒度數據集上的精度比較

圖4 不同樹數量設置下的GBDT方法精度

3.3 基于地面站點觀測的算法精度評價

為了評價提出算法的效果，本文一方面從地面站點測量數據提取短波反照率作為參考值，另一方面收集了GLASS V3版本反照率產品并提取各站點的短波反照率，同所提出算法生產的短波反照率進行比較。考慮到GLASS V3版本的生產時間到2013年，因此選擇精度評價的時間范圍是2011—2013年。圖5給出了本文算法和GLASS V3反照率產品在北美28個站點的驗證結果，圖中GLASS新算法反照率代表本文算法反照率，GLASS V3算法反照率代表GLASS V3版本反照率。從圖中可以看出，本文算法相比GLASS V3版反照率精度有所提升，兩者RMSE分別為0.257和0.268，本文算法相對提升4%；兩者偏差分別為-0.001 6和0.004 0，本文算法在絕對偏差上相對提升60%。

圖5 不同算法反照率與地面實測反照率的驗證結果

4 結束語

地球系統科學和氣候變化研究對長時間序列、空間連續且高質量的地表反照率有著廣泛的應用需求，目前的對地觀測可以支持1980年至今的長時間序列全球衛星反照率和近年高時空分辨率高精度全球衛星反照率產品的生成。本文提出了基于GBDT算法的地表反照率估算方法，并在分格網和格網組合兩個尺度分別構建GBDT模型進行反照率估算。分格網的GBDT模型算法相比MLR精度有顯著提升，在RMSE小于0.02的格網統計中，GBDT方法格網占比94.5%，明顯高于MLR方法格網占比(69.6%)。在網格組合下的GBDT模型算法相比MLR精度同樣有顯著提升，統計的MLR方法平均RMSE為0.017到0.02，GBDT方法平均RMSE為0.009到0.013。同時，研究發現，提升GBDT模型的復雜度(決策樹樹木的數量)能小幅提升算法精度。基于地面站點觀測數據進一步的評價可知，相比于GLASS V3反照率產品，本文方法在RMSE和絕對偏差上相對提升4%和60%。研究表明，新型機器學習方法在優化遙感經驗或半經驗模型中具有重要潛力。