劉 堃，何祺勝，荊琛琳，李金陽，陳 麗

(河海大學(xué)地球科學(xué)與工程學(xué)院,江蘇 南京 211100)

蒸散量(ET)作為地表能量平衡和水文循環(huán)的重要組成部分[1]，決定著地表生態(tài)系統(tǒng)的水熱傳輸，其準(zhǔn)確的測定和估算，對于區(qū)域水資源評定、水分利用效率、干旱預(yù)測等均具有十分重要的價(jià)值。經(jīng)過20多年的發(fā)展，渦動(dòng)相關(guān)儀已經(jīng)廣泛應(yīng)用于地表和大氣的能量交換觀測[2-3]。渦動(dòng)相關(guān)儀的觀測數(shù)據(jù)通常是以30 min為一個(gè)周期，采集1 d和1 a的通量數(shù)據(jù)[4]。然而，在觀測過程中，由于各種原因(降雨、儀器故障、人為誤操作等)部分觀測數(shù)據(jù)會缺失，研究表明一年中有17%～50%的觀測數(shù)據(jù)會缺失和被剔除[5]。因此，如何建立有效的數(shù)據(jù)插補(bǔ)方法來形成完整的通量數(shù)據(jù)集成為當(dāng)前亟待解決的問題[4]。國內(nèi)外學(xué)者針對渦動(dòng)相關(guān)儀缺失的蒸散量提出了多種數(shù)據(jù)插補(bǔ)方法[6]，主要有非線性回歸、動(dòng)態(tài)線性回歸、查找表、晝夜平均、卡爾曼濾波、人工神經(jīng)網(wǎng)絡(luò)(ANN)方法等。

目前，機(jī)器學(xué)習(xí)算法已經(jīng)廣泛應(yīng)用于遙感研究。對于蒸散量的插補(bǔ)，國內(nèi)外學(xué)者主要采用ANN方法進(jìn)行插補(bǔ)[7-9]，并取得了較好的精度。但是使用其他機(jī)器學(xué)習(xí)算法(例如決策樹、隨機(jī)森林等)進(jìn)行插補(bǔ)的研究卻鮮有報(bào)道。此外，土壤濕度作為蒸散發(fā)的重要因子，對許多生態(tài)系統(tǒng)的蒸散發(fā)都有不同程度的影響[10]，而有關(guān)土壤水分要素的參與對基于機(jī)器學(xué)習(xí)的蒸散量插補(bǔ)結(jié)果的影響研究很少。

本文以典型干旱半干旱區(qū)濕地、農(nóng)田、草地、林地生態(tài)系統(tǒng)為研究對象，研究作為輸入變量的氣象因子與蒸散量觀測值之間的相關(guān)性，分析多元線性回歸(MLR)、決策樹(CART)、隨機(jī)森林(RF)、支持向量回歸(SVR)、BP人工神經(jīng)網(wǎng)絡(luò)(BPANN)、深度學(xué)習(xí)算法(DL)對缺失蒸散量的插補(bǔ)效果，以及土壤水分的參與對機(jī)器學(xué)習(xí)算法插補(bǔ)精度的影響，以期能為不同生態(tài)系統(tǒng)渦動(dòng)相關(guān)儀觀測蒸散量的插補(bǔ)提供理論支持。

1 數(shù)據(jù)和方法

1.1 數(shù)據(jù)選擇

