基于集成學(xué)習(xí)方法的CLDAS土壤濕度降尺度研究

韓慧敏 沈潤(rùn)平 黃安奇 狄文麗

0 引言

土壤水分參與地-氣水分和能量交換，從而對(duì)作物生長(zhǎng)、流域水文過(guò)程和氣候變化產(chǎn)生重要影響[1-3]，準(zhǔn)確獲取土壤濕度時(shí)空變化分布信息具有重要意義.近年來(lái)，多套基于陸面模式發(fā)展起來(lái)的土壤濕度陸面同化系統(tǒng)，包括全球陸面數(shù)據(jù)同化系統(tǒng)(GLDAS,0.25°×0.25°)[4]、北美陸面數(shù)據(jù)同化系統(tǒng)(NLDAS,0.125°×0.125°)[5]和中國(guó)氣象局(CMA)陸面數(shù)據(jù)同化系統(tǒng)(CLDAS,0.062 5°×0.062 5°)[6].土壤濕度產(chǎn)品已被廣泛用于陸面分析和氣象業(yè)務(wù)，但是空間分辨率較低，極大地限制了其進(jìn)一步應(yīng)用，特別是在需要更高空間分辨率的精確農(nóng)業(yè)管理和干旱監(jiān)測(cè)領(lǐng)域.

為獲得更高空間分辨率數(shù)據(jù)，許多學(xué)者對(duì)土壤濕度的降尺度方法開展了研究，目前主要有基于衛(wèi)星遙感數(shù)據(jù)融合的方法、基于地理信息數(shù)據(jù)的方法以及基于模型的方法[7].基于衛(wèi)星遙感數(shù)據(jù)融合的方法主要使用遙感數(shù)據(jù)，而不依賴站點(diǎn)資料，適用于區(qū)域大尺度的土壤濕度降尺度研究，但受衛(wèi)星觀測(cè)時(shí)間限制和云層覆蓋的影響[8-9].基于地理信息數(shù)據(jù)的方法主要是利用地形、土壤質(zhì)地和植被覆蓋等參數(shù)，建立地理信息與土壤濕度之間的關(guān)系，獲得高分辨率土壤濕度，但需要大量的實(shí)地?cái)?shù)據(jù)來(lái)構(gòu)建地統(tǒng)計(jì)或分形插值模型，從而限制了在較大尺度區(qū)域內(nèi)的應(yīng)用[10].基于模型的降尺度方法主要包括數(shù)理統(tǒng)計(jì)模型(如基于地統(tǒng)計(jì)學(xué)、多重分形或小波)和陸面模型，該類方法應(yīng)用時(shí)需考慮模型關(guān)系的時(shí)空普適性，以及大量站點(diǎn)數(shù)據(jù)的輸入[11].隨著計(jì)算機(jī)性能的提高和人工智能技術(shù)的發(fā)展，機(jī)器學(xué)習(xí)方法被引入土壤濕度降尺度研究中，它能在缺乏連續(xù)數(shù)據(jù)的情況下，建立土壤濕度與陸表參數(shù)之間的關(guān)系，且具有較強(qiáng)的非線性問(wèn)題學(xué)習(xí)能力和整合多源數(shù)據(jù)的靈活性，成為提高土壤濕度空間分辨率的有效技術(shù).但單一機(jī)器學(xué)習(xí)方法易表現(xiàn)出對(duì)非線性及表征空間大的數(shù)據(jù)性能的不足，且易產(chǎn)生過(guò)擬合[12]，難以全面考慮土壤濕度變化特征，導(dǎo)致估算精度不高、模型魯棒性低等問(wèn)題.而集成學(xué)習(xí)方法能結(jié)合多種學(xué)習(xí)器的優(yōu)勢(shì)，具有更高的模型準(zhǔn)確性、魯棒性和整體歸納能力，目前在土壤濕度降尺度研究中還鮮有應(yīng)用.

CLDAS是國(guó)家氣象信息中心研發(fā)的中國(guó)氣象局陸面數(shù)據(jù)同化系統(tǒng)(CMA Land Data Assimilation System,CLDAS)，其土壤濕度產(chǎn)品是利用我國(guó)多種資料融合和同化獲得的大氣強(qiáng)迫數(shù)據(jù)，驅(qū)動(dòng)多種陸面模式模擬得到[13].該產(chǎn)品時(shí)空連續(xù)、不受天氣影響，在中國(guó)區(qū)域的表現(xiàn)優(yōu)于國(guó)際同類產(chǎn)品[14]，目前空間分辨率只有6 km.為此，本文利用梯度提升機(jī)、深度前饋神經(jīng)網(wǎng)絡(luò)、隨機(jī)森林以及Stacking集成學(xué)習(xí)方法，以華北地區(qū)為例，開展了對(duì)CLDAS土壤濕度產(chǎn)品進(jìn)行降尺度研究，將其空間分辨率降尺度至1 km，以獲得高時(shí)空分辨率連續(xù)的土壤濕度估算.

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

1.1 研究區(qū)概況

研究區(qū)位于我國(guó)華北地區(qū)，空間范圍為110°21′～122°43′E，31°23′～41°36′N，包括北京、天津、河北、河南和山東5個(gè)省(市)，占地約54萬(wàn)km2，是我國(guó)最主要的糧食產(chǎn)區(qū)，耕地面積大(面積占比71%)(圖1)，歷史上多次遭受重大干旱，土壤濕度能夠直接表征地表水分狀態(tài)，是關(guān)鍵地表參數(shù).

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

1.2.1 CLDAS土壤濕度數(shù)據(jù)

