拋石輸運(yùn)行為的連續(xù)時(shí)間隨機(jī)行走模型應(yīng)用

王茂枚，李志鵬，趙 鋼*，魯程鵬，孫洪廣

(1.江蘇省水利科學(xué)研究院，南京 210017； 2.河海大學(xué) 水文水資源學(xué)院，南京 210098； 3.河海大學(xué) 力學(xué)與材料學(xué)院，南京 211100)

拋石護(hù)岸是河道治理的傳統(tǒng)方法之一，作為一種在國內(nèi)外廣泛使用的防護(hù)技術(shù)，其工程造價(jià)小，實(shí)施方法簡便，被廣泛的應(yīng)用于各類河道治理工程中[1]。然而，不同地域的各類河流自然條件差異巨大，在工程上廣泛使用的經(jīng)驗(yàn)公式往往與實(shí)際觀測(cè)產(chǎn)生較大的出入[2]，因此研究拋石的輸運(yùn)行為在河道治理、工程應(yīng)用上成為了迫切需要解決的問題。目前國內(nèi)外對(duì)拋石輸運(yùn)行為的研究方法主要以實(shí)驗(yàn)和傳統(tǒng)水沙數(shù)值模型為主[3-4]。然而研究發(fā)現(xiàn)，河流中泥沙及拋石的輸運(yùn)行為往往是隨機(jī)反常的，即使同樣的水沙條件，河床中礫石的啟動(dòng)行為也是一個(gè)概率事件[5]。在拋石的輸運(yùn)過程中，由于受到河床及流場(chǎng)演變的影響，使得拋石的輸運(yùn)行為往往是非馬爾科夫性質(zhì)的，因此基于菲克定律的經(jīng)典方程在描述該類輸運(yùn)行為時(shí)往往達(dá)不到預(yù)期效果。這類具有非馬爾可夫性質(zhì)的輸運(yùn)行為被定義為“反常輸運(yùn)”[6]。

近年來，“反常輸運(yùn)”的研究受到國內(nèi)外的廣泛關(guān)注。“反常輸運(yùn)”在金融、醫(yī)學(xué)、半導(dǎo)體以及地質(zhì)領(lǐng)域上都得到了成功的應(yīng)用[7-8]。目前該方向研究以兩類方法為主，一類是以確定性方程為基礎(chǔ)的，例如：“分?jǐn)?shù)階導(dǎo)數(shù)”、“Hausdorff導(dǎo)數(shù)”[9]；而另一類方法以統(tǒng)計(jì)理論為基礎(chǔ)，基于事件發(fā)生的統(tǒng)計(jì)性質(zhì)得到輸運(yùn)行為的隨機(jī)運(yùn)動(dòng)方程[10]。本文所應(yīng)用的連續(xù)時(shí)間隨機(jī)行走理論就是基于隨機(jī)運(yùn)動(dòng)方程描述拋石的輸運(yùn)行為的。

連續(xù)時(shí)間隨機(jī)行走模型是研究“反常輸運(yùn)”行為最直觀的工具，在連續(xù)時(shí)間隨機(jī)行走理論框架中，一次成功的行走被分解為兩部分：(1)隨機(jī)的空間跳躍間隔；(2)兩次成功跳躍之間的等待時(shí)間；兩者分布由各自的概率分布密度支配[6，11]。本文應(yīng)用連續(xù)時(shí)間隨機(jī)行走模型，通過長尾的等待時(shí)間及跳躍步長分布模擬拋石在河道中的“反常輸運(yùn)”行為，并與實(shí)驗(yàn)室尺度下的實(shí)驗(yàn)觀測(cè)數(shù)據(jù)對(duì)比進(jìn)一步驗(yàn)證了模型對(duì)于拋石輸運(yùn)行為預(yù)測(cè)的有效性。

1 連續(xù)時(shí)間隨機(jī)行走理論

