河相关系的随机微分方程建模与研究

宋晓龙，钟德钰，王光谦

（1.清华大学 水沙科学与水利水电工程国家重点实验室，北京 100084；2.天津大学 水利工程仿真与安全国家重点实验室，天津 300072）

1 研究背景

河流生态系统，是陆地与海洋联系的纽带，在生物圈的物质循环中起着主要作用。随着近代以来全球化的推进，借助基础河流理论和数字分析手段的发展，河流资源开发力度逐步加大，同时与河流保护的矛盾问题也日益凸显。在制定河流保护和水资源开发计划之前，应该最大限度地掌握河流相关信息，以做到合理科学地规划河流。在中国等发展中国家，河流动力学和河流管理始终是相关基础研究领域的重要课题，对实现资源的可持续利用意义重大。不过研究目标需要分阶段逐步实现，河流系统的非线性和突变性等特征始终是摆在专家学者面前的一道难题。

在很长一段时间内，线性模型在研究河湾演化规律中扮演重要的角色，揭示了一些河湾现象的主要特征［1-4］，例如:河道几何特征如河宽和弯曲率等的局部波动可以引起河道中纵向沙坝和横向上河道平面的大尺度变化；小尺度自由沙坝在小宽深比条件下倾向于向下游移动，而在大宽深比条件下有可能向上游移动；在小扰动下，河流可在长期的发展过程中达到准稳定状态，等等。虽然各种研究层出不穷，但远远称不上完美。事实上，在洪水影响下，河流通常呈现出季节性和不规则波动。在极端天气增多增强的今天，突发强降雨等事件无疑增加了我们的疑虑:是否河流系统可以与时时变化的水沙输入保持平衡，达到所谓的“动态平衡状态”；用一个单独的造床流量是否可以代表影响河道形态的各种非恒定因素［5］；是否可以另外认为河道形态是一种历史扰动的综合反映，而不是同外界时时保持平衡的一种准稳定状态［6］；文献［7-9］中有一些对气候变化下的水文响应特征研究，通过自然实验或卫星遥感分析了极端降雨对河道的潜在影响，但至今缺少对随机扰动下的河流特征关系的机制性研究。

近几十年来，众多物理学家转向现代系统理论拓展本领域的研究。其方法与以往经典科学的显著不同之处是，把自己研究的对象不再看作某个单独事物，而是看作一个整体系统，要求做整体综合性的研究和把握。特别是，表示线性对象的输入与输出间关系的传递函数法获得了广泛应用。我们关注的自然河流系统，具有高度的规律性和交互性，其非平衡态特征正获得越来越多的关注，而简单的确定性方程并不能满足相关研究的要求。1980年代以来，越来越多的研究基于扩展的传递函数进行河流演化过程的追踪和识别，典型案例包括:用传递函数噪声模型预测产汇流过程［10］，用随机微分方程代替确定性方程来数值模拟水质量因子的动态概率演化过程［11］，研究河岸带宽度和植被密度随流量变化的随机响应［12］，研究土壤水平衡动态模型在随机降雨下的随机变化［13］，研究地表径流中悬浮颗粒的随机运动［14］等等。这些方法为河流特征研究和流域管理提供了新的思路:将水流输入项视为随机源，通过概率化的研究来有效处理气候突变等外界影响下的河流响应问题。

河相关系是一种描述河流形态要素特征关系的经典工具，文献［15-18］中有大量的关于其非线性特征的研究。然而问题是，大多研究只是将焦点聚焦在利用实测数据来求解其幂函数的系数和指数本身，而很少考虑过该函数关系在随机事件的影响下可能发生怎样的变化。因此，本文将利用随机微分方程的方法，对河相关系的随机扰动提出新的见解，推动相关基础理论的发展，有效提高河流管理水平。

2 河相关系的随机微分方程模型

