杜曉彤, 方秀琴, 汪 偉, 郭曉萌, 袁 玲

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

土壤濕度(Soil Moisture，SM)信息對(duì)于天氣和氣候預(yù)測(cè)、水文、干旱預(yù)報(bào)和流域管理至關(guān)重要[1]，它在調(diào)節(jié)陸地表面和大氣之間的水、能量和磁通量方面也意義非凡[2]。在地表，SM是與大氣相互作用的關(guān)鍵因素，在地下，尤其是在干旱和半干旱地區(qū)的根系層SM是植被根系水分吸收和蒸騰過(guò)程的主要控制因素之一[3]，特別是在作物生長(zhǎng)階段，根系層SM與農(nóng)業(yè)干旱密切相關(guān)。此外，生態(tài)系統(tǒng)水和能源的收支也都高度依賴于根系層SM的空間分布[4-5]。因此，對(duì)根系層SM的預(yù)測(cè)可以為農(nóng)業(yè)規(guī)劃、干旱診斷等提供依據(jù)[6]。

由于土壤性質(zhì)、植被等條件的異質(zhì)性，土壤濕度存在非常大的時(shí)空變異性。傳統(tǒng)站點(diǎn)觀測(cè)到的SM精度高，但是站點(diǎn)數(shù)量有限，且分布不均勻，無(wú)法進(jìn)行大面積的連續(xù)觀測(cè)，因此站點(diǎn)觀測(cè)無(wú)法獲取大面積的SM連續(xù)觀測(cè)信息。衛(wèi)星遙感是獲取SM的一種可行的方式，它們可以對(duì)地表SM狀況進(jìn)行連續(xù)、大規(guī)模的監(jiān)測(cè)，然而，它們只是提供了關(guān)于表土幾厘米的SM信息，不能準(zhǔn)確地表示根系層水分的變化[6]，而大多數(shù)生態(tài)系統(tǒng)的水分和能力收支高度依賴根系層SM的空間分布[7]。因此，近年來(lái)很多學(xué)者運(yùn)用不同方法對(duì)根系層SM估算進(jìn)行研究，主要有同化法[8-10]、指數(shù)濾波器法[11-13,15]、土壤水分分析關(guān)系(Soil Moisture Analytical Relationship，SMAR)法[1,6-7,14]。通過(guò)數(shù)據(jù)同化[8]將地表土壤水分觀測(cè)數(shù)據(jù)集成到水文模型中，已被證明是一種預(yù)測(cè)根區(qū)SM的有效方法，且應(yīng)用較廣泛。數(shù)據(jù)同化方法通過(guò)更新水文模型中的土壤水分狀態(tài)來(lái)估計(jì)剖面土壤水分，可以合理地反映土壤水分的時(shí)空分布，其中集成卡爾曼濾波(EnKF)是使用最多的同化方法[8-10]。但是在數(shù)據(jù)同化過(guò)程中，需要具備建模和觀測(cè)不確定性的先驗(yàn)知識(shí)，而且基于EnKF的方法不僅計(jì)算量大，而且對(duì)非線性系統(tǒng)有很強(qiáng)的局限性。Wagner等[15]在1999年就建議使用exp(-t/T)形式的指數(shù)濾波器從地表測(cè)量的時(shí)間序列中捕捉地下土壤濕度的動(dòng)態(tài)，T是特征時(shí)間長(zhǎng)度或衰退常數(shù)，它控制近地表水分含量在一段時(shí)間內(nèi)移動(dòng)和停留在根系層的多少。這種濾波器最大的優(yōu)點(diǎn)在于它的簡(jiǎn)單性，因?yàn)樗皇褂靡粋€(gè)參數(shù)，此外推導(dǎo)出的土壤水分指數(shù)(Soil Water Index，SWI)僅依賴于地面觀測(cè)數(shù)據(jù)。很多研究采用指數(shù)濾波器的方法進(jìn)行了根系層SM的估算，均獲得了良好的效果(Albergel[11]，F(xiàn)ord[12]，Gao[13])。但是，指數(shù)濾波器法中未發(fā)現(xiàn)T參數(shù)與土壤性質(zhì)之間的關(guān)系，限制了指數(shù)濾波器在衛(wèi)星遙感數(shù)據(jù)方面的應(yīng)用[1,4]。2014年Manfreda等[1]將Wagner等的模型擴(kuò)展為兩層基于物理的入滲模型，即為SMAR模型，提出了一種新的估算根區(qū)SM的方法。SMAR是由半干旱環(huán)境中一個(gè)簡(jiǎn)化的土壤水分平衡方程推導(dǎo)出來(lái)的，該方程提供了根區(qū)與地表土壤水分之間的關(guān)系，并且具有有限的物理一致性參數(shù)，也揭示了上述關(guān)系及其在利用衛(wèi)星遙感反演土壤水分方面的可能應(yīng)用。相關(guān)學(xué)者[1,6-7，14]應(yīng)用該模型估算根系層的土壤濕度，大都局限在站點(diǎn)尺度的監(jiān)測(cè)，可以很好地反映出站點(diǎn)根系層SM的動(dòng)態(tài)，但是對(duì)于區(qū)域面上分布的應(yīng)用研究較少。