連續(xù)時(shí)間隨機(jī)行走模型中粒子的跳躍過程由空間跳躍距離和兩次成功跳躍之間的等待時(shí)間組成，兩者由共同的聯(lián)合概率密度分布決定[6，8]。在實(shí)際研究中發(fā)現(xiàn)跳躍步長與等待時(shí)間往往是不相關(guān)的，即跳躍步長和等待時(shí)間分別由各自的分布密度決定。研究中λ(x)，ω(t)表示跳躍步長與等待時(shí)間的概率密度函數(shù)，則聯(lián)合概率密度函數(shù)可表示為

ψ(x,t)=λ(x)·ω(t)

(1)

連續(xù)時(shí)間隨機(jī)行走模型的廣義主方程可描述為

(2)

式中：η(x,t)表示粒子t時(shí)刻到達(dá)x處的概率密度，粒子的初始分布為P0(x)。

粒子的概率密度函數(shù)P(x,t)滿足方程

(3)

將廣義主方程(2)代入(3)式,最終得到連續(xù)時(shí)間隨機(jī)行走的概率方程為

(4)

在連續(xù)時(shí)間隨機(jī)行走框架中，空間和時(shí)間都被定義為連續(xù)變量，因此對(duì)方程(4)進(jìn)行Fourier-Laplace變換，得

(5)

傅里葉—拉普拉斯域內(nèi)的方程(5)被稱為Montroll-Weiss方程，該方程中通過跳躍步長和等待時(shí)間的概率密度函數(shù)的性質(zhì)可以區(qū)分不同輸運(yùn)行為，當(dāng)跳躍步長的概率密度分布二階矩發(fā)散而等待時(shí)間的概率密度分布一階矩收斂時(shí)，輸運(yùn)行為定義為超擴(kuò)散輸運(yùn)；跳躍步長的概率密度分布二階矩收斂而等待時(shí)間概率密度分布一階矩發(fā)散時(shí)，輸運(yùn)行為定義為次擴(kuò)散輸運(yùn)；當(dāng)兩者皆發(fā)散，即同時(shí)具有超擴(kuò)散輸運(yùn)行為和次擴(kuò)散輸運(yùn)行為，兩者皆收斂則為正常輸運(yùn)行為。

2 蒙特卡羅模擬方法

蒙特卡羅方法，又稱隨機(jī)模擬方法，利用重復(fù)隨機(jī)的抽樣方法，求得統(tǒng)計(jì)量的數(shù)字特征。當(dāng)重復(fù)次數(shù)足夠多時(shí)，蒙特卡羅方法得到的數(shù)值解逼近真實(shí)解。本研究中由于方程(4)在轉(zhuǎn)換到時(shí)域上具有較大的困難，因此采用蒙特卡羅方法求得其數(shù)值結(jié)果，該方法是典型的拉格朗日法，可以追蹤粒子運(yùn)動(dòng)過程中任意時(shí)刻的位置，完整呈現(xiàn)粒子運(yùn)動(dòng)軌跡。

經(jīng)典的對(duì)流擴(kuò)散方程(ADE)表征流體的質(zhì)量傳輸規(guī)律，并能描述滿足菲克定律的流體中粒子輸運(yùn)行為[9]。但實(shí)際觀測(cè)中發(fā)現(xiàn)很多情況下菲克定律會(huì)被打破，例如含水層的異質(zhì)性、河床結(jié)構(gòu)的非均勻性都會(huì)導(dǎo)致ADE方程無法準(zhǔn)確描述該類環(huán)境中粒子的輸運(yùn)行為。這類打破經(jīng)典菲克定律的輸運(yùn)行為被稱為“反常輸運(yùn)”。連續(xù)時(shí)間隨機(jī)行走理論描述輸運(yùn)行為時(shí)，當(dāng)概率密度分布滿足中心極限定理時(shí)，即跳躍步長分布為正態(tài)分布，等待時(shí)間分布為指數(shù)分布時(shí)，方程(5)將回歸到ADE方程。在本研究中，為了反映河床結(jié)構(gòu)的異質(zhì)性導(dǎo)致的拋石“反常輸運(yùn)”行為，采用 穩(wěn)定分布作為跳躍步長分布，Mittag-Leffler分布作為等待時(shí)間分布。