2.1 理论模型在河床演变学中，平滩水位是指河流滩地（唇）高程齐平，河道发生漫滩的临界水位。水位上升到平滩水位时，由于相应水流的流速大，输沙能力高，造床能力强，因此对应的所谓“平滩流量”成为反映河道排洪输沙能力和设计稳定河道几何形态的重要指标。经典研究普遍认为平滩流量可以由一个具体数值代表，无论来水来沙是否发生趋势性变化，都会随着河床的冲淤演变而随时间变化。如吴保生［19］基于河流自动调整的基本原理，根据河床受到扰动后调整速率与其当前状态和平衡状态之间差值成正比的基本规律，建立了一种平滩流量的滞后响应模型:

式中:Qt为t时刻的平滩流量；β为平滩流量的变化速率，是一个与水流能量大小和河床边界可动性有关的参数；Qe为相对某一给定水沙条件下，河床调整达到相对平衡时的平滩流量，可由式（2）表示:

式中:Qf为汛期平均流量；ξf=CtfQf，称作来沙系数；Ctf为汛期平均悬移质含沙量；K、b、c分别为待定系数和指数。

在一定的来沙水沙条件下，河床将不断调整其比降和断面形态，力求挟沙能力与上游来水来沙相适应。在此调整过程中河床所能达到的最适宜状态，或者说平衡状态，我们称为河相关系。它是通过某一特征流量，如平滩流量式（1），将河床不同断面的比降、河宽、水深、流速等特征联系在一起，用以描述冲积河床横断面形态以及河床的总轮廓［20］。如马卡维耶夫公式和Leopold［21］等都是从稳定断面形态开始，得到一组描述天然河流河床形态的经验关系式:

式中:St、Bt、Dt、Ut分别为平滩流量下t时刻某河段平均比降及入口断面的河宽、平均水深、平均流速；下标不同的4种α和m为待定系数和指数。

根据式（3）的幂函数形式，我们可以用一概化变量X来代表河相关系特征参量:

幂函数形式的河相关系满足非线性系统进行线性化处理的条件，可对式（4）求导得:

将式（1）、（2）带入式（5），并转换为微分形式，有:

如果取平滩流量滞后响应模型的多步递推模式［19］来表示Qt的通解:

式中:Q0为t=0时Qt的值，即初始平滩流量；Δt为时段长度；n为迭代时段数；i为时段编号。

根据微分方程式（6）得到一组描述河相关系的确定性的迭代解。

河流系统是一个自然结构、生态环境和经济社会相互耦合的开放系统，由于水体的流动性，系统与外界不断进行物质和能量的交换以及信息的传递，从而表现出明显的耗散性和协同性［22］。河流系统不能完全消除外界影响而独立存在，因此普通微分方程在真实反映河流演化过程的效果要大打折扣。考虑到河床形态因子受各种因素的影响，比如气候或人为环境扰动，引起水沙条件、河床边界条件等产生不确定性，因此可参照金融领域的Black-Shloes模型等经典的随机扩散模型［23-25］，假设河相关系的动态行为具有随机性，在稳态模型式（6）的基础上增加噪声模型来对其进行分析:

经典随机微分方程的噪声项通常由高斯白噪声表示。然而其缺点是高斯白噪声只能表现连续性强、较稳定的随机扰动，而不能控制突发情况的发生。实际的白噪声就像白色光包含所有可见光谱一样，也包含各种类型的噪声，依据概率分布形式不同，可主要分为高斯白噪声和非高斯白噪声。其中，常用泊松噪声（也称泊松跳）来模拟的非高斯噪声项，起着重要的传递突发信息的作用。其他领域的随机研究中已经发现了泊松跳的重要价值，并进行了大量的推广应用［26-28］。此外，传统的线性思维模式可能并不符合现实情况:天然河流系统就像分子生物系统，并不是纯粹的机械系统，而是一个生命系统，其可表现出某些自发性和特殊的复杂性。研究表明，真实的系统输入具有长相关和自相似性特征，而分子布朗运动是描述这种有机随机过程的典型工具。因此，我们有必要将分形理论也应用进随机微分方程以更好地刻画河相关系的演变特征。考虑到环境变化可能给河相关系带来较稳定的扰动，或类似脉冲的随机冲击，而这些随机过程不管是发生频率、发生时间或发生强度都是随机的，因此笔者使用三种流行的噪声模型，包括高斯白噪声、泊松噪声和分数白噪声，进行三种组合，评估随机微分方程（8）描述河相关系随机演变特征的适用性。详情如下:

2.1.1 简单的高斯白噪声模型

式中:{Wt:0≤t≤T}为标准维纳过程；m和σ分别为漂移项和扩散项参数；X0为初值。

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

2.1.2 高斯白噪声+泊松噪声组合模型

式中:Vj为跳跃强度；Nj(λj)为两种相互独立的泊松过程；λj为跳跃事件发生频率；j=u*，d*分别表示上跳和下跳；跳跃强度的对数分布规律(Y=ln(V))可由下面的双指数分布函数表示:

式中:ηu＞1，

这里用一个混合方法来模拟跳跃发生的时间点(τ1，τ2，…) 。其中，时间差Ri+1=τi+1-τi可根据参数为1λ的指数分布律计算，而λ可根据历史数据进行估计［29］。Xt∗从一个跳跃点到下一个跳跃点之间按几何布朗运动方式演化。假设一时间点t位于两跳跃点τ1，τi+1之间(τ1＜t＜τi+1) ，带跳随机微分方程（11）的离散解可表示为:

2.1.3 分数白噪声+泊松噪声组合模型

该模型与模型2的不同之处在于，含有Hurst参数H的{WtH，t≥0}是一个分数布朗运动（fBm），它是一个连续中心化高斯过程，有协方差函数:

当H=1/2时，fBm是标准布朗运动，模型3退化为模型2。

假设跳跃时间点为τ1，τ2，…，方程式（14）可以用如下改进的欧拉方法进行离散求解:

其中，为满足“准确与效率”的双重标准，使用一种快速的一维fBm生成器［30］模拟分数白噪声WtH，而用Yin［31］的方法处理fBm的差值:

2.2 模型参数的非参数估计法基于短期序列离散测量值进行参数估计是随机微分方程相关研究中的一项重要工作。国内外文献中已有大量方法，例如最大似然估计法（MLE）、广义矩估计法、模拟矩估计法和MCMC方法等，每一种都有自己的独特特征。Sørensen［32］经过详细调查和比较，证明MLE估计法由于其连续性、渐进正态性、渐进有效性等特征，在弱正则条件下表现最好。并且，MLE在从扩散项中分离跳跃项方面也具有优越性。因此，针对本研究提出的随机微分方程，下面将使用一种非参数估计方法，充分发挥MLE法的优越性，并排除传统参数估计法带来的限制，特别是重点将放在转移密度函数的有效构建上。

假设{xi，i=0，1，…，N} 为实测值序列；θ为待定参数集向量（包括确定项河道参数及噪声项参数）；Pθ(ti，xi； (ti-1，xi-1)，θ)是从xi-1到xi的转移密度函数；那么θ的极大似然函数就为下面根据如下步骤采用蒙特卡洛模拟方法对L()θ进行估计:

（1）假设时间范围{Ft:t∈[0 ，T]}。实测值时间间隔:考虑时间区间为[ ]ti-1，ti，把该区间插入M个长度为的子区间，那么以Xti-1为初值，利用Euler-Maruyama算法，就可以通过随机微分方程的计算得到Xti的近似值。重复R次，可以得到R维向量

其中，核函数取标准正态形式，带宽hpi取流行的Silverman［33］形式:

其中si表示的标准差。

对于带跳随机微分方程的情况，假设lr=ln(Xr)-ln(x)(r=1，2，…，R)，用m，n分别表iii-1ud示时间间隔Δ内发生上跳和下跳情况的数目，可以用如下4种条件概率密度函数的泊松加权和来表示区间[ti-1，ti]内的非条件概率密度函数，进而得到转移密度函数:

（3）对每一个xi重复以上步骤，就可以得到极大似然函数。经过对数变换，我们可以用Matlab中的fmincon函数求解如下负对数极大似然函数的极小值，从而得到需要的参数估计值:

3 实例分析

3.1 研究区域黄河下游高村至孙口河段（如图1）经常出现过流能力明显小于其上下游河段的驼峰现象，且在高村附近河底高程抬升幅度最大［34］。笔者也加入众多对此驼峰河段重点关注的学者［34-38］行列，在此选取1952年到2013年间黄河下游高村至孙口段的深泓线纵比降，以及高村站河流横断面形态要素（河宽、水深、流速）资料［39］作为原始数据，对本研究提出的随机微分方程模型进行非参数估计以及演变趋势研究。

根据吴保生方法［19］，利用高村站历年平滩流量、汛期平均流量、来沙系数等资料（表1），通过最小二乘法得到流量公式（1）中的系数和指数K、b、c和β分别为30.26、-0.47、0.42和0.346。根据该公式计算高村站平滩流量，并与实测值比较，如图2所示，得出结论符合较好。

图1 黄河下游平面示意图

图2 高村站的平滩流量实测值与计算值对比

表1 高村站历年汛期平均流量、来沙系数、实测及计算平滩流量

3.2 模型参数估计对于未知参数集θ，采用2.2节介绍方法进行参数估计。本例中，时间区间为［1952，1982］，初值取1952年河相关系各特征数据，时间间隔Δ=1，步长取在单次估值试算中，由于Matlab中的fmincon函数对初值依赖性较强，导致结果失真严重。经过多次尝试确定的解决方法是:先给模型一个经验初值进行计算，之后每次计算前先计算已得结果的均值作为新的初值，经过多次重复计算取全体估值的平均值为最终结果。当求解带跳模型的跳跃发生频率λj(j=u，d)时，可根据其历史数据计算［29］，例如:在离散时间区间［0，T］内，将极端事件（X≥xu或X≤xd）发生次数定义为nj，则λj=njT。对于本例中的河相关系特征量，使用1000个样本重复计算5000次，从而得到估计结果。此外为便于对比分析，应用最小二乘法估计计算了确定性方程（6）的参数m。

3.3 模拟计算基于以上的模型参数估计结果，可以通过多次重复计算（5000次），近似得到3种模型条件下全时段的河相关系概率分布演化结果。其中特别是，在一给定跳跃强度条件下，通过最大化随机分布均值与实测值的相关系数，从而可近似估计出最优跳跃发生时间点。如图3—图6为不同条件下的河相关系概率演化结果与实测数据的平面对比关系图。

从图3—图6可以看出，模拟的概率稳定分布图与实测值趋势总体相同，在一定程度上证明了随机微分方程描述河相关系动态演化趋势的有效性。然而对于模型1来说，问题是，在1990年左右随机平均解大幅度偏离真实值；似乎跟预想的一样，高斯白噪声驱动的随机过程连续性有余，突变性不足，局部模拟值过大或过小，这无疑凸显了模型对突发情况模拟能力的不足。实际动态河流系统可能需要更加复杂的噪声模型来刻画，现在的白噪声还不足以满足要求。就像有足够的离心力才能使得汽车平衡保持在赛道，白噪声随机方程还需要进一步加强。

图3 3种模型条件下比降的概率分布演化与实测值对比

图4 3种模型条件下河宽的概率分布演化与实测值对比

包含高斯白噪声和对数正态分布泊松噪声的随机微分方程最早由Merton［40］用于研究股票期权波动。考虑到气候条件突变或人为扰动可能给河相关系带来类似脉冲的随机影响，而这些随机过程不管是发生频率、发生时间或发生强度都是随机的，因此笔者考虑在高斯白噪声模型基础上，增加泊松跳来对河相关系进行研究。在该部分，河相关系的演化过程由几何布朗运动表示的连续性过程与泊松跳表示的非连续性过程组成。可以看出，跟图3—图6中各子图（a）相比，各子图（b）的计算结果有了大幅度提升，特别是，平面概率稳定带明显从平缓变得更加精致化，随机平均解也更加贴近真实解的变化。这证明带有泊松跳的随机微分方程能显著提高模拟水平，是基础随机理论的一个重要扩展。唯一的遗憾是图3—图6的（b）中局部突出的变量大大超出了自定义的有效区域解范围。某些时刻，较大（或较小）的河相关系特征量会发生过大的突变。因此带跳随机方程还不完美。除了突发性，可能还需要增加考虑非线性。

