瞿 波,愛(ài)迪生·保爾, 董艷梅

(1.南通大學(xué)理學(xué)院,江蘇 南通 226007; 2.英國(guó)卡迪歐數(shù)值公司,英國(guó) 愛(ài)丁堡 EH33 1EH;3.昆明理工大學(xué)城市學(xué)院,云南 昆明 650051)

?

開(kāi)放海域污染物擴(kuò)散的分形模擬

瞿 波,愛(ài)迪生·保爾2, 董艷梅3

(1.南通大學(xué)理學(xué)院,江蘇 南通 226007; 2.英國(guó)卡迪歐數(shù)值公司,英國(guó) 愛(ài)丁堡 EH33 1EH;3.昆明理工大學(xué)城市學(xué)院,云南 昆明 650051)

通過(guò)對(duì)傳統(tǒng)分?jǐn)?shù)布朗運(yùn)動(dòng)模型以及改進(jìn)了的FBMINC模型的分析,與粒子跟蹤模型結(jié)合,提出應(yīng)用范圍更廣的加速分?jǐn)?shù)布朗運(yùn)動(dòng)(AFBM)粒子跟蹤模型,并用于模擬海洋表面漂浮污染團(tuán)的擴(kuò)散。模擬結(jié)果表明,與傳統(tǒng)模型的模擬結(jié)果比較,加速分?jǐn)?shù)布朗運(yùn)動(dòng)粒子跟蹤模型更準(zhǔn)確可靠;由于豪斯特指數(shù)的取值范圍得到擴(kuò)展,加速分?jǐn)?shù)布朗運(yùn)動(dòng)粒子跟蹤模型的應(yīng)用范圍更廣泛。

擴(kuò)散；加速分?jǐn)?shù)布朗運(yùn)動(dòng);粒子追蹤模型;分形;開(kāi)放海域

研究海洋表面的浮標(biāo)軌跡及污染物的擴(kuò)散方法有很多，如蘇京志等[1]建立了以海表地轉(zhuǎn)流為主要回歸自變量的回歸模型來(lái)模擬浮標(biāo)漂流軌跡;郭民權(quán)等[2]在Daniel 等[3]提供的方法的基礎(chǔ)上,通過(guò)追蹤海表漂移物軌跡和目標(biāo)大小,用拉格朗日算法,建立了海表目標(biāo)物漂移的預(yù)測(cè)模型,用樹(shù)枝模擬的方法研究了浮子漂移在平潭近岸海域的軌跡;王勝正等[4]結(jié)合海表目標(biāo)漂移物的形狀和統(tǒng)計(jì)模型,模擬了海表漂移物的軌跡。海表漂浮污染物的傳播和風(fēng)速有關(guān)；陳麗萍等[5]發(fā)現(xiàn),隨著風(fēng)速的增大,水氣交界面附近的速度梯度會(huì)增大,從而增大揮發(fā)性污染團(tuán)的揮發(fā);韓龍喜等[6]在研究連云港近岸海域污染物隨潮流輸移和濃度增量的變化規(guī)律時(shí),采用了二維數(shù)學(xué)模型模擬研究海域的流場(chǎng)和污染物輸移軌跡,發(fā)現(xiàn)污染物的濃度增量與污染物的排放量成正比,與水體流速、流量成反比;武國(guó)正等[7]研究了分形維數(shù)算法在水體營(yíng)養(yǎng)狀況評(píng)價(jià)中的適用性。本文用分形中的分?jǐn)?shù)布朗運(yùn)動(dòng)來(lái)研究流體中污染物的傳播。