本研究中我們將拋石的運(yùn)動(dòng)行為劃分為沿水流方向和垂直于水流方向的運(yùn)動(dòng)。因此，運(yùn)動(dòng)過程的時(shí)間計(jì)量由下式確定

tn+1=tn+τβ
t0=0

(6)

式中：τβ為滿足Mittag-Leffler分布的隨機(jī)數(shù)。

沿水流方向的位移由流場(chǎng)產(chǎn)生的確定位移和粒子相互碰撞的噪聲項(xiàng)確定，表述為

(7)

式中：v表示粒子由對(duì)流產(chǎn)生的速度，可由實(shí)驗(yàn)獲得。

垂直于水流方向的位移僅僅由粒子的相互碰撞產(chǎn)生，由下式確定

(8)

式中：ξα1和ξα2是滿足α穩(wěn)定分布的隨機(jī)數(shù)。研究中隨機(jī)數(shù)的生成分別通過以下公式生成

其中滿足跳躍步長的α穩(wěn)定分布隨機(jī)數(shù)采用式(9)生成[12]

(9)

滿足Mittag-Leffler分布的等待時(shí)間隨機(jī)數(shù)采用式(10)生成[13]

(10)

式中：u1,u2∈(0，1)為均勻分布的隨機(jī)數(shù)，γx，γt是隨機(jī)數(shù)產(chǎn)生的尺度系數(shù)。

圖1驗(yàn)證了產(chǎn)生隨機(jī)數(shù)的概率密度分布以及不同參數(shù)下α穩(wěn)定分布和Mittag-Leffler分布的變化特征。通過圖1-a與1-c隨機(jī)數(shù)概率密度與實(shí)際分布對(duì)比，隨機(jī)數(shù)的生成很好的滿足了所需概率密度分布。進(jìn)一步分析不同參數(shù)下分布的概率密度特性，發(fā)現(xiàn)α穩(wěn)定分布在參數(shù)α=2時(shí)回歸到正態(tài)分布，并隨著α的減小分布逐漸變“寬”，說明隨著參數(shù)α的減小，粒子的跳躍步長愈加發(fā)散，超擴(kuò)散輸運(yùn)行為越強(qiáng)烈；Mittag-Leffler分布具有和α穩(wěn)定分布相同的性質(zhì)，當(dāng)參數(shù)β=1時(shí)，Mittag-Leffler分布回歸到指數(shù)分布，當(dāng)參數(shù)β變小，等待時(shí)間的分布更加發(fā)散，表明其次擴(kuò)散行為越強(qiáng)烈。

1-a 穩(wěn)定分布隨機(jī)數(shù)概率密度驗(yàn)證1-b 不同參數(shù)下穩(wěn)定分布的狀態(tài)1-c Mittag-Leffler 分布隨機(jī)數(shù)概率密度驗(yàn)證1-d 不同參數(shù)下Mittag-Leffler分布的狀態(tài)圖1 跳躍步長與等待時(shí)間隨機(jī)數(shù)概率密度驗(yàn)證及不同參數(shù)下概率密度函數(shù)的差異Fig.1 Verification of the probability of random numbers for jump steps and waiting times, and the difference in probability density function under various parameters

圖2 Martin[2012] 實(shí)驗(yàn)結(jié)果及連續(xù)時(shí)間隨機(jī)行走模型和ADE模型預(yù)測(cè)效果比較Fig.2 Martin [2012] experimental results and comparison of predictive effects between continuous time random walk model and ADE

3 實(shí)驗(yàn)結(jié)果預(yù)測(cè)

