基于遙感數(shù)據(jù)同化的土壤含鹽量估算方法

張智韜 黃小魚(yú) 陳欽達(dá) 張珺銳 臺(tái) 翔 韓 佳

(1.西北農(nóng)林科技大學(xué)旱區(qū)農(nóng)業(yè)水土工程教育部重點(diǎn)實(shí)驗(yàn)室， 陜西楊凌 712100；2.西北農(nóng)林科技大學(xué)水利與建筑工程學(xué)院， 陜西楊凌 712100)

0 引言

土壤鹽漬化是導(dǎo)致干旱半干旱地區(qū)土地退化的重要成因，土地退化會(huì)引起農(nóng)作物減產(chǎn)，進(jìn)而影響國(guó)家和社會(huì)的可持續(xù)性發(fā)展[1]。因此，準(zhǔn)確連續(xù)地監(jiān)測(cè)土壤含鹽量對(duì)灌區(qū)的農(nóng)業(yè)和經(jīng)濟(jì)發(fā)展有著十分重要的現(xiàn)實(shí)意義。

衛(wèi)星遙感已經(jīng)成為監(jiān)測(cè)土壤含鹽量的重要手段[2]，遙感反演具有快速、準(zhǔn)確地大面積監(jiān)測(cè)土壤含鹽量的優(yōu)點(diǎn)[3]，但是其在時(shí)間上不連續(xù)，只能反映衛(wèi)星過(guò)境時(shí)刻的瞬時(shí)狀態(tài)，并且無(wú)法揭示土壤鹽分運(yùn)移的內(nèi)在機(jī)制以及其與氣象、土壤、植被等因素的關(guān)系。鹽分運(yùn)移模型通過(guò)輸入氣象數(shù)據(jù)、土壤參數(shù)等，在時(shí)間序列上能連續(xù)模擬單點(diǎn)尺度的土壤含鹽量[4]，但是其受模型輸入數(shù)據(jù)的不確定性限制，導(dǎo)致模型誤差在時(shí)間上不斷累積，因而模擬精度不高。數(shù)據(jù)同化能夠把遙感信息和鹽分運(yùn)移模型各自的優(yōu)勢(shì)相結(jié)合，能較為準(zhǔn)確地連續(xù)監(jiān)測(cè)區(qū)域土壤含鹽量的時(shí)空變化情況，是精準(zhǔn)農(nóng)業(yè)的發(fā)展趨勢(shì)。

數(shù)據(jù)同化就是通過(guò)引入時(shí)間序列的觀測(cè)數(shù)據(jù)來(lái)校正模型模擬過(guò)程，從而提高模型模擬精度[5]。陳鶴等[6]構(gòu)建集合卡爾曼濾波數(shù)據(jù)同化方案，將其輸入到HELP模型中，有效地提高了土壤水分的模擬精度。王澤人等[7]以集合卡爾曼濾波為基礎(chǔ)，構(gòu)建太湖葉綠素a質(zhì)量濃度同化試驗(yàn)系統(tǒng)，提高了太湖藍(lán)藻水華的預(yù)測(cè)精度。上述研究是基于實(shí)測(cè)數(shù)據(jù)在單點(diǎn)尺度的數(shù)據(jù)同化，而未基于遙感數(shù)據(jù)開(kāi)展區(qū)域尺度的數(shù)據(jù)同化研究，因此同化結(jié)果具有一定的局限性。文獻(xiàn)[8-11]采用同化算法把遙感數(shù)據(jù)同化到模型中，開(kāi)展區(qū)域尺度數(shù)據(jù)同化研究，多應(yīng)用在產(chǎn)量估計(jì)、土壤水分等方面，關(guān)于區(qū)域尺度的土壤含鹽量同化研究較少。

本文以高分一號(hào)衛(wèi)星影像為數(shù)據(jù)源，結(jié)合研究區(qū)不同深度時(shí)間序列的實(shí)測(cè)土壤含鹽量，以嶺回歸法遙感反演的土壤含鹽量模型作為觀測(cè)算子，采用集合卡爾曼濾波同化算法將遙感反演值引入到HYDRUS-1D模型中，對(duì)不同深度的土壤鹽分運(yùn)移過(guò)程進(jìn)行同化模擬，通過(guò)對(duì)集合數(shù)進(jìn)行敏感性分析，確定單點(diǎn)尺度數(shù)據(jù)同化方案的最佳集合數(shù)，從而得到最佳的單點(diǎn)尺度監(jiān)測(cè)土壤含鹽量數(shù)據(jù)同化方案，并由此擴(kuò)展到區(qū)域尺度監(jiān)測(cè)土壤含鹽量的數(shù)據(jù)同化研究，采用同化遙感數(shù)據(jù)為連續(xù)準(zhǔn)確地監(jiān)測(cè)土壤含鹽量時(shí)空變化提供了可行方案。

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)位于內(nèi)蒙古河套灌區(qū)沙壕渠灌域(40°52′～41°00′N(xiāo)，107°05′～107°10′E)，形狀類(lèi)似于狹長(zhǎng)型的倒三角形，整體地勢(shì)呈南高北低。沙壕渠灌域地處干旱半干旱地區(qū)，屬于典型的溫帶大陸性氣候，灌域內(nèi)常見(jiàn)作物有小麥、葵花、玉米等，灌溉用水多數(shù)來(lái)自于過(guò)境的黃河水，年平均降水量約為140 mm，年平均蒸發(fā)量約為2 000 mm，降水少、蒸發(fā)大的特點(diǎn)以及長(zhǎng)期的不合理灌溉方式，使得研究區(qū)的土壤鹽漬化問(wèn)題日趨嚴(yán)重。