由于半干旱區(qū)的生態(tài)環(huán)境較為脆弱，具有對(duì)人為擾動(dòng)和氣候變化響應(yīng)敏感性的特征，在全球變化和水文研究中具有重要意義，因此本文選取位于半干旱區(qū)的老哈河流域，將遙感表層SM數(shù)據(jù)與適用于半干旱環(huán)境的SMAR方法相結(jié)合來(lái)估算根系層SM信息，獲取老哈河流域根系層SM的時(shí)空分布。本文探究SMAR模型在半干旱流域的空間應(yīng)用性能，為在更大范圍內(nèi)估算根系層SM提供支撐；另外，本研究的成果可以為區(qū)域干旱監(jiān)測(cè)和改善水文模擬提供有效的數(shù)據(jù)源。

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

1.1 研究區(qū)概況

考慮到半干旱區(qū)在全球變化和水文研究的重要意義以及遙感數(shù)據(jù)在該區(qū)域比較完備(全年有值天數(shù)達(dá)40%左右，遠(yuǎn)高于濕潤(rùn)地區(qū))等因素，故選取半干旱區(qū)的老哈河流域作為研究區(qū)，其地理范圍位于北緯41°—43°、東經(jīng)117°—120°，基于遙感土壤濕度數(shù)據(jù)估算根系層土壤濕度時(shí)空分布。這里以興隆坡水文站控制的集水區(qū)作為老哈河流域研究區(qū)，位于西遼河上游，流域地跨河北、遼寧兩省和內(nèi)蒙古自治區(qū)，海拔427～2 054 m，面積約為2.1萬(wàn)km2。

1.2 土壤濕度數(shù)據(jù)

研究區(qū)采用的土壤濕度數(shù)據(jù)是1993—2006年的4—9月共84個(gè)月的農(nóng)氣站土壤濕度觀測(cè)資料和CCI土壤濕度產(chǎn)品。農(nóng)氣站實(shí)測(cè)的SM數(shù)據(jù)來(lái)源于中國(guó)農(nóng)作物氣象數(shù)據(jù)網(wǎng)的中國(guó)農(nóng)作物生長(zhǎng)發(fā)育和農(nóng)田土壤濕度旬值數(shù)據(jù)集，包括10 cm，20 cm，50 cm，70 cm和100 cm多個(gè)根區(qū)剖面不同深度土壤相對(duì)濕度的逐旬觀測(cè)值，原始數(shù)據(jù)是相對(duì)土壤濕度(%)。本研究采用10 cm深度處的相對(duì)飽和值作為表層SM，而在10 cm以下根部剖面上相對(duì)飽和度的均值作為根系層SM?？紤]到數(shù)據(jù)的有效性，選取1993—2006年的逐年4月、5月、6月、7月、8月、9月實(shí)測(cè)的逐旬?dāng)?shù)據(jù)，分別求取表層和根系層不同作物逐月各旬的均值作為月平均，再剔除個(gè)別缺測(cè)月份最終作為各農(nóng)氣站多年逐月實(shí)測(cè)的相對(duì)土壤濕度數(shù)據(jù)。由于SMAR模型采用的簡(jiǎn)化土壤水分平衡方程適用于干旱半干旱區(qū)，因此本文根據(jù)中國(guó)300～500 mm年等降雨量線選擇了82個(gè)農(nóng)氣站點(diǎn)的數(shù)據(jù)(圖1)。

圖1 年降雨量300～500 mm農(nóng)氣站分布圖

歐空局CCI土壤濕度數(shù)據(jù)集來(lái)源于：https:∥www.esa-soilmoisture-cci.org，本研究使用的是ESA CCI_SM v04.4產(chǎn)品，涵蓋了從1978—2018年6月30日的(近)40 a記錄，包括主動(dòng)微波遙感數(shù)據(jù)集、被動(dòng)微波遙感數(shù)據(jù)集和主被動(dòng)合成的微波遙感數(shù)據(jù)集，這里使用的是主被動(dòng)合成的微波遙感數(shù)據(jù)集。主動(dòng)微波遙感數(shù)據(jù)集是由維也納大學(xué)根據(jù)ERS-1，ERS-2和METOP-A和METOP-B的C波段散射計(jì)觀測(cè)結(jié)果獲得，被動(dòng)微波遙感數(shù)據(jù)集是根據(jù)Nimbus 7 SMMR，DMSP SSM/I，TRMM TMI，Aqua AMSR-E，GCOM-W1 AMSR2和SMOS的微波觀測(cè)結(jié)果生成[16]。提供的數(shù)據(jù)存儲(chǔ)格式為NetCDF-4，空間分辨率為0.25°，時(shí)間分辨率為日，原始單位是體積含水量(m3/m3)。CCI SM選取與農(nóng)氣站觀測(cè)數(shù)據(jù)相同的時(shí)間長(zhǎng)度，將逐日的數(shù)據(jù)在ArcGIS中對(duì)應(yīng)處理成月均值，最終提取出逐年逐月逐像元的數(shù)據(jù)。

