抛石输运行为的连续时间随机行走模型应用

2019-11-19王茂枚李志鹏鲁程鹏孙洪广

水道港口 2019年5期

王茂枚，李志鹏，赵钢*，鲁程鹏，孙洪广

(1.江苏省水利科学研究院，南京 210017； 2.河海大学水文水资源学院，南京 210098； 3.河海大学力学与材料学院，南京 211100)

抛石护岸是河道治理的传统方法之一，作为一种在国内外广泛使用的防护技术，其工程造价小，实施方法简便，被广泛的应用于各类河道治理工程中[1]。然而，不同地域的各类河流自然条件差异巨大，在工程上广泛使用的经验公式往往与实际观测产生较大的出入[2]，因此研究抛石的输运行为在河道治理、工程应用上成为了迫切需要解决的问题。目前国内外对抛石输运行为的研究方法主要以实验和传统水沙数值模型为主[3-4]。然而研究发现，河流中泥沙及抛石的输运行为往往是随机反常的，即使同样的水沙条件，河床中砾石的启动行为也是一个概率事件[5]。在抛石的输运过程中，由于受到河床及流场演变的影响，使得抛石的输运行为往往是非马尔科夫性质的，因此基于菲克定律的经典方程在描述该类输运行为时往往达不到预期效果。这类具有非马尔可夫性质的输运行为被定义为“反常输运”[6]。

近年来，“反常输运”的研究受到国内外的广泛关注。“反常输运”在金融、医学、半导体以及地质领域上都得到了成功的应用[7-8]。目前该方向研究以两类方法为主，一类是以确定性方程为基础的，例如：“分数阶导数”、“Hausdorff导数”[9]；而另一类方法以统计理论为基础，基于事件发生的统计性质得到输运行为的随机运动方程[10]。本文所应用的连续时间随机行走理论就是基于随机运动方程描述抛石的输运行为的。

连续时间随机行走模型是研究“反常输运”行为最直观的工具，在连续时间随机行走理论框架中，一次成功的行走被分解为两部分：(1)随机的空间跳跃间隔；(2)两次成功跳跃之间的等待时间；两者分布由各自的概率分布密度支配[6，11]。本文应用连续时间随机行走模型，通过长尾的等待时间及跳跃步长分布模拟抛石在河道中的“反常输运”行为，并与实验室尺度下的实验观测数据对比进一步验证了模型对于抛石输运行为预测的有效性。

1 连续时间随机行走理论

连续时间随机行走模型中粒子的跳跃过程由空间跳跃距离和两次成功跳跃之间的等待时间组成，两者由共同的联合概率密度分布决定[6，8]。在实际研究中发现跳跃步长与等待时间往往是不相关的，即跳跃步长和等待时间分别由各自的分布密度决定。研究中λ(x)，ω(t)表示跳跃步长与等待时间的概率密度函数，则联合概率密度函数可表示为

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

(1)

连续时间随机行走模型的广义主方程可描述为

(2)

式中：η(x,t)表示粒子t时刻到达x处的概率密度，粒子的初始分布为P0(x)。

粒子的概率密度函数P(x,t)满足方程

(3)

将广义主方程(2)代入(3)式,最终得到连续时间随机行走的概率方程为

(4)

在连续时间随机行走框架中，空间和时间都被定义为连续变量，因此对方程(4)进行Fourier-Laplace变换，得

(5)

傅里叶—拉普拉斯域内的方程(5)被称为Montroll-Weiss方程，该方程中通过跳跃步长和等待时间的概率密度函数的性质可以区分不同输运行为，当跳跃步长的概率密度分布二阶矩发散而等待时间的概率密度分布一阶矩收敛时，输运行为定义为超扩散输运；跳跃步长的概率密度分布二阶矩收敛而等待时间概率密度分布一阶矩发散时，输运行为定义为次扩散输运；当两者皆发散，即同时具有超扩散输运行为和次扩散输运行为，两者皆收敛则为正常输运行为。

2 蒙特卡罗模拟方法

蒙特卡罗方法，又称随机模拟方法，利用重复随机的抽样方法，求得统计量的数字特征。当重复次数足够多时，蒙特卡罗方法得到的数值解逼近真实解。本研究中由于方程(4)在转换到时域上具有较大的困难，因此采用蒙特卡罗方法求得其数值结果，该方法是典型的拉格朗日法，可以追踪粒子运动过程中任意时刻的位置，完整呈现粒子运动轨迹。