1.2 田間實(shí)測(cè)數(shù)據(jù)

本研究在沙壕渠灌域內(nèi)共設(shè)置120個(gè)采樣點(diǎn)，研究區(qū)及采樣點(diǎn)分布如圖1所示。采用五點(diǎn)法對(duì)每個(gè)采樣點(diǎn)深度0～20 cm、20～40 cm、40～60 cm進(jìn)行分層采樣，用手持GPS儀器記錄該采樣點(diǎn)的具體位置信息，將土樣放置到已編號(hào)的鋁盒中，并帶回實(shí)驗(yàn)室。采用干燥法測(cè)定土壤含水率。把土樣研磨配制成土水比為1∶5溶液，經(jīng)過(guò)攪拌、靜置、沉淀、過(guò)濾等操作后，并使用電導(dǎo)率儀(DDS-307A型，上海佑科儀器公司)測(cè)定土壤電導(dǎo)率，取每個(gè)采樣點(diǎn)的5個(gè)電導(dǎo)率平均值作為該采樣點(diǎn)的電導(dǎo)率，并通過(guò)經(jīng)驗(yàn)公式[12]計(jì)算土壤含鹽量(SSC，%)。采樣時(shí)間分別為2019年4月15日、5月1日、5月15日、6月1日、6月15日、7月1日、7月15日、8月1日、8月15日、9月1日和9月15日，共11次。

圖1 研究區(qū)及采樣點(diǎn)分布圖Fig.1 Map of study area and distribution of sampling points

1.3 遙感反演土壤含鹽量

遙感數(shù)據(jù)來(lái)自于中國(guó)衛(wèi)星資源應(yīng)用中心(http:∥www.cresda.com/CN/)的高分一號(hào)衛(wèi)星。遙感數(shù)據(jù)預(yù)處理包括幾何精校正、輻射校正、大氣校正、裁剪等過(guò)程，并提取波段反射率，進(jìn)而計(jì)算光譜指數(shù)。本研究采用4個(gè)波段反射率組合的光譜指數(shù)，共計(jì)24個(gè)，分別為SI-T、BI、SI、SI1、SI2、SI3[13]、S1、S2、S3、S5、S6、NDSI[14]等12個(gè)鹽分指數(shù)和CRSI[15]、SR、NDVI、EVI、DVI、TVI、SRVI[16]、MSAVI[17]、ARVI[18]、NDWI[19]、NDGI[20]、GRVI[21]等12個(gè)植被指數(shù)。在模型構(gòu)建過(guò)程中，鹽分指數(shù)和植被指數(shù)分別經(jīng)過(guò)灰度關(guān)聯(lián)法篩選，選取排序前4的鹽分指數(shù)和植被指數(shù)的光譜指數(shù)組合作為輸入反演模型的自變量，以不同深度的實(shí)測(cè)土壤含鹽量為因變量，把所有采樣點(diǎn)的土壤含鹽量按照從小到大排序，每隔2個(gè)樣本選取出1個(gè)樣本作為驗(yàn)證集，其余則為建模集，分別對(duì)建模集和驗(yàn)證集構(gòu)建基于不同深度土壤含鹽量的嶺回歸法反演模型，得到不同深度的土壤含鹽量遙感反演值。衛(wèi)星遙感影像獲取的時(shí)間分別為2019年5月1日、5月13日、5月30日、6月28日、7月1日、7月30日、8月7日、8月15日、9月1日和9月18日。

1.4 HYDRUS-1D模型