黑河流域位于中國西北干旱半干旱地區(qū)(97.1°E～102.1°E，37.7°N～42.7°N)，是中國第二大內(nèi)陸河流域，屬于典型的溫帶大陸性干旱氣候，氣候干燥，降水稀少。黑河流域上游祁連山是高寒半干旱區(qū)，中、下游分別為河西走廊和額濟(jì)納平原干旱區(qū)，東西和南北差異特征顯著，生態(tài)環(huán)境比較脆弱。本文主要研究黑河流域濕地、草地、農(nóng)田、林地(檉柳、胡楊林、混合林)等生態(tài)系統(tǒng)。研究所用數(shù)據(jù)主要來自6個(gè)站點(diǎn)(表1)，這6個(gè)站點(diǎn)均建有自動(dòng)氣象站和渦動(dòng)相關(guān)儀，具有蒸散量、顯熱通量、氣象因子等的長期連續(xù)觀測數(shù)據(jù)(來源于黑河生態(tài)水文遙感試驗(yàn)(HiWATER)[11-12]，可以從寒區(qū)旱區(qū)科學(xué)數(shù)據(jù)中心申請獲得(http://westdc.westgis.ac.cn))。每個(gè)站點(diǎn)渦動(dòng)相關(guān)儀輸出的都是30 min蒸散量均值，自動(dòng)氣象站輸出的為10 min氣象因子均值；由于儀器損壞、質(zhì)量控制等原因，實(shí)際蒸散量輸出值都有不同程度的缺失，數(shù)據(jù)缺失狀況見表1。

表1 站點(diǎn)基本信息

采用每個(gè)站點(diǎn)30 min的有效氣象因子和蒸散量觀測值研究機(jī)器學(xué)習(xí)算法的插補(bǔ)效果，將氣象因子作為輸入變量，蒸散量觀測值作為響應(yīng)變量。草地、農(nóng)田、林地生態(tài)系統(tǒng)用的氣象因子為凈輻射、溫度、土壤熱通量、相對濕度、風(fēng)速、土壤體積含水率(表層2 cm深度)。對于濕地生態(tài)系統(tǒng)，因?yàn)槿鄙賹?shí)測土壤體積含水量，同時(shí)考慮到濕地常年積水，土壤體積含水量對蒸散發(fā)的影響低于其他氣象因子[13]，因此選用凈輻射、溫度、土壤熱通量、相對濕度和風(fēng)速作為氣象因子。

1.2 數(shù)據(jù)預(yù)處理

由于輸入變量各不相同，也不在同一個(gè)數(shù)量級上，為了消除各變量數(shù)量級差異對插補(bǔ)結(jié)果的影響，同時(shí)也為了加快人工神經(jīng)網(wǎng)絡(luò)和深度學(xué)習(xí)的訓(xùn)練和收斂速度，需要對輸入變量做歸一化處理。歸一化公式如下：

(1)

式中：Y—— 輸入變量歸一化值；x—— 輸入變量；xmin、xmax—— 輸入變量的最小值和最大值。

1.3 6種機(jī)器學(xué)習(xí)算法

本文采用的機(jī)器學(xué)習(xí)算法包括MLR、CART、RF、SVR、BPANN、DL。

a. 多元線性回歸。該算法是指2個(gè)或2個(gè)以上的自變量組成線性表達(dá)式對應(yīng)一個(gè)因變量，通常用在自變量和因變量具有線性相關(guān)的條件下。

b. 決策樹。決策樹是以平方誤差最小化準(zhǔn)則來不斷地分裂節(jié)點(diǎn)，從而遞歸地構(gòu)建二叉決策樹，并最終以葉子節(jié)點(diǎn)的均值作為預(yù)測值。經(jīng)反復(fù)試錯(cuò)，該算法的參數(shù)中，樹的最大深度和節(jié)點(diǎn)樣本數(shù)分別設(shè)置為20和11，其他參數(shù)保持缺省值不變。

c. 隨機(jī)森林。隨機(jī)森林通過集成學(xué)習(xí)的思想將多棵決策樹集成[14-15]，再通過取所有樹預(yù)測均值得到RF的預(yù)測值 。其不易發(fā)生過擬合，并具有很好的魯棒性。經(jīng)反復(fù)試驗(yàn)，該算法參數(shù)中，決策樹個(gè)數(shù)、每顆樹的最大深度和節(jié)點(diǎn)樣本數(shù)分別設(shè)置為100、30和1，其他參數(shù)保持缺省值不變。

d. 支持向量回歸。支持向量回歸是一種建立在統(tǒng)計(jì)學(xué)習(xí)理論的機(jī)器學(xué)習(xí)方法，通過結(jié)構(gòu)風(fēng)險(xiǎn)最小原理取得全局的最優(yōu)解[16-18]。SVR速度快，泛化能力強(qiáng)，能夠解決小樣本和高維輸入空間的問題。經(jīng)測試，SVR選用RBF作為核函數(shù)，懲罰系數(shù)C設(shè)置為200，核函數(shù)參數(shù)G設(shè)置為25。