1.3 空間數(shù)據(jù)集

本文采用的NASA GLDAS V2.0數(shù)據(jù)集是GLDAS中NOAH陸面模式計(jì)算的蒸散發(fā)產(chǎn)品，數(shù)據(jù)來(lái)源：http:∥disc.sci.gsfc.nasa.gov/hydrology/data-holdings，數(shù)據(jù)格式是NetCDF，空間分辨率是0.25°，時(shí)間分辨率是月。

歸一化植被指數(shù)NDVI數(shù)據(jù)采用的是GIMMS NDVI3 g V1數(shù)據(jù)集，來(lái)源于：https:∥ecocast.arc.nasa.gov/data/pub/gimms/3 g.v1/，原始數(shù)據(jù)格式為NetCDF，每年有兩個(gè)nc4文件，每個(gè)nc4文件包含6個(gè)月NDVI數(shù)據(jù)，時(shí)間分辨率為15 d，空間分辨率為0.083°，共12景。本研究通過(guò)格式轉(zhuǎn)換、最大值合成等預(yù)處理過(guò)程，使得每一個(gè)像元為該月最大NDVI值，從而獲得研究區(qū)逐年逐月的NDVI數(shù)據(jù)集。

1.4 土壤屬性數(shù)據(jù)

本文搜集農(nóng)氣站以及研究區(qū)逐像元的土壤屬性數(shù)據(jù)來(lái)源于《中國(guó)土種志》[17-19]。中國(guó)土種志以大區(qū)為單位，分別介紹了各大區(qū)主要土種在土壤分類系統(tǒng)中的歸屬、分布、面積、主要性狀、典型剖面、生產(chǎn)性能和理化分析數(shù)據(jù)。本研究主要整理了82個(gè)農(nóng)氣站對(duì)應(yīng)的典型剖面的土壤屬性組成，包括粗砂(2～0.2 mm)、細(xì)砂(0.2～0.02 mm)、粉砂(0.02～0.002 mm)、黏粒(<0.002 mm)以及各深度對(duì)應(yīng)的質(zhì)地名稱，將表層深度以下的各屬性均值作為站點(diǎn)根系層土壤的各屬性數(shù)據(jù)。

2 研究方法

本文目標(biāo)是獲得研究區(qū)長(zhǎng)時(shí)間序列根系層土壤濕度的空間分布。首先采用CDF匹配方法以農(nóng)氣站實(shí)測(cè)數(shù)據(jù)為基準(zhǔn)對(duì)CCI土壤濕度產(chǎn)品進(jìn)行校正，以支持基于物理基礎(chǔ)的SMAR模型。其次以中國(guó)300～500 mm等降雨量線內(nèi)的部分農(nóng)氣站實(shí)測(cè)的土壤表層SM作為SMAR模型的輸入，將實(shí)測(cè)根系層相對(duì)SM與模型模擬的根系層相對(duì)SM的RMSE最小作為目標(biāo)函數(shù)，利用MATLAB的遺傳算法工具箱進(jìn)行模型參數(shù)優(yōu)化，得到每個(gè)站點(diǎn)最優(yōu)的參數(shù)組合。然后采用多元統(tǒng)計(jì)回歸分析法結(jié)合土壤物理性質(zhì)、NDVI以及年均ET建立估算SMAR參數(shù)的多元線性回歸模型。最后以校正后的CCI土壤濕度作為模型輸入數(shù)據(jù)，運(yùn)用SMAR模型估算研究區(qū)逐像元根系層的土壤相對(duì)濕度。

2.1 CDF匹配