HYDRUS-1D模型適用于模擬一維條件下飽和-非飽和多孔介質(zhì)的水分、能量、溶質(zhì)運(yùn)移以及根系吸水等[22]。本研究主要使用HYDRUS-1D模型的水分運(yùn)移模塊和溶質(zhì)運(yùn)移模塊，不考慮根系吸水模塊。在開(kāi)展土壤鹽分運(yùn)移模擬時(shí)，需要輸入的數(shù)據(jù)包括氣象數(shù)據(jù)、地下水位、初始條件及模型相關(guān)參數(shù)等。氣象數(shù)據(jù)來(lái)自于中國(guó)氣象數(shù)據(jù)網(wǎng)(http:∥data.cma.cn/)，包括降水量、蒸發(fā)量、氣溫、風(fēng)速、日照時(shí)數(shù)等氣象要素。選取河套灌區(qū)內(nèi)的杭錦后旗、臨河、五原、磴口、烏拉特前旗等5個(gè)氣象站，對(duì)這些站點(diǎn)的逐日氣象要素進(jìn)行克里金插值。根據(jù)氣象數(shù)據(jù)，使用Penman-Monteith公式計(jì)算潛在蒸散發(fā)量[23]。模型的土壤水力特征參數(shù)、溶質(zhì)運(yùn)移參數(shù)等的初始值基于已有的研究選定[24-25]，通過(guò)比較不同深度的實(shí)測(cè)土壤含鹽量與模擬值的偏差，不斷地對(duì)模型參數(shù)值進(jìn)行調(diào)整，從而確定HYDRUS-1D模型主要參數(shù)的具體值。選取大氣邊界作為水分運(yùn)移模塊的上邊界，變壓力水頭作為水分運(yùn)移模塊的下邊界；溶質(zhì)運(yùn)移模塊的上、下邊界條件均為濃度邊界。初始條件為初始時(shí)刻實(shí)測(cè)的土壤含水率和土壤電導(dǎo)率。

1.5 集合卡爾曼濾波數(shù)據(jù)同化方案

1.5.1集合卡爾曼濾波原理

集合卡爾曼濾波作為順序數(shù)據(jù)同化算法之一，它以蒙特卡羅方法為基礎(chǔ)，采用一組符合高斯分布的隨機(jī)變量集合，通過(guò)模型向前積分，計(jì)算下一時(shí)刻狀態(tài)量的統(tǒng)計(jì)特征[26]。它很好地解決了卡爾曼濾波在實(shí)際應(yīng)用中誤差協(xié)方差矩陣的估計(jì)，克服了卡爾曼濾波線性化的問(wèn)題[27]。

(1)進(jìn)行初始化。假設(shè)有N個(gè)集合數(shù)，在t時(shí)刻，添加高斯白噪聲，從而對(duì)土壤含鹽量進(jìn)行初始擾動(dòng)，得到符合高斯分布的狀態(tài)變量集合Xi,0(i=0，1，…，N)。

(2)把t時(shí)刻的土壤含鹽量集合輸入到HYDRUS-1D模型中，得到t+1時(shí)刻的土壤含鹽量預(yù)測(cè)值集合，計(jì)算式為

(1)

Mt,t+1——模型算子，本研究為HYDRUS-1D模型

ωi,t——模型誤差，服從期望為0、方差為Qt的高斯分布

Qt——t時(shí)刻模型誤差擾動(dòng)所對(duì)應(yīng)的誤差協(xié)方差矩陣

(3)計(jì)算t+1時(shí)刻土壤含鹽量預(yù)測(cè)值集合的均值和誤差協(xié)方差矩陣，計(jì)算式為

(2)

(3)

(4)計(jì)算t+1時(shí)刻卡爾曼增益，計(jì)算式為

(4)

(5)

(6)

式中Kt+1——t+1時(shí)刻的卡爾曼增益

Ht+1——t+1時(shí)刻的觀測(cè)算子

Rt+1——t+1時(shí)刻觀測(cè)誤差擾動(dòng)所對(duì)應(yīng)的誤差協(xié)方差矩陣

(5)計(jì)算t+1時(shí)刻土壤含鹽量分析值集合的均值和誤差協(xié)方差矩陣，計(jì)算式為

(7)

(8)

(9)

νi,t+1——觀測(cè)誤差，服從期望為0、方差為Rt+1的高斯分布

1.5.2不同深度土壤含鹽量集合卡爾曼濾波數(shù)據(jù)同化方案

數(shù)據(jù)同化系統(tǒng)由觀測(cè)算子、模型算子和同化算法3部分組成，在本研究中，觀測(cè)算子為時(shí)間序列上基于不同深度土壤含鹽量的嶺回歸法反演模型，模型算子為HYDRUS-1D模型，模擬時(shí)間序列上不同深度土壤鹽分運(yùn)移，同化算法則為耦合觀測(cè)算子(遙感數(shù)據(jù))和模型算子(模型模擬數(shù)據(jù))的集合卡爾曼濾波，在考慮集合數(shù)、觀測(cè)誤差以及模型誤差等同化條件的基礎(chǔ)上，通過(guò)同化遙感數(shù)據(jù)到HYDRUS-1D模型中，開(kāi)展單點(diǎn)尺度不同深度土壤含鹽量數(shù)據(jù)同化方案，進(jìn)而將單點(diǎn)同化方案擴(kuò)展到不同深度的區(qū)域尺度土壤含鹽量數(shù)據(jù)同化方案。其流程如圖2所示。

圖2 數(shù)據(jù)同化流程圖Fig.2 Flow diagram of data assimilation

(1)開(kāi)始進(jìn)行HYDRUS-1D模型模擬，得到t時(shí)刻不同深度土壤含鹽量模擬值。

