宋曉龍，鐘德鈺，王光謙

（1.清華大學(xué) 水沙科學(xué)與水利水電工程國家重點實驗室，北京 100084；2.天津大學(xué) 水利工程仿真與安全國家重點實驗室，天津 300072）

1 研究背景

河流生態(tài)系統(tǒng)，是陸地與海洋聯(lián)系的紐帶，在生物圈的物質(zhì)循環(huán)中起著主要作用。隨著近代以來全球化的推進(jìn)，借助基礎(chǔ)河流理論和數(shù)字分析手段的發(fā)展，河流資源開發(fā)力度逐步加大，同時與河流保護(hù)的矛盾問題也日益凸顯。在制定河流保護(hù)和水資源開發(fā)計劃之前，應(yīng)該最大限度地掌握河流相關(guān)信息，以做到合理科學(xué)地規(guī)劃河流。在中國等發(fā)展中國家，河流動力學(xué)和河流管理始終是相關(guān)基礎(chǔ)研究領(lǐng)域的重要課題，對實現(xiàn)資源的可持續(xù)利用意義重大。不過研究目標(biāo)需要分階段逐步實現(xiàn)，河流系統(tǒng)的非線性和突變性等特征始終是擺在專家學(xué)者面前的一道難題。

在很長一段時間內(nèi)，線性模型在研究河灣演化規(guī)律中扮演重要的角色，揭示了一些河灣現(xiàn)象的主要特征［1-4］，例如:河道幾何特征如河寬和彎曲率等的局部波動可以引起河道中縱向沙壩和橫向上河道平面的大尺度變化；小尺度自由沙壩在小寬深比條件下傾向于向下游移動，而在大寬深比條件下有可能向上游移動；在小擾動下，河流可在長期的發(fā)展過程中達(dá)到準(zhǔn)穩(wěn)定狀態(tài)，等等。雖然各種研究層出不窮，但遠(yuǎn)遠(yuǎn)稱不上完美。事實上，在洪水影響下，河流通常呈現(xiàn)出季節(jié)性和不規(guī)則波動。在極端天氣增多增強(qiáng)的今天，突發(fā)強(qiáng)降雨等事件無疑增加了我們的疑慮:是否河流系統(tǒng)可以與時時變化的水沙輸入保持平衡，達(dá)到所謂的“動態(tài)平衡狀態(tài)”；用一個單獨的造床流量是否可以代表影響河道形態(tài)的各種非恒定因素［5］；是否可以另外認(rèn)為河道形態(tài)是一種歷史擾動的綜合反映，而不是同外界時時保持平衡的一種準(zhǔn)穩(wěn)定狀態(tài)［6］；文獻(xiàn)［7-9］中有一些對氣候變化下的水文響應(yīng)特征研究，通過自然實驗或衛(wèi)星遙感分析了極端降雨對河道的潛在影響，但至今缺少對隨機(jī)擾動下的河流特征關(guān)系的機(jī)制性研究。

近幾十年來，眾多物理學(xué)家轉(zhuǎn)向現(xiàn)代系統(tǒng)理論拓展本領(lǐng)域的研究。其方法與以往經(jīng)典科學(xué)的顯著不同之處是，把自己研究的對象不再看作某個單獨事物，而是看作一個整體系統(tǒng)，要求做整體綜合性的研究和把握。特別是，表示線性對象的輸入與輸出間關(guān)系的傳遞函數(shù)法獲得了廣泛應(yīng)用。我們關(guān)注的自然河流系統(tǒng)，具有高度的規(guī)律性和交互性，其非平衡態(tài)特征正獲得越來越多的關(guān)注，而簡單的確定性方程并不能滿足相關(guān)研究的要求。1980年代以來，越來越多的研究基于擴(kuò)展的傳遞函數(shù)進(jìn)行河流演化過程的追蹤和識別，典型案例包括:用傳遞函數(shù)噪聲模型預(yù)測產(chǎn)匯流過程［10］，用隨機(jī)微分方程代替確定性方程來數(shù)值模擬水質(zhì)量因子的動態(tài)概率演化過程［11］，研究河岸帶寬度和植被密度隨流量變化的隨機(jī)響應(yīng)［12］，研究土壤水平衡動態(tài)模型在隨機(jī)降雨下的隨機(jī)變化［13］，研究地表徑流中懸浮顆粒的隨機(jī)運動［14］等等。這些方法為河流特征研究和流域管理提供了新的思路:將水流輸入項視為隨機(jī)源，通過概率化的研究來有效處理氣候突變等外界影響下的河流響應(yīng)問題。