e. BP人工神經(jīng)網(wǎng)絡(luò)。BP人工神經(jīng)網(wǎng)絡(luò)是一種學(xué)習(xí)過程為信號正向傳播和誤差反向傳播的前饋神經(jīng)網(wǎng)絡(luò)，是目前應(yīng)用最廣泛的一種神經(jīng)網(wǎng)絡(luò)[19]。目前，國內(nèi)外主要采用3層神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)來插補(bǔ)蒸散量[7,9]，且從理論上說，單隱層BP網(wǎng)絡(luò)已經(jīng)有足夠的映射或逼近能力。因此，本文采用3層網(wǎng)絡(luò)結(jié)構(gòu)，分別由輸入層、隱藏層和輸出層組成。每層的激活函數(shù)均設(shè)為relu函數(shù),經(jīng)實(shí)際測試，把隱藏層單元個(gè)數(shù)設(shè)為200，此時(shí)插補(bǔ)結(jié)果最佳。

f. 深度學(xué)習(xí)。深度學(xué)習(xí)是在傳統(tǒng)人工神經(jīng)網(wǎng)絡(luò)的基礎(chǔ)上添加更多的隱藏層，以此建立復(fù)雜的非線性網(wǎng)絡(luò)結(jié)構(gòu)，從而在少量有限的樣本中挖掘出數(shù)據(jù)的本質(zhì)特征，用較少的參數(shù)建立復(fù)雜的函數(shù)。本文采用增加BP人工神經(jīng)網(wǎng)絡(luò)隱藏層數(shù)量的深度學(xué)習(xí)算法，把隱藏層設(shè)為4層，每層70個(gè)神經(jīng)元，每層的激活函數(shù)設(shè)為relu函數(shù)，并且為了防止過擬合，把參數(shù)L2(正則化項(xiàng))的值設(shè)為0.01。

上述算法中，除了深度學(xué)習(xí)是基于Python語言的keras模塊建立的模型，其他算法均是基于Python語言的skicit-learn模塊建立的模型。

1.4 插補(bǔ)模型的建立與性能評價(jià)

綜合考慮數(shù)據(jù)質(zhì)量和氣象數(shù)據(jù)的缺失程度，分別選取草地(2016年全年)，農(nóng)田、林地(均為2015年全年),濕地(2015-01-01至2015-09-25)生態(tài)系統(tǒng)中的10 min氣象數(shù)據(jù)和30 min蒸散量觀測值作為研究數(shù)據(jù)。因?yàn)闅庀髷?shù)據(jù)缺失很少，因此采用線性插值法補(bǔ)充完整[20]，同時(shí)為了與30 min蒸散量相對應(yīng)，把10 min氣象數(shù)據(jù)取均值處理為30 min氣象數(shù)據(jù)。選用各生態(tài)系統(tǒng)30 min蒸散量未缺失部分和對應(yīng)的30 min氣象數(shù)據(jù)作為訓(xùn)練和測試樣本數(shù)據(jù)的來源，其中訓(xùn)練樣本用于插補(bǔ)模型的訓(xùn)練，測試樣本作為人工制造的缺失值，通過訓(xùn)練好的模型在測試樣本上的模擬值與測試樣本的觀測值(真值)進(jìn)行比較，研究插補(bǔ)的效果。選用決定系數(shù)R2、均方根誤差RMSE、絕對平均誤差MAE評價(jià)測試樣本模擬值與觀測值之間的相關(guān)性與離散程度。

2 結(jié)果與分析

2.1 蒸散量與氣象因子的相關(guān)性分析

