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

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

0 引言

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

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

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

1 數(shù)據(jù)

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

為了評估算法效果，利用地面站點(diǎn)觀測的反照率同衛(wèi)星反演反照率進(jìn)行分析評價(jià)，地面參考反照率數(shù)據(jù)來自北美AmeriFlux、SURFRAD和ARM等觀測網(wǎng)絡(luò)的28個(gè)站點(diǎn)，站點(diǎn)信息如表1所示。

表1 反照率地面驗(yàn)證站點(diǎn)信息

2 方法

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

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

(1)

(2)

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

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

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

1)模型的構(gòu)建策略。根據(jù)模型建立的粒度(即訓(xùn)練數(shù)據(jù)集的粒度)，模型的構(gòu)建包括以下三種尺度。

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

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

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

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

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

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

3 結(jié)果與討論

分別基于訓(xùn)練數(shù)據(jù)集采用MLR方法和GBDT方法建立模型，使用預(yù)測數(shù)據(jù)集的輸入估算反照率，并以預(yù)測數(shù)據(jù)集的標(biāo)記反照率為參考計(jì)算RMSE和R2，評價(jià)兩種方法的精度。

3.1 分格網(wǎng)訓(xùn)練結(jié)果評價(jià)

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

圖1 太陽主平面下角度格網(wǎng)的精度分布比較

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

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

3.2 單類不同粒度格網(wǎng)組合訓(xùn)練結(jié)果評價(jià)

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

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

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

圖4 不同樹數(shù)量設(shè)置下的GBDT方法精度

3.3 基于地面站點(diǎn)觀測的算法精度評價(jià)

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

圖5 不同算法反照率與地面實(shí)測反照率的驗(yàn)證結(jié)果

4 結(jié)束語

地球系統(tǒng)科學(xué)和氣候變化研究對長時(shí)間序列、空間連續(xù)且高質(zhì)量的地表反照率有著廣泛的應(yīng)用需求，目前的對地觀測可以支持1980年至今的長時(shí)間序列全球衛(wèi)星反照率和近年高時(shí)空分辨率高精度全球衛(wèi)星反照率產(chǎn)品的生成。本文提出了基于GBDT算法的地表反照率估算方法，并在分格網(wǎng)和格網(wǎng)組合兩個(gè)尺度分別構(gòu)建GBDT模型進(jìn)行反照率估算。分格網(wǎng)的GBDT模型算法相比MLR精度有顯著提升，在RMSE小于0.02的格網(wǎng)統(tǒng)計(jì)中，GBDT方法格網(wǎng)占比94.5%，明顯高于MLR方法格網(wǎng)占比(69.6%)。在網(wǎng)格組合下的GBDT模型算法相比MLR精度同樣有顯著提升，統(tǒng)計(jì)的MLR方法平均RMSE為0.017到0.02，GBDT方法平均RMSE為0.009到0.013。同時(shí)，研究發(fā)現(xiàn)，提升GBDT模型的復(fù)雜度(決策樹樹木的數(shù)量)能小幅提升算法精度。基于地面站點(diǎn)觀測數(shù)據(jù)進(jìn)一步的評價(jià)可知，相比于GLASS V3反照率產(chǎn)品，本文方法在RMSE和絕對偏差上相對提升4%和60%。研究表明，新型機(jī)器學(xué)習(xí)方法在優(yōu)化遙感經(jīng)驗(yàn)或半經(jīng)驗(yàn)?zāi)Ｐ椭芯哂兄匾獫摿Α?/p>