河相關(guān)系是一種描述河流形態(tài)要素特征關(guān)系的經(jīng)典工具，文獻(xiàn)［15-18］中有大量的關(guān)于其非線性特征的研究。然而問題是，大多研究只是將焦點聚焦在利用實測數(shù)據(jù)來求解其冪函數(shù)的系數(shù)和指數(shù)本身，而很少考慮過該函數(shù)關(guān)系在隨機(jī)事件的影響下可能發(fā)生怎樣的變化。因此，本文將利用隨機(jī)微分方程的方法，對河相關(guān)系的隨機(jī)擾動提出新的見解，推動相關(guān)基礎(chǔ)理論的發(fā)展，有效提高河流管理水平。

2 河相關(guān)系的隨機(jī)微分方程模型

2.1 理論模型在河床演變學(xué)中，平灘水位是指河流灘地（唇）高程齊平，河道發(fā)生漫灘的臨界水位。水位上升到平灘水位時，由于相應(yīng)水流的流速大，輸沙能力高，造床能力強(qiáng)，因此對應(yīng)的所謂“平灘流量”成為反映河道排洪輸沙能力和設(shè)計穩(wěn)定河道幾何形態(tài)的重要指標(biāo)。經(jīng)典研究普遍認(rèn)為平灘流量可以由一個具體數(shù)值代表，無論來水來沙是否發(fā)生趨勢性變化，都會隨著河床的沖淤演變而隨時間變化。如吳保生［19］基于河流自動調(diào)整的基本原理，根據(jù)河床受到擾動后調(diào)整速率與其當(dāng)前狀態(tài)和平衡狀態(tài)之間差值成正比的基本規(guī)律，建立了一種平灘流量的滯后響應(yīng)模型:

式中:Qt為t時刻的平灘流量；β為平灘流量的變化速率，是一個與水流能量大小和河床邊界可動性有關(guān)的參數(shù)；Qe為相對某一給定水沙條件下，河床調(diào)整達(dá)到相對平衡時的平灘流量，可由式（2）表示:

式中:Qf為汛期平均流量；ξf=CtfQf，稱作來沙系數(shù)；Ctf為汛期平均懸移質(zhì)含沙量；K、b、c分別為待定系數(shù)和指數(shù)。

在一定的來沙水沙條件下，河床將不斷調(diào)整其比降和斷面形態(tài)，力求挾沙能力與上游來水來沙相適應(yīng)。在此調(diào)整過程中河床所能達(dá)到的最適宜狀態(tài)，或者說平衡狀態(tài)，我們稱為河相關(guān)系。它是通過某一特征流量，如平灘流量式（1），將河床不同斷面的比降、河寬、水深、流速等特征聯(lián)系在一起，用以描述沖積河床橫斷面形態(tài)以及河床的總輪廓［20］。如馬卡維耶夫公式和Leopold［21］等都是從穩(wěn)定斷面形態(tài)開始，得到一組描述天然河流河床形態(tài)的經(jīng)驗關(guān)系式:

式中:St、Bt、Dt、Ut分別為平灘流量下t時刻某河段平均比降及入口斷面的河寬、平均水深、平均流速；下標(biāo)不同的4種α和m為待定系數(shù)和指數(shù)。

根據(jù)式（3）的冪函數(shù)形式，我們可以用一概化變量X來代表河相關(guān)系特征參量:

冪函數(shù)形式的河相關(guān)系滿足非線性系統(tǒng)進(jìn)行線性化處理的條件，可對式（4）求導(dǎo)得:

將式（1）、（2）帶入式（5），并轉(zhuǎn)換為微分形式，有:

如果取平灘流量滯后響應(yīng)模型的多步遞推模式［19］來表示Qt的通解:

式中:Q0為t=0時Qt的值，即初始平灘流量；Δt為時段長度；n為迭代時段數(shù)；i為時段編號。

根據(jù)微分方程式（6）得到一組描述河相關(guān)系的確定性的迭代解。

河流系統(tǒng)是一個自然結(jié)構(gòu)、生態(tài)環(huán)境和經(jīng)濟(jì)社會相互耦合的開放系統(tǒng)，由于水體的流動性，系統(tǒng)與外界不斷進(jìn)行物質(zhì)和能量的交換以及信息的傳遞，從而表現(xiàn)出明顯的耗散性和協(xié)同性［22］。河流系統(tǒng)不能完全消除外界影響而獨立存在，因此普通微分方程在真實反映河流演化過程的效果要大打折扣。考慮到河床形態(tài)因子受各種因素的影響，比如氣候或人為環(huán)境擾動，引起水沙條件、河床邊界條件等產(chǎn)生不確定性，因此可參照金融領(lǐng)域的Black-Shloes模型等經(jīng)典的隨機(jī)擴(kuò)散模型［23-25］，假設(shè)河相關(guān)系的動態(tài)行為具有隨機(jī)性，在穩(wěn)態(tài)模型式（6）的基礎(chǔ)上增加噪聲模型來對其進(jìn)行分析:

經(jīng)典隨機(jī)微分方程的噪聲項通常由高斯白噪聲表示。然而其缺點是高斯白噪聲只能表現(xiàn)連續(xù)性強(qiáng)、較穩(wěn)定的隨機(jī)擾動，而不能控制突發(fā)情況的發(fā)生。實際的白噪聲就像白色光包含所有可見光譜一樣，也包含各種類型的噪聲，依據(jù)概率分布形式不同，可主要分為高斯白噪聲和非高斯白噪聲。其中，常用泊松噪聲（也稱泊松跳）來模擬的非高斯噪聲項，起著重要的傳遞突發(fā)信息的作用。其他領(lǐng)域的隨機(jī)研究中已經(jīng)發(fā)現(xiàn)了泊松跳的重要價值，并進(jìn)行了大量的推廣應(yīng)用［26-28］。此外，傳統(tǒng)的線性思維模式可能并不符合現(xiàn)實情況:天然河流系統(tǒng)就像分子生物系統(tǒng)，并不是純粹的機(jī)械系統(tǒng)，而是一個生命系統(tǒng)，其可表現(xiàn)出某些自發(fā)性和特殊的復(fù)雜性。研究表明，真實的系統(tǒng)輸入具有長相關(guān)和自相似性特征，而分子布朗運動是描述這種有機(jī)隨機(jī)過程的典型工具。因此，我們有必要將分形理論也應(yīng)用進(jìn)隨機(jī)微分方程以更好地刻畫河相關(guān)系的演變特征?？紤]到環(huán)境變化可能給河相關(guān)系帶來較穩(wěn)定的擾動，或類似脈沖的隨機(jī)沖擊，而這些隨機(jī)過程不管是發(fā)生頻率、發(fā)生時間或發(fā)生強(qiáng)度都是隨機(jī)的，因此筆者使用三種流行的噪聲模型，包括高斯白噪聲、泊松噪聲和分?jǐn)?shù)白噪聲，進(jìn)行三種組合，評估隨機(jī)微分方程（8）描述河相關(guān)系隨機(jī)演變特征的適用性。詳情如下:

2.1.1 簡單的高斯白噪聲模型

式中:{Wt:0≤t≤T}為標(biāo)準(zhǔn)維納過程；m和σ分別為漂移項和擴(kuò)散項參數(shù)；X0為初值。

方程式（9）的解可表示為:

2.1.2 高斯白噪聲+泊松噪聲組合模型

式中:Vj為跳躍強(qiáng)度；Nj(λj)為兩種相互獨立的泊松過程；λj為跳躍事件發(fā)生頻率；j=u*，d*分別表示上跳和下跳；跳躍強(qiáng)度的對數(shù)分布規(guī)律(Y=ln(V))可由下面的雙指數(shù)分布函數(shù)表示:

式中:ηu＞1，