由于遙感獲取的CCI土壤濕度產(chǎn)品和農(nóng)氣站實(shí)測(cè)的地表土壤濕度數(shù)據(jù)存在系統(tǒng)誤差，導(dǎo)致這兩組數(shù)據(jù)存在不一致性，因此為了使遙感數(shù)據(jù)與農(nóng)氣站實(shí)測(cè)數(shù)據(jù)更好地匹配，采用累積概率分布函數(shù)(CDF)[16]重調(diào)整方法，基于實(shí)測(cè)土壤濕度數(shù)據(jù)對(duì)遙感數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化，使其具有相同的取值范圍和累積概率分布特征。CDF重調(diào)整過(guò)程是以像元點(diǎn)為基礎(chǔ)進(jìn)行時(shí)間序列統(tǒng)計(jì)并調(diào)整。在實(shí)踐中，農(nóng)氣站的分布數(shù)量有限，且不均勻，因此需要通過(guò)對(duì)站點(diǎn)數(shù)據(jù)進(jìn)行空間插值獲取所研究區(qū)域SM的整體時(shí)空分布狀況。本文選擇老哈河流域周圍的赤峰等8個(gè)站點(diǎn)的1993—2006年4—9月的有效數(shù)據(jù)進(jìn)行克里金空間插值[16]，以此獲得老哈河流域表層實(shí)測(cè)SM的時(shí)空分布，再在此基礎(chǔ)上利用CDF匹配技術(shù)將CCI土壤濕度產(chǎn)品按照與實(shí)測(cè)數(shù)據(jù)相同的時(shí)間序列進(jìn)行重調(diào)整，建立研究區(qū)逐像元的CCI土壤濕度產(chǎn)品的改進(jìn)模型。

CDF匹配技術(shù)[20]是基于分段線性匹配，研究中，調(diào)整后的土壤濕度表達(dá)式如下:

CDFsite(x′)=CDFCCI(x)

(1)

式中：CDFsite和CDFCCI分別表示站點(diǎn)和CCI土壤濕度的累積分布函數(shù)；x表示調(diào)整之前的CCI土壤濕度。

2.2 SMAR模型

Manfreda等[1]假定土壤由兩層組成，第一層是幾厘米深的表層(合理的假設(shè)是衛(wèi)星傳感器觀測(cè)深度可以代表5—10 cm的表層動(dòng)態(tài))，而第二層位于第一層之下，延伸到相當(dāng)于植被的生根深度。兩層土壤最重要的水分交換是入滲，橫向流動(dòng)和毛細(xì)管上升等過(guò)程則可以忽略不計(jì)[14]。Manfreda等[1]定義了一個(gè)土壤水分平衡方程，它是由半干旱環(huán)境中一個(gè)簡(jiǎn)化的土壤水分平衡方程推導(dǎo)出來(lái)的，在這個(gè)方程中入滲期不是表示為降雨的函數(shù)，而是表層中土壤含水量的函數(shù)，分別用“1”和“2”區(qū)分上下層的變量和參數(shù)。只有當(dāng)土壤水分超過(guò)田間持水量時(shí)，表層的水分通量才能被認(rèn)為是顯著的。他們提出了描述從頂層到下層的瞬時(shí)入滲通量的方程：

(2)

第二層和深層土壤水分平衡受入滲和水土流失兩個(gè)因素控制，水土流失主要又包括蒸發(fā)蒸騰損失和滲透。定義土壤第二層的“有效”土壤相對(duì)飽和度為：

x2=(s2-sw2)/(1-sw2)

(3)

土壤貯水量為：

ω0=(1-sw2)n2Zr2

(4)

則土壤水分平衡可以用式(5) 來(lái)表述：

(5)

式中:s2[-]表示第二層的相對(duì)飽和度；sw2[-]表示凋萎點(diǎn)的相對(duì)飽和度；n2[-]表示第二層的土壤孔隙度；Zr2[L]表示第二層的土層深度；V2[LT-1]表示蒸散和滲透損失的土壤水分損失系數(shù)；x2[-]表示第二層土壤的“有效”相對(duì)飽和度。此方程不考慮土壤水分高值時(shí)土壤流失函數(shù)的高度非線性，這種簡(jiǎn)化，再加上滲透不考慮飽和效應(yīng)，意味該方法在潮濕環(huán)境中有一定的局限性,更適用于干旱半干旱環(huán)境。

上述方程可以用歸一化系數(shù)和簡(jiǎn)化為：

(6)

“擴(kuò)散系數(shù)”是土壤表層總可能儲(chǔ)水量(n1Zr1)與總根區(qū)儲(chǔ)水量((1-sw2)n2Zr2)之比確定，“水分流失系數(shù)”控制根系層水分損失，是蒸散發(fā)ET與排水損失之和與總根區(qū)儲(chǔ)水量之比確定，這些參數(shù)的取值可直接與兩土層的深度和水土流失系數(shù)之比有關(guān)。因此，土壤水分平衡方程即為式(7)：

(7)

假設(shè)相對(duì)飽和度x2(t)初始條件等于0，可以推導(dǎo)出這種線性微分方程的解析解：

(8)

在實(shí)際具體應(yīng)用中，可能還需要離散的形式見(jiàn)公式(9)：

(9)

擴(kuò)展式中，假設(shè)Δt=(tj-t(j-1))，根據(jù)表層土壤濕度的時(shí)間序列可以推導(dǎo)出下層土壤濕度的公式(10)：

x2(tj)=x2(tj-1)e-a(tj-tj-1)+by(tj)(tj-tj-1)

(10)

也可以寫(xiě)成s2的函數(shù)：

s2(tj)=[sw2+(s2(tj-1)-sw2)e-a(tj-tj-1)]+[(1-sw2)by(tj)(tj-tj-1)]