中國(guó)氣象局陸面數(shù)據(jù)同化系統(tǒng)(CLDAS-V2.0)土壤濕度資料是基于數(shù)據(jù)融合和同化技術(shù)，利用氣溫、氣壓、濕度、風(fēng)速、降水輻射數(shù)據(jù)和初始場(chǎng)信息，驅(qū)動(dòng)CLM和Noah-MP陸面模式集合模擬而獲得[15].CLDAS-V2.0土壤濕度數(shù)據(jù)垂直分為5層，分別為0～5、0～10、10～40、40～100、100～200 cm，單位為m3/m3，產(chǎn)品覆蓋東亞區(qū)域(60°～160°E，0°～65°N)，空間分辨率為0.062 5°，時(shí)間分辨率為逐小時(shí).將每24 h的土壤濕度值求平均，獲得逐日土壤濕度，研究利用2019年4月1日至2019年10月31日0～10 cm土層土壤濕度數(shù)據(jù)(中國(guó)氣象數(shù)據(jù)服務(wù)中心:http:∥data.cma.cn)開展研究.

1.2.2 遙感數(shù)據(jù)

采用Terra與Aqua衛(wèi)星的MODIS數(shù)據(jù)(來(lái)源于美國(guó)國(guó)家航空航天局(http:∥ladsweb.modaps.eosdis.nasa.gov))獲得的高分辨率陸面地表參量.時(shí)間范圍為2019年4—10月.

1)地表溫度數(shù)據(jù)：來(lái)自MODIS的MOD11A1和MYD11A1 V6產(chǎn)品，空間分辨率為1 km，時(shí)間分辨率為1 d.該數(shù)據(jù)采用廣義分裂窗算法反演得到，同時(shí)它能夠有效地消除大氣的影響[16].為保持土壤濕度數(shù)據(jù)與地表溫度數(shù)據(jù)的時(shí)間一致性，從數(shù)據(jù)集中提取白天地表溫度數(shù)據(jù)和夜晚地表溫度數(shù)據(jù)，經(jīng)過(guò)投影、系數(shù)轉(zhuǎn)換，并對(duì)同日次的白天和夜晚地表溫度數(shù)據(jù)進(jìn)行平均合成，得到逐日平均地表溫度數(shù)據(jù).

2)地表反照率數(shù)據(jù)：來(lái)自MODIS BRDF/ALBEDO系列MCD43A3產(chǎn)品，空間分辨率為500 m，時(shí)間分辨率為1 d，該數(shù)據(jù)是經(jīng)過(guò)雙向反射函數(shù)(BRDF)模型計(jì)算修正的反照率產(chǎn)品[17].本數(shù)據(jù)集包括7個(gè)窄波段和3個(gè)寬波段的黑空和白空反照率.白空反照率是漫射-半球反照率，反映了陰天條件下的地表反射狀況；黑空反照率能夠較為準(zhǔn)確地反映正午時(shí)刻地球表面對(duì)太陽(yáng)直射光線的反射情況[18].由于白空反照率和黑空反照率的平均值差異較小，且高度相關(guān)，實(shí)際地表反照率可選擇白空、黑空反照率的平均值計(jì)算獲得[19].

3)地表反射率數(shù)據(jù)：來(lái)自MODIS的MOD09A1和MYD09A1產(chǎn)品，空間分辨率為500 m，時(shí)間分辨率為8 d，共包含7個(gè)波段，其對(duì)應(yīng)的波長(zhǎng)范圍如表1所示.該數(shù)據(jù)是低觀測(cè)角度條件下，受云、云陰影及氣溶膠等影響最小的8 d日數(shù)據(jù)合成產(chǎn)品[20].歸一化差異水體指數(shù)(NDWI)計(jì)算公式[21]如下：

(1)

其中，ρ(858 nm)和ρ(1 240 nm)分別對(duì)應(yīng)反射率數(shù)據(jù)的2、5波段.

表1 反射率數(shù)據(jù)的波長(zhǎng)范圍

1.2.3 土壤質(zhì)地?cái)?shù)據(jù)

土壤質(zhì)地?cái)?shù)據(jù)來(lái)源于中國(guó)科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心(http:∥www.resdc.cn)，投影為Albers正軸等面積雙標(biāo)準(zhǔn)緯線圓錐投影，空間分辨率為1 km.該數(shù)據(jù)是根據(jù)1∶100萬(wàn)土壤類型圖和第二次土壤普查獲取到的土壤剖面數(shù)據(jù)編輯制作而成的，依據(jù)砂粒、粉粒、黏粒含量進(jìn)行土壤質(zhì)地劃分，將數(shù)據(jù)分為Sand(砂土)、Silt(粉砂土)與Clay(黏土)3大類，每一類數(shù)據(jù)均通過(guò)百分比來(lái)反映不同質(zhì)地顆粒的含量.

1.2.4 地形數(shù)據(jù)

數(shù)字高程模型(Digital Elevation Model,DEM)數(shù)據(jù)選用SRTM DEM數(shù)據(jù)，來(lái)源于中國(guó)科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心(http:∥www.resdc.cn).該數(shù)據(jù)是基于雷達(dá)測(cè)圖技術(shù)通過(guò)美國(guó)“奮進(jìn)”號(hào)航天飛機(jī)獲得的，涵蓋了60°N～56°S間陸地地表80%面積范圍，經(jīng)NASA噴氣推進(jìn)實(shí)驗(yàn)室處理完成[22].研究采用的數(shù)據(jù)是基于版本4.1的DEM數(shù)據(jù)，經(jīng)重采樣生成1 km全國(guó)數(shù)據(jù).