這里用一個混合方法來模擬跳躍發(fā)生的時間點(τ1，τ2，…) 。其中，時間差Ri+1=τi+1-τi可根據(jù)參數(shù)為1λ的指數(shù)分布律計算，而λ可根據(jù)歷史數(shù)據(jù)進(jìn)行估計［29］。Xt?從一個跳躍點到下一個跳躍點之間按幾何布朗運動方式演化。假設(shè)一時間點t位于兩跳躍點τ1，τi+1之間(τ1＜t＜τi+1) ，帶跳隨機(jī)微分方程（11）的離散解可表示為:

2.1.3 分?jǐn)?shù)白噪聲+泊松噪聲組合模型

該模型與模型2的不同之處在于，含有Hurst參數(shù)H的{WtH，t≥0}是一個分?jǐn)?shù)布朗運動（fBm），它是一個連續(xù)中心化高斯過程，有協(xié)方差函數(shù):

當(dāng)H=1/2時，fBm是標(biāo)準(zhǔn)布朗運動，模型3退化為模型2。

假設(shè)跳躍時間點為τ1，τ2，…，方程式（14）可以用如下改進(jìn)的歐拉方法進(jìn)行離散求解:

其中，為滿足“準(zhǔn)確與效率”的雙重標(biāo)準(zhǔn)，使用一種快速的一維fBm生成器［30］模擬分?jǐn)?shù)白噪聲WtH，而用Yin［31］的方法處理fBm的差值:

2.2 模型參數(shù)的非參數(shù)估計法基于短期序列離散測量值進(jìn)行參數(shù)估計是隨機(jī)微分方程相關(guān)研究中的一項重要工作。國內(nèi)外文獻(xiàn)中已有大量方法，例如最大似然估計法（MLE）、廣義矩估計法、模擬矩估計法和MCMC方法等，每一種都有自己的獨特特征。S?rensen［32］經(jīng)過詳細(xì)調(diào)查和比較，證明MLE估計法由于其連續(xù)性、漸進(jìn)正態(tài)性、漸進(jìn)有效性等特征，在弱正則條件下表現(xiàn)最好。并且，MLE在從擴(kuò)散項中分離跳躍項方面也具有優(yōu)越性。因此，針對本研究提出的隨機(jī)微分方程，下面將使用一種非參數(shù)估計方法，充分發(fā)揮MLE法的優(yōu)越性，并排除傳統(tǒng)參數(shù)估計法帶來的限制，特別是重點將放在轉(zhuǎn)移密度函數(shù)的有效構(gòu)建上。

假設(shè){xi，i=0，1，…，N} 為實測值序列；θ為待定參數(shù)集向量（包括確定項河道參數(shù)及噪聲項參數(shù)）；Pθ(ti，xi； (ti-1，xi-1)，θ)是從xi-1到xi的轉(zhuǎn)移密度函數(shù)；那么θ的極大似然函數(shù)就為下面根據(jù)如下步驟采用蒙特卡洛模擬方法對L()θ進(jìn)行估計:

（1）假設(shè)時間范圍{Ft:t∈[0 ，T]}。實測值時間間隔:考慮時間區(qū)間為[ ]ti-1，ti，把該區(qū)間插入M個長度為的子區(qū)間，那么以Xti-1為初值，利用Euler-Maruyama算法，就可以通過隨機(jī)微分方程的計算得到Xti的近似值。重復(fù)R次，可以得到R維向量

其中，核函數(shù)取標(biāo)準(zhǔn)正態(tài)形式，帶寬hpi取流行的Silverman［33］形式:

其中si表示的標(biāo)準(zhǔn)差。

對于帶跳隨機(jī)微分方程的情況，假設(shè)lr=ln(Xr)-ln(x)(r=1，2，…，R)，用m，n分別表iii-1ud示時間間隔Δ內(nèi)發(fā)生上跳和下跳情況的數(shù)目，可以用如下4種條件概率密度函數(shù)的泊松加權(quán)和來表示區(qū)間[ti-1，ti]內(nèi)的非條件概率密度函數(shù)，進(jìn)而得到轉(zhuǎn)移密度函數(shù):