经典的对流扩散方程(ADE)表征流体的质量传输规律，并能描述满足菲克定律的流体中粒子输运行为[9]。但实际观测中发现很多情况下菲克定律会被打破，例如含水层的异质性、河床结构的非均匀性都会导致ADE方程无法准确描述该类环境中粒子的输运行为。这类打破经典菲克定律的输运行为被称为“反常输运”。连续时间随机行走理论描述输运行为时，当概率密度分布满足中心极限定理时，即跳跃步长分布为正态分布，等待时间分布为指数分布时，方程(5)将回归到ADE方程。在本研究中，为了反映河床结构的异质性导致的抛石“反常输运”行为，采用稳定分布作为跳跃步长分布，Mittag-Leffler分布作为等待时间分布。

本研究中我们将抛石的运动行为划分为沿水流方向和垂直于水流方向的运动。因此，运动过程的时间计量由下式确定

tn+1=tn+τβ
t0=0

(6)

式中：τβ为满足Mittag-Leffler分布的随机数。

沿水流方向的位移由流场产生的确定位移和粒子相互碰撞的噪声项确定，表述为

(7)

式中：v表示粒子由对流产生的速度，可由实验获得。

垂直于水流方向的位移仅仅由粒子的相互碰撞产生，由下式确定

(8)

式中：ξα1和ξα2是满足α稳定分布的随机数。研究中随机数的生成分别通过以下公式生成

其中满足跳跃步长的α稳定分布随机数采用式(9)生成[12]

(9)

满足Mittag-Leffler分布的等待时间随机数采用式(10)生成[13]

(10)

式中：u1,u2∈(0，1)为均匀分布的随机数，γx，γt是随机数产生的尺度系数。

图1验证了产生随机数的概率密度分布以及不同参数下α稳定分布和Mittag-Leffler分布的变化特征。通过图1-a与1-c随机数概率密度与实际分布对比，随机数的生成很好的满足了所需概率密度分布。进一步分析不同参数下分布的概率密度特性，发现α稳定分布在参数α=2时回归到正态分布，并随着α的减小分布逐渐变“宽”，说明随着参数α的减小，粒子的跳跃步长愈加发散，超扩散输运行为越强烈；Mittag-Leffler分布具有和α稳定分布相同的性质，当参数β=1时，Mittag-Leffler分布回归到指数分布，当参数β变小，等待时间的分布更加发散，表明其次扩散行为越强烈。

1-a 稳定分布随机数概率密度验证1-b 不同参数下稳定分布的状态1-c Mittag-Leffler 分布随机数概率密度验证1-d 不同参数下Mittag-Leffler分布的状态图1 跳跃步长与等待时间随机数概率密度验证及不同参数下概率密度函数的差异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] 实验结果及连续时间随机行走模型和ADE模型预测效果比较Fig.2 Martin [2012] experimental results and comparison of predictive effects between continuous time random walk model and ADE

3 实验结果预测

研究中为了验证连续时间随机行走模型在实际工况中的预测效果，本节对两组不同源条件下实验室尺度水槽实验结果进行预测。表1给出了两组实验数据下连续时间随机行走模型参数的匹配情况；比较发现，在实际运动中抛石的输运行为往往是次扩散的(均方位移指数小于1)，但运动过程中同时表现出超扩散行为(1<α<2)和次扩散行为( 0<β<1)。

3.1 瞬时源条件下实验验证及分析

图3 连续时间随机行走模型预测Martin[2012] 试验的时-空分布及均方位移Fig.3 Time-Space evolution and MSD using continuous time random walk model

此处，我们选择Martin于2012年在瞬时源条件下进行的实验[14]。实验结果表明，砾石在河床上的分布并不为对称形式的正态分布，反而表现出一定程度的空间右拖尾现象，这种现象可能是由于空间的异质性导致了砾石在河床上的运动出现了反常输运的现象。图2比较了连续时间随机行走模型与经典ADE方程在实验数据预测上的优劣，通过比较发现，连续时间随机行走模型很好的捕捉了砾石的拖尾现象，而ADE方程不能捕捉具有这类拖尾现象的行为。同时图3粒子的时-空分布及均方位移指数说明该实验中砾石在河床上的运动行为是次扩散的，即一部分砾石被水流及河床结构加速，但河床结构整体上对砾石运动有阻隔作用，故砾石的整体的运动效应表现为次扩散输运状态。

3.2 连续源实验验证及分析