(2)t+1時(shí)刻分2種情況：①當(dāng)t+1時(shí)刻不存在反演值時(shí)，t時(shí)刻的不同深度土壤含鹽量模擬值繼續(xù)輸入HYDRUS-1D模型中，得到t+1時(shí)刻的不同深度土壤含鹽量模擬值。②當(dāng)t+1時(shí)刻存在反演值時(shí)，開(kāi)始數(shù)據(jù)同化的計(jì)算流程。對(duì)t時(shí)刻不同深度土壤含鹽量模擬值進(jìn)行初始化，并輸入到HYDRUS-1D模型中，再加上模型誤差生成不同深度土壤含鹽量預(yù)測(cè)值集合；利用t+1時(shí)刻不同深度土壤含鹽量反演值加上觀測(cè)誤差生成不同深度土壤含鹽量觀測(cè)值集合；通過(guò)預(yù)測(cè)值集合和觀測(cè)值集合計(jì)算卡爾曼增益，更新得到不同深度土壤含鹽量分析值集合，并計(jì)算其均值，即為同化值。

進(jìn)入下一時(shí)刻，重復(fù)步驟(2)的計(jì)算流程，直至整個(gè)同化研究結(jié)束。

本研究以不同深度土壤含鹽量作為數(shù)據(jù)同化方案的狀態(tài)變量，以多時(shí)相嶺回歸法反演模型為觀測(cè)算子，通過(guò)集合卡爾曼濾波同化算法將其引入HYDRUS-1D模型中，HYDRUS-1D模型模擬時(shí)間為2019年4月15日至9月21日，共160 d，在2019年4月15日對(duì)HYDRUS-1D模型輸入初始時(shí)刻不同深度土壤含鹽量實(shí)測(cè)值，模擬時(shí)間序列上不同深度土壤鹽分運(yùn)移。模型運(yùn)行到2019年5月1日，即第1次出現(xiàn)遙感反演值時(shí)，開(kāi)始進(jìn)行數(shù)據(jù)同化，同化時(shí)間取為2019年5月1日至9月21日，共144 d。

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

本研究采用決定系數(shù)R2、相對(duì)誤差(RE)對(duì)遙感反演結(jié)果進(jìn)行評(píng)價(jià)；采用均方根誤差(RMSE)、同化效率系數(shù)(EFF)[28]和歸一化還原誤差指數(shù)(NER)[29]對(duì)單點(diǎn)尺度同化結(jié)果進(jìn)行評(píng)價(jià)；采用相關(guān)系數(shù)r和NER對(duì)區(qū)域尺度同化結(jié)果進(jìn)行評(píng)價(jià)。其中，EFF和NER的計(jì)算公式分別為

(10)

(11)

式中EFF——同化效率系數(shù)

NER——?dú)w一化還原誤差指數(shù)

Ot——t時(shí)刻的土壤含鹽量實(shí)測(cè)值

EFF的取值范圍為-∞～1，當(dāng)EFF大于0時(shí)，HYDRUS-1D模型模擬結(jié)果的精度提高，當(dāng)EFF小于0時(shí)，HYDRUS-1D模型模擬結(jié)果精度降低。EFF越接近1，精度越高。

NER的取值范圍為-∞～1，當(dāng)NER大于0時(shí)，表明改善了HYDRUS-1D模型模擬結(jié)果；當(dāng)NER小于0時(shí)，表明降低了HYDRUS-1D模型模擬結(jié)果。NER越接近1，精度越高。

2 結(jié)果與分析

2.1 基于光譜指數(shù)的土壤含鹽量估算

以2019年6月28日的高分一號(hào)衛(wèi)星遙感圖像為例，使用灰色系統(tǒng)把120個(gè)采樣點(diǎn)不同深度土壤含鹽量實(shí)測(cè)值分別與12個(gè)鹽分指數(shù)、12個(gè)植被指數(shù)進(jìn)行灰度關(guān)聯(lián)分析，它們的灰色關(guān)聯(lián)度和排序如表1、2所示。本研究分別選取關(guān)聯(lián)度排序前4的鹽分指數(shù)和植被指數(shù)構(gòu)成的光譜指數(shù)組合作為自變量，以不同深度的實(shí)測(cè)土壤含鹽量作為因變量，構(gòu)建基于嶺回歸法的不同深度土壤含鹽量反演模型。模型的建模和驗(yàn)證結(jié)果如表3所示。

表1 不同深度鹽分指數(shù)與土壤含鹽量的灰色關(guān)聯(lián)度Tab.1 Gray correlation degree between salinity index and soil salinity at different depths

表2 不同深度植被指數(shù)與土壤含鹽量的灰色關(guān)聯(lián)度Tab.2 Gray correlation degree between vegetation index and soil salinity at different depths

表3 基于不同深度土壤含鹽量的嶺回歸法模型Tab.3 Ridge regression models based on soil salinity at different depths

