河流岸边排放口允许排污负荷的伴随方程反演法

Inversion Method for Permitting Loadings of Pollutant from Lateral Effluents Based on Adjoint Equations

摘要

岸边侧向排放是河流接纳污水排放的主要方式，污染混合区控制的允许排污负荷是排放口管理的重要指标。伴随方程法在解决水力学反问题上表现出较大优势，但预测浓度和观测浓度误差的优化目标无法直接应用于允许排污负荷的反演，限制了伴随方程法在该问题上的应用。为此，提出一种伴随方程法误差源项的计算方法，进而构建河流岸边排放口允许排污负荷的反演模型。该方法在排放口水质模拟的基础上，利用伴随方程法，建立对流扩散方程边界条件的反演方法，进而分析河流岸边排放污染混合区的分布特点；以污染混合区长度为控制目标，设计伴随方程中优化目标的计算方法，实现基于污染混合区长度的排放口允许排放负荷的反演计算。将该方法应用于长江上游某典型支流排污口的允许排污负荷的反演，基于河段的水功能区划和污染混合区管理目标，得到排放口的允许排放负荷。模拟结果显示：正问题方法可以很好地预测岸边污染混合区的高锰酸盐指数（COD_Mn）分布，预测误差为16.7%，能满足实际工程中水质模拟的精度要求。所提出的伴随方程反演法能快速反演允许排污负荷，只需迭代6步即可使得污染混合区的长度达到1 m的精度，可以满足工程应用；继续迭代至18步，计算精度可达到0.01 m以下。该方法在实际问题中具有计算精度高、收敛速度快的特点，有广阔的工程应用前景，可为天然河流排放口的管理提供技术支撑。

Abstract

Lateral discharge serves as the primary pathway through which rivers receive sewage, and the permitted pollutant loadings, determined based on the pollutant mixing zone, represent critical parameters in discharge management. The adjoint equation method demonstrates substantial benefits in solving inverse problems in hydraulics. However, optimization objectives that rely on discrepancies between predicted and observed concentrations cannot be directly applied to determine the permissible loadings, limiting the application of the adjoint equation method to this issue. This study applies the adjoint equation method to derive both the control equation and boundary conditions specifically suited to lateral effluents utilizing the depth-averaged pollution transport equations for lateral discharges. Considering the narrow and elongated characteristics of the pollutant mixing zone in lateral discharges, a new formula for the error source term is introduced, with the length of the pollutant mixing zone defined as the primary objective. The adjustment value for lateral effluents is calculated by solving the adjoint equations and employing the BFGS optimization algorithm, which iteratively determines the permitted pollutant loadings from lateral discharges. The simulation of the forward problem establishes the foundation for solving the inverse problem. This research focuses on an outlet from a sewage treatment facility located in the upper reaches of the Yangtze River to evaluate the hydrodynamic and water quality model. The findings indicated that the water quality model accurately simulates the pollutant mixing zone, with the prediction error for the permanganate index (COD_Mn) maintained at 16.7%, meeting the precision requirements of water quality simulations in practical engineering. Following the accuracy verification in the forward problem, an experiment is conducted to evaluate the performance of the proposed inversion method. The inversion outcomes revealed that, after 18 iterations, the computational precision for the length of the pollutant mixing zone remains below 0.01 m despite two fluctuations during the convergence process due to inherent limitations of the BFGS method. In practical engineering applications, the required precision for controlling the mixing zone length is comparatively modest and is achieved within six iterations, reducing the error to 1 m. These results highlight the method's high computational accuracy and rapid convergence rate, providing valuable technical support for managing effluents in natural rivers.

Graphical abstract

关键词

河流岸边排放 / 允许排污负荷 / 伴随方程法 / 水质模拟 / 反问题 / BFGS算法

Key words

lateral effluents in rivers / permitting loadings of pollutant / adjoint equation / water quality simulation / inverse problems / bfgs method

本刊网刊

污水排入受纳水域后，会在排放口附近形成高浓度水域。污染混合区是指排放口附近污染物浓度超过管理目标浓度的过渡区域^[1‒2]。污染混合区设置的目的是利用水体的自净能力完成对污废水的最终净化，在水环境管理中是被允许的。然而，水体的自净能力有限，排放口污染负荷过高会对水环境造成极大破坏，因此国家地表水环境质量标准（GB 3838—2002）^[3]中规定：排放口所在水域划定的混合区，不得影响鱼类洄游通道及混合区外水域使用功能。由于不同水域生态环境的多样性，污染混合区的控制范围并没有统一的标准^[4‒5]，根据不同水域实际情况和控制标准，反演得到河道排放口的允许排污负荷，在河流的水质管理中极为重要。

