中農(nóng)大土壤水動(dòng)力學(xué)課件07土壤水分動(dòng)態(tài)模擬模型及數(shù)值解法

第7章土壤水分動(dòng)態(tài)模擬模型及數(shù)值解法

(Introductiontothemodelsofsoilwaterdynamicsandthenumericalsimulationmethod)土壤水分動(dòng)態(tài)模擬模型土壤水分動(dòng)態(tài)模擬模型的必要性模型的分類典型模型介紹土壤水分運(yùn)動(dòng)的數(shù)值解法定解問(wèn)題及常用數(shù)值計(jì)算方法Richards方程的差分離散邊界條件處理土壤水分運(yùn)動(dòng)參數(shù)取值差分方程的求解土壤水分動(dòng)態(tài)模擬模型的必要性土壤水分研究在水文學(xué)、土壤學(xué)、農(nóng)田水利學(xué)、生態(tài)學(xué)、環(huán)境科學(xué)等學(xué)科的重要作用，眾多學(xué)者對(duì)土壤水分變化規(guī)律進(jìn)行深入研究，有必要建立土壤水分動(dòng)態(tài)模擬模型，以模擬土壤水分在不同自然及人為條件下的動(dòng)態(tài)變化及分布規(guī)律，進(jìn)而對(duì)農(nóng)田及天然植被耗水過(guò)程、污染物遷移過(guò)程進(jìn)行模擬分析

土壤水分動(dòng)態(tài)模擬模型的類型研究尺度：土體、農(nóng)田、灌區(qū)模型原理：系統(tǒng)模型（經(jīng)驗(yàn)?zāi)Ｐ停?、概念模型（水量平衡模型）和機(jī)理模型（水動(dòng)力學(xué)模型）是否考慮土壤水分動(dòng)態(tài)變化的不確定性：確定性模型和隨機(jī)性模型研究尺度：土體、農(nóng)田、灌區(qū)其相應(yīng)的空間尺度范圍分別為1～10m2、102～106m2、>106m2土體尺度，一般不考慮氣象、土壤、作物等因素的空間變異性，模型和研究方法較多樣而對(duì)于區(qū)域尺度，則需要考慮以上因素的空間變化，所采用的方法主要是水量平衡分析土壤水分動(dòng)態(tài)模擬的系統(tǒng)模型、概念模型和機(jī)理模型(1)系統(tǒng)模型一般不著重考慮土壤水分動(dòng)態(tài)變化的機(jī)理，而是根據(jù)土壤水分變化與其主要影響因素的關(guān)系或土壤水分序列自身的動(dòng)態(tài)變化規(guī)律，建立各種各樣的經(jīng)驗(yàn)性模擬模型，主要包括：(a)統(tǒng)計(jì)回歸模型：通過(guò)對(duì)土壤水分與降水等影響因素的關(guān)系進(jìn)行分析，建立反映其相互關(guān)系的回歸模型進(jìn)行水分動(dòng)態(tài)模擬與預(yù)報(bào)。(b)時(shí)間序列分析模型：利用時(shí)間序列分析方法研究土壤水分序列的動(dòng)態(tài)變化特性，將土壤水分的變化分解為趨勢(shì)項(xiàng)、周期項(xiàng)和隨機(jī)項(xiàng)，各項(xiàng)疊加即可對(duì)土壤水分動(dòng)態(tài)變化進(jìn)行模擬。(c)指數(shù)消退模型：通過(guò)對(duì)土壤水分消退規(guī)律的分析可以建立土壤水分動(dòng)態(tài)模擬的指數(shù)消退模型。(d)人工神經(jīng)網(wǎng)絡(luò)（ANN）模型：利用ANN結(jié)構(gòu)來(lái)反映氣象、作物等因素對(duì)土壤水分變化的影響，建立土壤水分動(dòng)態(tài)模擬的ANN模型。(2)概念模型水量平衡模型是一種概念性模型，根據(jù)一定時(shí)段內(nèi)土壤水分的輸入和輸出來(lái)模擬作物根系層土壤水分的動(dòng)態(tài)變化。

農(nóng)田水量平衡的基本方程：ΔW=W2-W1=P+I-ET-Q-R