由于更適合對(duì)普遍的、各種不同程度擴(kuò)散的模擬[8-10],分?jǐn)?shù)布朗運(yùn)動(dòng)已經(jīng)越來(lái)越受到人們的青睞。Qu等[11]研究了分?jǐn)?shù)布朗運(yùn)動(dòng)的FBMINC 模型, 并與Mandelbrot的分?jǐn)?shù)布朗運(yùn)動(dòng)模型(FBM)進(jìn)行比較, 得出比FBM更精準(zhǔn)的結(jié)果。瞿波等[12]用FBMINC模型模擬了沿海海表面污染物的擴(kuò)散,并與傳統(tǒng)的布朗運(yùn)動(dòng)進(jìn)行了比較,發(fā)現(xiàn)控制擴(kuò)散程度的豪斯特指數(shù)H(00.5)。Sanderson等[14]發(fā)現(xiàn),衛(wèi)星跟蹤的海面浮標(biāo)的軌跡尺度在10～4 000 m之間,其分形維數(shù)Df在1.34左右。根據(jù)分形維數(shù)和豪斯特指數(shù)之間的關(guān)系(Df=1/H), 豪斯特指數(shù)H在0.8左右。Sanderson等[15]將擴(kuò)散范圍擴(kuò)展到了8～150 km尺度，指出對(duì)于寬闊的海面,漂流物的豪斯特指數(shù)可以超過(guò)1,并首次提出了加速分?jǐn)?shù)布朗運(yùn)動(dòng)。

分?jǐn)?shù)布朗運(yùn)動(dòng)因其靈活性,不僅被應(yīng)用于流體中[16], 而且被應(yīng)用于日常生活的很多分支, 包括金融、股市、醫(yī)療器械(腦電圖)、人體的毛細(xì)胞血管分布、巖石的裂縫等。而分?jǐn)?shù)布朗運(yùn)動(dòng)粒子追蹤模型在沿海、河流等區(qū)域有廣泛的應(yīng)用。對(duì)寬闊的海面,分?jǐn)?shù)布朗運(yùn)動(dòng)中的豪斯特指數(shù)的限制(0

1 分?jǐn)?shù)布朗運(yùn)動(dòng)模型

Mandelbrot的分?jǐn)?shù)布朗運(yùn)動(dòng)模型[17]為

(1)

式中:BH(t)為具有豪斯特指數(shù)H的分?jǐn)?shù)布朗運(yùn)動(dòng)在時(shí)刻t的軌跡；B(s)為布朗運(yùn)動(dòng)在時(shí)刻s的軌跡;Γ(x)為伽瑪函數(shù)。

式(1)中的時(shí)間是從負(fù)無(wú)窮大開(kāi)始的。如果限制時(shí)間從零點(diǎn)開(kāi)始, 得到改進(jìn)了的分?jǐn)?shù)布朗運(yùn)動(dòng)模型:

(2)

瞿波等[11-12]對(duì)FBM作了改進(jìn),得到FBMINC模型:

(i-j-1)H-1/2]R(j)+R(j-1)}

(3)

式中：i為時(shí)間；M為記憶長(zhǎng)度；R為隨機(jī)漫步。

FBMINC模型克服了FBM模型的一些弱點(diǎn),變得更精準(zhǔn)[11],已成功應(yīng)用于海灣和河流擴(kuò)散模擬中[9,12,16]。

一簇污染物隨時(shí)間擴(kuò)散的方差和時(shí)間t的關(guān)系為

(5)

式中σc(t)為標(biāo)準(zhǔn)差。

2 加速分?jǐn)?shù)布朗運(yùn)動(dòng)模型

Sanderson等[14]發(fā)現(xiàn)海洋表面浮標(biāo)的軌跡大多呈分?jǐn)?shù)布朗運(yùn)動(dòng),其豪斯特指數(shù)在0.73～0.84之間。海洋上的浮標(biāo)分?jǐn)?shù)維數(shù)可用結(jié)構(gòu)散布技術(shù)計(jì)算[18]。 圖1是一個(gè)海上浮標(biāo)的軌跡[18],其豪斯特指數(shù)是0.8。但是在寬闊的海面上,豪斯特指數(shù)可以超過(guò)1[14-15],需要新的方法來(lái)模擬這種情形的軌跡。

圖1 海上浮標(biāo)的移動(dòng)軌跡[15]

(6)

式中β為常數(shù),β≥ 0。如果β=0, 將變成分?jǐn)?shù)布朗運(yùn)動(dòng)(式(5))。因此,加速分?jǐn)?shù)布朗運(yùn)動(dòng)可以通過(guò)對(duì)時(shí)間的調(diào)整由分?jǐn)?shù)布朗運(yùn)動(dòng)得到。