為了確定所選氣象因子對于蒸散量插補(bǔ)的有效性，選取濕地、草地、農(nóng)田、林地的輸入氣象因子與對應(yīng)的蒸散量觀測值進(jìn)行相關(guān)性分析。結(jié)果(表2)表明，在置信度水平為0.05和0.01時(shí)，濕地、草地、農(nóng)田、林地的蒸散量均與所選氣象因子存在顯著相關(guān)關(guān)系。因此，可以選用凈輻射、溫度、土壤熱通量、相對濕度、風(fēng)速、土壤體積含水率6個(gè)氣象因子來建立插補(bǔ)模型。

2.2 蒸散量模型插補(bǔ)結(jié)果

采用不同機(jī)器學(xué)習(xí)算法計(jì)算濕地、草地、農(nóng)田、林地在有無土壤體積含水率情況下蒸散量的插補(bǔ)值，其與觀測值之間的R2、RMSE和MAE見表3。

表2 蒸散量與氣象因子間的相關(guān)性

注：*是在0.05水平上雙側(cè)顯著相關(guān)；**是在0.01水平上雙側(cè)顯著相關(guān)。

表3 6種機(jī)器學(xué)習(xí)模型蒸散量的插補(bǔ)結(jié)果

注：有無土壤體積含水率指輸入變量是否含有相應(yīng)的值。

整體來看，各生態(tài)系統(tǒng)最好插補(bǔ)結(jié)果的決定系數(shù)都比較高(R2>0.83)，其中草地和農(nóng)田模擬結(jié)果最好，其次是濕地、檉柳和混合林，胡楊林模擬結(jié)果最差。

由表3可知，MLR在每種生態(tài)系統(tǒng)的插補(bǔ)均表現(xiàn)最差，這表明氣象因子與蒸散量之間是復(fù)雜的非線性關(guān)系，不適合用簡單的線性模型。其次，CART要好于線性回歸的插補(bǔ)結(jié)果，而相比MLR、CART，RF、SVR、BPANN、DL在每種生態(tài)系統(tǒng)的結(jié)果均更好。在濕地、草地、農(nóng)田、胡楊林生態(tài)系統(tǒng)中，RF、SVR 、BPANN、DL算法的結(jié)果差異并不是很顯著。檉柳、混合林生態(tài)系統(tǒng)，SVR的結(jié)果要稍差于RF、BPANN、DL的結(jié)果，可能是因?yàn)闄f柳、混合林的實(shí)測氣象數(shù)據(jù)和蒸散量包含一些噪聲，而SVR對數(shù)據(jù)噪聲比較敏感[21]，導(dǎo)致噪聲帶來的錯(cuò)誤被擴(kuò)大，從而影響了插補(bǔ)結(jié)果。總之，對于干旱半干旱區(qū)的濕地、草地、農(nóng)田、林地生態(tài)系統(tǒng)，RF、BPANN、DL、SVR均可以得到較理想的插補(bǔ)結(jié)果，但是SVR穩(wěn)定性稍差于其他3種方法。

相比無土壤體積含水率情況，在土壤體積含水率參與的情況下，農(nóng)田、檉柳、胡楊林、混合林的R2、RMSE、MAE均有顯著提高，但草地相比卻并沒有顯著提高。在蒸發(fā)旺盛的生長季(3—10月)，草地的土壤體積含水率大部分都在30%以上(圖1)，相比其他生態(tài)系統(tǒng)，其下墊面水分比較充足，這導(dǎo)致土壤水分對蒸散發(fā)的影響相對較弱；同時(shí)草地生態(tài)系統(tǒng)中土壤水分與蒸散發(fā)的相關(guān)性最弱(表2)，這些因素表明，在半小時(shí)尺度下草地生態(tài)系統(tǒng)的土壤水分對蒸散發(fā)的影響程度相對較弱。因此，土壤水分的加入并未使草地的插補(bǔ)精度有明顯提高。

圖1 林地、農(nóng)田、草地的日土壤體積含水率變化Fig.1 Change of daily soil moisture in farmland, woodland and grassland