研究中為了驗(yàn)證連續(xù)時(shí)間隨機(jī)行走模型在實(shí)際工況中的預(yù)測(cè)效果，本節(jié)對(duì)兩組不同源條件下實(shí)驗(yàn)室尺度水槽實(shí)驗(yàn)結(jié)果進(jìn)行預(yù)測(cè)。表1給出了兩組實(shí)驗(yàn)數(shù)據(jù)下連續(xù)時(shí)間隨機(jī)行走模型參數(shù)的匹配情況；比較發(fā)現(xiàn)，在實(shí)際運(yùn)動(dòng)中拋石的輸運(yùn)行為往往是次擴(kuò)散的(均方位移指數(shù)小于1)，但運(yùn)動(dòng)過程中同時(shí)表現(xiàn)出超擴(kuò)散行為(1<α<2)和次擴(kuò)散行為( 0<β<1)。

3.1 瞬時(shí)源條件下實(shí)驗(yàn)驗(yàn)證及分析

圖3 連續(xù)時(shí)間隨機(jī)行走模型預(yù)測(cè)Martin[2012] 試驗(yàn)的時(shí)-空分布及均方位移Fig.3 Time-Space evolution and MSD using continuous time random walk model

此處，我們選擇Martin于2012年在瞬時(shí)源條件下進(jìn)行的實(shí)驗(yàn)[14]。實(shí)驗(yàn)結(jié)果表明，礫石在河床上的分布并不為對(duì)稱形式的正態(tài)分布，反而表現(xiàn)出一定程度的空間右拖尾現(xiàn)象，這種現(xiàn)象可能是由于空間的異質(zhì)性導(dǎo)致了礫石在河床上的運(yùn)動(dòng)出現(xiàn)了反常輸運(yùn)的現(xiàn)象。圖2比較了連續(xù)時(shí)間隨機(jī)行走模型與經(jīng)典ADE方程在實(shí)驗(yàn)數(shù)據(jù)預(yù)測(cè)上的優(yōu)劣，通過比較發(fā)現(xiàn)，連續(xù)時(shí)間隨機(jī)行走模型很好的捕捉了礫石的拖尾現(xiàn)象，而ADE方程不能捕捉具有這類拖尾現(xiàn)象的行為。同時(shí)圖3粒子的時(shí)-空分布及均方位移指數(shù)說明該實(shí)驗(yàn)中礫石在河床上的運(yùn)動(dòng)行為是次擴(kuò)散的，即一部分礫石被水流及河床結(jié)構(gòu)加速，但河床結(jié)構(gòu)整體上對(duì)礫石運(yùn)動(dòng)有阻隔作用，故礫石的整體的運(yùn)動(dòng)效應(yīng)表現(xiàn)為次擴(kuò)散輸運(yùn)狀態(tài)。

3.2 連續(xù)源實(shí)驗(yàn)驗(yàn)證及分析

為了驗(yàn)證連續(xù)時(shí)間隨機(jī)行走理論在連續(xù)源條件下以及短時(shí)對(duì)長時(shí)的預(yù)測(cè)效果，我們?cè)谶@里選取了Chang and Yen于2002年進(jìn)行的連續(xù)源實(shí)驗(yàn)[15]。圖4對(duì)比了連續(xù)時(shí)間隨機(jī)行走模型和經(jīng)典ADE方程的預(yù)測(cè)結(jié)果，此處我們采用初始時(shí)刻t=7.5 s作為模型參數(shù)的濾定參考項(xiàng)，之后各個(gè)時(shí)刻均采用初始時(shí)刻的參數(shù)。相比于經(jīng)典ADE方程結(jié)果，連續(xù)時(shí)間隨機(jī)行走模型在短時(shí)間濾定參數(shù)進(jìn)而預(yù)測(cè)長時(shí)結(jié)果時(shí)具有更好的效果。進(jìn)一步分析發(fā)現(xiàn)，ADE方程在短時(shí)預(yù)測(cè)長時(shí)的效果上，隨著時(shí)間的增加其預(yù)測(cè)誤差會(huì)逐漸增大，圖4-a明顯表現(xiàn)出ADE方程在描述空間尾部效應(yīng)時(shí)，其曲線尾部的衰減以指數(shù)形式衰減，而連續(xù)時(shí)間隨機(jī)行走模型尾部以冪律衰減，這是兩類模型在后期預(yù)測(cè)效果產(chǎn)生巨大差異的重要原因。

