朱悅璐 吳奇俞

摘要：針對無資料區(qū)土體含水率數(shù)據(jù)難以獲取的問題，提出了一種基于衛(wèi)星反演-相空間重構(gòu)-非飽和入滲計算的組合方案，以研究區(qū)110 d土體表層含水率為基礎(chǔ)，預(yù)測未來100 d無資料時段土體表層及內(nèi)部含水率分布規(guī)律。計算結(jié)果表明：研究區(qū)含水率時間序列具備混沌特征，可由一維時間序列拓?fù)錇橐粋€嵌入維數(shù)m=5，遲滯τ=10的相空間，由該相空間預(yù)測的土體表層含水率在驗證期最小相對誤差為0.7%，最大相對誤差為2.4%，在預(yù)測期最小相對誤差為2.2%，最大相對誤差為8.3%，均滿足工程需求，因此將其用于后續(xù)非飽和入滲計算的邊界條件是真實(shí)有效的。該方案具有動力學(xué)特性和物理力學(xué)意義，可為無資料地區(qū)土體含水率估計借鑒。

關(guān)鍵詞：土體含水率； 非飽和入滲； 相空間重構(gòu)； 混沌理論； Richards方程； 非飽和土

中圖法分類號： TU43

文獻(xiàn)標(biāo)志碼： A

DOI：10.16232/j.cnki.1001-4179.2024.04.028

0引 言

基于事物已有信息準(zhǔn)確預(yù)測未來發(fā)展趨勢，一直以來都是科學(xué)探索最前沿的目標(biāo)之一［1-3］。在土壤水領(lǐng)域，預(yù)測土體表層及內(nèi)部含水率變化規(guī)律，是進(jìn)一步探索陸氣水量交換［4］、地下水動力學(xué)［5］、非飽和土力學(xué)［6］等很多相關(guān)課題的基礎(chǔ)，在工程中具有重要意義。

現(xiàn)階段，一種獲取長時間序列某區(qū)域土體含水率的有效手段，是通過大尺度衛(wèi)星遙感技術(shù)（例如高分5號光譜衛(wèi)星［7］、GF-1多光譜衛(wèi)星［8］、Sentinel -2衛(wèi)星等［9］），反演逐日土體表層含水率數(shù)據(jù)，經(jīng)整合處理后，最終生成高分辨率含水率數(shù)據(jù)庫。但受傳感器或軌道限制，現(xiàn)有的深空衛(wèi)星探測技術(shù)均存在衛(wèi)星掃描盲區(qū)（固定時段盲區(qū)或固定位置盲區(qū)），因此無法避免數(shù)據(jù)庫中出現(xiàn)資料不連續(xù)的時段，即無資料區(qū)，當(dāng)研究區(qū)無資料時段過多時，數(shù)據(jù)庫可信度就會下降，從而大大降低預(yù)測的精度。

現(xiàn)有針對無資料區(qū)含水率的數(shù)據(jù)修復(fù)，可分為野外單點(diǎn)實(shí)測和室內(nèi)數(shù)值模擬。傳統(tǒng)的野外監(jiān)測雖然精度較高，但同時亦具有成本高、耗時長、數(shù)據(jù)同步性差等缺點(diǎn)；而室內(nèi)數(shù)值模擬，大多通過統(tǒng)計手段插補(bǔ)無資料區(qū)表層土體含水率數(shù)據(jù)，再基于機(jī)器學(xué)習(xí)構(gòu)建不同深度的土體內(nèi)部含水率模型［10-11］，這導(dǎo)致預(yù)測結(jié)果大多僅具備統(tǒng)計意義而缺乏物理力學(xué)依據(jù)。

基于此，本次研究以衛(wèi)星反演生成的土體表層含水率數(shù)據(jù)庫為基礎(chǔ)，通過相空間重構(gòu)預(yù)測方法［12-14］，將一維含水率時間序列拓?fù)錇楦呔S相空間，在相空間內(nèi)對無資料區(qū)土體表層含水率進(jìn)行非線性動力學(xué)預(yù)測，并將預(yù)測結(jié)果代入非飽和入滲控制方程，最終構(gòu)建土體內(nèi)部含水率變化規(guī)律。相比于傳統(tǒng)研究，該方案土體表層含水率預(yù)測結(jié)果是基于實(shí)測數(shù)據(jù)的動力學(xué)反演，具備物理意義；土體內(nèi)部含水率預(yù)測結(jié)果是基于非飽和入滲控制方程迭代計算，具備力學(xué)意義。

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

