韓慧敏 沈潤平 黃安奇 狄文麗

0 引言

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

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

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

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

1.1 研究區(qū)概況

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

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

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

中國氣象局陸面數(shù)據(jù)同化系統(tǒng)(CLDAS-V2.0)土壤濕度資料是基于數(shù)據(jù)融合和同化技術，利用氣溫、氣壓、濕度、風速、降水輻射數(shù)據(jù)和初始場信息，驅動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°，時間分辨率為逐小時.將每24 h的土壤濕度值求平均，獲得逐日土壤濕度，研究利用2019年4月1日至2019年10月31日0～10 cm土層土壤濕度數(shù)據(jù)(中國氣象數(shù)據(jù)服務中心:http:∥data.cma.cn)開展研究.

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

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

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

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

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

(1)

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

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

1.2.3 土壤質地數(shù)據(jù)

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

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

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

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

1.2.5 站點觀測數(shù)據(jù)

站點數(shù)據(jù)來源于國家氣象信息中心資料服務室的2019年逐小時觀測資料[23].該站點土壤濕度觀測在垂直方向上分為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.利用頻域反射技術(Frequency Domain Reflectometry,FDR)來測定土壤體積含水量[24].研究參考韓帥[25]提出的方法對土壤濕度觀測數(shù)據(jù)進行質量控制，得到有效站點223個，采用隨機方法，抽出185個作為訓練站點數(shù)據(jù)用于建模，38個站點數(shù)據(jù)用于驗證(圖1).

2 研究方法

2.1 降尺度因子的選擇

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

2.2 梯度提升機(GBM)

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

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

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

3)迭代n=1,2,…,N(N為樣本數(shù))，k=1,2,3,…,K(K為基學習器的個數(shù))，計算rkn:

(2)

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

4)對rkn建立基學習器模型Tk(x)，對梯度提升進行更新：

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

(3)

5)得到強學習器：

(4)

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

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

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

圖2 梯度提升機土壤濕度降尺度結構示意Fig.2 Downscaling of soil moisture based on gradient boosting machine

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

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

y=f(x,θ),

(5)

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

yh1=ReLU(W1X+b1),

(6)

yh2=ReLU(W2yh1+b2),

(7)

y=softmax(W3yh2+b3).

(8)

模型訓練及預測過程如下：

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

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

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

(9)

式中，xi為一列自變量中的第i個值，mean(x)和std(x)分別是x所在列自變量的均值和標準差.

2.4 隨機森林(RF)

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

SSMO=fRF(C)+ε,

(10)

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

(11)

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

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

(12)

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

2.5 Stacking集成學習方法

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

圖4 Stacking集成學習結構示意Fig.4 Structure of Stacking ensemble learning

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

(13)

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

(14)

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

(15)

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

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

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

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

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

圖5 Stacking集成學習方法土壤濕度降尺度結構示意Fig.5 Downscaling of soil moisture based on Stacking ensemble learning

2.6 評價指標

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

(16)

(17)

(18)

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

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

(19)

3 結果與分析

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

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

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

3.2 不同方法降尺度結果空間分布

為了比較不同方法降尺度效果，研究比較了梯度提升機(GBM)、深度前饋神經(jīng)網(wǎng)絡(DFNN)、隨機森林(RF)和Stacking集成學習等4種方法，因土壤溫度小于0 ℃時，觀測儀器異常導致缺少可靠的觀測數(shù)據(jù)，研究僅分析2019年4—10月的CLDAS土壤濕度產(chǎn)品降尺度效果.從降尺度結果和原土壤濕度日均值分布(圖6)可以看出，總體上，華北地區(qū)的南部和沿海區(qū)域土壤濕度較高，中部和北部土壤濕度較低，平均土壤濕度均達到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之間，并且降尺度后土壤濕度的空間分布細節(jié)更加豐富，這與降尺度過程中融合了對土壤濕度影響較大的高分辨率地表溫度、地表反照率、地形等地表變量數(shù)據(jù)有關.

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

利用站點觀測數(shù)據(jù)逐日土壤濕度評估表明(圖7和圖8)：不同降尺度方法預測的降尺度結果與站點觀測土壤濕度之間存在顯著的相關性，相關系數(shù)介于0.699 6～0.756 8，高于原土壤濕度的相關系數(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，說明降尺度后土壤濕度精度得到提高，相比于原土壤濕度，降尺度結果更接近于實測值.

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

3.4 不同方法降尺度結果時間序列及誤差分析

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

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

圖7 CLDAS逐日土壤濕度精度站點評估散點圖Fig.7 Scatter plot of CLDAS daily soil moisture and observations

從土壤濕度與觀測值的相關系數(shù)和誤差分析(圖10)來看，與原土壤濕度相比，4種降尺度方法降尺度結果相關系數(shù)均有所提升，其中，Stacking集成學習土壤濕度相關系數(shù)提高較大，平均提高0.13，其次為RF土壤濕度，平均提高0.12.與原土壤濕度相比，4種不同降尺度方法降尺度結果的均方根誤差均有明顯的降低，其中，Stacking集成學習土壤濕度和RF土壤濕度比原土壤濕度的均方根誤差，分別平均降低了0.012 1和0.011 7 m3/m3，更接近于觀測值.4種方法降尺度結果的偏差絕大多數(shù)日次在0以下，即降尺度結果存在一定程度的低估現(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集成學習方法為-0.005 2 m3/m3，4種降尺度結果均改善了原土壤濕度的高估問題，其中，以Stacking集成學習方法最優(yōu).

圖8 不同方法降尺度結果站點精度評估散點圖 a.GBM;b.DFNN;c.RF;d.Stacking集成學習Fig.8 Scatter plots of downscaled soil moistures and observations a.GBM;b.DFNN;c.RF;d.Stacking ensemble learning

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

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

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

4 結論與討論

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

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

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

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

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

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