分数布朗运动的研究最早可追溯到1940年代。Mandelbrot and Aizenman［41］命名并大大推动了分形理论的发展。如今其已广泛应用于自然科学和工程科学的广泛领域，例如水力学中的湍流研究。因此，笔者建立了分数-带跳随机微分方程。整体来看，图3—图6中各子图（c）比相应子图（b）表现更好。对于河宽和流速算例，是由于其Hurst参数（0.3490和0.3014）距离0.5较远，因此促成了显著地时增扩散现象。换句话说，河宽和流速过程时间序列是显著负相关的，在相邻时间点之间长期发生高低值转换现象。同时，对于比降和水深，该模型（3）计算结果与带跳模型（2）计算结果十分接近，均成功预测出了局部跳跃点，但也漏掉了一些。还有令人遗憾的是，仍然有大量散点远离概率图中的核心有效区间。笔者怀疑可能是由于真实的泊松项要远比我们这里的简单概化方法更复杂，其跳跃过大或过小都不合适。可能有效的改进手段包括在跳跃强度或跳跃频率中增加缩放因子，或者将本文随机模型中的互不相关噪声改为有相关关系的噪声。

4 模型讨论与应用

4.1 模拟比较从图3—图6可以看出，应用随机微分方程，传统的确定性的河相关系解进化成为一种概率稳定带，从而为河流控制和河流管理提供了一种重要的新思路。显然，如果外轮廓令人满意，模型好坏评判准则主要是与有效区域解，或者，随机平均解与测量值相关度（如，相对误差-δ；Nash-Sutcliff系数-NSE）等有关。也就是说，假如我们定义离散值在概率分布平面演化图中的发生概率为有效概率，内限值为有效稳定宽度；那么可以说，有效概率越大，有效稳定宽度越小，随机平均解与测量值的相对误差越小、Nash-Sutcliff系数越大，模型整体评价越高。所以，考虑到以上因子对模型调查和比较的重要意义，在图7—图8中基于上文的概率分布演化结果，给出相对有效稳定宽度、相对误差、Nash-Sutcliff系数的演化和比较结果图。结合图3—图6，接下来我们对各河相关系特征量给出针对性的模型比较和最优建议。

图5 3种模型条件下水深的概率分布演化与实测值对比

图6 3种模型条件下流速的概率分布演化与实测值对比

图7 基于随机微分方程计算的河相关系有效概率稳定宽度时变演化

图8 基于随机微分方程计算的河相关系随机平均解与实测值比较

（1）比降:图3—图6表明泊松跳能够大幅度提高系统动态演化的模拟效果。简单高斯白噪声模型理应被抛弃。同时根据图7—图8可以判断，分数-泊松模型由于较小的有效稳定宽度、较大的有效概率和NSE系数，要胜过简单带跳随机模型。

（2）河宽:与比降算例类似，根据图3—图6，我们可以直观发现后两种泊松模型胜过简单白噪声模型的地方。然后，再分析图7—图8可以发现，简单带跳随机模型并没有大幅度落后分数-泊松模型，其有效概率稳定宽度在后半段大部分时间，和NSE系数在全程的表现，均领先对手。分数-泊松模型的优势是相对误差-δ较小、随机均解的波动情况模拟较好。综合来看，两种泊松模型应该结合在一起作为河流管理的参考资源。

（3）水深:与以上两算例不同，图3—图6和图7—图8均显示出高斯白噪声+泊松噪声模型在模拟水深动态演进过程中的优势，换句话说，水深主要受到突发随机事件的线性影响。