1.1全球衛(wèi)星反演數(shù)據(jù)庫

本次研究以SMAP衛(wèi)星（Soil Moisture Active Passive）與AMSR-E（Advanced Microwave Scanning Radiometer）系列傳感器所獲的2002～2020年土體表層0～5 cm逐日全球高精度、長系列含水率（NNsm）為基礎(chǔ)數(shù)據(jù)集，其空間分辨率精度為0.3°×0.3°經(jīng)緯度，數(shù)據(jù)來源于“國家青藏高原科學(xué)數(shù)據(jù)中心”（http：∥data.tpdc.ac.cn）。

1.2本研究計算數(shù)據(jù)集

研究區(qū)地理中心坐標(biāo)為N34°22′55″，E107°17′01″，位于中國陜西省寶雞市遠(yuǎn)郊，人類活動影響因素較弱。經(jīng)資料分析，該地區(qū)受氣象條件、云層遮擋等因素影響，存在數(shù)據(jù)間斷性缺失，符合無資料區(qū)特征。具體計算中，選擇無資料時段前具有一段完整數(shù)據(jù)的時段，作為基礎(chǔ)數(shù)據(jù)集。若基礎(chǔ)數(shù)據(jù)集選擇過短，會導(dǎo)致相空間重構(gòu)不準(zhǔn)，而基礎(chǔ)數(shù)據(jù)集選擇過長，則會導(dǎo)致計算復(fù)雜度呈幾何倍數(shù)增長。基于以上原因，本次研究選取2020年第103～212天連續(xù)110組數(shù)據(jù)作為計算序列x（t），213～227天15組數(shù)據(jù)作為驗證序列，如表1所列。前人研究成果表明，這一長度的序列可滿足相空間重構(gòu)相關(guān)要求［15］。

2研究方法

本次研究預(yù)測無資料區(qū)土體剖面含水率的具體步驟為：① 通過相空間重構(gòu)，將一維含水率時間序列θ（t）拓?fù)錇楦呔S相空間，在相空間內(nèi)對土體表層含水率展開預(yù)測；② 將預(yù)測結(jié)果作為未來第一日土體邊界含水率θ01，結(jié)合測定的初始含水率θi1，代入入滲控制方程，計算第一日入滲曲線θ（z，t）；③ 當(dāng)t=24 h時，對應(yīng)的含水率曲線即為第二日土體初始含水率θi2；再結(jié)合相空間內(nèi)逐步預(yù)測的第二日表層含水率θ02，作為第三日計算基礎(chǔ)；④ 重復(fù)以上步驟，即可建立未來時期無資料地區(qū)土體剖面含水率預(yù)測模型。迭代計算示意如圖1所示。

2.1一維時間序列相空間重構(gòu)

對于一組一維時間序列x（t），t=1，2，…n，Packard等［16］提出一種相空間重構(gòu)方法，若該序列具有混沌特性，則可通過引入嵌入維數(shù)m和遲滯時間τ，將一維時間序列擴(kuò)展成多維相空間，以此為基礎(chǔ)預(yù)測原時間序列。一個標(biāo)準(zhǔn)的一維時間序列相空間重構(gòu)如式（1）所示。

當(dāng)m、τ均?。?，1）時，按照式（1）重構(gòu)規(guī)則，上述一維時間序列演化為一個二維相空間，其坐標(biāo)依次為（17.3，19.2），（19.2，16.7），（1.67，12.4）…，相點(diǎn)在二維相空間（x，y）的移動軌跡如圖3所示。

同理，當(dāng)m、τ均?。?，2）時，上述二維相空間進(jìn)一步變化為三維相空間，其坐標(biāo)依次為（17.3，19.2，16.7），（19.2，16.7，12.4），（16.7，12.4，18.3）…，相點(diǎn)在三維相空間（x，y，z）的移動軌跡如圖4所示。