設(shè)基于分?jǐn)?shù)布朗運(yùn)動(dòng)所需要的時(shí)間是t, 要成為加速分?jǐn)?shù)布朗運(yùn)動(dòng),所需要調(diào)整的時(shí)間為t′=tγ, 這里γ>1。根據(jù)原來(lái)的分?jǐn)?shù)布朗運(yùn)動(dòng)有

(7)

式中Xp為加速分?jǐn)?shù)布朗運(yùn)動(dòng)在時(shí)間t的軌跡。由于Xp(t)的方差滿(mǎn)足Var(Xp(t))∝t′2H,有

t′2H=(tγ)2H=t2H+β

(8)

(9)

所以,在加速分?jǐn)?shù)布朗運(yùn)動(dòng)模型中,時(shí)間步長(zhǎng)、記憶M和時(shí)間間隔Δt都需根據(jù)式(9)重新調(diào)整。

可以用lgσc-lgt圖中的斜率來(lái)計(jì)算豪斯特指數(shù)[16, 18]。 當(dāng)運(yùn)動(dòng)是加速分?jǐn)?shù)布朗運(yùn)動(dòng)時(shí),lgσc-lgt的斜率不再限于1以?xún)?nèi)，此時(shí)用Hη來(lái)表示此斜率,并稱(chēng)之為表面豪斯特指數(shù)。

假設(shè) lgσc-lgt曲線的斜率為Hη，則2H+β=2Hη,從而有

β=2(Hη-H)

(10)

因此, 對(duì)加速分?jǐn)?shù)布朗運(yùn)動(dòng)模型, 一組粒子擴(kuò)散時(shí)的標(biāo)準(zhǔn)差ση和時(shí)間的關(guān)系式如下:

ση=(2Dft)Hη

(11)

式中Df為擴(kuò)散系數(shù)。兩邊取對(duì)數(shù), 得:

lgση=Hηlgt+Hηlg(2Df)

(12)

這里, lgση-lgt曲線的斜率即為表面豪斯特指數(shù)Hη,其值可以大于1, 截距是

A=Hηlg(2Df)

(13)

可得Df=0.5×10A/Hη

(14)

(15)

將式(9)代入可得:

(16)

由此就可以模擬污染物在加速分?jǐn)?shù)布朗運(yùn)動(dòng)時(shí)的軌跡。一般的水利工程中,擴(kuò)散系數(shù)在0.01～0.02之間[16]。對(duì)海面上漂浮的污染物軌跡來(lái)說(shuō),可以取H=0.8,Df=0.01，然后計(jì)算出表面豪斯特指數(shù)Hη。此時(shí)Hη>1。

3 加速分?jǐn)?shù)布朗運(yùn)動(dòng)的粒子跟蹤模型

流體中污染團(tuán)的傳播軌跡由分流和擴(kuò)散兩個(gè)運(yùn)動(dòng)的組合來(lái)完成，其位置變化為:

Δxi=U(xi,yi)Δt+ΔBHx(xi,yi)

(17)

Δyi=V(xi,yi)Δt+ΔBHy(xi,yi)

(18)

式中：Δxi、Δyi分別為污染團(tuán)在兩個(gè)坐標(biāo)方向上的位置變化；U(xi,yi)、V(xi,yi)為在時(shí)刻t沿兩個(gè)相互垂直方向上的水流速度；ΔBHx(xi,yi)、ΔBHy(xi,yi)為在U、V方向上的擴(kuò)散位移, 可由分?jǐn)?shù)布朗運(yùn)動(dòng)FBMINC模型計(jì)算得到。在寬闊的海面上模擬污染團(tuán)的擴(kuò)散,將會(huì)用到加速分?jǐn)?shù)布朗運(yùn)動(dòng)。圖2為海面上漂浮的一個(gè)污染團(tuán)在3個(gè)時(shí)刻的形狀和位置(圖中O1B和O2C分別代表時(shí)刻1和時(shí)刻2的x方向，O2D代表時(shí)刻1至?xí)r刻2污染團(tuán)的擴(kuò)散方向，O2E代表時(shí)刻2流速方向的垂直方向)。

圖2 沿?cái)U(kuò)散方向的運(yùn)動(dòng)