將以上遙感數(shù)據(jù)、土壤質(zhì)地?cái)?shù)據(jù)以及地形數(shù)據(jù)投影統(tǒng)一轉(zhuǎn)換至經(jīng)緯度投影，利用華北地區(qū)行政邊界裁剪等預(yù)處理，經(jīng)雙線性插值方法，對(duì)數(shù)據(jù)重采樣，分別生成分辨率6 km和1 km的兩種數(shù)據(jù)，6 km數(shù)據(jù)和CLDAS 6 km土壤濕度數(shù)據(jù)相匹配，用來(lái)訓(xùn)練降尺度模型，1 km數(shù)據(jù)作為降尺度模型輸入數(shù)據(jù)，用于估算高分辨土壤濕度.

1.2.5 站點(diǎn)觀測(cè)數(shù)據(jù)

站點(diǎn)數(shù)據(jù)來(lái)源于國(guó)家氣象信息中心資料服務(wù)室的2019年逐小時(shí)觀測(cè)資料[23].該站點(diǎn)土壤濕度觀測(cè)在垂直方向上分為8層：0～10 cm、10～20 cm、20～30 cm、30～40 cm、40～50 cm、50～60 cm、60～80 cm、80～100 cm.利用頻域反射技術(shù)(Frequency Domain Reflectometry,FDR)來(lái)測(cè)定土壤體積含水量[24].研究參考韓帥[25]提出的方法對(duì)土壤濕度觀測(cè)數(shù)據(jù)進(jìn)行質(zhì)量控制，得到有效站點(diǎn)223個(gè)，采用隨機(jī)方法，抽出185個(gè)作為訓(xùn)練站點(diǎn)數(shù)據(jù)用于建模，38個(gè)站點(diǎn)數(shù)據(jù)用于驗(yàn)證(圖1).

2 研究方法

2.1 降尺度因子的選擇

華北地區(qū)的土壤濕度在時(shí)間尺度上具有明顯的季節(jié)性，受到夏季風(fēng)的影響，夏季土壤濕度高，冬季土壤濕度低.在地域上中部區(qū)域地勢(shì)較平坦，土地利用類型以農(nóng)業(yè)用地為主，土壤濕度高；北部與西南部地形復(fù)雜，土地利用類型以草地和灌木叢為主，涵養(yǎng)水源的作用較差，土壤濕度低；南部和沿海地區(qū)粉砂粒和黏粒含量高、砂粒含量低，土壤濕度高.總體上除降水外，本地區(qū)土壤濕度變化受到溫度、地形、植被、土壤等因素的共同影響[26-27].地表溫度為監(jiān)測(cè)和降尺度土壤濕度中最重要的變量，土壤表層溫度發(fā)生變化，其內(nèi)在因素是土壤熱慣量，且隨著土壤濕度的增大，土壤熱慣量增大，因此，地表溫度與土壤濕度有密切的相關(guān)性[28].高程和坡度是影響土壤濕度空間分布的關(guān)鍵因素[26]，華北地區(qū)高程為-23～2 539 m，高程變化大，坡度為0°～22°，其中，94%的區(qū)域坡度為0°～5°，坡度變化小，所以引入高程作為降尺度因子之一.在可見(jiàn)光和近紅外區(qū)域，土壤濕度與植被光譜響應(yīng)之間存在顯著關(guān)系，研究表明，短波紅外區(qū)域?qū)ν寥罎穸鹊谋O(jiān)測(cè)效果更好[29-31].同時(shí)，短波紅外波段是植物葉片吸收水分的區(qū)域，植被反射率與葉片含水量呈負(fù)相關(guān)[32].因此，選擇基于短波紅外波段的歸一化差異水體指數(shù)，作為降尺度因子之一.土壤異質(zhì)性通過(guò)土壤質(zhì)地的變化，包括土壤顆粒和孔隙分布的變化，影響土壤濕度的分布.此外，地表反照率受到土壤顏色的影響，從而影響植被稀疏土壤的蒸發(fā)效率，進(jìn)而影響土壤濕度[33].因此，參照前人研究[34-38]，選用更高分辨率地表溫度、高程、歸一化差異水體指數(shù)、土壤質(zhì)地和地表反照率等對(duì)土壤水分較為敏感的因子指標(biāo)，作為降尺度因子來(lái)開展研究.

2.2 梯度提升機(jī)(GBM)

梯度提升機(jī)(Gradient Boosting Machine,GBM)是一種可以解決分類、回歸和重要性排序問(wèn)題的機(jī)器學(xué)習(xí)模型，是Boosting算法的典型代表.它遵循了集成學(xué)習(xí)的一種思想，即分多個(gè)階段迭代訓(xùn)練一系列可疊加的基學(xué)習(xí)器模型，在迭代的推進(jìn)過(guò)程中不斷進(jìn)行優(yōu)化和提升，使每一次新的迭代都是為了減少上一次迭代的殘差，使模型沿著殘差減少最快的方向進(jìn)行，由此產(chǎn)生一系列弱分類器，每個(gè)弱分類器都是一棵二叉樹，最終將這些弱分類器組合形成能使損失函數(shù)達(dá)到極小的模型[39].GBM對(duì)異常值和不平衡數(shù)據(jù)具有魯棒性，確保了高效的性能.其算法過(guò)程[40]如下：

1)輸入訓(xùn)練集數(shù)據(jù)T={(x1,y1),(x2,y2),…,(xn,yn)}，損失函數(shù)L(y,f(x)).

2)初始化模型f0(x)=arg min∑L(yk,c)，c為常量，是用來(lái)估計(jì)損失函數(shù)最小化的常數(shù)值.

3)迭代n=1,2,…,N(N為樣本數(shù))，k=1,2,3,…,K(K為基學(xué)習(xí)器的個(gè)數(shù))，計(jì)算rkn:

(2)