4-a t=7.5 s4-b t=22.5 s4-c t=37.5 s

4-d t=52.5 s4-e t=67.5 s4-f t=82.5 s圖4 Chang and Yen[2002] 實(shí)驗(yàn)結(jié)果及連續(xù)時(shí)間隨機(jī)行走模型和ADE模型預(yù)測(cè)效果比較Fig.4 Chang and Yen[2002] experimental results and comparison of predictive effects between continuous time random walk model and ADE model

表1 連續(xù)時(shí)間隨機(jī)行走模型參數(shù)及實(shí)驗(yàn)數(shù)據(jù)對(duì)比Tab.1 Continuous time random walk model parameters and experimental data

4 模擬結(jié)果分析

為了表征不同輸運(yùn)行為在拋石輸運(yùn)過程中的影響，本節(jié)對(duì)比了不同位置、不同拋石方式下，拋石的運(yùn)動(dòng)行為。分析中參數(shù)采用無量綱化處理,粒子的平均速度統(tǒng)一為v=0.5。

4.1 上游中心點(diǎn)源釋放

粒子于上游中心釋放，對(duì)比了正常輸運(yùn)行為、次擴(kuò)散輸運(yùn)行為以及超擴(kuò)散輸運(yùn)行為下拋石的空間分布。圖5比較了次擴(kuò)散-正常擴(kuò)散-超擴(kuò)散的輸運(yùn)行為差異，隨著超擴(kuò)散行為的增強(qiáng)，拋石的空間離散性逐漸增大。同時(shí)由于次擴(kuò)散產(chǎn)生的長程等待時(shí)間，在次擴(kuò)散狀態(tài)下相當(dāng)一部分拋石停留在離點(diǎn)源位置較近的區(qū)域。圖中X方向是粒子具有對(duì)流速度的方向，箭頭表示水流方向，Y方向是粒子僅具有隨機(jī)彌散位移的方向。

5-a 正常輸運(yùn)行為下粒子不同時(shí)刻的空間分布5-b 次擴(kuò)散輸運(yùn)行為下粒子不同時(shí)刻的空間分布5-c 超擴(kuò)散輸運(yùn)行為下粒子不同時(shí)刻的空間分布圖5 沿中心釋放拋石Fig.5 Release of riprap along the center

4.2 上游岸邊點(diǎn)源釋放

將拋石的釋放位置放置于岸邊，對(duì)比正常擴(kuò)散-次擴(kuò)散-超擴(kuò)散三種輸運(yùn)行為下拋石的空間分布。比較圖6得出，在超擴(kuò)散輸運(yùn)狀態(tài)下，拋石會(huì)在最短時(shí)間到達(dá)對(duì)岸，這是由于跳躍步長概率密度分布的二階矩發(fā)散導(dǎo)致拋石在短時(shí)間產(chǎn)生長程的輸運(yùn)行為；隨著次擴(kuò)散輸運(yùn)行為的增強(qiáng)，拋石到達(dá)對(duì)岸的時(shí)間逐漸增加。

6-a 正常輸運(yùn)行為下粒子不同時(shí)刻的空間分布6-b 次擴(kuò)散輸運(yùn)行為下粒子不同時(shí)刻的空間分布6-c 超擴(kuò)散輸運(yùn)行為下粒子不同時(shí)刻的空間分布圖6 沿上游河岸釋放拋石Fig.6 Release of riprap along the upper reaches of the riverside

4.3 上游沿水流方向釋放拋石以及同時(shí)釋放不同粒徑拋石