(11)

式中：左半部分計(jì)算給定的時(shí)間步長(zhǎng)中相對(duì)水分損失，水分損失函數(shù)隨s2線性減小，直至達(dá)到sw2，此時(shí)水分損失為零，蒸散和排水損失假定在凋萎點(diǎn)值或以下可忽略不計(jì)，右半部分表示從表面到根區(qū)擴(kuò)散的水分含量。Manfreda等[1]提出的方法是一種代表了兩個(gè)狀態(tài)變量之間土壤水分分析關(guān)系式，需要確定的4個(gè)參數(shù)分別是：sw2，a，b，sc1，所有參數(shù)均可從土壤質(zhì)地分?jǐn)?shù)、土層深度和土壤水分損失等估計(jì)。

2.3 遺傳算法

為了分別得到最優(yōu)的模型參數(shù)組合，本文使用MATLAB的遺傳算法工具箱GUI進(jìn)行模型參數(shù)的優(yōu)化估計(jì)[14,21]。MATLAB遺傳算法工具，即圖形用戶界面GUI，它使我們可以使用遺傳算法而不用命令行的方式。遺傳算法[22]是一種基于自然選擇、生物進(jìn)化過(guò)程來(lái)求解問(wèn)題的方法。在每一步中，遺傳算法隨機(jī)地從當(dāng)前種群中選擇若干個(gè)體作為父輩，并且使用它們產(chǎn)生下一代的子種群。在連續(xù)若干代后，種群朝著優(yōu)化解的方向進(jìn)化，最終得到最優(yōu)解組合。本文以各站點(diǎn)表層SM數(shù)據(jù)模擬得到的根系層相對(duì)SM與農(nóng)氣站實(shí)際測(cè)量得到的根系層相對(duì)SM均值之間的RMSE最小為優(yōu)化準(zhǔn)則進(jìn)行優(yōu)化處理，假設(shè)每個(gè)參數(shù)的上、下邊界分別相當(dāng)于最大和最小可能值(例如，0≤sc1≤1)。參數(shù)估計(jì)時(shí)，第二層土壤相對(duì)飽和度s2初始值設(shè)置為各農(nóng)氣站實(shí)測(cè)的逐旬田間持水量相對(duì)飽和值的均值。應(yīng)該注意的是，SMAR模型可能產(chǎn)生高于1的值，這些值需自動(dòng)設(shè)置為1。

3 結(jié)果與分析

3.1 CCI土壤濕度產(chǎn)品的CDF重調(diào)整

CDF重調(diào)整過(guò)程是以像元點(diǎn)為基礎(chǔ)進(jìn)行時(shí)間序列統(tǒng)計(jì)并調(diào)整，本文以研究區(qū)赤峰站所在像元進(jìn)行示例分析，CDF重調(diào)整過(guò)程為：

(1) 分別繪制農(nóng)氣站實(shí)測(cè)表層和遙感SM數(shù)據(jù)的CDF曲線，以0%，5%，10%，20%，30%，40%，50%，60%，70%，80%，90%，95%，100%這13個(gè)節(jié)點(diǎn)得到曲線上對(duì)應(yīng)點(diǎn)的坐標(biāo)，將曲線劃分成12段，見(jiàn)圖2A、B)所示。

(2) 根據(jù)上一步得到各節(jié)點(diǎn)的CDF值建立CCI土壤濕度與農(nóng)氣站表層實(shí)測(cè)SM的分段線性回歸曲線，見(jiàn)圖2C所示。對(duì)每個(gè)部分進(jìn)行分段線性回歸，從而得到12組線性回歸方程。

(3) 利用得到的分段線性回歸方程，進(jìn)行CDF匹配調(diào)整，從而得到農(nóng)氣站實(shí)測(cè)土壤濕度、CCI土壤濕度和調(diào)整后的土壤濕度(圖2D)，圖中實(shí)線表示農(nóng)氣站實(shí)測(cè)土壤濕度，粗虛線表示原始CCI土壤濕度，細(xì)虛線表示調(diào)整后的土壤濕度數(shù)據(jù)。

圖2 實(shí)測(cè)土壤濕度的CDF曲線(A)、CCI土壤濕度的CDF曲線(B)、CCI與實(shí)測(cè)土壤濕度分段線性回歸曲線(C)、CDF重調(diào)整前后比較(D)