由表3可知，同一深度上，嶺回歸法模型的建模集和驗(yàn)證集的決定系數(shù)R2和相對(duì)誤差RE相差不大，并且隨著深度的增加，模型精度隨之降低。0～20 cm模型反演精度最高，其次是20～40 cm，模型反演精度最差的是40～60 cm。在3個(gè)深度上，嶺回歸法模型的建模集和驗(yàn)證集的R2均達(dá)到0.64以上，擬合精度較高，建模集和驗(yàn)證集的RE介于0.14～0.22之間，表明使用嶺回歸法模型得到的遙感反演值與實(shí)測(cè)值的誤差相對(duì)較小，模型反演精度較高。

使用嶺回歸法模型得到其余時(shí)刻的不同深度土壤含鹽量遙感反演值，并把時(shí)間序列上的不同深度土壤含鹽量遙感反演模型作為數(shù)據(jù)同化方案的觀測(cè)算子，開(kāi)展不同深度的土壤含鹽量數(shù)據(jù)同化研究。

2.2 單點(diǎn)數(shù)據(jù)同化方案

2.2.1集合數(shù)的敏感性分析

基于不同深度土壤含鹽量的集合卡爾曼濾波數(shù)據(jù)同化方案，集合數(shù)是影響數(shù)據(jù)同化模擬精度的重要因素。本研究設(shè)置的集合數(shù)分別為10、25、50、75、100，對(duì)不同深度的土壤鹽分運(yùn)移進(jìn)行同化模擬，分析集合數(shù)對(duì)同化結(jié)果的影響，如圖3所示。

圖3 不同集合數(shù)對(duì)同化結(jié)果的影響Fig.3 Effects of different ensemble numbers on assimilation results

由圖3a可知，隨著集合數(shù)的增大，每個(gè)深度土壤含鹽量同化值的均方根誤差RMSE都隨之減小。當(dāng)集合數(shù)小于等于50時(shí)，RMSE降低速度較大；當(dāng)集合數(shù)大于50時(shí)，RMSE降低速度相對(duì)較小。當(dāng)集合數(shù)為10和25時(shí)，RMSE均高于0.01%；當(dāng)集合數(shù)為50時(shí)，RMSE介于0.006%～0.011%之間；當(dāng)集合數(shù)為50、75和100時(shí)，其RMSE相差不大。由圖3b可知，每個(gè)深度土壤含鹽量同化值的同化效率系數(shù)EFF都大于0，說(shuō)明在不同的集合數(shù)條件下，均具有同化效果。各深度的EFF都隨著集合數(shù)增加而增加，增加的速度在減小，當(dāng)集合數(shù)為10時(shí)，同化效果最差，各深度的EFF都小于0.55；當(dāng)集合數(shù)為25時(shí)，同化效果相對(duì)較差；當(dāng)集合數(shù)為50時(shí)，各深度的EFF均達(dá)到0.8以上，表明經(jīng)過(guò)數(shù)據(jù)同化后，提高了HYDRUS-1D模型模擬精度，同化效果較好；當(dāng)集合數(shù)為75和100時(shí)，各深度的EFF與集合數(shù)為50時(shí)的EFF相差不大。由圖3c可知，每個(gè)深度土壤含鹽量同化值的歸一化還原誤差指數(shù)NER均大于0，說(shuō)明在不同的集合數(shù)條件下，土壤含鹽量模擬精度均得以提高。各深度的NER變化趨勢(shì)表現(xiàn)為隨著集合數(shù)的增加而增加，增加的趨勢(shì)在減緩，當(dāng)集合數(shù)為10和25時(shí)，各深度的NER均低于0.6，表明同化方案對(duì)改善模型模擬精度的能力較差；隨著集合數(shù)的增加，當(dāng)集合數(shù)為50時(shí)，NER介于0.61～0.73之間；當(dāng)集合數(shù)為75和100時(shí)，與集合數(shù)為50時(shí)的NER相差不大。綜上所述，為了提高同化運(yùn)算效率和滿足同化精度的要求，本研究的最佳集合數(shù)為50。

由圖3可知，當(dāng)集合數(shù)相同時(shí)，各深度同化值的RMSE均表現(xiàn)為隨著深度的增加而增大，而各深度同化值的EFF和NER都表現(xiàn)為隨著深度的增加而減小。在3個(gè)深度上，0～20 cm同化值的RMSE最小，EFF最大，NER最大，0～20 cm同化效果最好；20～40 cm次之；40～60 cm最差，其RMSE最大，EFF最小，NER最小。故同化精度是隨著深度的增加而降低。

2.2.2單點(diǎn)同化結(jié)果分析