（3）對每一個xi重復(fù)以上步驟，就可以得到極大似然函數(shù)。經(jīng)過對數(shù)變換，我們可以用Matlab中的fmincon函數(shù)求解如下負(fù)對數(shù)極大似然函數(shù)的極小值，從而得到需要的參數(shù)估計值:

3 實例分析

3.1 研究區(qū)域黃河下游高村至孫口河段（如圖1）經(jīng)常出現(xiàn)過流能力明顯小于其上下游河段的駝峰現(xiàn)象，且在高村附近河底高程抬升幅度最大［34］。筆者也加入眾多對此駝峰河段重點關(guān)注的學(xué)者［34-38］行列，在此選取1952年到2013年間黃河下游高村至孫口段的深泓線縱比降，以及高村站河流橫斷面形態(tài)要素（河寬、水深、流速）資料［39］作為原始數(shù)據(jù)，對本研究提出的隨機(jī)微分方程模型進(jìn)行非參數(shù)估計以及演變趨勢研究。

根據(jù)吳保生方法［19］，利用高村站歷年平灘流量、汛期平均流量、來沙系數(shù)等資料（表1），通過最小二乘法得到流量公式（1）中的系數(shù)和指數(shù)K、b、c和β分別為30.26、-0.47、0.42和0.346。根據(jù)該公式計算高村站平灘流量，并與實測值比較，如圖2所示，得出結(jié)論符合較好。

圖1 黃河下游平面示意圖

圖2 高村站的平灘流量實測值與計算值對比

表1 高村站歷年汛期平均流量、來沙系數(shù)、實測及計算平灘流量

3.2 模型參數(shù)估計對于未知參數(shù)集θ，采用2.2節(jié)介紹方法進(jìn)行參數(shù)估計。本例中，時間區(qū)間為［1952，1982］，初值取1952年河相關(guān)系各特征數(shù)據(jù)，時間間隔Δ=1，步長取在單次估值試算中，由于Matlab中的fmincon函數(shù)對初值依賴性較強(qiáng)，導(dǎo)致結(jié)果失真嚴(yán)重。經(jīng)過多次嘗試確定的解決方法是:先給模型一個經(jīng)驗初值進(jìn)行計算，之后每次計算前先計算已得結(jié)果的均值作為新的初值，經(jīng)過多次重復(fù)計算取全體估值的平均值為最終結(jié)果。當(dāng)求解帶跳模型的跳躍發(fā)生頻率λj(j=u，d)時，可根據(jù)其歷史數(shù)據(jù)計算［29］，例如:在離散時間區(qū)間［0，T］內(nèi)，將極端事件（X≥xu或X≤xd）發(fā)生次數(shù)定義為nj，則λj=njT。對于本例中的河相關(guān)系特征量，使用1000個樣本重復(fù)計算5000次，從而得到估計結(jié)果。此外為便于對比分析，應(yīng)用最小二乘法估計計算了確定性方程（6）的參數(shù)m。

3.3 模擬計算基于以上的模型參數(shù)估計結(jié)果，可以通過多次重復(fù)計算（5000次），近似得到3種模型條件下全時段的河相關(guān)系概率分布演化結(jié)果。其中特別是，在一給定跳躍強(qiáng)度條件下，通過最大化隨機(jī)分布均值與實測值的相關(guān)系數(shù)，從而可近似估計出最優(yōu)跳躍發(fā)生時間點。如圖3—圖6為不同條件下的河相關(guān)系概率演化結(jié)果與實測數(shù)據(jù)的平面對比關(guān)系圖。

從圖3—圖6可以看出，模擬的概率穩(wěn)定分布圖與實測值趨勢總體相同，在一定程度上證明了隨機(jī)微分方程描述河相關(guān)系動態(tài)演化趨勢的有效性。然而對于模型1來說，問題是，在1990年左右隨機(jī)平均解大幅度偏離真實值；似乎跟預(yù)想的一樣，高斯白噪聲驅(qū)動的隨機(jī)過程連續(xù)性有余，突變性不足，局部模擬值過大或過小，這無疑凸顯了模型對突發(fā)情況模擬能力的不足。實際動態(tài)河流系統(tǒng)可能需要更加復(fù)雜的噪聲模型來刻畫，現(xiàn)在的白噪聲還不足以滿足要求。就像有足夠的離心力才能使得汽車平衡保持在賽道，白噪聲隨機(jī)方程還需要進(jìn)一步加強(qiáng)。