以內(nèi)蒙古赤峰站所在像元為例，已知1993—2006年逐年4—9月的站點(diǎn)實(shí)測(cè)土壤濕度值及CCI遙感土壤濕度，將實(shí)測(cè)數(shù)據(jù)轉(zhuǎn)化成以體積含水量為單位，進(jìn)行CDF重調(diào)整，從而獲得調(diào)整后的土壤濕度數(shù)據(jù)時(shí)間序列曲線(圖3)，從圖3可以看出，原始CCI SM整體上小于實(shí)測(cè)的SM，但是兩套數(shù)據(jù)的趨勢(shì)具有良好的一致性，調(diào)整后的數(shù)據(jù)普遍向上抬升，在數(shù)值和形態(tài)表現(xiàn)上都更加靠近基準(zhǔn)的SM數(shù)據(jù)。對(duì)該站點(diǎn)所在像元CDF重調(diào)整前后的CCI SM與實(shí)測(cè)數(shù)據(jù)之間的誤差進(jìn)行統(tǒng)計(jì)分析(表1)，可以看出CCI SM經(jīng)重調(diào)整后各項(xiàng)誤差指數(shù)都得到降低，尤其是MAE和RMSE減少12%以上，表明調(diào)整后的數(shù)據(jù)更加接近于實(shí)測(cè)的SM數(shù)據(jù)。綜上，赤峰站所在像元的實(shí)測(cè)土壤數(shù)據(jù)經(jīng)過(guò)CDF重調(diào)整后精度得到明顯改善，整個(gè)研究區(qū)像元調(diào)整方法類似，誤差也得到減小，精度得到提高。

表1 赤峰站所在像元CDF調(diào)整前后誤差分析

3.2 農(nóng)氣站點(diǎn)SMAR模型的應(yīng)用及參數(shù)優(yōu)化

本文首先根據(jù)中國(guó)農(nóng)作物氣象數(shù)據(jù)網(wǎng)的中國(guó)農(nóng)作物生長(zhǎng)發(fā)育和農(nóng)田土壤濕度旬值數(shù)據(jù)集整理了分布在中國(guó)年降雨量線300～500 mm內(nèi)的82個(gè)站點(diǎn)不同剖面深度實(shí)測(cè)相對(duì)SM,將農(nóng)氣站表層10 cm的相對(duì)SM的時(shí)間序列作為SMAR模型的輸入數(shù)據(jù)，采用MATLAB的遺傳算法優(yōu)化由表層相對(duì)SM模擬作為輸入得到的根系層相對(duì)SM與其實(shí)測(cè)的10 cm以下的根系層相對(duì)SM均值之間的RMSE，從而得到每個(gè)農(nóng)氣站SMAR模型最優(yōu)的參數(shù)組合，為后續(xù)SMAR參數(shù)的統(tǒng)計(jì)分析奠定基礎(chǔ)。該遺傳算法能夠在大多數(shù)站點(diǎn)位置找到每個(gè)參數(shù)的近似全局最小值。這里以老哈河流域周圍的包括赤峰等6個(gè)站點(diǎn)為例運(yùn)用SMAR模型估算站點(diǎn)尺度根系層的相對(duì)土壤濕度。

圖3 赤峰站所在像元CDF匹配調(diào)整前后的土壤濕度時(shí)間序列曲線

圖4和表2給出了SMAR模型在研究區(qū)及周邊幾個(gè)農(nóng)氣站14 a的初步應(yīng)用，可以看到該方法能很好地捕捉下層土壤水分的動(dòng)態(tài)變化，與實(shí)測(cè)根系層的相對(duì)SM趨勢(shì)變化具有良好的一致性，并接近于每個(gè)站點(diǎn)的觀察結(jié)果。RMSE和R分別表示農(nóng)氣站實(shí)測(cè)10 cm以下根系層的相對(duì)SM均值與SMAR模型模擬的根系層的相關(guān)關(guān)系。根據(jù)MATLAB遺傳算法[14]對(duì)站點(diǎn)的參數(shù)進(jìn)行校準(zhǔn)，在均方根誤差方面表現(xiàn)較好，而在相關(guān)性方面各站點(diǎn)表現(xiàn)好壞存在差異，R最低的奈曼旗站為0.471，最高的阜新站為0.895。

表2中的參數(shù)值和誤差值都與站點(diǎn)實(shí)測(cè)值存在密切關(guān)系，產(chǎn)生誤差的一部分原因是SMAR模型估算根系層相對(duì)土壤濕度是具有時(shí)間序列的，但在實(shí)際運(yùn)行時(shí)由于部分月份可能因?yàn)閿?shù)據(jù)缺測(cè)而被剔除，造成前后時(shí)間間斷，時(shí)間序列的中斷也會(huì)對(duì)后續(xù)模型估算產(chǎn)生影響。從模型運(yùn)行來(lái)看，根區(qū)田間持水量的初始值對(duì)于模型的預(yù)測(cè)影響是可以忽略不計(jì)的。另外各站點(diǎn)氣候條件不同或區(qū)域的特定空間異質(zhì)性也導(dǎo)致了它們SMAR模型參數(shù)值的不同。

表2 6個(gè)農(nóng)氣站點(diǎn)SMAR模型的參數(shù)組合