式中,rkn為損失函數(shù)負(fù)梯度在當(dāng)前模型的值，若損失函數(shù)已達(dá)到最小值，則進(jìn)行步驟5)，否則進(jìn)行下一步.

4)對(duì)rkn建立基學(xué)習(xí)器模型Tk(x)，對(duì)梯度提升進(jìn)行更新：

fm(x)=fm-1(x)+Tk(x).

(3)

5)得到強(qiáng)學(xué)習(xí)器：

(4)

梯度提升機(jī)方法基于陸地表面變量和土壤濕度之間的統(tǒng)計(jì)關(guān)系，降尺度過(guò)程主要涉及兩個(gè)階段：

1)訓(xùn)練.基于CLDAS土壤濕度數(shù)據(jù),與土壤濕度數(shù)據(jù)空間分辨率保持一致的陸表變量數(shù)據(jù)和站點(diǎn)數(shù)據(jù)建立梯度提升機(jī)回歸模型.

2)預(yù)測(cè).將高分辨率陸表變量數(shù)據(jù)輸入第1階段建立的回歸模型，以生成高分辨率土壤濕度數(shù)據(jù)(圖2).

圖2 梯度提升機(jī)土壤濕度降尺度結(jié)構(gòu)示意Fig.2 Downscaling of soil moisture based on gradient boosting machine

2.3 深度前饋神經(jīng)網(wǎng)絡(luò)(DFNN)

本研究采用深度前饋神經(jīng)網(wǎng)絡(luò)(Deep Feedforward Neural Network,DFNN)作為深度學(xué)習(xí)算法.在回歸任務(wù)中，該模型可以從大量變量中提取高級(jí)特征，以實(shí)現(xiàn)較高的預(yù)測(cè)精度[41].每一層的神經(jīng)元可以接收前一層神經(jīng)元的信號(hào)，并產(chǎn)生信號(hào)輸出到下一層.第0層叫做輸入層，最后一層叫做輸出層，其他中間層叫做隱藏層(圖3).這種網(wǎng)絡(luò)模型在各層之間具有全連接的神經(jīng)元結(jié)構(gòu)，其中隱藏層根據(jù)模型的復(fù)雜程度可以設(shè)計(jì)成任意數(shù)量的多層，各層之間的連接表示特征的權(quán)重，其中信息沒(méi)有反饋的從左向右傳輸.連接輸入和輸出的每層神經(jīng)元結(jié)構(gòu)具有以下映射關(guān)系[42]：

y=f(x,θ),

(5)

其中x和y分別代表輸入和輸出，θ表示已知輸入和期望輸出值之間映射的最優(yōu)參數(shù)解.為了避免深層網(wǎng)絡(luò)反向傳播可能帶來(lái)的梯度消失和梯度爆炸的問(wèn)題，各層的激活函數(shù)采用ReLU函數(shù).輸出層概率分布計(jì)算采用softmax函數(shù).假設(shè)輸入樣本表示為：X=(X1,X2,X3,X4,X5)，兩個(gè)隱藏層h1，h2的維度分別為H1，H2，則兩個(gè)隱藏層的輸出分別如式(6)、式(7)所示，輸出層的輸出如式(8)所示.W1,W2,W3為各連接層的連接權(quán)重矩陣，b1,b2,b3為各層偏移量.模型訓(xùn)練優(yōu)化的參數(shù)集合為θ={W1,W2,W3,b1,b2,b3}.

yh1=ReLU(W1X+b1),

(6)

yh2=ReLU(W2yh1+b2),

(7)

y=softmax(W3yh2+b3).

(8)

模型訓(xùn)練及預(yù)測(cè)過(guò)程如下：

1)變量因子標(biāo)準(zhǔn)化及輸入.本研究使用標(biāo)準(zhǔn)差標(biāo)準(zhǔn)化方法，如式(9)所示，將處理好的標(biāo)準(zhǔn)化6 km空間分辨率自變量與因變量因子輸入DFNN模型.

2)模型調(diào)參.需要通過(guò)調(diào)整模型參數(shù)以達(dá)到最好的訓(xùn)練效果，主要參數(shù)包括隱藏層層數(shù)L、隱藏層每層神經(jīng)元數(shù)量N和訓(xùn)練迭代次數(shù)epochs.

3)土壤濕度預(yù)測(cè).得到最優(yōu)訓(xùn)練模型后，將1 km空間分辨率的自變量因子輸入訓(xùn)練好的模型，得到1 km空間分辨率預(yù)測(cè)土壤濕度(圖3).

(9)

式中，xi為一列自變量中的第i個(gè)值，mean(x)和std(x)分別是x所在列自變量的均值和標(biāo)準(zhǔn)差.

2.4 隨機(jī)森林(RF)

隨機(jī)森林是Bagging算法的典型代表，作為一種增強(qiáng)型決策樹模型，可用于分類、回歸等任務(wù)[43].此外，模型對(duì)異常值不敏感，在樣本和變量的隨機(jī)化訓(xùn)練階段表現(xiàn)出良好的性能.與其他機(jī)器學(xué)習(xí)方法相比，隨機(jī)森林模型被廣泛應(yīng)用于微波土壤濕度產(chǎn)品降尺度[44].隨機(jī)森林模型的主要思想是基于回歸樹建立輸入變量與輸出土壤濕度之間的非線性函數(shù)[45]：

SSMO=fRF(C)+ε,

(10)

C=(LST,Albedo,DEM,soiltexture,NDWI,CLDAS-SSM),

(11)

其中:SSMO表示訓(xùn)練階段的實(shí)測(cè)土壤濕度值；C為輸入向量，表示輸入變量，包括地表溫度(LST)、地表反照率(Albedo)、高程(DEM)、土壤質(zhì)地(soil texture)、歸一化差異水體指數(shù)(NDWI)和CLDAS土壤濕度(CLDAS-SSM)；fRF是一個(gè)非線性函數(shù)，在輸入變量和輸出SSMO之間建立關(guān)系.

