基于混沌理論的非飽和土含水率預測

朱悅璐 吳奇俞

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

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

中圖法分類號： TU43

文獻標志碼： A

DOI：10.16232/j.cnki.1001-4179.2024.04.028

0引 言

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

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

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

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

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）為基礎數(shù)據(jù)集，其空間分辨率精度為0.3°×0.3°經(jīng)緯度，數(shù)據(jù)來源于“國家青藏高原科學數(shù)據(jù)中心”（http：∥data.tpdc.ac.cn）。

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

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

2研究方法

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

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

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

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

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

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

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

2.2.1遲滯τ計算

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

2.2.2混沌特性判斷

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

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

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

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

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

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

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

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

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

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

3模型計算結(jié)果

3.1混沌方案預測表層土壤含水率

3.1.1遲滯計算

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

3.3討 論

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

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

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

4結(jié) 論

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

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

參考文獻：

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

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

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

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

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

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

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

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

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

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

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

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

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

［14］周翠英，陳恒，朱鳳賢.邊坡演化的非線性時間序列多元混沌判別［J］.地球科學-中國地質(zhì)大學學報，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］羅哲賢，馬鏡嫻.混沌動力學及其在干旱預測中的應用［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］.北京：北京大學出版社，2001.

［20］高俊杰.混沌時間序列預測研究及應用［D］.上海：上海交通大學，2013.

［21］朱悅璐，陳磊.基于最小作用原理的Richards方程變分解［J］.巖土力學，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］雷志棟，楊詩秀，謝森傳.土壤水動力學［M］.北京：清華大學出版社，1988.

（編輯：黃文晉）