圖3 3種模型條件下比降的概率分布演化與實測值對比

圖4 3種模型條件下河寬的概率分布演化與實測值對比

包含高斯白噪聲和對數(shù)正態(tài)分布泊松噪聲的隨機(jī)微分方程最早由Merton［40］用于研究股票期權(quán)波動?？紤]到氣候條件突變或人為擾動可能給河相關(guān)系帶來類似脈沖的隨機(jī)影響，而這些隨機(jī)過程不管是發(fā)生頻率、發(fā)生時間或發(fā)生強(qiáng)度都是隨機(jī)的，因此筆者考慮在高斯白噪聲模型基礎(chǔ)上，增加泊松跳來對河相關(guān)系進(jìn)行研究。在該部分，河相關(guān)系的演化過程由幾何布朗運動表示的連續(xù)性過程與泊松跳表示的非連續(xù)性過程組成?？梢钥闯?，跟圖3—圖6中各子圖（a）相比，各子圖（b）的計算結(jié)果有了大幅度提升，特別是，平面概率穩(wěn)定帶明顯從平緩變得更加精致化，隨機(jī)平均解也更加貼近真實解的變化。這證明帶有泊松跳的隨機(jī)微分方程能顯著提高模擬水平，是基礎(chǔ)隨機(jī)理論的一個重要擴(kuò)展。唯一的遺憾是圖3—圖6的（b）中局部突出的變量大大超出了自定義的有效區(qū)域解范圍。某些時刻，較大（或較小）的河相關(guān)系特征量會發(fā)生過大的突變。因此帶跳隨機(jī)方程還不完美。除了突發(fā)性，可能還需要增加考慮非線性。

分?jǐn)?shù)布朗運動的研究最早可追溯到1940年代。Mandelbrot and Aizenman［41］命名并大大推動了分形理論的發(fā)展。如今其已廣泛應(yīng)用于自然科學(xué)和工程科學(xué)的廣泛領(lǐng)域，例如水力學(xué)中的湍流研究。因此，筆者建立了分?jǐn)?shù)-帶跳隨機(jī)微分方程。整體來看，圖3—圖6中各子圖（c）比相應(yīng)子圖（b）表現(xiàn)更好。對于河寬和流速算例，是由于其Hurst參數(shù)（0.3490和0.3014）距離0.5較遠(yuǎn)，因此促成了顯著地時增擴(kuò)散現(xiàn)象。換句話說，河寬和流速過程時間序列是顯著負(fù)相關(guān)的，在相鄰時間點之間長期發(fā)生高低值轉(zhuǎn)換現(xiàn)象。同時，對于比降和水深，該模型（3）計算結(jié)果與帶跳模型（2）計算結(jié)果十分接近，均成功預(yù)測出了局部跳躍點，但也漏掉了一些。還有令人遺憾的是，仍然有大量散點遠(yuǎn)離概率圖中的核心有效區(qū)間。筆者懷疑可能是由于真實的泊松項要遠(yuǎn)比我們這里的簡單概化方法更復(fù)雜，其跳躍過大或過小都不合適?？赡苡行У母倪M(jìn)手段包括在跳躍強(qiáng)度或跳躍頻率中增加縮放因子，或者將本文隨機(jī)模型中的互不相關(guān)噪聲改為有相關(guān)關(guān)系的噪聲。

4 模型討論與應(yīng)用