研究中對(duì)比了以線源方式釋放拋石及同時(shí)釋放不同粒徑拋石的情況。由于河床上由不同粒徑的砂石組成的聚集體對(duì)一些拋石有阻隔作用，而另一些拋石有概率進(jìn)入河床上的裂隙結(jié)構(gòu)，這樣的裂隙結(jié)構(gòu)使進(jìn)入其中的拋石快速向下游運(yùn)動(dòng)，導(dǎo)致實(shí)際運(yùn)動(dòng)中拋石同時(shí)具有超擴(kuò)散和次擴(kuò)散輸運(yùn)行為。為表征這類輸運(yùn)行為，研究中參數(shù)設(shè)置為α=1.6，β=0.6。圖7-a展示了以線源釋放后，拋石的空間分布較為均勻，對(duì)水流的阻隔效果最好。圖7-b顯示次擴(kuò)散行為較強(qiáng)的粗粒徑拋石更易停留在離源較近的區(qū)域，而超擴(kuò)散行為更強(qiáng)的細(xì)粒徑拋石會(huì)更快向下游運(yùn)動(dòng)。

7-a 垂直流速方向上游連續(xù)釋放拋石的空間分布7-b 不同粒徑拋石同時(shí)釋放后的空間分布圖7 連續(xù)源釋放拋石Fig.7 Release of riprap under continuous source

4.4 不同輸運(yùn)行為的MSD及時(shí)間-空間演化分析

物理上通常將輸運(yùn)行為用均方位移來表征，其中滿足菲克定律的正常輸運(yùn)行為描述為

〈x2〉∝t

(11)

然而，當(dāng)輸運(yùn)行為的均方位移滿足下式時(shí)，則被定義為反常輸運(yùn)行為

〈x2〉∝tγ(γ≠1)

(12)

式中：0<γ<1時(shí)為次擴(kuò)散輸運(yùn)行為，當(dāng)γ>1時(shí)為超擴(kuò)散輸運(yùn)行為。

8-a 次擴(kuò)散輸運(yùn)行為下時(shí)-空演化及均方位移對(duì)比8-b 正常輸運(yùn)行為下時(shí)-空分布及均方位移8-c 超擴(kuò)散輸運(yùn)行為下時(shí)-空分布及均方位移圖8 三種擴(kuò)散行為下粒子的空間演化圖及均方位移圖Fig.8 The space evolution and mean square displacement of particles under three diffusion behaviors

通過圖8-a、8-b、8-c比較可以發(fā)現(xiàn)隨著次擴(kuò)散輸運(yùn)狀態(tài)向正常輸運(yùn)狀態(tài)再向超擴(kuò)散輸運(yùn)狀態(tài)轉(zhuǎn)變，粒子的空間分布逐漸變寬，隨著次擴(kuò)散向超擴(kuò)散現(xiàn)象的過渡，粒子由中心集中逐漸向空間均勻分布過渡。

5 結(jié)論

本研究通過連續(xù)時(shí)間隨機(jī)行走模型結(jié)合蒙特卡羅方法對(duì)拋石的輸運(yùn)行為進(jìn)行研究，并與室內(nèi)水槽實(shí)驗(yàn)結(jié)果進(jìn)行比較，得出了以下2點(diǎn)結(jié)論：

(1)連續(xù)時(shí)間隨機(jī)行走模型通過“寬”分布的等待時(shí)間和跳躍步長從微觀粒子隨機(jī)運(yùn)動(dòng)的角度解釋了拋石反常輸運(yùn)的物理機(jī)理。

(2)通過與室內(nèi)試驗(yàn)比較得出，連續(xù)時(shí)間隨機(jī)行走模型在捕捉拋石反常輸運(yùn)行為時(shí)效果比經(jīng)典的對(duì)流擴(kuò)散方程好；拋石在河道中的輸運(yùn)行為同時(shí)存在超擴(kuò)散輸運(yùn)行為和次擴(kuò)散輸運(yùn)行為。