為了研究插補(bǔ)模型的泛化能力，選用草地(2015年)、農(nóng)田(2016年)、檉柳(2016年)、胡楊林(2014年)、濕地(2016年)生態(tài)系統(tǒng)的氣象數(shù)據(jù)和蒸散量觀測值(氣象因子和時(shí)間尺度與訓(xùn)練模型保持一致)作為測試樣本，用已訓(xùn)練好的插補(bǔ)模型計(jì)算測試樣本的模擬值，并與測試樣本的真值進(jìn)行比較。從整體結(jié)果看，無論土壤水分是否參與，各生態(tài)系統(tǒng)在其他年份的插補(bǔ)精度均有一定程度的降低。其中草地插補(bǔ)精度降低的幅度最小(R2降低0.02～0.03)，最好模擬精度為R2=0.9～0.91，RMSE=31.817～32.446 W/m2，MAE=17.232～18.627 W/m2，而濕地、農(nóng)田、檉柳、胡楊林、混合林插補(bǔ)精度降低的幅度較大(R2降低0.1～0.2)。

2.3 插補(bǔ)效果分析

為了檢驗(yàn)機(jī)器學(xué)習(xí)算法對蒸散量的插補(bǔ)效果，以農(nóng)田為例，采用人工神經(jīng)網(wǎng)絡(luò)算法，分別選用有土壤體積含水率和無土壤體積含水率2種數(shù)據(jù)集，對2015年缺失蒸散量進(jìn)行插補(bǔ)，得到一年連續(xù)的蒸散量，并以6月8日為例分析插補(bǔ)效果(圖2中虛線框代表原始觀測值缺失部分)，同時(shí)將插補(bǔ)完整的蒸散量與凈輻射進(jìn)行比較和分析，結(jié)果見圖2～3。

圖2 蒸散量插補(bǔ)效果Fig.2 Gap-filling effect of latent heat flux

圖3 插補(bǔ)完整蒸散量與凈輻射日變化對比Fig.3 Comparison of interpolation completed latent heat flux to net radiation

由圖2可以看到有土壤水分參與的結(jié)果要優(yōu)于無土壤水分參與的結(jié)果，更貼近實(shí)測曲線的變化趨勢。凈輻射作為蒸散發(fā)的能量來源，是影響蒸散量變化的主要因子，從蒸散量與凈輻射的變化趨勢可對插補(bǔ)效果做定性分析。從圖3可以看到，兩者變化趨勢比較接近，有較好的相關(guān)關(guān)系。通過上述分析，表明插補(bǔ)結(jié)果比較合理可信。

3 結(jié) 語

MLR在各個(gè)生態(tài)系統(tǒng)的蒸散量插補(bǔ)精度均最差(R2=0.6～0.7)，CART次之(R2=0.78～0.9)，而RF、SVR、BPANN、DL的插補(bǔ)精度較高(R2=0.83～0.93)；在檉柳、混合林生態(tài)系統(tǒng)，相比RF、BPANN、DL，SVR的插補(bǔ)精度稍低，表明其穩(wěn)定性偏差。此外，土壤水分參與插補(bǔ)要比土壤水分不參與插補(bǔ)可以獲得更高的精度(R2提高了0.01～0.06)。以已建立的插補(bǔ)模型去插補(bǔ)其他年份的蒸散量，發(fā)現(xiàn)其精度有不同程度的下降。

綜合考慮模型的插補(bǔ)精度和穩(wěn)定性, 對干旱半干旱區(qū)濕地、草地、農(nóng)田、林地生態(tài)系統(tǒng),隨機(jī)森林(RF)、BP人工神經(jīng)網(wǎng)絡(luò)(BPANN)、深度學(xué)習(xí)(DL)更適合用來插補(bǔ)，可以取得較好的結(jié)果，同時(shí)加入土壤水分可以在一定程度上提升機(jī)器學(xué)習(xí)模型插補(bǔ)的精度。

本研究雖較全面地分析了MLR、CART、RF、SVR、BPANN、DL對干旱半干旱區(qū)生態(tài)系統(tǒng)缺失蒸散量的插補(bǔ)效果、但這些算法在濕潤區(qū)、干旱區(qū)等其他環(huán)境的生態(tài)系統(tǒng)是否適用，還需要進(jìn)一步的研究和驗(yàn)證。