4.1 模擬比較從圖3—圖6可以看出，應(yīng)用隨機(jī)微分方程，傳統(tǒng)的確定性的河相關(guān)系解進(jìn)化成為一種概率穩(wěn)定帶，從而為河流控制和河流管理提供了一種重要的新思路。顯然，如果外輪廓令人滿意，模型好壞評判準(zhǔn)則主要是與有效區(qū)域解，或者，隨機(jī)平均解與測量值相關(guān)度（如，相對誤差-δ；Nash-Sutcliff系數(shù)-NSE）等有關(guān)。也就是說，假如我們定義離散值在概率分布平面演化圖中的發(fā)生概率為有效概率，內(nèi)限值為有效穩(wěn)定寬度；那么可以說，有效概率越大，有效穩(wěn)定寬度越小，隨機(jī)平均解與測量值的相對誤差越小、Nash-Sutcliff系數(shù)越大，模型整體評價越高。所以，考慮到以上因子對模型調(diào)查和比較的重要意義，在圖7—圖8中基于上文的概率分布演化結(jié)果，給出相對有效穩(wěn)定寬度、相對誤差、Nash-Sutcliff系數(shù)的演化和比較結(jié)果圖。結(jié)合圖3—圖6，接下來我們對各河相關(guān)系特征量給出針對性的模型比較和最優(yōu)建議。

圖5 3種模型條件下水深的概率分布演化與實測值對比

圖6 3種模型條件下流速的概率分布演化與實測值對比

圖7 基于隨機(jī)微分方程計算的河相關(guān)系有效概率穩(wěn)定寬度時變演化

圖8 基于隨機(jī)微分方程計算的河相關(guān)系隨機(jī)平均解與實測值比較

（1）比降:圖3—圖6表明泊松跳能夠大幅度提高系統(tǒng)動態(tài)演化的模擬效果。簡單高斯白噪聲模型理應(yīng)被拋棄。同時根據(jù)圖7—圖8可以判斷，分?jǐn)?shù)-泊松模型由于較小的有效穩(wěn)定寬度、較大的有效概率和NSE系數(shù)，要勝過簡單帶跳隨機(jī)模型。

（2）河寬:與比降算例類似，根據(jù)圖3—圖6，我們可以直觀發(fā)現(xiàn)后兩種泊松模型勝過簡單白噪聲模型的地方。然后，再分析圖7—圖8可以發(fā)現(xiàn)，簡單帶跳隨機(jī)模型并沒有大幅度落后分?jǐn)?shù)-泊松模型，其有效概率穩(wěn)定寬度在后半段大部分時間，和NSE系數(shù)在全程的表現(xiàn)，均領(lǐng)先對手。分?jǐn)?shù)-泊松模型的優(yōu)勢是相對誤差-δ較小、隨機(jī)均解的波動情況模擬較好。綜合來看，兩種泊松模型應(yīng)該結(jié)合在一起作為河流管理的參考資源。

（3）水深:與以上兩算例不同，圖3—圖6和圖7—圖8均顯示出高斯白噪聲+泊松噪聲模型在模擬水深動態(tài)演進(jìn)過程中的優(yōu)勢，換句話說，水深主要受到突發(fā)隨機(jī)事件的線性影響。

（4）流速:從圖3—圖6和圖7—圖8可以看出，分?jǐn)?shù)-泊松模型要優(yōu)越于高斯-泊松模型，特別是在隨機(jī)均解波動情況的模擬、相對誤差和NSE系數(shù)上的表現(xiàn)。顯然分?jǐn)?shù)-泊松模型更令人信服。

總結(jié)來看，分?jǐn)?shù)-泊松模型較適合模擬比降和流速的演進(jìn)，兩種泊松模型最好結(jié)合起來共同管控河寬，而簡單高斯-泊松模型對于探索水深特征已經(jīng)足夠。

4.2 系統(tǒng)應(yīng)用系統(tǒng)穩(wěn)定性的評估始終都是理解河流的必要工作。有一些很好的方法已經(jīng)用于從整體上描述河流系統(tǒng)，包括能夠區(qū)分河型的河床穩(wěn)定指數(shù)（Zw）［42］（游蕩型:Zw＜5； 分汊型:5＜Zw＜15；彎曲型:Zw＞15），代表河岸侵蝕強(qiáng)度的寬深比系數(shù)（ζ），以及表示水流強(qiáng)度的河流功率（Ω）（式（24））。接下來，我們將使用分?jǐn)?shù)-泊松隨機(jī)微分方程對沿高村站向下游的Zw、ζ、Ω進(jìn)行概率分布演化的統(tǒng)計計算，結(jié)果如圖9所示。