英國(guó)諾森伯蘭海洋所在1995—1996年在英國(guó)諾森伯蘭海洋投放了幾個(gè)燃料污染團(tuán)， 英國(guó)沃靈夫特水利公司(HRWallingford,HR)用傳統(tǒng)的布朗運(yùn)動(dòng)模型對(duì)此做了模擬,筆者用加速分?jǐn)?shù)布朗運(yùn)動(dòng)粒子擴(kuò)散模型對(duì)相同的數(shù)據(jù)也進(jìn)行了模擬,結(jié)果表明筆者的模擬略勝英國(guó)沃靈夫特水利公司的模擬[16]。

圖3是對(duì)在英國(guó)諾森伯蘭海洋Horden 區(qū)域的污染物W (W1, W2, W3, W4, W5, W6代表污染物擴(kuò)散的不同時(shí)刻的等高圖)的模擬結(jié)果。本文用一簇2 000個(gè)粒子組成的粒子云(P=2 000)來(lái)模擬一個(gè)污染物體。每個(gè)粒子都遵循式(17)(18)描述的運(yùn)動(dòng)軌跡(分流加分散,而分散部分用了加速分?jǐn)?shù)布朗運(yùn)動(dòng))。對(duì)此粒子云,可以計(jì)算出平面的等高線，而每個(gè)時(shí)刻的質(zhì)量濃度由外而內(nèi)分別為0.01 μg/L、0.1 μg/L、1 μg/L、 10 μg/L、 100 μg/L。 當(dāng)?shù)雀呔€圖擴(kuò)散開(kāi)去時(shí),最大的質(zhì)量濃度從100 μg/L降低到10 μg/L,或甚至1 μg/L。先計(jì)算出沿x、y、r(x,y是坐標(biāo)軸方向,r是沿漩渦半徑的方向)、U、V方向上的粒子束擴(kuò)散運(yùn)動(dòng)的標(biāo)準(zhǔn)差(σx,σy,σr,σU,σV)。表1是沿這5個(gè)方向標(biāo)準(zhǔn)差的計(jì)算結(jié)果。 根據(jù)式(12), 可以計(jì)算出5個(gè)方向表面豪斯特指數(shù)Hη分別為1.065 4、1.090 4、1.081 6、1.232 6、0.683 3，可見(jiàn)，對(duì)污染物W,Hη的值在0.383 3～1.232 6之間,有4個(gè)方向的Hη都超過(guò)了1。

圖3 FBMINC模型的粒子云擴(kuò)散軌跡

對(duì)模擬結(jié)果和原始觀測(cè)數(shù)據(jù)以及英國(guó)Wollingford水利公司的模擬結(jié)果作了比較，比較的方法用了兩種:一種是比較標(biāo)準(zhǔn)差,另一種是比較濃度的水平輪廓面積的大小。結(jié)果表明，標(biāo)準(zhǔn)差的比較是最穩(wěn)定的方法。沿U,V方向的標(biāo)準(zhǔn)差的計(jì)算公式如下:

表1 污染團(tuán)擴(kuò)散過(guò)程的lgt 和lgσ

(19)

(20)

表2為污染團(tuán)W的模擬結(jié)果在U、V方向上的不同時(shí)刻污染物擴(kuò)散的標(biāo)準(zhǔn)差和相應(yīng)的觀測(cè)值的比較，結(jié)果顯示,模擬結(jié)果和觀測(cè)值的差距不大。模擬的結(jié)果在W1和W4低于觀測(cè)值,而在W2、W3 高于觀測(cè)值,相對(duì)誤差在±30%之內(nèi)。因?yàn)閃5、W6這兩個(gè)時(shí)刻的最大質(zhì)量濃度小于1 μg/L，本身的精確度不高，忽略而不作比較。

表2 U、 V方向模擬結(jié)果和觀測(cè)值的比較

對(duì)其他的污染物也作了比較,結(jié)果顯示,模擬和觀測(cè)結(jié)果的標(biāo)準(zhǔn)差之比在0.4～1.85之間?？紤]到在整個(gè)計(jì)算過(guò)程中始終選用了同一個(gè)H值(H=0.8),在此基礎(chǔ)上再計(jì)算表面豪斯特指數(shù)的值，因此認(rèn)為結(jié)果較為滿(mǎn)意。此外與HR的模擬面積以及觀察面積比較[16-17]， 結(jié)果表明,HR的模擬結(jié)果低估了通常的濃度水平,而本文的模擬結(jié)果高估了低濃度的水平。