水量平衡模型考慮了農(nóng)田水分收、支對(duì)土壤水分的影響，但對(duì)土壤水分運(yùn)動(dòng)考慮不夠深入。(3)水動(dòng)力學(xué)模型在Darcy定律和連續(xù)方程基礎(chǔ)上建立土壤水運(yùn)動(dòng)的基本方程——Richards方程，同時(shí)考慮土壤蒸發(fā)、作物蒸騰與根系吸水等界面過(guò)程，可建立農(nóng)田土壤水動(dòng)力學(xué)模型。根據(jù)考慮問(wèn)題的復(fù)雜程度，可以有單純的不考慮作物影響的SW模型，SPAC水分傳輸模型、SPAC水熱傳輸模型、GSPAC水熱傳輸模型、土壤凍融過(guò)程下的水熱傳輸模型以及土壤中的水熱鹽耦合運(yùn)移模型等等。與土壤水量平衡模型相比，偏微分方程形式的土壤水動(dòng)力學(xué)模型可以更細(xì)致地描述土壤水的運(yùn)動(dòng)與轉(zhuǎn)化，得到土壤水分動(dòng)態(tài)的時(shí)空分布規(guī)律。實(shí)際情況下，多采用數(shù)值方法來(lái)進(jìn)行求解土壤水分動(dòng)態(tài)模擬的確定性模型和隨機(jī)性模型（1）確定性模型影響農(nóng)田土壤水分動(dòng)態(tài)變化的因素（氣象、土壤等）在時(shí)間、空間上均有一定的隨機(jī)特性。如果不考慮上述因素的隨機(jī)性，則模型屬于確定性模型，目前常用的水量平衡模型和水動(dòng)力學(xué)模型多為確定性模型，適用于下墊面（土壤、作物等）均勻、氣象要素確定的情況。（2）隨機(jī)性模型考慮農(nóng)田土壤水分動(dòng)態(tài)變化的因素（氣象、土壤等）在時(shí)間、空間上的隨機(jī)特性，建立的模型。對(duì)于時(shí)域隨機(jī)性，首先用適當(dāng)?shù)碾S機(jī)過(guò)程模型來(lái)描述降水、騰發(fā)等的隨機(jī)變化特性。然后建立描述土壤水量平衡或水分運(yùn)動(dòng)的隨機(jī)微分（差分）方程模型或狀態(tài)空間模型，可以求解得到土壤水分動(dòng)態(tài)的概率分布。對(duì)于空間變異性有兩種處理方法，a.分布式模型方法，b.利用概率分布函數(shù)描述有關(guān)模型輸入及參數(shù)的空間變異特性。模型推導(dǎo)：土壤水分的減少是由蒸散發(fā)和深層滲漏造成的，而一般情況下下邊界水分通量比蒸散發(fā)量要小。在土壤水分脅迫條件下，蒸散發(fā)量與土壤含水量之間近似為線性關(guān)系。基于此，假設(shè)土壤水分消退階段水分消退率與貯水量W成正比，即

積分得到：W(t)＝W1

exp[－k(t－t1)]考慮降水及灌水的影響，土壤水分變化的遞推關(guān)系（以天為單位）可表示為Wt+1＝Wt·exp(－kΔt)＋P＋I(xiàn)典型模型1——土壤水分指數(shù)消退模型

土壤水分消退指數(shù)的確定方法模型的主要參數(shù)為土壤水分消退指數(shù)k。在無(wú)降水及灌水的時(shí)段內(nèi)，k可由土壤水分觀測(cè)資料推求：k＝ln(W1/W2)/(t2－t1)典型模型1——土壤水分指數(shù)消退模型

指數(shù)消退模型在冬小麥田間土壤水分動(dòng)態(tài)預(yù)報(bào)中的應(yīng)用利用北京永樂(lè)店試驗(yàn)站6試驗(yàn)小區(qū)冬小麥田間土壤水分觀測(cè)資料，計(jì)算出相應(yīng)的土壤水分消退指數(shù)(如圖)，擬合得到：土壤水分消退指數(shù)k動(dòng)態(tài)變化（冬小麥返青后）典型模型1——土壤水分指數(shù)消退模型

土壤水分動(dòng)態(tài)預(yù)報(bào)結(jié)果利用以上計(jì)算出不同時(shí)段的土壤水分消退指數(shù)k，并利用以上土壤水分動(dòng)態(tài)變化的遞推關(guān)系式即可進(jìn)行土壤水分動(dòng)態(tài)預(yù)報(bào)不同小區(qū)1m土層貯水量（W）實(shí)測(cè)值與預(yù)報(bào)值的比較

典型模型1——土壤水分指數(shù)消退模型