圖9 基于分?jǐn)?shù)-泊松隨機(jī)微分方程計算的河床穩(wěn)定指數(shù)、寬深比系數(shù)和河流功率隨機(jī)演化圖

從圖9（a）中可以看出，河床穩(wěn)定指數(shù)從1986年左右開始大幅下滑，在此期間，其隨機(jī)均值幾乎跌破Zw=5的游蕩閾值線，直到接近2000年才又逐漸增大并回到分汊區(qū)間（5＜Zw＜15）。與此同時，高村站的水流強(qiáng)度經(jīng)歷緩慢增加過程（圖9（b）），且在2000年之后經(jīng)受的沖深作用大于展寬作用（圖9（c））。其原因，經(jīng)過查找文獻(xiàn)，我們可以了解到，是由于上游小浪底水庫的運行化解了下游長期少水多沙的負(fù)擔(dān)。這也與相關(guān)研究［35］做出的結(jié)論，認(rèn)為駝峰現(xiàn)象與水流流量大小成反比相吻合。整個圖9顯示出，河流穩(wěn)定性的低谷已經(jīng)在1999年被挽救，然而2008年左右出現(xiàn)的危機(jī)依然值得關(guān)注。并且，概率分布穩(wěn)定寬度與隨機(jī)均值間的正相關(guān)關(guān)系，也提醒我們未來需要提高對黃河下游河道演變的警戒水平。此外，考慮到測量解并不能完全反應(yīng)真實變化過程，例如2000年左右的寬深比系數(shù)（ζ），我們也因此可以看出隨機(jī)微分方程的優(yōu)越性。換句話說，對河相關(guān)系特征量個別單獨的測量值并不能幫助我們很好地了解河流系統(tǒng)性的變化，而在僅有少量汛期流量和懸移質(zhì)沙量的有限條件下，我們可以通過本研究的隨機(jī)方法表示的，例如寬深比系數(shù)（ζ），來進(jìn)行河流的系統(tǒng)性評估。通過Zw、ζ、Ω的隨機(jī)均值與平灘流量相關(guān)系數(shù)的比較（表2），我們可以得知流量主導(dǎo)的黃河下游適合用隨機(jī)狀態(tài)空間上的河流功率（Ω）來系統(tǒng)表征。近年我們收集的最新資料顯示，受水利與水土保持工程攔截、上游水庫調(diào)節(jié)、沿河用水增加及氣候變化等多重因素影響，黃河干流水沙形勢發(fā)生了明顯變化，表現(xiàn)為來水來沙量的持續(xù)減小。在新的水沙形勢和河道狀況下，河道自身過流能力顯著增強(qiáng)，同時中小流量出現(xiàn)的機(jī)率增加。這些都降低了洪水漫灘的機(jī)會。而長期相關(guān)的影響還需要進(jìn)行系統(tǒng)性地測量、識別和評估。本研究具有重大的理論和現(xiàn)實意義。

表2 Zw、ζ、Ω的隨機(jī)平均解與平灘流量的相關(guān)系數(shù)

5 結(jié)論

河流是一個由流量驅(qū)動，以河相關(guān)系為重要表現(xiàn)形式的動態(tài)系統(tǒng)。本文通過在確定性方程基礎(chǔ)上加入隨機(jī)噪聲模型（包括簡單高斯白噪聲模型、高斯-泊松模型和分?jǐn)?shù)-泊松模型），成功地建立并應(yīng)用隨機(jī)微分方程研究了黃河下游高村-孫口段的河相關(guān)系隨時間演化的概率分布動態(tài)演化過程，為河流基礎(chǔ)理論做出了重要貢獻(xiàn)，同時為相關(guān)河流決策提供了一個新思路。本方法可以幫助我們在有限的條件下系統(tǒng)性地掌握河流變化特征。通過模型比較，本文認(rèn)定同時考慮非線性和突發(fā)性擾動特征的分?jǐn)?shù)-泊松擴(kuò)散模型是適宜研究河相關(guān)系隨機(jī)演化特征的較優(yōu)模型。未來將把該方法進(jìn)行更多的推廣應(yīng)用和進(jìn)一步的完善，包括避免使用簡化的噪聲模型，而使用互相相關(guān)的更真實的復(fù)雜噪聲模型等。