在本研究的回歸任務(wù)中，首先在訓(xùn)練期間內(nèi)建立若干棵決策樹，每棵決策樹由bootstrap樣本建立，其中訓(xùn)練輸入數(shù)據(jù)約占總樣本的2/3，其余(1/3)的樣本用于驗(yàn)證每棵樹.為了進(jìn)一步提高隨機(jī)森林模型的泛化能力，通過(guò)對(duì)許多獨(dú)立回歸樹的結(jié)果求算術(shù)平均來(lái)生成最終模型預(yù)測(cè)值，模型最終結(jié)果表示為

(12)

式中，p(SSMO|C)為最終預(yù)測(cè)結(jié)果，m是回歸樹的數(shù)量，Pi(SSMO|C)表示第i棵樹的預(yù)測(cè)結(jié)果.

2.5 Stacking集成學(xué)習(xí)方法

集成學(xué)習(xí)的優(yōu)勢(shì)在于將多個(gè)基學(xué)習(xí)器的結(jié)果進(jìn)行優(yōu)化組合輸出，以獲得比任意基學(xué)習(xí)器更好的結(jié)果[46].集成學(xué)習(xí)主要包括并行化集成的Bagging、序列化集成的Boosting以及堆疊式集成的Stacking等.以Bagging和Boosting為代表的集成方法可以對(duì)訓(xùn)練效果差的樣本賦以較高的權(quán)重進(jìn)行二次學(xué)習(xí)，提高組合預(yù)測(cè)的泛化能力.然而該類方法只能集成同類決策樹模型，難以融合其他模型的優(yōu)勢(shì)特性，不同算法間數(shù)據(jù)觀測(cè)的差異性難以體現(xiàn)[47].Stacking是一種分層模型集成框架(圖4)，首先調(diào)用不同類型的學(xué)習(xí)器對(duì)數(shù)據(jù)集進(jìn)行訓(xùn)練學(xué)習(xí)，將各學(xué)習(xí)器得到的訓(xùn)練結(jié)果組成一個(gè)新的訓(xùn)練樣例，作為元學(xué)習(xí)器的輸入，最終第2層模型中元學(xué)習(xí)器綜合多個(gè)基學(xué)習(xí)器的輸出特征，作出最后的決策[48].因此，Stacking增加了模型的準(zhǔn)確性、魯棒性和整體歸納能力.

圖4 Stacking集成學(xué)習(xí)結(jié)構(gòu)示意Fig.4 Structure of Stacking ensemble learning

廣義線性模型(Generalized Linear Models,GLM) 作為一般線性模型的擴(kuò)展，基本思想是通過(guò)概率分布函數(shù)模擬非線性過(guò)程，它具有清晰的變量權(quán)重結(jié)構(gòu)，能夠模擬非線性的響應(yīng)關(guān)系，模型不會(huì)出現(xiàn)明顯的過(guò)擬合，對(duì)每個(gè)基學(xué)習(xí)器的效應(yīng)產(chǎn)生清晰的認(rèn)識(shí)[49].GLM模型主要通過(guò)連接函數(shù)，建立響應(yīng)變量的數(shù)學(xué)期望值，及其與代表線性組合的預(yù)測(cè)變量間的關(guān)系.一個(gè)廣義線性模型包括隨機(jī)成分、系統(tǒng)成分和連接函數(shù)三部分：

(13)

ni=β0+β1xi1+β2xi2+…+βpxip,

(14)

E(Yi)=μi=g-1(ni).

(15)

隨機(jī)成分即因變量的概率分布，其中a(φi)，b(θi)，c(yi,φi)為已知的函數(shù)；系統(tǒng)成分即自變量的線性組合.連接函數(shù)建立了隨機(jī)成分與系統(tǒng)成分之間的特定關(guān)系.

本文基于Stacking集成學(xué)習(xí)的土壤濕度降尺度模型框架(圖5)，算法步驟如下：

1)輸入原始數(shù)據(jù)集T，即包括CLDAS土壤濕度數(shù)據(jù)、與土壤濕度數(shù)據(jù)空間分辨率保持一致的陸表變量數(shù)據(jù)和站點(diǎn)數(shù)據(jù)，并按照3∶1的比例隨機(jī)劃分訓(xùn)練集T1和測(cè)試集T2，T=T1∪T2，T1∩T2=?.

2)學(xué)習(xí)并生成新的數(shù)據(jù)集.第1層包含3種基學(xué)習(xí)器：梯度提升機(jī)、深度前饋神經(jīng)網(wǎng)絡(luò)和隨機(jī)森林，采用K折交叉驗(yàn)證來(lái)訓(xùn)練第1層模型，3種模型擴(kuò)展之后生成第2層訓(xùn)練集T′1.在基學(xué)習(xí)器模型進(jìn)行K折交叉驗(yàn)證過(guò)程中，對(duì)測(cè)試集K次計(jì)算結(jié)果求平均，3種模型擴(kuò)展后生成第2層測(cè)試集T′2.T′1和T′2構(gòu)成新的數(shù)據(jù)集T′.

3)將得到的T′1用于訓(xùn)練第2層元學(xué)習(xí)器廣義線性模型，并用T′2驗(yàn)證模型性能.訓(xùn)練得到最終土壤濕度降尺度模型.

圖5 Stacking集成學(xué)習(xí)方法土壤濕度降尺度結(jié)構(gòu)示意Fig.5 Downscaling of soil moisture based on Stacking ensemble learning