在更高的維度下，相點(diǎn)分布及移動軌跡已不能用坐標(biāo)系直觀描述，但可通過數(shù)學(xué)計算求解。當(dāng)系統(tǒng)坐標(biāo)確定、演化軌跡已知時（例如圖中箭頭方向），即可根據(jù)演化軌跡對系統(tǒng)未來發(fā)展做出預(yù)測，顯然不同的嵌入維數(shù)m和遲滯τ所構(gòu)成的相空間包含的信息不同，因此如何選擇合理的重構(gòu)參數(shù)十分重要，以下小節(jié)將對這一內(nèi)容具體討論。

2.2相空間參數(shù)及混沌特性判斷

2.2.1遲滯τ計算

對于延遲時間的選擇，常用的手段為式（2）所示的自相關(guān)系數(shù)法，通過計算序列在不同延遲時間k=1，2…下的自相關(guān)系數(shù)rk，當(dāng)自相關(guān)系數(shù)rk第一次接近于0時，所對應(yīng)的延遲時間k即定義為遲滯τ。

2.2.2混沌特性判斷

只有具備混沌特性的序列，相空間重構(gòu)才有意義。最常見的混沌特性判別法為G-P算法［18］，即識別時間序列是否具有飽和關(guān)聯(lián)維數(shù)D2（m），具體計算方案為：任取一m（實(shí)際操作中一般從m=2向后連續(xù)取自然數(shù)），將m、τ代入式（1）生成一個相空間，在該相空間中，人為設(shè)定一個r值且滿足r>0。

由式（3）、（4）及Heaviside函數(shù)定義可知，Cmr的含義為相空間內(nèi)兩個相點(diǎn)距離小于r的概率。求得Cmr后，計算lnCmr和lnr，二者是一組確定的常數(shù)。由小到大重復(fù)設(shè)定不同的r值代入（3）式計算，在lnCr-lnr坐標(biāo)系中擬合曲線，此時曲線的直線段斜率即為關(guān)聯(lián)維數(shù)D2。

重新選擇不同的m，重復(fù)以上步驟，計算D2，若隨著m的增加，D2-m擬合曲線有水平漸近線，則表明原時間序列具有混沌特性；若隨著m增加，D2一直增加（或減?。?，擬合曲線無漸近線，則認(rèn)為原時間序列為一隨機(jī)序列，不具備混沌特性。

2.2.3嵌入維數(shù)m計算

若經(jīng)過判斷，時間序列具有混沌特性，則可取D2不再增加時所對應(yīng)的m為重構(gòu)相空間的嵌入維數(shù)。當(dāng)最終參數(shù)m、τ確定后，即將原一維時間序列x（t）擴(kuò)展成一個m維相空間，此時可在該相空間內(nèi)，對原序列進(jìn)行預(yù)測。

2.3相空間重構(gòu)數(shù)據(jù)預(yù)測

本文采用局部相點(diǎn)平均法對生成的相空間預(yù)測［20］，在式（1）所表達(dá)的相空間中，最后一個相點(diǎn)（最后一列）Xn的下一步演化結(jié)果顯然受到臨近相點(diǎn)的影響。因此，可尋找與Xn最接近的若干個相點(diǎn)，以這些相點(diǎn)下一步演化結(jié)果的平均值作為Xn的演化結(jié)果，記作Xn+1；新相點(diǎn)對應(yīng)向量最后一個分量xn+1即為原時間序列在t+1時刻的預(yù)測值。局部相點(diǎn)平均預(yù)測法原理如圖5（a）所示，算法如圖5（b）所示。圖中與Xn距離最近的3個相點(diǎn)分別為Xp、Xq、Xl，它們下一步演化為Xp+1、Xq+1、Xl+1，演化后相點(diǎn)平均值即為Xn的演化結(jié)果Xn+1。

在實(shí)際計算中，可采用式（4）計算相空間中每一列與最后一列的距離，并找出距離最短的若干列，用這些列對應(yīng)的下一列取平均值，即為新相點(diǎn)Xn+1，該新相點(diǎn)即為預(yù)測結(jié)果。以圖5（b）所示為例，新相點(diǎn)計算公式為

2.4非飽和入滲曲線迭代計算