根據(jù)上述研究，集合數(shù)取為50，觀測(cè)誤差和模型誤差分別設(shè)置為0.01、0.01，對(duì)單點(diǎn)尺度不同深度土壤含鹽量開(kāi)展集合卡爾曼濾波同化研究。基于不同深度的土壤含鹽量實(shí)測(cè)數(shù)據(jù)(簡(jiǎn)稱(chēng)實(shí)測(cè)值)，采用HYDRUS-1D模型模擬得到土壤含鹽量模擬值(簡(jiǎn)稱(chēng)模擬值)、采用嶺回歸法模型反演得到土壤含鹽量遙感反演值(簡(jiǎn)稱(chēng)反演值)、采用集合卡爾曼濾波同化算法得到土壤含鹽量同化值(簡(jiǎn)稱(chēng)同化值)，比較不同深度土壤含鹽量的模擬值、反演值、同化值和實(shí)測(cè)值，如圖4所示，并進(jìn)行誤差分析，如表4所示。

圖4 不同深度土壤含鹽量的模擬值、反演值、同化值與實(shí)測(cè)值的比較Fig.4 Comparisons of simulated values, inversion values, assimilation values and measured values of soil salinity at different depths

由圖4可知，在同一深度上，模擬值與實(shí)測(cè)值的誤差最大，其次是反演值，同化值最接近實(shí)測(cè)值；在不同深度上，模擬值、反演值和同化值與實(shí)測(cè)值的誤差都隨著深度的增加而增加。由表4可知：①?gòu)耐`差的角度分析，在同一深度上，與模擬值、反演值相比，同化值與實(shí)測(cè)值的誤差最小，其在3個(gè)深度的RMSE均低于0.011%，反演值次之，反演值的RMSE與同化值的RMSE相差不大，模擬值與實(shí)測(cè)值的誤差最大，其在3個(gè)深度的RMSE為0.024%～0.028%；在3個(gè)深度上，同化值、反演值和模擬值的誤差隨著深度的增加而增加，表明無(wú)論是模型模擬、遙感反演還是同化模擬對(duì)土壤含鹽量的作用都是表層優(yōu)于深層。②就同化效率系數(shù)而言，在同一深度上，同化值和反演值的EFF大于0.65，并且同化值的EFF更接近于1，表明引入反演值到HYDRUS-1D模型中進(jìn)行同化模擬，能較好地改善模型模擬精度，使同化值更接近實(shí)測(cè)值；在3個(gè)深度上，同化值和反演值的EFF都為正值，同化值的EFF大于反演值的EFF，它們的EFF隨著深度的增加而降低，表明引入遙感數(shù)據(jù)對(duì)0～20 cm土壤含鹽量的同化作用優(yōu)于20～40 cm、40～60 cm。③從提高模型精度程度來(lái)看，在同一深度上，同化值和反演值的NER均大于0.4，并且同化值大于反演值，表明同化遙感數(shù)據(jù)，使得模型的模擬軌跡更加接近實(shí)測(cè)值，減小了模擬值與實(shí)測(cè)值的誤差，提高了HYDRUS-1D模型的模擬精度；在3個(gè)深度上，同化值和反演值的NER都為正值，隨著深度的增加，同化值和反演值的NER逐漸減少，表明同化遙感數(shù)據(jù)對(duì)每個(gè)深度土壤含鹽量的模擬精度均能提高，并且對(duì)于每個(gè)深度的土壤含鹽量，同化作用隨著深度的增加而減小。

表4 不同深度土壤含鹽量的模擬值、反演值、同化值與實(shí)測(cè)值的誤差分析Tab.4 Error analysis of simulated values, inversion values, assimilation values and measured values of soil salinity at different depths

綜上所述，在同一深度上，同化值與實(shí)測(cè)值最為接近，反演值次之，模擬值和實(shí)測(cè)值相差最大。在3個(gè)深度上，同化值、反演值和模擬值的精度都是隨著深度的增加而降低。在單點(diǎn)尺度上，開(kāi)展不同深度的集合卡爾曼濾波數(shù)據(jù)同化方案，結(jié)合了遙感反演和模型模擬各自的優(yōu)勢(shì)，通過(guò)集合卡爾曼濾波同化算法把時(shí)間序列上的不同深度土壤含鹽量反演值應(yīng)用到HYDRUS-1D模型，提高了土壤剖面模擬精度，得到了較為準(zhǔn)確的同化值。并且引入遙感數(shù)據(jù)對(duì)各深度土壤含鹽量的同化作用隨著深度的增加而減小，0～20 cm土壤含鹽量的同化模擬優(yōu)于20～40 cm，40～60 cm最差。

2.3 區(qū)域數(shù)據(jù)同化方案

本研究開(kāi)展區(qū)域尺度的不同深度土壤含鹽量集合卡爾曼濾波數(shù)據(jù)同化模擬，在不同采樣時(shí)間，基于區(qū)域尺度的實(shí)測(cè)值、反演值、模擬值以及同化值，通過(guò)r和NER計(jì)算它們與實(shí)測(cè)值之間的誤差，以2019年6月15日為例，由表5可知，在3個(gè)深度上，模擬值的精度隨著深度的增加而減小，表明HYDRUS-1D模型對(duì)于0～20 cm的土壤含鹽量模擬精度較高，而20～40 cm、40～60 cm的土壤含鹽量模擬精度相對(duì)較差。在3個(gè)深度上，反演值的精度隨著深度的增加而減少，表明嶺回歸法模型對(duì)0～20 cm的土壤含鹽量反演精度最高，20～40 cm反演精度次之，40～60 cm深度最差，但是20～40 cm、40～60 cm反演值的r和NER均較大，它們的模型反演精度較高。在3個(gè)深度上，同化值的同化精度隨著深度的增加而降低，表明區(qū)域尺度的集合卡爾曼濾波數(shù)據(jù)同化方案對(duì)于0～20 cm土壤含鹽量的同化效果最好，高于0～20 cm的反演值和模擬值。與0～20 cm相比，20～40 cm和40～60 cm同化值的r和NER相對(duì)較差，但仍?xún)?yōu)于同一深度的反演值和模擬值。它們的r和NER都是隨著深度的增加而減小，表明0～20 cm同化值、反演值、模擬值的精度均為最高，它們?cè)?0～40 cm精度相對(duì)較差，在40～60 cm精度最差。在同一深度上，同化值與實(shí)測(cè)值的誤差最小，其r為0.94～0.97，NER為0.61～0.69，反演值的精度次之，其與同化值的r相差較小，其N(xiāo)ER為0.31～0.35，模擬值與實(shí)測(cè)值的誤差最大，其r為0.65～0.72。在3個(gè)深度上，反演值和同化值的NER均大于0，并且同化值遠(yuǎn)高于反演值，表明區(qū)域尺度的集合卡爾曼濾波數(shù)據(jù)同化方案把遙感數(shù)據(jù)引入HYDRUS-1D模型中，能較好地監(jiān)測(cè)區(qū)域尺度每個(gè)深度的土壤含鹽量。

表5 區(qū)域尺度不同深度土壤含鹽量模擬值、反演值、同化值與實(shí)測(cè)值的誤差分析Tab.5 Error analysis of simulated values, inversion values, assimilation values and measured values of soil salinity at different depths in regional scale

故基于集合卡爾曼濾波同化算法開(kāi)展區(qū)域尺度的數(shù)據(jù)同化研究，使得在區(qū)域尺度上，同化值優(yōu)于反演值和模擬值，能更好地模擬區(qū)域尺度時(shí)間序列上不同深度土壤含鹽量時(shí)空變化，從而更加接近實(shí)際情況。

以2019年6月15日為例，不同深度土壤含鹽量同化圖如圖5所示。

圖5 不同深度土壤含鹽量同化圖Fig.5 Assimilation maps of soil salinity at different depths

3 討論

在單點(diǎn)尺度上，與不同深度的模擬值和反演值相比，同化值與實(shí)測(cè)值更為接近，這是因?yàn)榛诩峡柭鼮V波同化算法把遙感數(shù)據(jù)同化到HYDRUS-1D模型中，提高了土壤剖面模型模擬精度，得到了更為準(zhǔn)確的同化值，使得同化值與實(shí)測(cè)值的誤差最小。同化精度隨著深度的增加而減小，0～20 cm同化效果最好，20～40 cm次之，40～60 cm最差。本研究探究了集合數(shù)對(duì)單點(diǎn)尺度集合卡爾曼濾波數(shù)據(jù)同化方案的敏感性，結(jié)果表明集合數(shù)對(duì)同化模擬過(guò)程較為敏感，在各深度上均表現(xiàn)為隨著集合數(shù)的增加，RMSE先大幅度降低，而EFF和NER則是先大幅度增加，然后達(dá)到相對(duì)穩(wěn)定，故最佳的集合數(shù)為50，既能保證同化精度，又能提高同化運(yùn)算效率。此外，還可以從同化頻率[30]、觀測(cè)誤差[31]和引入不同深度的實(shí)測(cè)數(shù)據(jù)[32]等方面對(duì)數(shù)據(jù)同化方案進(jìn)行敏感性分析，找到對(duì)數(shù)據(jù)同化方案的較為敏感性因素。

在區(qū)域尺度上，各深度的同化值與實(shí)測(cè)值的r均達(dá)到0.94以上，高于反演值和模擬值，并且各深度同化值的NER均在0.61以上，皆為正值，表明本研究能較好地應(yīng)用于區(qū)域尺度的同化模擬，得到較為準(zhǔn)確的時(shí)間序列區(qū)域尺度土壤含鹽量時(shí)空變化信息。在區(qū)域尺度的集合卡爾曼濾波數(shù)據(jù)同化方案中，同化精度最好的深度為0～20 cm，其次是20～40 cm，40～60 cm最差，與單點(diǎn)尺度數(shù)據(jù)同化方案的結(jié)論一致。這是因?yàn)殡S著深度的增加，HYDRUS-1D模型的輸入數(shù)據(jù)等對(duì)20～40 cm、40～60 cm土壤含鹽量模擬的作用較小，導(dǎo)致它們與實(shí)測(cè)值的誤差隨著深度的增加而增加，并且?guī)X回歸法模型的反演精度隨著深度的增加，也隨之降低，在二者作用下，同化精度表現(xiàn)為隨著深度的增加而降低。

本研究通過(guò)集合卡爾曼濾波同化遙感反演和模型模擬的土壤含鹽量，開(kāi)展區(qū)域尺度的不同深度土壤含鹽量數(shù)據(jù)同化研究，有效地提高了時(shí)間序列的區(qū)域尺度土壤剖面含鹽量模擬精度。一方面模型模擬彌補(bǔ)了遙感數(shù)據(jù)的不連續(xù)性問(wèn)題，另一方面同化遙感數(shù)據(jù)可以提高模型模擬精度，但是同化結(jié)果仍存在誤差，主要來(lái)源于遙感反演模型、HYDRUS-1D模型數(shù)據(jù)的不確定性和同化算法帶來(lái)的誤差。本研究使用的遙感反演模型為統(tǒng)計(jì)回歸模型，未來(lái)研究可以采用機(jī)器學(xué)習(xí)模型提高遙感監(jiān)測(cè)的精度[33]，把遙感數(shù)據(jù)引入到HYDRUS-1D模型中，對(duì)區(qū)域尺度的各深度土壤含鹽量進(jìn)行同化模擬，從而提高同化精度。在HYDRUS-1D模型中，土壤水力特征參數(shù)、溶質(zhì)運(yùn)移參數(shù)等都會(huì)對(duì)同化結(jié)果產(chǎn)生較大的影響。因此，為提高同化精度，不僅需提高反演模型精度，也需提高模型參數(shù)的精度，從而更為準(zhǔn)確地監(jiān)測(cè)區(qū)域尺度各深度土壤含鹽量的時(shí)空變化，使同化結(jié)果更具有可靠性。

本研究只考慮了土壤含鹽量作為狀態(tài)變量，在開(kāi)展區(qū)域尺度的同化模擬時(shí)，尚未考慮其他狀態(tài)變量隨時(shí)間變化的特性是否會(huì)對(duì)同化結(jié)果造成一定影響。因此，今后研究可以探討多狀態(tài)變量的數(shù)據(jù)同化方案，如同時(shí)同化狀態(tài)變量和模型參數(shù)[34-35]，進(jìn)一步準(zhǔn)確地提高區(qū)域尺度土壤含鹽量的同化精度，或是將來(lái)基于遙感數(shù)據(jù)與鹽分運(yùn)移模型的數(shù)據(jù)同化研究重點(diǎn)。此外，單一的同化算法仍存在一定的局限性，結(jié)合多種不同類(lèi)型的同化算法，如集合-變分?jǐn)?shù)據(jù)同化算法[36-37]，發(fā)展基于遙感數(shù)據(jù)的耦合多種同化算法的數(shù)據(jù)同化方案，或是未來(lái)利用數(shù)據(jù)同化方法獲取大面積土壤含鹽量時(shí)空信息的一個(gè)重要研究方向。并且，植被覆蓋和人工作業(yè)等影響含鹽量因素也是本研究考慮不足的地方。在植被覆蓋階段，遙感反演受植被覆蓋度、土壤質(zhì)地、作物種類(lèi)等多種因素的影響，從而導(dǎo)致同化精度的降低。對(duì)于人工作業(yè)，如灌溉等，對(duì)土壤含鹽量造成淋洗，同時(shí)它具有局部、粗化效應(yīng)，從而對(duì)同化結(jié)果造成一定的影響。上述不足有待進(jìn)一步深入研究。

4 結(jié)論

(1)在進(jìn)行遙感反演時(shí)，各深度嶺回歸法模型的R2均達(dá)到0.64以上，RE在0.14以上，反演精度較好，能得到比較準(zhǔn)確的不同深度土壤含鹽量遙感反演值。

(2)基于單點(diǎn)尺度的數(shù)據(jù)同化方案最佳集合數(shù)為50，不同深度土壤含鹽量同化值的RMSE均小于0.011%、EFF均大于0.8、NER均大于0.6，優(yōu)于反演值和模擬值，表明單點(diǎn)尺度集合卡爾曼濾波數(shù)據(jù)同化方案能提高土壤剖面含鹽量模型模擬精度。

(3)在區(qū)域尺度的不同深度土壤含鹽量集合卡爾曼濾波數(shù)據(jù)同化方案中，不同深度土壤含鹽量同化值與實(shí)測(cè)值的相關(guān)性較強(qiáng)，r均為0.94以上，NER均為0.61以上，能較好地模擬在時(shí)間序列上區(qū)域尺度的不同深度土壤含鹽量時(shí)空變化趨勢(shì)。

(4)集合卡爾曼濾波數(shù)據(jù)同化方案在單點(diǎn)尺度和區(qū)域尺度上，對(duì)3個(gè)深度的土壤含鹽量進(jìn)行數(shù)據(jù)同化，均取得較好的同化效果。對(duì)比3種不同深度同化結(jié)果可知，同化效果最佳的是0～20 cm，其次是20～40 cm，40～60 cm最差。