2.6 評(píng)價(jià)指標(biāo)

本研究采用相關(guān)系數(shù)(R)、偏差(Bias，其量值記為B)和均方根誤差(RMSE)3個(gè)指標(biāo)定量地分析原始CLDAS土壤濕度和降尺度土壤濕度，計(jì)算公式[50]如下：

(16)

(17)

(18)

精度評(píng)估過(guò)程中，采用雙線性內(nèi)插法，將CLDAS土壤濕度和不同方法降尺度結(jié)果內(nèi)插到站點(diǎn)，然后與站點(diǎn)觀測(cè)數(shù)據(jù)進(jìn)行分析，計(jì)算相關(guān)系數(shù)、均方根誤差和偏差.

采用決定系數(shù)(R2)和均方根誤差(RMSE)對(duì)深度學(xué)習(xí)模型預(yù)測(cè)結(jié)果進(jìn)行精度判定.均方根誤差如式(18)所示，決定系數(shù)計(jì)算公式[51]如下：

(19)

3 結(jié)果與分析

3.1 深度前饋神經(jīng)網(wǎng)絡(luò)參數(shù)調(diào)優(yōu)

研究使用決定系數(shù)R2和均方根誤差RMSE來(lái)表征深度前饋神經(jīng)網(wǎng)絡(luò)的擬合效果.通過(guò)調(diào)整隱藏層層數(shù)L、神經(jīng)元數(shù)量N和迭代次數(shù)epochs，得到模型訓(xùn)練集和測(cè)試集的決定系數(shù)R2和均方根誤差RMSE.當(dāng)模型包含過(guò)多參數(shù)時(shí)，為了避免模型結(jié)果過(guò)擬合，提高模型的泛化能力，有必要同時(shí)考慮訓(xùn)練集和測(cè)試集具有最高的決定系數(shù)和最低的均方根誤差.結(jié)果表明(表2)，固定神經(jīng)元數(shù)量N=400，逐漸增加隱藏層層數(shù)L或迭代次數(shù)epochs時(shí)，總體上，訓(xùn)練集和測(cè)試集的決定系數(shù)逐漸增大，均方根誤差逐漸降低.當(dāng)隱藏層層數(shù)L=5，迭代次數(shù)epochs=600時(shí)，訓(xùn)練集和測(cè)試集的決定系數(shù)分別為0.936和0.830，均方根誤差分別為0.018和0.030 1 m3/m3，此時(shí)再逐漸增加隱藏層層數(shù)L或迭代次數(shù)epochs，訓(xùn)練集和測(cè)試集的決定系數(shù)逐漸降低，均方根誤差逐漸增大.因此，在本研究中最終選擇以下數(shù)值作為模型的初始輸入?yún)?shù)：隱藏層層數(shù)L=5，神經(jīng)元數(shù)量N=400，迭代次數(shù)epochs=600.

表2 模型參數(shù)調(diào)整結(jié)果

3.2 不同方法降尺度結(jié)果空間分布

為了比較不同方法降尺度效果，研究比較了梯度提升機(jī)(GBM)、深度前饋神經(jīng)網(wǎng)絡(luò)(DFNN)、隨機(jī)森林(RF)和Stacking集成學(xué)習(xí)等4種方法，因土壤溫度小于0 ℃時(shí)，觀測(cè)儀器異常導(dǎo)致缺少可靠的觀測(cè)數(shù)據(jù)，研究?jī)H分析2019年4—10月的CLDAS土壤濕度產(chǎn)品降尺度效果.從降尺度結(jié)果和原土壤濕度日均值分布(圖6)可以看出，總體上，華北地區(qū)的南部和沿海區(qū)域土壤濕度較高，中部和北部土壤濕度較低，平均土壤濕度均達(dá)到0.2 m3/m3以上，降尺度前后產(chǎn)品均較好地反映出此變化規(guī)律.但降尺度前土壤濕度日均值的平均值為0.226 m3/m3，降尺度后土壤濕度有所降低，日均值的平均值分別為：GBM(0.208 m3/m3)>DFNN (0.207 m3/m3)>RF(0.207 m3/m3)>Stacking (0.206 m3/m3).特別是華北地區(qū)的南部和沿海區(qū)域降尺度后降低明顯，降尺度前土壤濕度日均值大于0.25 m3/m3，降尺度后介于0.2～0.3 m3/m3之間，并且降尺度后土壤濕度的空間分布細(xì)節(jié)更加豐富，這與降尺度過(guò)程中融合了對(duì)土壤濕度影響較大的高分辨率地表溫度、地表反照率、地形等地表變量數(shù)據(jù)有關(guān).

3.3 不同方法降尺度效果精度分析

利用站點(diǎn)觀測(cè)數(shù)據(jù)逐日土壤濕度評(píng)估表明(圖7和圖8)：不同降尺度方法預(yù)測(cè)的降尺度結(jié)果與站點(diǎn)觀測(cè)土壤濕度之間存在顯著的相關(guān)性，相關(guān)系數(shù)介于0.699 6～0.756 8，高于原土壤濕度的相關(guān)系數(shù)0.651 1；降尺度后均方根誤差介于0.050 5～0.055 3 m3/m3之間，低于原土壤濕度均方根誤差0.062 3 m3/m3；降尺度后偏差介于-0.008 1～-0.005 2 m3/m3之間，比原土壤濕度偏差0.023 9 m3/m3更接近于0，說(shuō)明降尺度后土壤濕度精度得到提高，相比于原土壤濕度，降尺度結(jié)果更接近于實(shí)測(cè)值.