預(yù)測得到的研究區(qū)土體表層含水率數(shù)據(jù)，即可作非飽和入滲迭代初始條件代入Richards方程進(jìn)行計算。該方程是研究非飽和土入滲的重要工具之一，在一定程度上可以揭示非飽和土入滲的主要特征，在非飽和土計算領(lǐng)域應(yīng)用廣泛［21］。本文僅針對第一類邊界條件展開探討，其控制方程及邊界條件、初始條件如下［22］：

3模型計算結(jié)果

3.1混沌方案預(yù)測表層土壤含水率

3.1.1遲滯計算

以式（2）計算表1含水率序列各階自相關(guān)系數(shù)，計算結(jié)果如圖6（a）所示，選取自相關(guān)系數(shù)第一次接近0時對應(yīng)的τ作為系統(tǒng)遲滯，由圖中可以看出，各階自相關(guān)系數(shù)隨時間t呈現(xiàn)波動狀態(tài)，當(dāng)t=10時，r10=0.134 1第一次接近于0，故本系統(tǒng)遲滯取τ=10。

3.1.2混沌特性判定及嵌入維數(shù)的選擇

本文選用嵌入維數(shù)m=2～10，應(yīng)用G-P算法進(jìn)行試算，代入式（3）～（4），在假設(shè)不同的m下變換r值，生成一組lnCr-lnr曲線，如圖7所示。

從圖中可以看出，m取不同值時，各曲線均存在直線段（圖中方框區(qū)域），且隨著m的增加，直線段傾角不斷增大，但均未超過斜率為tan 64°的漸近線，為更清晰表征這一特征，繪制直線段斜率D2與嵌入維數(shù)m關(guān)系曲線如圖6（b）。從圖中可以看出，D2-m曲線有水平漸近線D2=2.068，這表明關(guān)聯(lián)維數(shù)D2不會隨著m的增大一直增大，序列具有混沌特性，因此可取D2趨于穩(wěn)定時對應(yīng)的m=5為系統(tǒng)的最終嵌入維數(shù)。

3.1.3相空間重構(gòu)及預(yù)測

以含水率時間序列x（t）、τ=10、m=5，應(yīng)用式（1）重構(gòu)相空間，重構(gòu)的含水率相空間共有n-m-1τ=110-5-1×10=70個相點(diǎn)，每個相點(diǎn)有m=5個分維。即該相空間是一個5行70列矩陣。限于篇幅所限，本文僅按照重構(gòu)規(guī)則，羅列第1，15，30，50，70號相點(diǎn)，以供讀者檢驗。

上述5個影響相點(diǎn)經(jīng)過一步演化后的相點(diǎn)為X59、X69、X63、X30、X49，如式（10）所示。

表1數(shù)據(jù)庫可知，這一區(qū)域含水率在2020年第213天（即t=111時刻），含水率為0.136 7，當(dāng)取本文研究精度小數(shù)點(diǎn)后3位時，預(yù)測值與實(shí)測值在千分位上精確一致。采用同樣方法，預(yù)測未來15 d含水率數(shù)據(jù)，并與第213～227天的驗證序列進(jìn)行對比，對每組數(shù)據(jù)，采用工程中常用的相對誤差公式進(jìn)行計算，其結(jié)果如圖8所示。

可以看出，在驗證期的15組數(shù)據(jù)中，預(yù)測值與實(shí)測值最大相對誤差為2.4%（第225 d），最小相對誤差為0.7%（第213，214，215 d），對于一般工程，可以認(rèn)為該精度可以滿足工程需求。重復(fù)以上步驟，計算未來100 d土體表層含水率數(shù)據(jù)，如圖9所示。

由圖9可以看出，在預(yù)測期內(nèi)，大部分時段實(shí)測數(shù)據(jù)缺失，但在第12天、27天、45天、78天、94天仍有衛(wèi)星反演數(shù)據(jù)，將這些數(shù)據(jù)作為實(shí)測值應(yīng)用相對誤差公式與預(yù)測值進(jìn)行對比，其整理結(jié)果如表2所示。