4 結(jié) 語(yǔ)

本文在以往的分?jǐn)?shù)布朗運(yùn)動(dòng)的基礎(chǔ)上,介紹了加速分?jǐn)?shù)布朗運(yùn)動(dòng)模型的粒子追蹤模型，并利用AFBM粒子追蹤模型進(jìn)行模擬。根據(jù)英國(guó)森伯蘭海洋所提供的污染團(tuán)隨時(shí)間擴(kuò)散的水平面等高線圖, 將其轉(zhuǎn)化成網(wǎng)格點(diǎn)數(shù)據(jù),然后計(jì)算粒子束運(yùn)動(dòng)的標(biāo)準(zhǔn)差。對(duì)AFBM模型作了推導(dǎo),主要思路是在FBMINC模型的基礎(chǔ)上,通過(guò)調(diào)節(jié)時(shí)間參數(shù),從而調(diào)節(jié)豪斯特指數(shù)和擴(kuò)散系數(shù), 以達(dá)到加速運(yùn)動(dòng)的目的。本文用AFBM粒子跟蹤模型對(duì)2 000個(gè)粒子云進(jìn)行擴(kuò)散追蹤,計(jì)算出粒子云的標(biāo)準(zhǔn)差,并和實(shí)際數(shù)據(jù)(等高圖)進(jìn)行了比較。模擬的結(jié)果和觀測(cè)值的差距不大。盡管本文的模擬結(jié)果高估了低濃度的水平, 而HR的模擬結(jié)果低估了通常的濃度水平,總的來(lái)說(shuō),本文的模擬效果略勝一籌。

加速分?jǐn)?shù)布朗運(yùn)動(dòng)粒子跟蹤模型的表面豪斯特指數(shù)Hη具有更大的變化范圍(Hη可以大于1)，顯示出在模擬開(kāi)闊海面污染物擴(kuò)散方面的優(yōu)越性。本文介紹的方法方便簡(jiǎn)捷、實(shí)用，并可以推廣到三維的情形。

[ 1 ] 蘇京志,王東曉,陳舉,等.利用回歸模型模擬衛(wèi)星跟蹤海洋漂流浮標(biāo)軌跡[J].地球科學(xué)進(jìn)展,2005,20(6): 607-617.(SU Jingzhi,WANG Dongxiao,CHENG Ju,et al.Modeling the trajectories of satellite tracked drifters with regression models[J].Advances in Earth Science,2005,20(6): 607-617.(in Chinese)).

[ 2 ] 郭民權(quán),曾銀東,李雪丁,等.平潭近岸海域浮子漂移軌跡及其數(shù)值模擬[J].應(yīng)用海洋學(xué)學(xué)報(bào),2014,33(4): 449-454.(GUO Minquan,ZENG Yindong,LI Xueding,et al.Buoy drift trajectory and numerical simulation in Pingtan coastal waters[J].Journal of Applied Oceanography,2014,33(4): 449-454.(in Chinese)).

[ 3 ] DANIEL P,JAN G,CABIOC H F,et al.Drift modeling of cargo containers[J].Spill Science & Technology Bulletin,2002,7(5/6): 279-288.

[ 4 ] WANG Shengzheng,NIE Haobin,SHI Chaojian.A drifting trajectory prediction model based on object shape and stochastic motion features[J].Journal of Hydrodynamics,2014,26(6): 951-959.

[ 5 ] 陳麗萍,程璟濤,蔣軍成,等.風(fēng)速引起的湍流對(duì)揮發(fā)性污染物傳質(zhì)的影響[J].河海大學(xué)學(xué)報(bào) (自然科學(xué)版),2012,40(6)：610-614.(CHEN Liping,CHENG Jingtao,JIANG Juncheng,et al.Effect of turbulence induced by wind velocity on volatile pollutants mass transfer[J].Journal of Hohai University (Natural Sciences),2012,40(6): 610-614.(in Chinese))