为了验证连续时间随机行走理论在连续源条件下以及短时对长时的预测效果，我们在这里选取了Chang and Yen于2002年进行的连续源实验[15]。图4对比了连续时间随机行走模型和经典ADE方程的预测结果，此处我们采用初始时刻t=7.5 s作为模型参数的滤定参考项，之后各个时刻均采用初始时刻的参数。相比于经典ADE方程结果，连续时间随机行走模型在短时间滤定参数进而预测长时结果时具有更好的效果。进一步分析发现，ADE方程在短时预测长时的效果上，随着时间的增加其预测误差会逐渐增大，图4-a明显表现出ADE方程在描述空间尾部效应时，其曲线尾部的衰减以指数形式衰减，而连续时间随机行走模型尾部以幂律衰减，这是两类模型在后期预测效果产生巨大差异的重要原因。

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] 实验结果及连续时间随机行走模型和ADE模型预测效果比较Fig.4 Chang and Yen[2002] experimental results and comparison of predictive effects between continuous time random walk model and ADE model

表1 连续时间随机行走模型参数及实验数据对比Tab.1 Continuous time random walk model parameters and experimental data

4 模拟结果分析

为了表征不同输运行为在抛石输运过程中的影响，本节对比了不同位置、不同抛石方式下，抛石的运动行为。分析中参数采用无量纲化处理,粒子的平均速度统一为v=0.5。

4.1 上游中心点源释放

粒子于上游中心释放，对比了正常输运行为、次扩散输运行为以及超扩散输运行为下抛石的空间分布。图5比较了次扩散-正常扩散-超扩散的输运行为差异，随着超扩散行为的增强，抛石的空间离散性逐渐增大。同时由于次扩散产生的长程等待时间，在次扩散状态下相当一部分抛石停留在离点源位置较近的区域。图中X方向是粒子具有对流速度的方向，箭头表示水流方向，Y方向是粒子仅具有随机弥散位移的方向。

5-a 正常输运行为下粒子不同时刻的空间分布5-b 次扩散输运行为下粒子不同时刻的空间分布5-c 超扩散输运行为下粒子不同时刻的空间分布图5 沿中心释放抛石Fig.5 Release of riprap along the center

4.2 上游岸边点源释放

将抛石的释放位置放置于岸边，对比正常扩散-次扩散-超扩散三种输运行为下抛石的空间分布。比较图6得出，在超扩散输运状态下，抛石会在最短时间到达对岸，这是由于跳跃步长概率密度分布的二阶矩发散导致抛石在短时间产生长程的输运行为；随着次扩散输运行为的增强，抛石到达对岸的时间逐渐增加。

6-a 正常输运行为下粒子不同时刻的空间分布6-b 次扩散输运行为下粒子不同时刻的空间分布6-c 超扩散输运行为下粒子不同时刻的空间分布图6 沿上游河岸释放抛石Fig.6 Release of riprap along the upper reaches of the riverside

4.3 上游沿水流方向释放抛石以及同时释放不同粒径抛石

研究中对比了以线源方式释放抛石及同时释放不同粒径抛石的情况。由于河床上由不同粒径的砂石组成的聚集体对一些抛石有阻隔作用，而另一些抛石有概率进入河床上的裂隙结构，这样的裂隙结构使进入其中的抛石快速向下游运动，导致实际运动中抛石同时具有超扩散和次扩散输运行为。为表征这类输运行为，研究中参数设置为α=1.6，β=0.6。图7-a展示了以线源释放后，抛石的空间分布较为均匀，对水流的阻隔效果最好。图7-b显示次扩散行为较强的粗粒径抛石更易停留在离源较近的区域，而超扩散行为更强的细粒径抛石会更快向下游运动。

7-a 垂直流速方向上游连续释放抛石的空间分布7-b 不同粒径抛石同时释放后的空间分布图7 连续源释放抛石Fig.7 Release of riprap under continuous source

4.4 不同输运行为的MSD及时间-空间演化分析

物理上通常将输运行为用均方位移来表征，其中满足菲克定律的正常输运行为描述为

〈x2〉∝t

(11)

然而，当输运行为的均方位移满足下式时，则被定义为反常输运行为

〈x2〉∝tγ(γ≠1)

(12)

式中：0<γ<1时为次扩散输运行为，当γ>1时为超扩散输运行为。

8-a 次扩散输运行为下时-空演化及均方位移对比8-b 正常输运行为下时-空分布及均方位移8-c 超扩散输运行为下时-空分布及均方位移图8 三种扩散行为下粒子的空间演化图及均方位移图Fig.8 The space evolution and mean square displacement of particles under three diffusion behaviors

通过图8-a、8-b、8-c比较可以发现随着次扩散输运状态向正常输运状态再向超扩散输运状态转变，粒子的空间分布逐渐变宽，随着次扩散向超扩散现象的过渡，粒子由中心集中逐渐向空间均匀分布过渡。