一维圣维南方程的反问题研究与计算方法

一维圣维南方程的反问题研究与计算方法

董文军，杨则燊

(1.天津大学理学院数学系)

摘要：作者对一维圣维南方程中曼宁糙率的参数辩识问题进行了研究与计算。使用最小二乘逼近的思想建立了相应最优模型的目标函数。通过Fréchet微分的概念和构造协态方程来进一步确定目标函数的下降方向，再用牛顿-辛普森迭代法来求解上述的最优模型。最后本文给出了一个计算实例。

关键词：曼宁糙率；参数辨识；拟牛顿算法；下降方向

作者简介：董文军(1969-)，男，山东文登人，博士，副教授，从事河口海岸泥沙数学模型研究。

数学模型已经成为分析明渠河道中非恒定水流运动的一个强有力的工具，而将此模型应用到实际问题并想获得精确的结果，就必须确定控制微分方程中的参数，但是从实际观点来说这些参数往往不是可测量的。其中曼宁糙率n在河道水流及其冲淤变化计算中占有重要作用，同时其影响因素又很多并且复杂。因此探求n值的变化规律长期以来一直是水力学中的重要课题[１]，目前常用河道的实测的水文资料推求n值。虽然国内已有许多学者在理论上和实验方面作了大量工作，得到了一些计算公式和方法[2～4]，但是这些计算方法还是属于半经验性的，往往与实际偏差较大，数学模型中的问题可以分为两类：正问题和反问题。在正问题中有关模型的所有参数都是已知的，模型的目的就是根据给定的初、边值条件确定一个相应的结果。反问题可归类于参数辩识的方法，它是一个数学过程，在这个过程中，根据计算域中测点观测值，通过对所建立的目标函数求最小的过程来确定控制微分方程中的参数，因此在参数辩识中往往含有数学规划和最优控制等技术。

关于河流水质模型辨识的研究，过去已展开了一些工作，有关辨识问题的数学理论和计方法也日趋完善[５]。目前国外已有了很多成果。Yih[６]等利用数值计算方法，给出了辨识顺河流方向扩散系数的例子。Slater和Durrer[7]使用最小二乘和线性规划对水库模拟的调整记性进行了研究。Yeh和Tauxe[８]使用拟线性的方法讨论了蓄水扩散性的参数辩识进行了计算。Yoon和Yeh[９]利用有限元方法计算了不均匀介质的参数辩识。Becker和Yeh[10]构造了一个关于明渠非恒定流参数辩识的“影响系数算法”。本文利用最小二乘逼近思想建立了关于一维圣维南方程中曼宁糙率的辩识模型，利用Fréchet微分的概念和构造协态方程确定了目标泛函的下降方向，并使用拟牛顿法[11]制定了计算n值的计算算法。并利用此算法给出了一个计算实例。

1 问题定义

反映明渠非恒定流运动的控制方程是圣维南方程，写成矢量形式如下：

(1)

式中：矢量=(h(x,t),u(x,t))T ；h(x,t)为沿明渠的水深；u(x,t)为流速；x 为沿河流方向的位移；t 为时间；q(x,t)为源汇项；g 为重力加速度；J 0为底坡；J f 为能坡，其值可由经验公式得到：

J f =n 2u 2R -4/3 (2) 式中：R 为水力半径；n 为曼宁糙率，为待定的参数。可以证明，对于任何的n ，

方程(1)都存在唯一的解y,我们用y=y(n)=y(x,t ；n)表示y 对n 的依赖关系。 本文的最终目标，就是要根据流场中的某些测点处各水力要素的观测值，建立相应目标泛函，用最优化方法来确定n 的分布规律，这一过程就是前面所说的反问题，为此假设在x=L 1，…，L m

处有下列观测值：

根据上述边界观测值，我们可以确

定

，使得式(1)在上述测点处的对应解y(Li,t ；n)尽量拟合z i (t),即：

(3) 式中：

为容许参数集，Ω=[0，L]×[0，T]，

V=L 2 (0，T ），L 2(Ω)=

表示空间V 中的范数，

(w,z)V 表示V 中二元素w 与z 间的内积，且有

为解决优化问题经常的病态性，代替它可以考虑目标泛函：

(4)

常数α＞0为给定系数。通过求目标函数J(n)的最小值来确定n 值。

2 基本算法

显然这是一个无约束最优化问题，可以用收敛速度快、易于计算的拟牛顿

(Quasi Newton Method)算法[11]解决此问题。具体步骤如下：(1)

给定初始

并确定J(n)的下降方向g

0=-J(n

)。令B

=1。(2)由式(1)解得(x,t；

n 0)，由式(4)根据测点x=L

1

,…，L

m

处的计算值y(L

i

,t；n

)和观测值z

i

(t)(i=1，

2，…，m)计算J(n

0)。(3)若J(n

)<λ

1

或｜g

｜<λ

2

(λ

1

,λ

2

为某特定数)，则n

→n*，

结束；否则转(4)。(4)置搜索方向d

n =g

/B

，一维搜索求α使得J(n

+ad

n

)=min

｛J(n

0+α′d

n

)｜α′≥0｝。(5)修正n

为n

1

=n

+ad

n

U

ad

，计算g

1

=-J(n

1

)。由

BFGS公式计算B

１=g

1

-g

/n

1

-n

。(6)n

1

→n

,g

1

→g

，转(2)。在计算中若｜

J(n

1)-J(n

)｜<λ

1

，且‖n1-n0‖<λ2，即两次迭代无显著变化；也可结束。在上述

算法中关键是要确定目标函数(4)的下降方向-J(n)。

3 下降方向的确定

及增量η，可选择充分小的ε，使得，由内积运算并经如下计算有：

(5)

式中=(n)η为y关于n沿方向η的Fréchet微分，且由式(1)可导出它满足下式

(6

) 构造对应式(6)的协态方程：