对于中国内陆水域，岸边排放是普遍使用的排污方式^[6]。根据污染混合区的管理要求或者排放口附近的浓度场得到岸边排污口允许排污负荷属于环境水力学中的反问题，其传统解法主要有解析法^[7]、试错法^[8]和优化方法^[9‒10]等。但由于不同工况下公式概化方法存在差异，导致解析法不具备普适性，且概化过程导致计算结果精度较低；而试错法工作量大、计算效率低，只适用于简单条件下的反演。因此，对这两种方法的研究和应用已逐渐减少。优化方法是指根据模拟值与观测值之间的差距不断调整模型，得到最优解的方法。其自动化程度高、精度高、对边界条件及设计条件的适应能力强，在计算机技术高速发展的时代，已成为解决环境水力学反问题主流的方法之一^[11‒12]。Zhang等^[13]基于1维非稳态水质方程，利用遗传算法进行寻优，建立的污染源识别模型误差不超过5%。刘晓东等^[14]针对非恒定水流条件，基于Nelder‒Mead单纯形算法，提出河道1维污染源自动识别的FDM‒SIM方法。岸边排放污染混合区的模拟是2维问题，计算量相对较大，反问题计算更为复杂，目前相关的文献较少^[4]，亟需提出一种高效率、高精度的计算方法。

伴随方程法是一种基于变分法的反方法，源于气象和海洋预报领域的数据同化研究，近年来许多学者开始研究将伴随方法应用于环境水力学反问题的可行性。胡煜等^[15]基于1维对流扩散方程，利用伴随方法反演了1维河道的污染源，平均误差小于15%。杨中华等^[16]利用伴随同化法，开展了1维河流中4个污染源源强反演的数值试验，取得了较好的演算效果。Gordillo等^[17]基于梯度下降的伴随方程法，仅利用下游测站上少量的观测数据，重建了污染负荷上游边界条件。Gordillo等^[18]将伴随方程法应用于存在汊点的渠道、河网上，在可以接受的迭代次数下，可以较好地反演出入口的污染物浓度。Jing等^[19]基于伴随方程法，建立了一种估算污染负荷释放质量的模型，在瞬时源和连续源的测试工况下，实测值与反演值的相对误差控制在10%左右。Ding等^[20]利用变分原理，推导了非线性圣维南方程的伴随方程法模型，用于明渠水流的优化控制。

常用的伴随方程法将浓度观测值作为反问题的优化目标，将正问题的预测值和观察值的误差作为伴随方程的源项，进而通过伴随方程的求解将预测误差反向传播，作为正问题参数调整的基础。而排放口允许排污负荷反演的优化目标是污染混合区的范围，缺乏对应的浓度数据，现有的方法不能直接应用，需要对伴随方程法进行改进。因此，本文在岸边排放口模拟得到充分验证的基础上，设计一种伴随方程法误差源项的计算方法，进而提出岸边侧向排放口允许排污染负荷的伴随方程反演法，并在实际河流的污染排放口开展应用。

1 计算模型和方法

河流岸边污染排放问题中，侧向入流相对于河流来流的动量较小，仅仅在排污口很小的范围内产生影响，对污染混合区的范围影响不大^[21]。本文采用文献[22]深度平均的2维水动力模型求解流场，在此基础上，开展水质正问题和反问题的求解。

1.1 平面2维水质模型

根据文献[4]的研究，对于污染混合区的模拟预测，输移扩散等的物理作用占主导地位，而降解、挥发等生化反应具次要作用，相关影响不到10%。因此，描述河流岸边排放的深度平均污染物输移扩散方程可以表示为：

\frac{\partial h{U}_{i}C}{\partial x_i}=\frac{\partial }{\partial x_i}\left(h{D}_{\mathrm{t}}\frac{\partial C}{\partial x_i}\right)（1）

式中：x_i （i=1，2）为水平坐标；C为污染物指标，这里代表典型指标高锰酸盐指数（COD_Mn）；h为水深；U_i 为深度平均流速；D_t为水平混合系数。对于深度平均的2维模型，水平混合系数呈现各向异性的特性：在水流的横向方向，横向扩散起主导作用；在水流的纵向方向，纵向离散的作用远远大于横向扩散，占主导地位。然而，对于恒定排放问题，输移在水流纵向方向的作用远大于各种混合作用^[23]，也就是说，模拟结果对纵向混合系数并不敏感。为计算方便，这里采用了各向同性的水平扩散系数，并采用紊动黏性系数求解：

D_{\mathrm{t}}=\frac{v_{\mathrm{t}}}{{\sigma }_{\mathrm{c}}}（2）

式中：v_t为紊动黏性系数，由水动力模型求出；σ_c为Prandtl‒Schmidt浓度数，一般取值为1.0。

式（1）不包括时间项，所描述的模型为恒定模型。模型中，流速和水深由水动力模型给出。边界条件分为两类，在河流上游和排放口给定固定的值，在河流下游和岸边界给定零梯度条件，具体表示为：