綜上，由老哈河流域及周邊實(shí)測(cè)站點(diǎn)的根系層土壤濕度估算可以看出，SMAR模型估算半干旱的站點(diǎn)根系層SM具有良好的效果，為區(qū)域根系層SM的估算奠定了基礎(chǔ)。

3.3 SMAR參數(shù)回歸

考慮到SMAR模型參數(shù)的物理意義與土壤和植被的關(guān)系，因此可以將土壤性質(zhì)、NDVI和ET[3,6]用于預(yù)測(cè)SMAR參數(shù)(表3)，所構(gòu)建的回歸關(guān)系與歐洲相關(guān)學(xué)者[15]運(yùn)用土壤性質(zhì)作為變量建立的pedo—轉(zhuǎn)換函數(shù)具有類似的預(yù)測(cè)能力。

本文先根據(jù)各農(nóng)氣站的質(zhì)地編號(hào)結(jié)合《中國(guó)土種志》[17-19]整理了各站點(diǎn)土壤不同剖面深度下的機(jī)械組成(包括粗砂、細(xì)砂、粉砂、黏土含量)，剔除未知土壤類型的站點(diǎn)以及參數(shù)異常的站點(diǎn)后剩余38個(gè)站點(diǎn)用于多元回歸分析。用這些農(nóng)氣站點(diǎn)土壤物理性質(zhì)、NDVI和年平均ET數(shù)據(jù)(這里提取的是逐年4—9月的均值為年平均)作為變量，建立了估計(jì)SMAR參數(shù)的回歸方程，并通過(guò)了SMAR參數(shù)與土壤性質(zhì)之間(p<0.05)的雙尾t檢驗(yàn)。

站點(diǎn)根系層凋萎點(diǎn)水平與土壤性質(zhì)的相關(guān)性最強(qiáng)，與根區(qū)細(xì)砂含量呈明顯的負(fù)相關(guān)關(guān)系(雙尾t檢驗(yàn)，p<0.001)(圖5A)，從模型(表3)中也可以看出細(xì)砂含量對(duì)sw2有明顯的作用。年均ET是SMAR各參數(shù)的一個(gè)重要且具有統(tǒng)計(jì)學(xué)意義的變量，因?yàn)橹参锼治蘸捅砻嬲舭l(fā)會(huì)影響近地表和根系層水分狀況，與各參數(shù)都密切相關(guān)。其中擴(kuò)散系數(shù)b與年平均ET呈顯著的正線性相關(guān)(雙尾t檢驗(yàn)，p=0.008)，這表明在b的回歸中，應(yīng)加入ET，通過(guò)近地表擴(kuò)散率與表層田間持水量的聯(lián)系，也應(yīng)將其考慮到田間持水量sc1的回歸中。水分流失系數(shù)a反映的是蒸發(fā)蒸騰和排水損失之和，從圖5C中可以看出與年平均ET呈正相關(guān)關(guān)系(雙尾t檢驗(yàn)，p=0.038)，圖5D可以看出與上層粉砂的含量也呈正相關(guān)(雙尾t檢驗(yàn)，p=0.04)，建立的模型見(jiàn)表3。黏土含量和孔隙度與其他變量相比影響不大，在模型中也沒(méi)有表現(xiàn)出來(lái)。NDVI能夠檢測(cè)植被的生長(zhǎng)狀態(tài)，可以反映植物冠層的背景影響，如土壤等，因此也可以用來(lái)估計(jì)SMAR模型的相關(guān)參數(shù)。從表3中可以看出擴(kuò)散系數(shù)與NDVI呈負(fù)相關(guān)關(guān)系，這里的NDVI采取的是各站點(diǎn)所有時(shí)間序列的平均值。綜上，表3所示的多元線性回歸函數(shù)均能良好地估算SMAR模型的4個(gè)參數(shù)，并都通過(guò)了雙尾t，p<0.05的顯著性檢驗(yàn)，這為區(qū)域逐像元模型的應(yīng)用提供了支撐。

圖4 6個(gè)站點(diǎn)實(shí)測(cè)的表層、根系層土壤濕度和SMAR模擬得到的土壤濕度時(shí)間序列曲線

表3 SMAR模型各參數(shù)的多元線性回歸函數(shù)

參數(shù)多元線性回歸函數(shù)R2p值sw2-1.994**-細(xì)砂2(0.347)*+ln(AET)(0.462)***0.497<0.001a-0.303+粉砂1(0.085)*+ln(AET)(0.057)*0.2340.004b-2.137**+粉砂1(0.317)*-NDVI(0.496)*+ln(AET)(0.397)**0.2820.002sc1-1.596*+粉砂2(0.448)*+ln(AET)(0.34)*0.2890.001

注：AET：年均蒸散發(fā)；1：表示土壤表層；2：表示土壤根系層

3.4 SMAR模型在流域空間內(nèi)的應(yīng)用