對(duì)比4種不同的降尺度方法結(jié)果的相關(guān)系數(shù)表明(圖8)，Stacking集成學(xué)習(xí)土壤濕度和RF土壤濕度相關(guān)系數(shù)較高，分別為0.756 8和0.740 2，其次為DFNN和GBM土壤濕度，分別為0.7155和0.699 6，且Stacking集成學(xué)習(xí)土壤濕度和RF土壤濕度均方根誤差較低，分別為0.050 5 m3/m3和0.050 9 m3/m3，其次為GBM土壤濕度和DFNN土壤濕度，分別為0.055 3 m3/m3和0.055 4 m3/m3.4種不同降尺度方法結(jié)果偏差均在0以下，存在一定程度上的低估，但絕對(duì)偏差小于CLDAS產(chǎn)品，一定程度上改善了其高估現(xiàn)象，以Stacking集成學(xué)習(xí)方法低估程度最小，土壤濕度絕對(duì)偏差為0.005 2 m3/m3.因此，相對(duì)來(lái)說(shuō)，Stacking集成學(xué)習(xí)方法優(yōu)于其他3種方法，其相關(guān)系數(shù)最高，均方根誤差和偏差相對(duì)較小，這與Stacking集成學(xué)習(xí)方法能夠更好地挖掘出輸入變量和土壤濕度的相關(guān)性，提升模型擬合效果有關(guān).

3.4 不同方法降尺度結(jié)果時(shí)間序列及誤差分析

圖6 不同方法結(jié)果土壤濕度日均值空間分布 a.CLDAS;b.GBM;c.DFNN;d.RF;e.Stacking集成學(xué)習(xí)Fig.6 Spatial distribution of original and downscaled daily average soil moisture a.CLDAS;b.GBM;c.DFNN;d.RF;e.Stacking ensemble learning

從不同方法降尺度結(jié)果0～10 cm土層土壤濕度時(shí)間序列來(lái)看(圖9)，4種降尺度方法的降尺度結(jié)果和原土壤濕度的變化趨勢(shì)與觀測(cè)值總體上相似，均能反映出土壤濕度隨時(shí)間變化的規(guī)律，但大多數(shù)日次原土壤濕度存在高估現(xiàn)象，在整個(gè)時(shí)間段內(nèi)，比觀測(cè)值高0.013 7 m3/m3.4種降尺度結(jié)果則存在一定程度上的低估，尤其在91～98、180～190、200～210和240～270日次，但整體上降尺度結(jié)果和觀測(cè)值的曲線更為接近.按其與觀測(cè)值偏離程度大小依次為：GBM>DFNN>RF>Stacking集成學(xué)習(xí)方法，其中，Stacking集成學(xué)習(xí)土壤濕度比觀測(cè)值低大約0.005 2 m3/m3，整體上，Stacking集成學(xué)習(xí)降尺度土壤濕度與觀測(cè)值曲線趨勢(shì)變化更相吻合，更接近站點(diǎn)觀測(cè)數(shù)據(jù).

圖7 CLDAS逐日土壤濕度精度站點(diǎn)評(píng)估散點(diǎn)圖Fig.7 Scatter plot of CLDAS daily soil moisture and observations

從土壤濕度與觀測(cè)值的相關(guān)系數(shù)和誤差分析(圖10)來(lái)看，與原土壤濕度相比，4種降尺度方法降尺度結(jié)果相關(guān)系數(shù)均有所提升，其中，Stacking集成學(xué)習(xí)土壤濕度相關(guān)系數(shù)提高較大，平均提高0.13，其次為RF土壤濕度，平均提高0.12.與原土壤濕度相比，4種不同降尺度方法降尺度結(jié)果的均方根誤差均有明顯的降低，其中，Stacking集成學(xué)習(xí)土壤濕度和RF土壤濕度比原土壤濕度的均方根誤差，分別平均降低了0.012 1和0.011 7 m3/m3，更接近于觀測(cè)值.4種方法降尺度結(jié)果的偏差絕大多數(shù)日次在0以下，即降尺度結(jié)果存在一定程度的低估現(xiàn)像.原土壤濕度偏差均值為0.023 9 m3/m3、GBM為-0.008 1 m3/m3、DFNN為-0.005 8 m3/m3、RF為-0.005 6 m3/m3、Stacking集成學(xué)習(xí)方法為-0.005 2 m3/m3，4種降尺度結(jié)果均改善了原土壤濕度的高估問(wèn)題，其中，以Stacking集成學(xué)習(xí)方法最優(yōu).

圖8 不同方法降尺度結(jié)果站點(diǎn)精度評(píng)估散點(diǎn)圖 a.GBM;b.DFNN;c.RF;d.Stacking集成學(xué)習(xí)Fig.8 Scatter plots of downscaled soil moistures and observations a.GBM;b.DFNN;c.RF;d.Stacking ensemble learning

圖9 不同方法降尺度結(jié)果0～10 cm土層土壤濕度時(shí)間序列Fig.9 Time series of 0-10 cm soil moisture downscaled by different methods

圖10 不同方法降尺度日均土壤濕度與觀測(cè)值的相關(guān)系數(shù)與誤差Fig.10 Correlation coefficients (up),RMSEs (middle),and biases (down) of downscaled daily average soil moistures compared with observations