模型特點(diǎn)：比較簡(jiǎn)單、參數(shù)較少，使用方便；主要局限性：模型中土壤消退系數(shù)地域、時(shí)域性較強(qiáng)。因此在應(yīng)用這一模型時(shí)，可先根據(jù)其它年份的試驗(yàn)資料或部分測(cè)點(diǎn)的試驗(yàn)資料推求土壤水分消退指數(shù)，再用于相似條件下的土壤水分動(dòng)態(tài)預(yù)報(bào)典型模型1——土壤水分指數(shù)消退模型

典型模型2——土壤水分模擬的BP網(wǎng)絡(luò)模型人工神經(jīng)網(wǎng)絡(luò)（artificialneuralnetwork,ANN）理論—

是人工智能研究的一個(gè)重要領(lǐng)域。ANN是模擬人腦神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)與功能特征的一種非線性信息并行處理系統(tǒng)，是由大量的處理單元（人工神經(jīng)元）相互連接而成的網(wǎng)絡(luò)。ANN具有自學(xué)習(xí)、自適應(yīng)、自組織等特性，目前已在系統(tǒng)預(yù)測(cè)、系統(tǒng)優(yōu)化、模式識(shí)別、數(shù)據(jù)挖掘等領(lǐng)域得到了廣泛的應(yīng)用。

BP(back-propagation)網(wǎng)絡(luò)模型—

是一種基于誤差反向傳播的多層前饋型網(wǎng)絡(luò)，是目前應(yīng)用最為廣泛的ANN模型之一。BP-ANN結(jié)構(gòu)：三層網(wǎng)絡(luò)結(jié)構(gòu)：輸入層：Xn隱含層：Hq輸出層：Ym誤差：典型模型2——土壤水分模擬的BP網(wǎng)絡(luò)模型BP網(wǎng)絡(luò)的訓(xùn)練（模型率定）：

輸入信息的正向傳播過(guò)程和誤差的反向傳播過(guò)程

。步驟：（1）隨機(jī)產(chǎn)生初始權(quán)重和閾值

（2)逐層處理計(jì)算得到網(wǎng)絡(luò)的輸出（3）計(jì)算輸出與實(shí)測(cè)值間的誤差（4）根據(jù)誤差調(diào)整權(quán)重和閾值W2-W1=P+I–(ET+Q)ΔtET=KsKcET0

Q=f(W)輸入：時(shí)間t時(shí)段長(zhǎng)Δt階段降水灌水總量PI階段平均氣溫T時(shí)段初貯水量W1輸出：時(shí)段末貯水量W2典型模型2——土壤水分模擬的BP網(wǎng)絡(luò)模型輸入輸出變量的確定：預(yù)測(cè)結(jié)果統(tǒng)計(jì)：典型模型2——土壤水分模擬的BP網(wǎng)絡(luò)模型C1：F4：典型模型2——土壤水分模擬的BP網(wǎng)絡(luò)模型ΔW=W2-W1=P+I-ET-Q-R騰發(fā)量：ET=KsETm，ETm=KcET0

ET0：Penman-Monteith公式Kc：Ks：典型模型3——水量平衡模型ΔW=W2-W1=P+I-ET-Q-R根系層底部水分交換量：簡(jiǎn)化法（近似認(rèn)為Q=0）零通量面法經(jīng)驗(yàn)方法：典型模型3——水量平衡模型參數(shù)：根據(jù)實(shí)測(cè)資料率定

Wj，n（土壤供水系數(shù)）

Kcm，tm，c（作物系數(shù)）