通過表2最后一列相對誤差變化趨勢可以看出，在預(yù)測期內(nèi)，預(yù)測誤差隨著日期延長，呈現(xiàn)波動上升趨勢，這表明預(yù)測精度有所下降，但盡管如此，最大相對誤差仍僅為8.3%（預(yù)測期第94 d），仍小于一般工程中含水率預(yù)測精度在10%以內(nèi)的要求。因此可以認(rèn)為，在預(yù)測區(qū)間，基于相空間重構(gòu)的含水率表層預(yù)測方案可用于工程實(shí)際，同時也可以認(rèn)為，其作為后續(xù)逐日非飽和入滲計算的邊界條件是真實(shí)有效的。

3.2土體含水率剖面迭代計算

本文計算中，土體初始含水率為0.037，以3.1小節(jié)所計算的t=111 d土體表層含水率預(yù)測值0.137作為起始時刻邊界含水率，其余參數(shù)見表3；計算第一日不同時段t不同深度z處的含水率，表3中僅羅列t=480 min時刻，0～80 cm深度處含水率計算數(shù)據(jù)以供驗證，其余數(shù)據(jù)類似，不再贅述。

將表3計算成果繪制成一簇曲線，如圖10（a）所示。當(dāng)t=1 440 min（t=24 h）時，所對應(yīng)的曲線即為第一天入滲結(jié)束時，0～80 cm處土體含水率分布，該含水率剖面分布，即可作為第二日土體初始含水率，而第二日邊界含水率為3.1小節(jié)應(yīng)用相空間重構(gòu)預(yù)測的第二日表層含水率，將二者繼續(xù)代入式（7），即可迭代計算出第二日非飽和入滲曲線，如圖10（b）所示。

重復(fù)上述步驟，即可獲得未來連續(xù)100 d土體非飽和含水率分布曲線，這些曲線即構(gòu)成土體含水率剖面分布模型；限于篇幅，第3～100天類型重復(fù)圖片不再羅列。通過圖集，即可方便獲取未來100 d內(nèi)任意時刻、任意位置的土體含水率分布。

3.3討 論

通過以上小節(jié)內(nèi)容可知，應(yīng)用本文方法，可以方便預(yù)測未來時段，無資料地區(qū)表層含水率分布和土體內(nèi)部含水率分布，但仍有兩個問題需要進(jìn)一步討論說明：

（1） 在混沌預(yù)測表層含水率的方案中，本次研究僅采用基礎(chǔ)的“局部相點(diǎn)法”進(jìn)行預(yù)測，其余例如基于Lyapunov指數(shù)預(yù)測法、基于關(guān)聯(lián)維數(shù)D2預(yù)測法等其他混沌預(yù)測方法，由于篇幅所限本文并未采用。

（2） 在非飽和入滲計算中，為更清晰展現(xiàn)本研究思路，文章僅采用最為廣泛的一種Richards方程解析解進(jìn)行迭代，并未使用例如Hydrus-1D/2D/3D系列數(shù)值計算軟件進(jìn)模擬，但這同樣并不妨礙本文的主旨和結(jié)論。

4結(jié) 論

本文提出了一種無資料區(qū)土體表層及內(nèi)部含水率的預(yù)測方案，結(jié)合相空間重構(gòu)理論與非飽和入滲控制方程，計算生成了未來100 d內(nèi)土體0～80 cm深度的含水率分布，預(yù)測結(jié)果表明，在驗證期和預(yù)測期，含水率預(yù)測值與實(shí)測值最大相對誤差分別為2.4%和8.3%，滿足工程精度，以該預(yù)測結(jié)果計算的土體剖面含水率模型真實(shí)可靠。

按照本文步驟，可以復(fù)現(xiàn)任意研究區(qū)土體含水率潛在分布狀態(tài)，對于估計無資料地區(qū)土體含水率在長時間序列下的變化規(guī)律，本方案具有計算簡單、無野外作業(yè)、成本低等優(yōu)勢，可為相關(guān)工程提供一定借鑒。

參考文獻(xiàn)：

［1］李彬楠，樊貴盛，申麗霞.基于BP神經(jīng)網(wǎng)絡(luò)方法的黃土水分特征曲線預(yù)測模型比選［J］.中國農(nóng)村水利水電，2023（1）：171-175.

［2］詹明強(qiáng)，陳波，劉庭赫，等.基于變權(quán)重組合預(yù)測模型的混凝土壩變形預(yù)測研究［J］.水電能源科學(xué)，2022，40（9）：115-119.