表3為4種不同降尺度方法的降尺度結(jié)果與觀測(cè)值在月尺度上的相關(guān)系數(shù)和誤差.4種方法降尺度結(jié)果與觀測(cè)值的相關(guān)系數(shù)達(dá)到0.45以上，其中，4月、6月、7月和8月相關(guān)系數(shù)較大，均達(dá)到0.7以上，9月較小，介于0.5～0.6之間，這可能與夏季降水頻繁，土壤水分具有較強(qiáng)的空間異質(zhì)性，增加了土壤濕度估算的不確定性有關(guān).整體來(lái)看，平均相關(guān)系數(shù)均大于0.65，其中，Stacking集成學(xué)習(xí)方法最高.除9月外，4種方法降尺度結(jié)果與觀測(cè)值在月尺度上的絕對(duì)偏差均小于0.01 m3/m3，其中，4月絕對(duì)偏差最小.在整個(gè)時(shí)間段內(nèi)，除RF土壤濕度在4月偏差為正值，其余各月偏差均為負(fù)值，各月土壤濕度估算值均低于觀測(cè)值.絕對(duì)偏差均值大小依次為：GBM>DFNN>RF>Stacking集成學(xué)習(xí).4種方法降尺度結(jié)果與觀測(cè)值在月尺度上的均方根誤差介于0.043～0.067 m3/m3之間，其中，9月均方根誤差較大，各方法降尺度結(jié)果的均方根誤差均大于0.059 m3/m3，4月均方根誤差較小，各方法降尺度結(jié)果的均方根誤差介于0.043～0.050 m3/m3之間，以Stacking集成學(xué)習(xí)方法最小.因此，Stacking集成學(xué)習(xí)方法在月尺度上優(yōu)于其他方法.

4 結(jié)論與討論

研究以華北地區(qū)為例，以地表溫度、地表反照率、土壤質(zhì)地、高程、歸一化差異水體指數(shù)以及站點(diǎn)數(shù)據(jù)作為建模數(shù)據(jù)，基于3種單一模型(梯度提升機(jī)、深度前饋神經(jīng)網(wǎng)絡(luò)和隨機(jī)森林)以及多模型Stacking集成學(xué)習(xí)的方法，開展了中國(guó)氣象局陸面數(shù)據(jù)同化系統(tǒng)(CLDAS-V2.0) 0～10 cm土層土壤濕度數(shù)據(jù)降尺度研究，使其空間分辨率從6 km降尺度至1 km，并以站點(diǎn)觀測(cè)數(shù)據(jù)對(duì)降尺度結(jié)果進(jìn)行精度分析，得到以下結(jié)論：

表3 不同方法降尺度月均土壤濕度與觀測(cè)值的相關(guān)系數(shù)與誤差

1)4種不同降尺度方法的降尺度結(jié)果和原土壤濕度在華北地區(qū)的空間分布具有相似規(guī)律，南部和沿海區(qū)域土壤濕度較高，中部和北部土壤濕度較低，平均土壤濕度均達(dá)到0.2 m3/m3以上.降尺度后土壤濕度日均值有所降低，特別是在華北地區(qū)的南部和沿海區(qū)域.

2)4種不同降尺度方法均有效提高了CLDAS土壤濕度產(chǎn)品的空間分辨率和精度，4種方法絕對(duì)偏差均小于CLDAS產(chǎn)品，一定程度上改善了高估現(xiàn)象.原土壤濕度與站點(diǎn)觀測(cè)數(shù)據(jù)的相關(guān)系數(shù)、均方根誤差和偏差分別為0.651 1、0.062 3 m3/m3和0.023 9 m3/m3，降尺度土壤濕度的精度高低依次是Stacking集成學(xué)習(xí)方法、隨機(jī)森林、深度前饋神經(jīng)網(wǎng)絡(luò)、梯度提升機(jī).Stacking集成學(xué)習(xí)降尺度方法估算精度最高，相關(guān)系數(shù)、均方根誤差和偏差分別為0.756 8、0.050 5 m3/m3和-0.005 2 m3/m3，梯度提升機(jī)估算效果較差，相關(guān)系數(shù)、均方根誤差和偏差分別為0.699 6、0.055 3 m3/m3和-0.008 1 m3/m3.

3)原土壤濕度和4種不同降尺度方法的降尺度結(jié)果均能較好地體現(xiàn)土壤濕度的日變化特征，但大多數(shù)日次原土壤濕度存在高估現(xiàn)象，4種降尺度結(jié)果存在一定程度上的低估.整體上，降尺度結(jié)果和觀測(cè)值的曲線更為接近，其中，Stacking集成學(xué)習(xí)方法最優(yōu)，與原土壤濕度相比，相關(guān)系數(shù)平均提高0.13，均方根誤差平均降低0.012 1 m3/m3，偏差均值為-0.005 2 m3/m3.月尺度結(jié)果中，4種不同降尺度方法的降尺度結(jié)果在9月相關(guān)系數(shù)均較小，均方根誤差和偏差較大，整體來(lái)看，與其他3種方法相比，Stacking集成學(xué)習(xí)土壤濕度在月尺度上相關(guān)系數(shù)較大，均方根誤差和偏差較低，估算的土壤濕度精度較高.

本文利用梯度提升機(jī)、深度前饋神經(jīng)網(wǎng)絡(luò)、隨機(jī)森林和Stacking集成學(xué)習(xí)方法，對(duì)CLDAS 0～10 cm土層土壤濕度產(chǎn)品開展了降尺度研究，使其空間分辨率從6 km降尺度至1 km，且精度有所提高.但本研究選取的土壤濕度數(shù)據(jù)時(shí)間范圍為2019年4—10月，即春夏秋三季，未考慮冬季降尺度模型的土壤濕度估算表現(xiàn)，主要是因?yàn)楫?dāng)前土壤濕度測(cè)量?jī)x器在土壤含有冰水混合物時(shí)，測(cè)量結(jié)果存在誤差與不確定性，冬季的降尺度結(jié)果有待進(jìn)一步研究.另外，與土壤濕度相關(guān)的降尺度因子數(shù)據(jù)的質(zhì)量對(duì)建模效果有較大影響，高質(zhì)量的數(shù)據(jù)會(huì)提高降尺度結(jié)果的精度.