a，d（下界面通量）典型模型3——水量平衡模型(a)A1小區(qū)（灌水4次）(b)B1小區(qū)（灌水3次）(c)C1小區(qū)（灌水2次）(d)D1小區(qū)（灌水2次）模擬結(jié)果：北京永樂(lè)店試驗(yàn)站1999年典型模型3——水量平衡模型永樂(lè)店試驗(yàn)站冬小麥返青后騰發(fā)量ET、土壤水利用量Su與降水灌溉總量的關(guān)系典型模型4——機(jī)理模型待續(xù)土壤水分運(yùn)動(dòng)的數(shù)值解法定解問(wèn)題：Richards方程邊界條件初始條件求解方法：解析法、半解析法：數(shù)值方法：解析法：求出解的函數(shù)表達(dá)式:精度高，但有很大局限性:簡(jiǎn)化、特殊情況。數(shù)值法：求用數(shù)值表示的有限個(gè)離散點(diǎn)和離散時(shí)段上的近似解:具有廣泛適用性；容易程序化，數(shù)值模擬軟件多；隨著計(jì)算機(jī)技術(shù)的提高和數(shù)值算法的改進(jìn)得以迅速發(fā)展。常用數(shù)值計(jì)算方法：全離散化方法：時(shí)間和空間進(jìn)行離散，利用微分方程推導(dǎo)出物理量在節(jié)點(diǎn)值的代數(shù)方程組，求解方程組得到微分方程在節(jié)點(diǎn)上的近似解有限差分法(FDM)：用差商替代導(dǎo)數(shù)控制容積法(CVM)：積分方程離散，守恒有限元(FEM)：借助基函數(shù)，直接求函數(shù)的近似解半離散方法：時(shí)間或空間離散有限解析法（FAM）線法(MOL)數(shù)值計(jì)算的收斂性與穩(wěn)定性收斂性：Δt、Δx→0時(shí)，數(shù)值方程解→微分方程的解穩(wěn)定性：計(jì)算誤差有界數(shù)值彌散問(wèn)題：對(duì)流占優(yōu)情況下容易出現(xiàn)數(shù)值彌散，利用特征線法等來(lái)克服有限元法的離散：有限元法變化形式較多，如：1）從權(quán)函數(shù)的選擇來(lái)說(shuō)，有配置法、矩量法、最小二乘法和伽遼金法；2）從計(jì)算單元網(wǎng)格的形狀來(lái)劃分，有三角形網(wǎng)格、四邊形網(wǎng)格和多邊形網(wǎng)格；3）從插值函數(shù)的精度來(lái)劃分，又分為線性插值函數(shù)和高次插值函數(shù)等。

fj有限元法－用簡(jiǎn)單的插值函數(shù)代替每個(gè)單元上未知函數(shù)分布，形成一系列代數(shù)方程組的方法

Galerkin法引入基函數(shù)ψi

并采用ψi作權(quán)函數(shù)，使得：

舉例—Galerkin法：對(duì)于：可用矩陣表示為：式中H、dH/dt分別為由h1，h2,…,hn和dh1/dt,dh2/dt,…dhn/dt組成的列矢量；D為導(dǎo)水矩陣，P為貯水矩陣。若采用隱式方法：通過(guò)求解大型稀疏矩陣方程，得到Hk+1有限差分法－利用Taylor級(jí)數(shù)用差商代替導(dǎo)數(shù)：前向差分：后向差分：中心差分：前差后差有限差分法的離散：i＝12。。。。。。。。Ij＝12

。

。

。

。

。

。

。Jfi,jΔx中心差一階導(dǎo)數(shù)的近似表達(dá)式：

若采用中心差的形式，二階導(dǎo)數(shù)如何表示？Richards方程的差分離散定解問(wèn)題：初始：地表瞬時(shí)濕潤(rùn)時(shí)間差分：顯式差分格式：α=0——不穩(wěn)定中心差分(Crank-Nicholson)格式：α=0.5—可能不穩(wěn)定隱式差分格式：α=1——穩(wěn)定參數(shù)取值：時(shí)段初(β=0)時(shí)段中(β=0.5)時(shí)段末(β=1)——效果較好，需要迭代邊界條件處理一類邊界：θ0=θs二類邊界（通量邊界）：直接差分半節(jié)點(diǎn)水量平衡虛擬節(jié)點(diǎn)土壤水分運(yùn)動(dòng)參數(shù)取值時(shí)間：顯式：用時(shí)段初含水率剖面計(jì)算參數(shù)預(yù)報(bào)校正：以時(shí)段初值求解方程，用新的含水率校正參數(shù)，再次求解迭代：半節(jié)點(diǎn)土壤水分運(yùn)動(dòng)參數(shù)：算術(shù)平均幾何平均調(diào)和平均…差分方程的求解對(duì)于非顯式差分，得到方程：代入邊界條件得到三對(duì)角方程組如下，用追趕法（Thomas算法）求解求解步驟求解土壤水分運(yùn)動(dòng)方程時(shí)，首先要給出有關(guān)土壤水分運(yùn)動(dòng)參數(shù)的表達(dá)式、初始條件與邊界條件，然后按一定的時(shí)間、空間步長(zhǎng)進(jìn)行離散。在每一離散時(shí)段內(nèi)，進(jìn)行以下迭代計(jì)算：(1)以時(shí)段初的土壤含水率作為時(shí)段末含水率的預(yù)報(bào)值。(2)計(jì)算相關(guān)的土壤水分運(yùn)動(dòng)參數(shù)，利用追趕法求解土壤水分運(yùn)動(dòng)的差分方程組（三對(duì)角方程組），得到時(shí)段末含水率。(3)