基于平均影響值算法和BP神經網絡的根區土壤濕度估算

袁玲， 方秀琴,2*， 郭曉萌， 楊露露， 張曉祥， 任立良

(1.河海大學水文與水資源學院， 南京 210098； 2.河海大學海岸災害及防護教育部重點實驗室， 南京 210024)

土壤濕度(soil moisture, SM)是陸地表面系統的關鍵水文和生態變量，而根區土壤濕度(root zone soil moisture, RZSM)不僅影響降雨在入滲水和地表徑流中的分配，是洪水災害發生的重要孕災環境，還影響著地表植被的水分吸收和蒸騰作用，進而影響植被的結構和組織[1]。此外，由于土壤含水量和土抗剪強度之間的內在關系，RZSM對土壤侵蝕的加速和產沙量的增加也有一定影響。因此，對RZSM的預測不僅可以為洪水風險分析過程提供依據，也是農業規劃和干旱監測工作中的重要環節[2]。

在過去的幾十年里，為了獲取準確的土壤水分信息，發展了多種技術和儀器。現場土壤傳感儀，如中子探頭、電磁傳感器和脈沖傳感器被廣泛應用于獲取監測站點的土壤濕度數據[3]，但是這些監測方法需要耗費大量人力物力，而且局限于站點，得到的土壤濕度資料存在空間上密度低、時間上不連續，進而存在大量缺測值的問題[4]。此外，由于土壤質地、植被、地形等地表特征具有較高的空間異質性，因而利用點或場尺度的觀測數據來反映土壤水分的空間分布具有相當大的局限性[5]。遙感技術因其能夠獲取大尺度連續區域的土壤水分資料而得到了廣泛的關注，取得了重要進展，楊丹陽等[6]基于多時相高分一號影像建立反演模型對實驗區土壤濕度進行了反演，達到了較高的精度，然而利用遙感技術獲得土壤水分信息限制在地表5 cm以內，無法獲取到根區的土壤濕度值。總體而言，現有的觀測技術不足以滿足RZSM獲取的需求，在許多未測量區域和觀測站密度較低的區域，RZSM觀測在時間和空間上都存在大量的空白，而遙感技術由于測量深度有限，無法填補根區土壤濕度的空白。

近年來，研究人員提出了多種利用地表信息反演根區土壤濕度的技術。在眾多方法中，將土壤-植被-大氣傳輸模型和地表模型與數據同化技術相結合被證明是一種有效的RZSM估算方法，然而基于物理過程的模型的參數優化需要大量觀測數據來確定，并且模型實現過程復雜，其中難免存在重要物理過程概化的缺陷。針對此問題，一些研究提出了基于數據驅動的RZSM估算技術。與傳統水文模型不同，基于數據驅動的估算方法不需要通過明確的模型結構來表示所涉及的水文過程，僅需保證與目標值相關的訓練數據的真實性和可靠性，劉立文等[7]研究發現土壤濕度通過各種水文過程與氣象因素存在高度相關性，因此可以通過地表氣象要素數據驅動來反演根區土壤濕度。

在各種類型的數據驅動方法中，人工神經網絡被廣泛應用于土壤濕度的估算。Lee等[8]利用研究區610 km2區域內的10個監測站點數據來訓練人工神經網絡，并對該區的一個獨立站點的SM進行預測，獲得了良好的估算結果；楊娜等[9]以甘肅省為例，對氣象要素與土壤濕度之間的作用關系進行了研究，并建立BP神經網絡完成了土壤濕度值的估算實驗，獲得了較好的預測結果；Meng等[10]提出了一種基于GF-3衛星數據和人工神經網絡的農業地區土壤水分反演方法，使用從高級積分方程模型生成的訓練樣本數據集上對神經網絡進行訓練和測試，得到了較好的估算性能。雖然這些研究都顯示了人工神經網絡方法在土壤水分估計中的潛力，但都局限于小范圍區域，鮮有關于模型在不同氣候和土壤條件區域的適用性評價。在前人的研究基礎上，現選擇中國境內4個典型土壤濕度長期監測臺站為研究對象，以地表信息為輸入，根區土壤濕度值為輸出構建BP神經網絡模型，分別對不同站點的RZSM進行預測，并基于平均影響值(mean impact value，MIV)算法計算輸入地表要素的重要性，探討不同氣候和土壤區域的根區土壤濕度的主要驅動因素和模型性能表現差異，對BP神經網絡利用地表特征信息計算根區土壤濕度的性能和區域適用性進行研究。