[ 6 ] 韓龍喜,易路,劉軍英,等.連云港近岸海域污染物輸移規(guī)律[J].河海大學(xué)學(xué)報(bào)(自然科學(xué)版),2011,39(3): 248-253.(HAN Longxi,YI lu,LIU Junying,et al.Transport laws for pollutants in coastal area of Lianyungang City[J].Journal of Hohai University (Natural Sciences)，2011,39(3): 248-253.(in Chinese))[ 7 ] 武國(guó)正,徐宗學(xué),李暢游.基于分形理論的水體富營(yíng)養(yǎng)狀況評(píng)價(jià)及其驗(yàn)證[J].水資源保護(hù),2012,28(4): 12-16.(WU Guozheng,XU Zongxue,LI Changyou.Water eutrophication assessment and validation based on fractal theory[J].Water Resources Protection,2012,28(4): 12-16.(in Chinese))

[ 8 ] OSBORNE A R,KIRWAN A D,PROVENZALE A,et al.Fractal drifter trajectories in the Kuroshio extension[J].Tellus,2010,41A(5): 416-435.

[ 9 ] QU B,ADDISON P S.Modelling flow trajectories using fractional Brownian motion[EB/OL].(2015-10-24) https://www.computer.org/web/csdl/index/-/csdl/proceedings/iwcfta/2010/4247/00/4247a420-abs.html

[10] LI M.On the long-range dependence of fractional Brownian motion,mathematical problems in engineering[M].New York.:Hindawi Publishing Corporation,2013.

[11] QU B,ADDISON P S.Development of FBMINC model for particle diffusion in fluids[J].International Journal of Sediment Research,2009,24(4): 439-454.

[12] 瞿波,保爾·愛(ài)迪生.流體中污染物擴(kuò)散的分形模擬[J].水利水電科技進(jìn)展,2009,29(6): 9-12.(QU B,ADDISON P S.Simulation of fractal pollutant dispersion in fluids[J].Advances in Science and Technology of Water Resauces,2009,29 (6): 9-12.(in Chinese))

[13] SANDERSON B G,PAL B K.Patch diffusion computed from Lagrangian data,with application to the Atlantic equatorial undercurrent[J].Atmosphere-Ocean,1990,28(4): 444-465.

[14] SANDERSON B G,GOULDING A.The fractal dimension of relative Lagrangian motion[J].Tellus,1990,42(5).Part 6: 550-556.

[15] SANDERSON B G,BOOTH D A.The fractal dimension of drifter trajectories and estimates of horizontal eddy-diffusivity[J].Tellus,1991,43(5): 334-349.

[16] 瞿波.分形幾何與流體[M].上海:上海社會(huì)科學(xué)院出版社,2016.

[17] MANDELBROT B B,van NESS J W.Fractional Brownian Motions,fractal noises and Applications[J].SIAM Review,1968,10(4): 422-437.

[18] ADDISON P S.Fractals and chaos: an illustrated course[M].London: Institute of Physics Publishing.1997.

Fractal simulation of diffusion of pollutants in open ocean

//QU Bo1, ADDISON Paul S2, DONG Yanmei3

(1.ScienceFaculty,NantongUniversity,Nantong226007,China; 2.CardioDigitalLtd.,EdinburghEH33 1EH,UK; 3.CityCollegeofKunmingUniversityofScienceandTechnology,Kunming650051,China)

The traditional fractional Brownian motion (FBM) model and the improved FBMINC model were analyzed. Based on these two models and the particle tracking model, an accelerated fractional Brownian motion (AFBM) particle tracking model was developed and used to simulate the diffusion of floating pollutants in oceans. The simulation results show that, compared with the traditional model, the AFBM particle tracking model is more reliable and has a wider scope of application, due to the wider range of the Hurst exponent.

diffusion; accelerated fractional Brownian motion; particle tracking model; fractal; open ocean

10.3880/j.issn.1006-7647.2016.06.004

國(guó)家自然科學(xué)基金 (41276097)

瞿波(1962—),女,副教授,博士,主要從事分形和環(huán)境數(shù)學(xué)模型研究。E-mail:qubo@ntu.edu.cn

O29

A

1006-7647(2016)06-0020-04

2015-10-24 編輯：鄭孝宇)