（4）流速:从图3—图6和图7—图8可以看出，分数-泊松模型要优越于高斯-泊松模型，特别是在随机均解波动情况的模拟、相对误差和NSE系数上的表现。显然分数-泊松模型更令人信服。

总结来看，分数-泊松模型较适合模拟比降和流速的演进，两种泊松模型最好结合起来共同管控河宽，而简单高斯-泊松模型对于探索水深特征已经足够。

4.2 系统应用系统稳定性的评估始终都是理解河流的必要工作。有一些很好的方法已经用于从整体上描述河流系统，包括能够区分河型的河床稳定指数（Zw）［42］（游荡型:Zw＜5； 分汊型:5＜Zw＜15；弯曲型:Zw＞15），代表河岸侵蚀强度的宽深比系数（ζ），以及表示水流强度的河流功率（Ω）（式（24））。接下来，我们将使用分数-泊松随机微分方程对沿高村站向下游的Zw、ζ、Ω进行概率分布演化的统计计算，结果如图9所示。

图9 基于分数-泊松随机微分方程计算的河床稳定指数、宽深比系数和河流功率随机演化图

从图9（a）中可以看出，河床稳定指数从1986年左右开始大幅下滑，在此期间，其随机均值几乎跌破Zw=5的游荡阈值线，直到接近2000年才又逐渐增大并回到分汊区间（5＜Zw＜15）。与此同时，高村站的水流强度经历缓慢增加过程（图9（b）），且在2000年之后经受的冲深作用大于展宽作用（图9（c））。其原因，经过查找文献，我们可以了解到，是由于上游小浪底水库的运行化解了下游长期少水多沙的负担。这也与相关研究［35］做出的结论，认为驼峰现象与水流流量大小成反比相吻合。整个图9显示出，河流稳定性的低谷已经在1999年被挽救，然而2008年左右出现的危机依然值得关注。并且，概率分布稳定宽度与随机均值间的正相关关系，也提醒我们未来需要提高对黄河下游河道演变的警戒水平。此外，考虑到测量解并不能完全反应真实变化过程，例如2000年左右的宽深比系数（ζ），我们也因此可以看出随机微分方程的优越性。换句话说，对河相关系特征量个别单独的测量值并不能帮助我们很好地了解河流系统性的变化，而在仅有少量汛期流量和悬移质沙量的有限条件下，我们可以通过本研究的随机方法表示的，例如宽深比系数（ζ），来进行河流的系统性评估。通过Zw、ζ、Ω的随机均值与平滩流量相关系数的比较（表2），我们可以得知流量主导的黄河下游适合用随机状态空间上的河流功率（Ω）来系统表征。近年我们收集的最新资料显示，受水利与水土保持工程拦截、上游水库调节、沿河用水增加及气候变化等多重因素影响，黄河干流水沙形势发生了明显变化，表现为来水来沙量的持续减小。在新的水沙形势和河道状况下，河道自身过流能力显著增强，同时中小流量出现的机率增加。这些都降低了洪水漫滩的机会。而长期相关的影响还需要进行系统性地测量、识别和评估。本研究具有重大的理论和现实意义。

表2 Zw、ζ、Ω的随机平均解与平滩流量的相关系数

5 结论

河流是一个由流量驱动，以河相关系为重要表现形式的动态系统。本文通过在确定性方程基础上加入随机噪声模型（包括简单高斯白噪声模型、高斯-泊松模型和分数-泊松模型），成功地建立并应用随机微分方程研究了黄河下游高村-孙口段的河相关系随时间演化的概率分布动态演化过程，为河流基础理论做出了重要贡献，同时为相关河流决策提供了一个新思路。本方法可以帮助我们在有限的条件下系统性地掌握河流变化特征。通过模型比较，本文认定同时考虑非线性和突发性扰动特征的分数-泊松扩散模型是适宜研究河相关系随机演化特征的较优模型。未来将把该方法进行更多的推广应用和进一步的完善，包括避免使用简化的噪声模型，而使用互相相关的更真实的复杂噪声模型等。