［3］吳浩，張金接，姜龍，等.基于灰色理論的邊坡變形時效特性預(yù)測研究［J］.中國水利水電科學(xué)研究院學(xué)報，2014，12（3）：270-275.

［4］劉藝朦，魏江峰，趙景文.WRF 4.0模式陸氣耦合對不同物理參數(shù)化方案的敏感性初探［J］.大氣科學(xué)學(xué)報，2023：46（3）：431-440.

［5］陸培榕，羅紈，賈忠華，等.HYDRUS數(shù)值模擬方法在土壤水動力學(xué)實(shí)驗教學(xué)中的應(yīng)用［J］.實(shí)驗技術(shù)與管理，2022，39（11）：173-176.

［6］謝定義.非飽和土力學(xué)［M］.北京：高等教育出版社，2015.

［7］曹粵，包妮沙，周斌，等.基于實(shí)測光譜和國產(chǎn)高分五號高光譜衛(wèi)星的鐵尾礦表層含水率遙感反演方法研究［J］.光譜學(xué)與光譜分析，2023，43（4）：1225-1233.

［8］姚一飛，王爽，張珺銳，等.基于GF -1衛(wèi)星遙感的河套灌區(qū)土壤含水率反演模型研究［J］.農(nóng)業(yè)機(jī)械學(xué)報，2022，53（9）：239-251.

［9］杜瑞麒，陳俊英，張智韜，等.Sentinel-2多光譜衛(wèi)星遙感反演植被覆蓋下的土壤鹽分變化［J］.農(nóng)業(yè)工程學(xué)報，2021，37（17）：107-115.

［10］ 冀榮華，李鑫，張舒蕾，等.基于時延神經(jīng)網(wǎng)絡(luò)的多深度土壤含水率預(yù)測［J］.農(nóng)業(yè)工程學(xué)報，2017，33（增1）：132-136.

［11］李旭強(qiáng)，鄭秀清，薛靜，等.基于BAS-BPNN模型的季節(jié)性凍融期土壤含水率預(yù)測［J］.節(jié)水灌溉，2022（10）：66-70.

［12］黃潤生，黃浩.混沌及其應(yīng)用（第二版）［M］.武漢：武漢大學(xué)出版社，2005.

［13］王利，岳聰，舒寶，等.基于混沌時間序列的黃土滑坡變形預(yù)測方法及應(yīng)用［J］.地球科學(xué)與環(huán)境學(xué)報，2021，43（5）：917-925.

［14］周翠英，陳恒，朱鳳賢.邊坡演化的非線性時間序列多元混沌判別［J］.地球科學(xué)-中國地質(zhì)大學(xué)學(xué)報，2008，33（3）：393-398.

［15］YANO M，HOMMA N，SAKAI M，et al.Phase-space reconstruction from observed time series using Lyapunov spectrum analysis［C］∥Proceeding of SICE Annual Conference 2002.Osaka（Japan）：SICE，2002，721-726.

［16］PACKARD N H，CRUTCHFIELD J P，F(xiàn)ARMER J D，et，al.Geometry from a time series［J］.Physical Review Letters，1980，45（9）：712-716.

［17］羅哲賢，馬鏡嫻.混沌動力學(xué)及其在干旱預(yù)測中的應(yīng)用［J］.甘肅氣象，1997，15（3）：1-5.

［18］GRASSBERGER P，PROCACCIA I.Measuring the strangeness of strange attractors［J］.Physica D-Nonlinear Phenomena，1983，9（1-2）：189-208.

［19］張恭慶，林源渠.泛函分析講義［M］.北京：北京大學(xué)出版社，2001.

［20］高俊杰.混沌時間序列預(yù)測研究及應(yīng)用［D］.上海：上海交通大學(xué)，2013.

［21］朱悅璐，陳磊.基于最小作用原理的Richards方程變分解［J］.巖土力學(xué)，2022，43（1）：119-126.

［22］ZHU Y L，XIAO Y T.Research on a key question in the parlange solution［J］.KSCE Journal of Civil Engineering，2021，26：472-781.

［23］雷志棟，楊詩秀，謝森傳.土壤水動力學(xué)［M］.北京：清華大學(xué)出版社，1988.

（編輯：黃文晉）