\left\{\begin{array}{l}C(x)=B,x\in \partial {\Omega }_{\mathrm{i}\mathrm{n}\mathrm{l}\mathrm{e}\mathrm{t}};\\ \frac{\partial C}{\partial {\mathit{n}}_i}(x)=\mathrm{0},x\in (\partial {\Omega }_{\mathrm{o}\mathrm{u}\mathrm{t}\mathrm{l}\mathrm{e}\mathrm{t}},\partial {\Omega }_{\mathrm{b}\mathrm{a}\mathrm{n}\mathrm{k}})\end{array}\right.（3）

式中：\partial {\Omega }_{\mathrm{i}\mathrm{n}\mathrm{l}\mathrm{e}\mathrm{t}}为进口边界，包括河流上游和排污口入流边界；B为河流上游入流或者排放的指标值；\partial {\Omega }_{\mathrm{o}\mathrm{u}\mathrm{t}\mathrm{l}\mathrm{e}\mathrm{t}}为出流边界，即河流下游边界；\partial {\Omega }_{\mathrm{b}\mathrm{a}\mathrm{n}\mathrm{k}}为岸边界，{\mathit{n}}_i为积分域边界的单位外法线向量。

1.2 边界条件的伴随方程法

岸边侧向排放的允许排污负荷反演属于边界条件的反演问题，一般要求污染混合区的面积、长度或者宽度不能超过规定目标值。污染混合区的定义中既没有规定混合区内的指标值分布，也没有规定混合区边界的形状，无法直接以污染物指标值作为目标值。这里，根据河流侧向排放污染混合区的特点，设定污染混合区的特征点，以特征点的管理目标指标值作为反问题的优化目标，具体表示为：

F(C(B))=\frac{\mathrm{1}}{\mathrm{2}}{\int }_{\Omega }(C-C^{\mathrm{o}\mathrm{b}\mathrm{j}}{)}^{\mathrm{2}}\mathrm{d}\Omega（4）

式中：C^{\mathrm{o}\mathrm{b}\mathrm{j}}为特征点的目标指标值；\Omega为计算域；C-C^{\mathrm{o}\mathrm{b}\mathrm{j}}为在指定的边界条件下，预测值和目标值的误差值。因此，针对排污口允许排污负荷进行反演即通过边界条件的优化，使得式（4）的误差最小。

根据伴随方程法的一般原理^[24]，将目标函数式（4）和水质模拟的控制方程式（1）相加，得到满足约束条件的Lagrange函数：

L(C,B,\lambda )=\frac{\mathrm{1}}{\mathrm{2}}{\int }_{\Omega }(C-C^{\mathrm{o}\mathrm{b}\mathrm{j}}{)}^{\mathrm{2}}\mathrm{d}\Omega +{\int }_{\Omega }\lambda E(C)\mathrm{d}\Omega（5）

式中：L(·)为Lagrange函数；λ为拉格朗日乘子；E（C）为水质模拟控制方程式（1）的变形，可表示为：

E(C)=\frac{\partial h{U}_{i}C}{\partial x_i}-\frac{\partial }{\partial x_i}\left(h{D}_{\mathrm{t}}\frac{\partial C}{\partial x_i}\right)（6）

为了得到Lagrange函数的最小值，对其进行泛函分析，得到Lagrange函数的微分算子：

\delta L(C,B,\lambda )={\int }_{\Omega }(C-C^{\mathrm{o}\mathrm{b}\mathrm{j}})\delta C+\lambda E(\delta C)\mathrm{d}\Omega =\mathrm{0}（7）

利用格林公式，代入相关的边界条件和物质守恒方程，对式（7）中的第2项进行变形处理得到：

\begin{array}{l}{\int }_{\Omega }\lambda E(\delta C)\mathrm{d}\Omega =\\ {\int }_{\Omega }\lambda \left[\frac{\partial h{U}_i\delta C}{\partial x_i}-\frac{\partial }{\partial x_i}\left(h{D}_{\mathrm{t}}\frac{\partial \delta C}{\partial x_i}\right)\right]\mathrm{d}\Omega =\\ {\int }_{\partial \Omega }-\lambda h{D}_{\mathrm{t}}\frac{\partial \delta C}{\partial x_i}{\mathit{n}}_i\mathrm{d}S+\\ {\int }_{\partial \Omega }\delta C\left(\lambda h{D}_{\mathrm{t}}{U}_i{\mathit{n}}_i+h{D}_{\mathrm{t}}\frac{\partial \lambda }{\partial x_i}{\mathit{n}}_i\right)\mathrm{d}S-\\ {\int }_{\Omega }\delta C\left[\frac{\partial h{U}_i\lambda }{\partial x_i}+\frac{\partial }{\partial x_i}\left(h{D}_{\mathrm{t}}\frac{\partial \lambda }{\partial x_i}\right)\right]\mathrm{d}\Omega \end{array}（8）

式中，S为计算域边界。对于式（3）中的边界条件，均有：

\frac{\partial \delta C}{\partial x_i}{\mathit{n}}_i=\mathrm{0}（9）

将式（8）代入式（7），Lagrange函数的微分算子具体表示为：

\begin{array}{l}\delta L\left(C,B,\lambda \right)=\\ {\int }_{\partial \Omega }\delta C\left(\lambda h{D}_{\mathrm{t}}{U}_i{\mathit{n}}_i+h{D}_{\mathrm{t}}\frac{\partial \lambda }{\partial x_i}{\mathit{n}}_i\right)\mathrm{d}S-\\ {\int }_{\Omega }\delta C\left[\frac{\partial h{U}_i\lambda }{\partial x_i}+\frac{\partial }{\partial x_i}\left(h{D}_{\mathrm{t}}\frac{\partial \lambda }{\partial x_i}\right)-\left(C-C^{\mathrm{o}\mathrm{b}\mathrm{j}}\right)\right]\mathrm{d}\Omega \end{array}（10）

假设式（10）中面积分和线积分都等于0，即可保障Lagrange函数的微分算子等于0，进而得到伴随方法的控制方程和边界条件：

\frac{\partial h{U}_i\lambda }{\partial x_i}+\frac{\partial }{\partial x_i}\left(h{D}_{\mathrm{t}}\frac{\partial \lambda }{\partial x_i}\right)-\left(C-C^{\mathrm{o}\mathrm{b}\mathrm{j}}\right)=\mathrm{0}（11）

\left\{\begin{array}{l}\lambda (x)=\mathrm{0},x\in \partial {\Omega }_{\mathrm{o}\mathrm{u}\mathrm{t}\mathrm{l}\mathrm{e}\mathrm{t}};\\ \frac{\partial \lambda }{\partial {\mathit{n}}_i}(x)=\mathrm{0},x\in (\partial {\Omega }_{\mathrm{i}\mathrm{n}\mathrm{l}\mathrm{e}\mathrm{t}},\partial {\Omega }_{\mathrm{b}\mathrm{a}\mathrm{n}\mathrm{k}})\end{array}\right.（12）

通过求解式（11）和（12）的伴随方程，可以将误差源项值反演到侧向排放口，进而利用优化算法（broyden fletcher goldfarb shanno，BFGS）^[25]，确定侧向排放口污染指标值的校正值。BFGS法是用逐次修改切线刚度矩阵的方法求新近似解的拟牛顿法，既吸收了适合求解大规模优化问题的共轭梯度法计算量小的优点，又吸收了牛顿法收敛速度快的优点，同时避免了牛顿法计算Hesse矩阵及其求逆过程，对高维问题来说，该方法比牛顿法的计算量小。从收敛速度、计算工作量、所需内存等各项指标综合衡量，BFGS法是求解中小规模无约束优化问题时最有效的算法之一^[26]。最后通过迭代得到收敛的侧向排放口的允许排污负荷。

1.3 误差源项的处理

污染混合区位于排放口附近，它的形状和范围与排放口附近的浓度分布密切相关。对于河流的岸边侧向排放口，污染物排入河流后，污染物在横向扩散的同时，主要随水流向下游输移；污染混合区呈狭长状，污染混合区的长度是该类混合区的控制指标^[4]。在大型河流侧向排污口的下游，污染物指标值沿着横向的分布，从同岸到对岸逐渐减小，污染混合区的长度往往为同岸岸边的长度。另外，在排放口附近，同侧岸边的浓度沿着流道方向呈先上升后下降的变化模式，在排污口处，浓度最大。基于河流岸边侧向排放污染混合区的特点，提出了伴随方程误差源项的计算方法。

图1为伴随方程误差源项的计算示意图。

图1中：S₁～S₄为4个特征点；虚线为采用给定的排放条件下预测的污染指标值分布，该曲线中控制目标值C^{\mathrm{c}\mathrm{t}\mathrm{l}}对应的特征点为S₁和S₂，对应的污染混合区长度为L=S₂-S₁；实线为假设的目标分布，该分布对应的污染混合区长度为L^{\mathrm{c}\mathrm{t}\mathrm{l}}=S_{\mathrm{4}}-S_{\mathrm{3}}；特征点S₁和S₂对应的指标平均值记为C^{\mathrm{o}\mathrm{b}\mathrm{j}}（这里假设S₁和S₂对应的指标值相差不大），利用梯度关系有：

C^{\mathrm{o}\mathrm{b}\mathrm{j}}-C^{\mathrm{c}\mathrm{t}\mathrm{l}}\approx \left(S_{\mathrm{1}}-S_{\mathrm{3}}\right)\frac{\partial C}{\partial S}\left|{}_{S_{\mathrm{1}}}\right.\approx \left(S_{\mathrm{2}}-S_{\mathrm{4}}\right)\frac{\partial C}{\partial S}\left|{}_{S_{\mathrm{2}}}\right.（13）

整理得到预测指标值和目标指标值误差的计算公式：

C-C^{\mathrm{o}\mathrm{b}\mathrm{j}}=C^{\mathrm{c}\mathrm{t}\mathrm{l}}-C^{\mathrm{o}\mathrm{b}\mathrm{j}}\approx \frac{L-L^{\mathrm{c}\mathrm{t}\mathrm{l}}}{\frac{\mathrm{1}}{\partial C/\partial S\left|{}_{S_{\mathrm{1}}}\right.}-\frac{\mathrm{1}}{\partial C/\partial S\left|{}_{S_{\mathrm{2}}}\right.}}（14）

由此可计算得到伴随方程式（11）的误差源项。

1.4 计算流程和程序实现

水质模拟控制方程的求解以及允许排污负荷的反演基于OpenFOAM开发完成。OpenFOAM采用C++语言编写，用于数值计算的开源软件工具包，能够用于计算流体力学和其他典型对流扩散方程的求解^[27]。本文采用OpenFOAM求解了水动力和水质模拟的控制方程，在此基础上实现了水质模拟反问题的求解。具体步骤如下：

1）读取计算网格，初始化模型参数。

2）读取流场数据U、h，或者调用水动力模块求解流场数据U、h。

3）设定排污口浓度条件，求解水质模拟的控制方程式（1），得到恒定的指标值C。

4）根据第1.3节方法，计算混合区的长度L以及预测指标值与目标指标值的误差；检查优化目标函数式（4），如果满足精度要求，则终止计算。

5）将误差源项代入伴随方程式（11），并采用边界条件式（12）进行数值求解伴随方程，得到反演参数\lambda。

6）根据排污口处反演参数\lambda的数值，采用BFGS优化算法，计算排污口指标边界的更新量，得到新的排污指标值。

7）重复步骤3）～6），直至收敛或者达到迭代的最大步数。

2 排污口水质模拟及验证

正问题模拟是反问题计算的基础。本文以长江上游某典型支流污水处理厂排污口为研究对象，以高锰酸盐指数（COD_Mn）为主要关注的水质指标，在水动力计算的基础上，模拟并分析排污口附近COD_Mn空间分布，通过与实测数据的对比，验证正问题水质模拟的可靠性。

2.1 网格划分和模型构建

排污口位于长江上游某典型支流右岸，上游河段弯曲，流向由右岸向左岸倾斜（图2（a））。枯水期河宽500~600 m，河道中泓线从河流左侧向右偏斜，中泓线水深最大超过24 m，最小约12 m（图2（b））。根据这些特点，河道模拟长度选取1 700 m，其中，排污口上游500 m，下游1 200 m。

排污口的几何尺度以及污染混合区的横向尺度远小于河面宽度，且排污口附近COD_Mn梯度较大，流速比较复杂，需要加密网格以保证计算精度；同时，该河段面积相对较大，远离排污口的水域的COD_Mn分布场和速度场均较为稳定，设定较大网格尺度以减小计算量。计算所用网格如图2（c）所示，网格总数为33 980，排污口附近的网格尺寸为1 m，其余区域的网格尺度为15 m，以兼顾局部精度和整体计算效率，确定数值稳定性与收敛性。水动力控制方程为2维深度平均浅水方程，采用有限体积法离散求解。

2.2 网格划分和模型构建

于2012年1月6日现场观测了研究河段以及排污口附近的流速和COD_Mn分布。根据监测，当天河段上游断面流量为3 635 m³/s，COD_Mn约为1.61 mg/L；下游断面水位为178.23 m；排污口宽度约6 m，污水量0.41 m³/s，COD_Mn排放值为4.65 mg/L。根据以上指标，设定模型的边界条件；同时，根据文献[28]的方法，综合糙率n计算公式如下：

n=n_{\mathrm{0}}+K/h（15）

式中，n₀为沙粒糙率，K为形状系数，二者均由实测值率定所得。针对本模拟河道，率定所得的河床糙率取值为0.03～0.04。

2.3 计算结果的验证和分析

2.3.1 水动力模拟

计算模型的水动力模块可以较为精确合理地反映水下地形和不规则复杂边界对水流的影响。在排污口附近及河道对岸收缩处均形成漩涡（图3），与现场观测情况相符合。水位从上游到下游逐渐降低，其中上游收缩段由于流速增大，水位降低情况明显。根据流速实际监测点位分布，从流场中提取Ｈ1～Ｈ4共4个断面（图4）的流速分布，流速计算值和监测值对比结果如图5所示，计算误差见表1。

从图5可看出，H4断面最靠近岸边监测点的流速接近于0，相对误差较大，这里没有进行比较（记为点1）。结果表明，模型模拟的流速在整体上与监测值吻合较好，总体误差不超过9.00%。目前应用较为广泛的2维水动力水质模型有EFDC模型、MIKE21等，许多学者将其用于河流湖库的水动力水质模拟，如：卢诚等^[29]基于EFDC模型构建了神定河的水动力水质模型，其中，河口断面模拟流量的平均误差为19.63%；李添雨等^[30]基于MIKE21模型对沙河水库的水动力水质情况进行模拟，水位相对误差为14.00%；杨彦龙等^[31]基于SIMPLE算法建立了平面2维水流模型，模拟了大渡河乐山段的水动力情况，流速的平均相对误差为5.30%。综上所述，本文的水动力模拟精度与前人的研究结果相比精度相当甚至更优，验证了水动力模型的可靠性。

2.3.2 水质模拟

在水动力模拟的基础上，计算了COD_Mn的空间分布，结果如图6所示。污染物从排污口排出后，受到主流流速的限制，横向扩散较小，主要沿着岸边向下游输移，水流漩涡的存在导致污染物向上游扩散一小段距离。因此，形成的污染混合区范围纵向长度远小于横向，且纵横比越往下游越大。4个断面COD_Mn计算值与监测值的对比如图7所示，计算误差见表1。

在第2.3.1节中提到的水动力水质模拟的相关文献中，EFDC模型模拟COD_Mn的最大误差大于30.00%^[29]，MIKE21模型水质模拟精度为12.64%~43.61%^[30]，可见本文的水质模型模拟的COD_Mn值在整体上与观测值吻合较好，在天然河流的水质模拟中具有一定优势。

3 允许排放负荷的反演

排污口允许排放负荷的大小与河流条件及污染混合区的控制范围相关。根据《全国水环境容量核定技术指南》^[32]，选取河流近10年最枯月平均流量（30Q10）作为入流的水文条件。2014年以后，研究区域上游先后有A、B两座大型水电站建成运行，本河段位于A电站下游约655 km，B电站下游约376 km。受这些水库运行调度的影响，该段的水文条件处于不断的变化过程中，最枯月的平均流量有增大的变化趋势。统计2015—2022年朱沱水文站的流量，2015年2月的月均流量为3 480 m³/s，2021年2月份的月均流量为3 650 m³/s，为相对较枯的月份，较前期枯水流量显著增加^[33]。本文采用3 480 m³/s作为研究河段入流的边界条件，另外，该河段位于C水库上游约636 km，水位条件主要受C水库的影响，保持178.23 m不变。上游入流的COD_Mn保持1.61 mg/L不变。

根据研究区域于2010年发布的《水功能区划分标准》^[34]，所研究河段的水功能目标为Ⅲ类水；因此采用《地表水环境质量标准（GB 3838—2002）》^[3]的Ⅲ类水质标准作为污染混合区边界的控制值，对于COD_Mn，控制标准为6 mg/L。由于受纳水域生态环境的多样性，污染混合区的控制范围没有统一的标准，美国环境保护署甚至建议不同河流的排污口应该采用不同的控制范围。美国Idaho州环保部门规定^[35]，Thompson Creek河的污染混合区的长度限制在排污口下游100 m以内，而Squaw Creek和Salmon河中污染混合区的长度限制在排污口下游50 m以内。文献[33]建议C水库污染混合区的控制长度为100 m，控制宽度为100 m或者河面宽度的1/10。本文所研究河段位于C水库库尾影响区内，污染混合区的控制范围采用长度不超过100 m的标准。综合考虑收敛速度和计算稳定性，迭代过程中松弛因子取0.8。

基于上述讨论，利用第2节所建立的伴随方程反演法，计算排污口的允许污染负荷，得到的值为2.55 g/s，计算的收敛过程如图8所示。从图8中可以看出，在迭代过程中，背景值保持不变，排污口的排污负荷、混合区的长度不断变化，经过18步的迭代计算，混合区的长度收敛到100 m，计算精度控制到0.01 m以下。在收敛过程中，混合区长度的偏差出现两次波动，均是收敛到小于1 m，但大于0.1 m的范围，从另外一个方向重新计算迭代，这是因为所采用的BFGS方法固有的限制，混合区的偏差出现正负的变化（即大于100 m和小于100 m的交替变化），计算方法失效。可以通过增大松弛因子的方式减小震荡，但同时会降低收敛的速度，因此在实际应用中需综合考虑计算稳定性和收敛速度。对于天然河流岸边排放的允许排放负荷问题，混合区长度的精度控制没有必要如图8所示的精细程度。从图8中还可以看出：经过3步迭代，偏差控制在10.00 m以内；经6步迭代，偏差下降到1.00 m以内，说明了本文所建立方法的高效性。

4 结论

针对河流水质管理中的排污口允许排污负荷问题，本文以污染混合区长度为控制目标，改进了传统伴随方程法的误差源项计算方法，提出了河流岸边排放口允许排污负荷的伴随方程反演法。计算方法应用于长江上游某典型支流污水处理厂排放口的模拟计算，计算结果表明水动力‒水质模拟精度较高，可以满足天然河流的模拟需求。在此基础上，讨论了该河段的水功能区划和污染混合区管理目标，计算得到了该排放口的允许排放负荷。本文提出的计算方法仅需6步迭代即可满足工程应用中混合区长度的精度控制要求（1.00 m）。

通过对伴随方程法的改进，该反演方法可以被成功应用于河流岸边排放口允许排污负荷的高效准确反演，为天然河流排放口的管理提供技术支撑，更好地指导水环境保护工作。

参考文献

[1] Wu Zhouhu, Jia Hongyu.Analytic method for pollutant mixing zone in rivers[J].Advances in Water Science,2009,20(4):544‒548. doi:10.3321/j.issn:1001-6791.2009.04.015

[2] 武周虎,贾洪玉.河流污染混合区的解析计算方法[J].水科学进展,2009,20(4):544‒548. doi:10.3321/j.issn:1001-6791.2009.04.015

[3] Wu Zhouhu, Wu Wen.Theoretical analysis of pollutant mixing zone considering lateral distribution of flow velocity and diffusion coefficient[J].Environmental Science and Pollution Research International,2019,26(30):30675‒30683. doi:10.1007/s11356-018-2746-z

[4] 国家环境保护总局,国家质量监督检验检疫总局. 地表水环境质量标准:GB 3838—2002 [S].北京:中国环境科学出版社,2002. doi:10.3969/j.issn.1009-6310.2007.02.025

[5] Liu Zhaowei, Chen Yongcan, Zhu Dejun,et al.Study on the method in determination of permitting loadings of pollutant from effluents in Three Gorges Reservoir[J].Journal of Hydroelectric Engineering,2009,28(2):101‒105.

[6] 刘昭伟,陈永灿,朱德军,等.三峡水库排污口允许负荷量的计算方法研究[J].水力发电学报,2009,28(2):101‒105.

[7] Rodríguez Benítez A J, García Gómez A, Álvarez Díaz C.Definition of mixing zones in rivers[J].Environmental Fluid Mechanics,2016,16:209‒244. doi:10.1007/s10652-015-9425-0

[8] Liu Zhaowei, Chen Yongcan, Fu Jian,et al.Computation and analysis of environmental capacity in Wanzhou reach of Three Gorges Reservoir[J].Journal of Hydroelectric Engineering,2006,25(4):51‒56.

[9] 刘昭伟,陈永灿,付健,等.三峡水库万州段岸边环境容量的计算与分析[J].水力发电学报,2006,25(4):51‒56.

[10] Han Longxi, Zhu Dangsheng, Yao Qi.Water environment capacity calculating method for shallow-broad rivers[J].Journal of Hehai University(Natural Sciences),2001,29(4):72‒75.

[11] 韩龙喜,朱党生,姚琪.宽浅型河道纳污能力计算方法[J].河海大学学报(自然科学版),2001,29(4):72‒75.

[12] Zhang Hongju, Gan Shengwei, Yuan Hongzhou,et al.Analysis of controlled concentration of water quality of rivers around Taihu Lake[J].Water Resources Protection,2012,28(6):8‒11. doi:10.3969/j.issn.1004-6933.2012.05.002

[13] 张红举,甘升伟,袁洪州,等.环太湖河流入湖水质控制浓度分析[J].水资源保护,2012,28(6):8‒11. doi:10.3969/j.issn.1004-6933.2012.05.002

[14] Chen Yuanhua, Wang Peng, Jiang Jiping,et al.Contaminant point source identification of rivers chemical spills based on correlation coefficients optimization method[J].China Environmental Science,2011,31(11):1802‒1807. doi:10.1080/00405000.2010.522047

[15] 陈媛华,王鹏,姜继平,等.基于相关系数优化法的河流突发污染源项识别[J].中国环境科学,2011,31(11):1802‒1807. doi:10.1080/00405000.2010.522047

[16] Jiang Jiping, Han Feng, Zheng Yi,et al.Inverse uncertainty characteristics of pollution source identification for river chemical spill incidents by stochastic analysis[J].Frontiers of Environmental Science & Engineering,2018,12(5):1‒16. doi:10.1007/s11783-018-1081-4

[17] Dong Fei, Liu Xiaobo, Peng Wenqi,et al.Calculation methods of water environmental capacity of surface waters:Review and prospect[J].Advances in Water Science,2014,25(3):451‒463.

[18] 董飞,刘晓波,彭文启,等.地表水水环境容量计算方法回顾与展望[J].水科学进展,2014,25(3):451‒463.

[19] Liu Xiaodong, Yao Qi, Xue Hongqin,et al.Advance in inverse problems of environmental hydraulics[J].Advances in Water Science,2009,20(6):885‒893.

[20] 刘晓东,姚琪,薛红琴,等.环境水力学反问题研究进展[J].水科学进展,2009,20(6):885‒893.

[21] Zhang Shouping, Xin Xiaokang.Pollutant source identification model for water pollution incidents in small straight rivers based on genetic algorithm[J].Applied Water Science,2017,7(4):1955‒1963. doi:10.1007/s13201-015-0374-z

[22] Liu Xiaodong, Chen Liqiang, Hua Zulin,et al.Inverse problem on pollutant source identification in 1D channel[J].Journal of Hydroelectric Engineering,2014,33(6):132‒135. doi:10.1088/1126-6708/2001/11/064

[23] 刘晓东,陈立强,华祖林,等.河道一维污染源识别反问题[J].水力发电学报,2014,33(6):132‒135. doi:10.1088/1126-6708/2001/11/064

[24] Hu Yu, Zhang Wenjun, Ren Huatang,et al.Adjoint method for source term inversion in reverse process of pollutant convection-diffusion[J].Advances in Science and Technology of Water Resources,2019,39(2):7‒11.

[25] 胡煜,张文俊,任华堂,等.污染物对流‒扩散逆过程源项反演的伴随方法[J].水利水电科技进展,2019,39(2):7‒11.

[26] Yang Zhonghua, Zhou Wugang, Bai Fengpeng,et al.Inversing parameters of 1‒D stream water quality model using adjoint data assimilation method[J].Journal of Basic Science and Engineering,2018,26(2):276‒284.

[27] 杨中华,周武刚,白凤朋,等.基于伴随同化法的一维河流水质模型参数反演[J].应用基础与工程科学学报,2018,26(2):276‒284.

[28] Gordillo G, Morales‒Hernández M, García‒Navarro P.A gradient-descent adjoint method for the reconstruction of boundary conditions in a river flow nitrification model[J].Environmental Science:Processes & Impacts,2020,22(2):381‒397. doi:10.1039/C9EM00500E

[29] Gordillo G, Morales‒Hernández M, García-Navarro P.Solute transport control at channel junctions using adjoint sensitivity[J].Mathematics,2022,10(1):93. doi:10.3390/math10010093

[30] Jing Pingfei, Yang Zhonghua, Zhou Wugang,et al.Inverse estimation of finite-duration source release mass in river pollution accidents based on adjoint equation method[J].Environmental Science and Pollution Research International,2020,27(13):14679‒14689. doi:10.1007/s11356-020-07841-1

[31] Ding Yan, Wang S S Y.Optimal control of open-channel flow using adjoint sensitivity analysis[J].Journal of Hydraulic Engineering,2006,132(11):1215‒1228. doi:10.1061/(asce)0733-9429(2006)132:11(1215)

[32] Liu Zhaowei, Chen Yongcan, Zhu Dejun.Study on the concentration distribution in a trapezoidal open-channel flow with a side discharge[J].Environmental Fluid Mechanics,2007,7(6):509‒517. doi:10.1007/s10652-007-9042-7

[33] Chen Yongcan, Wang Zhiyong, Liu Zhaowei,et al.1D-2D coupled numerical model for shallow-water flows[J].Journal of Hydraulic Engineering,2012,138(2):122‒132. doi:10.1061/(asce)hy.1943-7900.0000481

[34] Liu Z, Li Z, Chen Y.Examination of the effect of dispersion terms in the depth-averaged transport equation[M].Proceedings of the 6th International Symposium on Environmental Hydraulics,Athens:CRC Press,2010:1015‒1020. doi:10.1201/b10553-178

[35] Kanca F, Ismailov M I.The inverse problem of finding the time-dependent diffusion coefficient of the heat equation from integral overdetermination data[J].Inverse Problems in Science and Engineering,2012,20(4):463‒476. doi:10.1080/17415977.2011.629093

[36] Jourak M, Nezhadhosein S, Rahpeymaii F.A new version of augmented self-scaling BFGS method[J].Journal of Mathematical Modeling,2023,11(2):323‒342. doi:10.1007/s10288-023-00544-6

[37] Najib U, Jamilu S, Abdullah S.A derivative‐free scaling memoryless Broyden-Fletcher-Goldfarb-Shanno method for solving a system of monotone nonlinear equations[J].Numerical Linear Algebra with Applications,2021,28(5):e2374. doi:10.1002/nla.2374

[38] Weller H G, Tabor G, Jasak H,et al.A tensorial approach to computational continuum mechanics using object-oriented techniques[J].Computers in Physics,1998,12(6):620‒631. doi:10.1063/1.168744

[39] Xin Wenjie.Some problems in the numerical computation of 2-D tidal flow in estuaries and bays[J].Chinese Journal of Hydrodynamics,1993(3):348‒354.

[40] 辛文杰.河口、海湾平面潮流数值计算中的几个问题[J].水动力学研究与进展(A辑),1993(3):348‒354.

[41] Lu Cheng, An Kunda, Zhang Xiaotong,et al.Water quality simulation of Shending River based on EFDC model[J].Environmental Monitoring in China,2020,36(4):106‒114.

[42] 卢诚,安堃达,张晓彤,等.基于EFDC模型的神定河水质模拟[J].中国环境监测,2020,36(4):106‒114.

[43] Li Tianyu, Li Zhenhua, Huang Bingbin,et al.Simulation on water quantity and quality of Shahe Reservoir by MIKE21 model[J].Acta Scientiae Circumstantiae,2021,41(1):293‒300.

[44] 李添雨,李振华,黄炳彬,等.基于MIKE21模型的沙河水库水量水质响应模拟研究[J].环境科学学报,2021,41(1):293‒300.

[45] Yang Yanlong, Cheng Kaiyu, Shi Jiayue,et al.Application of 2D flow mathematical model to fish habitat design in braided river channel[J].Journal of Changjiang River Scientific Research Institute,2022,39(8):65‒70.

[46] 杨彦龙,程开宇,施家月,等.二维水流数学模型在多分汊河道鱼类栖息地设计中的应用[J].长江科学院院报,2022,39(8):65‒70.

[47] 中国环境规划院.全国水环境容量核定技术指南[R].北京:中国环境规划院,2003.

[48] 黄真理,李玉樑,陈永灿,等.三峡水库水质预测和环境容量计算[M].北京:中国水利水电出版社,2006.

[49] 中华人民共和国住房和城乡建设部. 水功能区划分标准:GB/T 50594—2010 [S].北京:中国计划出版社,2010.

[50] Idaho Department of Environmental Quality. Idaho mixing zone implementation guidance [S].Boise:Idaho Department of Environmental Quality,2017.

基金资助

四川省科技计划项目(2022YFG0120)

四川省重点研发计划项目(2022YFS0474)