1 實驗區與數據

1.1 實驗區概況

選擇分別位于中國境內東北黑土區、中南紅壤區、西南紫色土區和西北半干旱土區的四大整合農田土壤濕度監測場地為實驗區域，分別為遼寧沈陽(123.37°E, 41.52°N)、湖南桃源(111.47°E, 28.92°N)、四川鹽亭(105.46°E, 31.27°N)和新疆阿拉爾(80.85°E, 40.62°N)，站點土地利用類型和氣象條件情況如表1所示。

1.2 數據資料

1.2.1 土壤濕度數據

土壤濕度數據來自中國土壤數據庫(http://vdb3.soil.csdb.cn/)提供的站點土壤水分動態數據，該數據中土壤濕度用體積含水率來表示，即土壤中水分占有的體積和土壤總體積的比值。土壤濕度獲取方法為中子管法，數據時間步長為5 d，從每月1日開始，5 d一采樣，對每個監測站點在不同采樣地得到的土壤濕度取平均值作為該監測站點當日的土壤濕度值。5個站點提供10～150 cm每隔10 cm深度的土壤濕度監測值，綜合考慮數據的完整性，最終共選擇7層有效深度的土壤濕度數據為研究對象，對應深度和時間序列如下：鹽亭站第1層至第7層垂直深度分別為10、20、30、40、50、60、100 cm，時間范圍為2007年1月5日—2010年12月31日；沈陽站為10、20、30、40、50、70、90 cm，時間范圍為2004年4月30日—2011年10月4日；桃源站為10、20、30、40、50、70、90 cm，時間范圍為2004年1月11日—2011年12月20日；阿拉爾站為10、20、30、40、50、70、90 cm，時間范圍為2008年1月5日—2011年11月30日。其中，除沈陽站每年11月—次年3月存在數據缺測外，其余站點均為連續完整數據。將10 cm土壤濕度視為表層土壤濕度數據，其余層的土壤濕度視為對應不同深度的RZSM。

表1 實驗站點概況Table 1 Overview of experimental sites

對于SM時間序列中的少量缺失數據依據K鄰近插值算法(K-nearest neighbors，KNN)的基本思想進行插補，即在K個最近鄰對象中通過設定函數值來填充缺失值。研究設定K為10，使用10個最鄰近對象的加權平均值來填補缺失數據，得到4個野外臺站共1 310組完整時間序列的土壤濕度觀測值，其中鹽亭站287組，沈陽站306組，桃源站414組，阿拉爾站303組, 將每個站點的數據按照6∶2∶2的比例分為訓練集、驗證集和測試集，其中訓練集用于模型構建，驗證集用于調整模型的超參數，測試集用來評估模型最終的泛化能力。

1.2.2 氣象數據

氣象要素是土壤濕度計算中重要的地面特征變量，能在一定程度上反映土壤水分的變化情況，從中國氣象數據網(http://data.cma.cn)提供的觀測站點中提取與4個土壤濕度觀測站地理位置一致的站點，獲取對應站點的地面氣候資料日值數據集，并將其整理為模型所需的地表氣象輸入數據，包括日平均氣溫(daily mean temperature，TEM)、日累計降水量(daily cumulative precipitation，PRE)、日照時長(sunshine duration，SSD)、蒸發量(evaporation，EVP)、相對濕度(relative humidity，RHU)、風速(wind speed，WIN)6種氣象數據，并提取與土壤濕度數據的時間序列對應的各氣象因子的日值。

1.2.3 植被指數

近年來，歸一化植被指數(normalized difference vegetation index, NDVI)和歸一化差異紅外指數(normalized difference infrared index, NDII)被逐漸用于對土壤濕度的研究。Schnur等[11]利用美國一些土壤濕度站點實測數據和MODIS(moderate-resolution imaging spectroradiometer)的NDVI和EVI產品，研究發現植被指數與每個深度的土壤濕度都高度相關，當植被指數滯后5～10 d時，相關性達到最高，且將回歸模型應用到相似的遠處站點時，估算值能夠解釋觀測值的50%～88%；Sriwongsitanon等[12]發現在流域尺度上NDII與根區土壤濕度顯著相關，兩者的R2在旱季可以達到0.87；Joiner等[13]發現NDVI和NDII在幾天到幾周的時間內與根區土壤濕度有顯著的相關性；Castelli等[14]用Landsat提取的NDII序列與突尼斯干旱區的多個土壤濕度觀測站點數據進行分析，發現所有站點上兩者都有顯著相關(考慮整年數據的R2在0.62～0.67)，Castelli等[15]也將NDII成功用于非洲其他地區的根區土壤濕度估算研究；Pradhan[16]在美國西北部山區的研究證明NDVI能夠反映干旱/半干旱區的根區土壤水分的動態變化。Liu等[17]在中國華北小麥地的實驗結果表明太陽誘導葉綠素熒光(solar-induced chlorophyll fluorescence，SIF)和修正的NDVI都與根區土壤水分有顯著正相關。這些研究都證實植被的光譜指數(NDVI、NDII、SIF等)可能成為大尺度根區土壤濕度監測的有效途徑。

基于此，利用MODIS地表反射率產品MOD09GA，根據紅波段(620～670 nm)、近紅外波段(841～876 nm)和中紅外波段(1 628～1 652 nm)的平均反射率ρ，按式(1)和式(2)計算4個研究站點對應時間序列的NDVI和NDII值，即

NDVI=(ρ2-ρ1)/(ρ2+ρ1)

(1)

NDII=(ρ2-ρ6)/(ρ2+ρ6)

(2)

式中：ρ1、ρ2、ρ6分別為紅波段、近紅外波段和中紅外波段的平均反射率。

2 研究方法

2.1 輸入變量選擇

與根區土壤濕度相關的表層土壤濕度、氣象驅動數據和植被指數共9個模型作為輸入變量，為了使模型更加符合真實情況，考慮變量存在的時間滯后問題，將變量分為三組，其中，表層土壤濕度(SSM)、蒸發量(EVP)、相對濕度(RHU)、風速(WIN)、日照時常(SSD)為當日記錄值，反映日增量；日累計降水量(PRE)和日平均氣溫(TEM)為前a天的累計值，為模型提供對過去值的記憶；歸一化植被指數(NDVI)和植被水分指數(NDII)為滯后b天數據，這是由于植被生長狀況對土壤含水量的響應需要一定時間，且Pradhan[16]研究表明，NDVI和NDII與土壤水分在5～10 d時間滯后時表現出最為顯著的相關性。a和b的具體值通過比較每個站點變量與土壤濕度之間皮爾森相關系數值的大小及顯著性來確定，最終確定4個站點的輸入變量為：表層土壤濕度(SSM)、蒸發量(EVP)、相對濕度(RHU)、風速(WIN)、日照時常(SSD)、前45 d降水量累計值(PRE)、前45 d氣溫累計值(TEM)、滯后7 d NDVI值(NDVI)、滯后7 d NDII值(NDII)。

2.2 BP神經網絡

近20年來，人工神經網絡在水文學研究領域得到了普及與應用[18-19]，在多類神經網絡中，BP神經網絡因其較強的非線性映射能力而被廣泛使用。

2.2.1 模型結構和參數

采用三層BP神經網絡，包含輸入層、隱含層和輸出層。輸入層神經元個數為9，對應模型的9個輸入變量；輸出層神經元個數為1，對應站點當前垂直深度的根區土壤濕度值；隱含層神經元個數依據經驗公式[式(3)]并結合試錯法確定，公式為

(3)

式(3)中：M、P、N分別為輸入層、隱含層和輸出層神經元的數目；L∈[1，10]。

隱含層神經元個數的具體確定方法為：以每個站點20 cm深度土壤濕度值的估算為例，將隱含層神經元個數從4增加到14，計算測試集中預測值與目標值之間的均方根誤差(root mean squared error，RMSE)，結果如圖1所示。從圖1可以看出，RMSE并沒有大幅度上升或下降的趨勢，但總體來說，當隱含層神經元個數為6時，4組模型的綜合性能達到最優，故最終確定隱含層神經元個數為6。

圖1 均方根誤差RMSE與隱含層神經元個數的關系Fig.1 The relationship between the RMSE and the number of neurons in the hidden layer

以測試集預測精度為依據，設置網絡其他參數如下：學習率依據站點不同略有差異，鹽亭站為0.1、沈陽站為0.25、桃源站為0.45、阿拉爾站為0.15；最大迭代次數均為1 000次；激活函數為sigmoid函數；優化算法為梯度下降法。

2.2.2 模型性能評價指標

使用以下3個指標來反映估計值與目標值之間的擬合優度。

(1)決定系數R2。R2為相關系數的平方，反映輸入變量對結果的解釋程度，其值越大表示模型擬合優度越好。

(4)

(2)均方根誤差RMSE。RMSE為預測值與真值偏差的平方和與觀測次數比值的平方根，反映預測的精密度和模型的穩定性，其值越小表示模型仿真性能越好。

(5)

(3)平均絕對誤差MAE。MAE為所有預測值與真實值之間絕對誤差的平均值，反映預測值實際誤差的大小，其值越小表示模擬結果越接近真實值。

(6)

2.3 輸入變量重要性計算

為了計算9個輸入變量對4個站點不同深度土壤濕度預測結果的影響，更加直觀地反映地表各特征在根區土壤濕度預測中的重要程度，使用平均影響值(mean impact value，MIV)算法分別計算4個站點各輸入變量的影響值來反映其重要性。

平均影響值算法是Dombi等[20]提出的一種反映神經網絡中權重矩陣變化情況的方法，被廣泛應用于非線性算法中的變量排序，徐磊等[21]將BP神經網絡和MIV算法相結合，確定了公交場站污染物暴露水平的主要影響因素。

假設神經網絡共有l個輸入變量， MIV算法步驟如下。

步驟1依次將其中一個自變量Xi(i=1,2,…,l)的數值增加10%，其他l-1個變量不變，得到新的樣本后將其輸入訓練好的網絡進行仿真得到l組模擬結果Pi(i=1，2，…,l)。

步驟2依次將其中一個自變量Xi(i=1，2，…,l)的數值減小10%，其他l-1個變量不變，得到新的樣本后將其輸入訓練好的網絡進行仿真得到l組模擬結果Qi(i=1，2，…,l)。

步驟3增加自變量得到的模擬結果減去減小自變量的模擬結果求平均，得到變量Xi的MIV值。

3 結果與分析

3.1 模型性能

使用驗證集來評估模型性能，為了避免偶然誤差，提升模型的可信程度，對4個站點每個深度的模型各訓練500次，取500次R2、RMSE、MAE的平均值來衡量模型性能，結果如表2～表5所示。從站點來看，4個站點中桃源站的整體模型性能最好，而從土層深度來看，淺層土壤濕度的模擬精度高于要深層土壤。對4個站點的評估結果求平均(圖2)，結果如表6所示，結合表6和圖2可以看出，隨著土層深度的增加，模型的決定系數R2不斷減小，而均方根誤差(RMSE)和平均絕對誤差(MAE)則呈現明顯的增加趨勢。這與Pan等[5]的研究結果一致，即僅通過地表資料足以反演根區土壤水分狀況，但反演精度會隨著土壤深度的增加而減弱。

表2 鹽亭站不同土層深度的模型性能Table 2 Model performance in different depths at Yanting

表3 沈陽站不同土層深度的模型性能Table 3 Model performance in different depths at Shenyang

表4 桃源站不同土層深度的模型性能Table 4 Model performance in different depths at Taoyuan

表5 阿拉爾站不同土層深度的模型性能Table 5 Model performance in different depths at Alaer

表6 4個站點不同土層深度下模型評價指標均值Table 6 Averages of model evaluation indexes under different depths at four stations

圖2 4個站點不同土層深度下模型評價指標均值Fig.2 Averages of model evaluation indexes under different depths at four station

3.2 地表特征信息的重要性

使用平均影響值(MIV)算法分別計算4個站點各地表特征變量對RZSM計算模型的貢獻，結果如表7～表10所示，其中紅色背景代表變量對RZSM為正影響，且顏色越深影響值越大；綠色背景代表變量對RZSM為負影響，且顏色越深影響越大。

(1)表層土壤濕度(SSM)。在眾多地表特征變量中，無論對于哪個站點，SSM對于根區土壤水分計算模型來說都是最重要的輸入變量，這是因為，表層土壤水通過下滲進入深層土壤，是根區土壤水分的直接來源，也是RZSM計算的重要依據。

(2)降水(PRE)、歸一化植被指數(NDVI)和歸一化差異紅外指數(NDII)。這3個特征變量的影響值在4個站點多為正值，表明這3個地表特征在不同氣候和土壤區對不同深度的RZSM具有正面影響。這是因為，在自然狀態下，降水是土壤水分最重要的來源之一，降水量的多少直接影響土壤中水分含量，而地表植被對土壤濕度既存在正面影響也存在負面影響，當植被根系層對土壤水的保蓄作用大于冠層截留時，NDVI和NDII與土壤含水量之間表現為正相關。但從桃源和阿拉爾結果顯示，對于土地利用形式為水田和水澆地的兩個站點來說，PRE、NDVI和NDII指標的影響值都很小，這可能是由于外界因素的干擾，改變了土壤水分的主要來源途徑，使得降雨和植被的變化不能很好地反映土壤含水量的變化。

(3)蒸發量(EVP)。在沈陽和阿拉爾兩個位于溫帶地區的站點，EVP對RZSM表現為負影響，而在位于亞熱帶的兩個站點，EVP與淺層(20 cm)土壤水分存在較為明顯的負相關關系，而在深層(30 cm以下)對RZSM的影響不明顯。這可能是因為亞熱帶氣候雨熱同期的特點，當夏季地面蒸發量大的時候雨水也充沛，導致EVP這一地面特征并不能很好地反映深層根區土壤濕度的變化情況。

(4)日照(SSD)。在鹽亭、沈陽和桃源3個位于濕潤地區的站點，SSD對RZSM沒有表現出明顯的正影響或負影響，而在處于干旱區的阿拉爾站，SSD在淺層(50 cm及以上)對RZSM表現為明顯的負影響，這是因為在濕潤地區，土壤水分不僅具有多重來源，而且水分含量也受蒸發量、植被冠層截留等多重因素的影響，不確定性大，而在位于干旱區的阿拉爾站，土壤水分的主要來源不是自然降水，而是人工灌溉，日照時長會影響土壤水分的蒸發速度，對干旱區土壤濕度存在明顯的負影響。

(5)相對濕度(RHU)、風速(WIN)和氣溫(TEMP)。這3個特征變量對RZSM的影響在本研究中的4個實驗站點沒有表現出明顯的規律。

總的來說，根區土壤濕度值計算的主要依據是地表土壤濕度，多個地面特征之間的相互作用為輔助條件，且不同特征在不同氣候區和土壤類型區對RZSM的影響情況也不一致。

表7 鹽亭站不同土層深度變量的MIV值Table 7 MIV values of variables in different depths at Yanting

表8 沈陽站不同土層深度變量的MIV值Table 8 MIV values of variables in different depths at Shenyang

表9 桃源站不同土層深度變量的MIV值Table 9 MIV values of variables in different depths at Taoyuan

表10 阿拉爾站不同土層深度變量的MIV值Table 10 MIV values of variables in different depths at Alaer

3.3 模型應用

為了驗證模型的有效性，將訓練好的神經網絡應用到4個站點，可用時間序列的全部數據來預測不同深度的根區土壤濕度，并將其與觀測值進行比較，結果如圖3所示。

從整體上看，幾乎所有預測結果都能較好地反映根區土壤濕度值的變化趨勢，具體到預測精度上來說，桃源和阿拉爾兩個站點的預測結果與真實值最為接近，沈陽站模擬結果較好，而鹽亭站預測值與真實值之間的偏差較大。

同一站點不同深度的土壤濕度值略有差異但變化趨勢基本一致，說明以表層土壤水分數據為依據，結合氣象要素預測深層土壤濕度值這一方法是可行的。

4個站點的預測結果普遍存在高值偏高和低值偏低的問題，即模型對極值的預測還不夠準確。

模型對于淺層土壤濕度的估算性能普遍比深層土壤濕度要好，50 cm以上土層的擬合精度高于50 cm以下土層。

總的來說，構建BP神經網絡模型，依據表層土壤濕度和氣象要素來估算根區土壤濕度，在中國不同氣候區的4個站點都取得了較為理想的精度，證實了這一方法的可行性。

由于存在數據缺測，沈陽站缺失每年11月—次年3月數據值圖3 4個站點不同土層深度的模型預測值與真實值比較Fig.3 Comparison of predicted value with true value in different depths at the four sites

4 結論

選擇中國境內不同氣候區的4個土壤濕度監測站點為研究對象，構建了以地表信息為輸入，根區土壤濕度值為輸出的BP神經網絡模型，對不同站點的RZSM進行了預測，并基于MIV算法計算了輸入變量的重要性，對BP神經網絡的性能和區域適用性進行了探討。結論如下。

(1)在選擇合理參數的前提下，利用BP神經網絡模型來預測根區土壤濕度能到達較好的精度，在20～90 cm和20～100 cm深度，模型在驗證集上的平均R2值在4個站點分別為0.79、0.69、0.66、0.56、0.51和0.47；RMSE為1.91%、2.17%、2.51%、2.71%、2.82%和3.08%；MAE為1.44%、1.61%、1.75%、1.89%、2.04%和2.35%。說明構建的BP神經網絡模型能應用于中國不同氣候和土壤區的四大整合農田土壤濕度監測場根區土壤濕度的估算，且越接近表層，土壤水分估算模型的性能越好。

(2)從地表特征的重要性計算結果來看，SSM是模型最重要的輸入變量，此外，對位于不同氣候區的站點來說，氣象輸入要素對根區土壤濕度的影響都不相同。其中，PRE、NDVI和NDII 3個特征的影響值在4個站點多為正值；EVP在兩個溫帶地區的站點對RZSM表現為負影響，在位于亞熱帶的兩個站點，對淺層(20 cm)SM表現為負影響，而對深層(30 cm以下)SM的影響不明顯；SSD在3個位于濕潤地區的站點，沒有對RZSM表現出明顯的正影響或負影響，而在處于干旱區的站點，SSD對淺層(50 cm及以上)SM表現為明顯的負影響；RHU、WIN和TEMP對4個實驗站點RZSM的影響沒有表現出明顯的規律。

(3)從BP神經網絡的預測結果來看，該方法在所有站點幾乎都能較好地反映出根區土壤濕度值的變化趨勢，但在4個站點都存在極值預測不夠準確的問題，此外，模型對于淺層土壤濕度的估算性能普遍比深層土壤濕度要好，50 cm以上土層的擬合精度高于50 cm以下土層。