由于農(nóng)氣站點(diǎn)數(shù)量與規(guī)模的限制以及遙感SM只能反映表層的土壤水分狀況，因此本文將遙感SM產(chǎn)品與SMAR模型結(jié)合起來(lái)，進(jìn)行流域空間內(nèi)根系層土壤水分時(shí)空分布的估算。根據(jù)CCI遙感土壤濕度產(chǎn)品的空間分辨率0.25°×0.25°，提取覆蓋老哈河流域的所有像元，共77個(gè)。根據(jù)《中國(guó)土種志》逐像元整理出對(duì)應(yīng)的質(zhì)地類型和土壤機(jī)械組成數(shù)據(jù)，并對(duì)應(yīng)提取出各像元NDVI和年均ET數(shù)據(jù)，然后按照上文在農(nóng)氣站所使用的SMAR模型參數(shù)回歸的方法進(jìn)行參數(shù)回歸。將重調(diào)整后的逐像元遙感SM作為SMAR模型的最終輸入數(shù)據(jù)，再結(jié)合回歸分析得到的模型參數(shù)，最終可以估算出覆蓋流域范圍的根系層相對(duì)土壤濕度的分布狀況。以老哈河流域2006年4—9月為例，模擬出區(qū)域根系層相對(duì)土壤水分狀況(圖6)，從2006年的模擬圖中可以看出當(dāng)遙感土壤濕度作為SMAR輸入數(shù)據(jù)時(shí)，能夠較好地模擬出根系層的相對(duì)土壤水分狀態(tài)。

圖5 SMAR模型參數(shù)與變量的相關(guān)關(guān)系散點(diǎn)圖

圖6 老哈河流域遙感數(shù)據(jù)作為輸入模擬得到的2006年4—9月根系層土壤水分狀況

為了進(jìn)一步分析遙感SM數(shù)據(jù)作為SMAR模型輸入的結(jié)果，分別以表層實(shí)測(cè)和重調(diào)整后的遙感數(shù)據(jù)作為輸入，對(duì)覆蓋老哈河流域的矩形柵格區(qū)域進(jìn)行逐像元SMAR模擬，對(duì)兩種結(jié)果進(jìn)行相關(guān)關(guān)系和誤差分析(圖7)?？梢钥闯觯瑑烧叩南嚓P(guān)系數(shù)總體較好，R平均值為0.692，經(jīng)過(guò)SPSS雙變量相關(guān)分析得到，97%的像元通過(guò)p<0.05的雙尾t顯著性檢驗(yàn)。從4種誤差分析來(lái)看，誤差總體均小于0.1，其中平均絕對(duì)誤差MAE表現(xiàn)最好，其余3種相差不大?？傮w來(lái)說(shuō)，使用遙感數(shù)據(jù)與SMAR模型結(jié)合，能夠較好地模擬出區(qū)域根系層的土壤濕度，為更大尺度預(yù)測(cè)根系層土壤水分提供了可行性參考，并為干旱預(yù)報(bào)和其他水文模擬業(yè)務(wù)提供了參考。

圖7 老哈河流域逐像元SAMR模型估算的根系層相對(duì)土壤濕度相關(guān)分析

4 結(jié)論與展望

根系層土壤濕度是陸地—大氣相互作用中的一個(gè)重要變量，尤其對(duì)干旱半干旱區(qū)土壤濕度的研究具有重要意義。本研究介紹了將遙感數(shù)據(jù)與適用于半干旱環(huán)境的SMAR模型結(jié)合有效估算區(qū)域根系層土壤濕度的方法，選取半干旱區(qū)的老哈河流域作為研究對(duì)象，SMAR模型參數(shù)是在站點(diǎn)遺傳算法優(yōu)化的基礎(chǔ)上，與土壤物理屬性、NDVI和ET建立多元線性回歸函數(shù)得到，并通過(guò)雙尾t，p<0.05的顯著性檢驗(yàn)。再將累積概率分布函數(shù)(CDF)重調(diào)整后的遙感數(shù)據(jù)應(yīng)用到SMAR模型中，有效獲得了老哈河流域根系層土壤濕度的時(shí)空分布，從而探究出SMAR模型在半干旱流域空間上具有良好的應(yīng)用性能，可以為大范圍估算根系層SM提供支撐，并且這種模型預(yù)測(cè)能力也可以為更大范圍內(nèi)估計(jì)干旱風(fēng)險(xiǎn)提供關(guān)鍵的資源。

未來(lái)SMAR模型參數(shù)估算可以加上盡可能多的可用物理信息以及變量獲得更良好的結(jié)果，比如葉面積指數(shù)(LAI)、增強(qiáng)植被指數(shù)(EVI)等，使得SMAR模型與遙感數(shù)據(jù)結(jié)合更有效地進(jìn)行大范圍的空間操作，也能更好地運(yùn)用于干旱半干旱區(qū)農(nóng)業(yè)規(guī)劃、干旱監(jiān)測(cè)和其他水文模型。