中国城镇化的快速发展，导致地面硬化、径流系数加大、汇流时间缩短，加之城镇区域局地气候的影响以及单位面积上社会经济资产的增多，使得“小雨大水，小水大灾，城市看海”的情况频现，尤其在城市化水平较高的平原河网区。作为平原河网水流运动模拟的重要参数—糙率，是体现河床组成、河道形态、岸壁特征、水位高低等因素共同作用下水流运动阻力大小的综合性参数，对水力计算结果影响较大。

目前，难以对河道糙率进行直接量测，常采用查表法、经验公式法、试错法进行估算或率定，但是这些方法存在明显的不足^[1–2]：1）糙率选取标准不够明确，人为因素影响较大，结果因人而异；2）当河网规模较大时，糙率率定工作量大，对技术人员要求高；3）率定结果的可靠性难以评估，准确性得不到保证。因此，河道糙率反演理论研究引起了国内外学者的广泛重视。

Becker等^[3]采用结合敏感度分析的单纯形法进行了河道糙率的反演研究；Wasantha Lal^[4]建立了基于奇异值分解的河道糙率反演模型，并提出了糙率分辨率的概念；Atanov等^[5]采用拉格拉日算子变分法对规则棱柱形河道的糙率进行了反演研究；Khatibi等^[6]采用高斯–牛顿优化方法进行河道糙率的反演研究，分析了准则函数和样本数量对反演结果的影响；Ramesh等^[7]建立了基于顺序二次规划算法的河道糙率反演模型，分析了不同噪声水平下反演模型的率定效果；董文军等^[8]利用Frechet微分概念和构造相应的协态方程确定目标函数的下降方向，并用牛顿–辛普森迭代法求解糙率优化反演模型；程伟平等^[9]基于广义逆和Backus-Gilbert反演理论，构造了用于河道糙率反演的牛顿–广义逆方法和自然逆方法，并推导了两种算法反演解估计的分辨率矩阵和单位协方差矩阵；Ding等^[10]建立了基于误差最小化分析的有限记忆准牛顿法，并在浅水河道中进行糙率的反演分析；Bilgil等^[11]将人工神经网络法应用于河道糙率的优化率定，得到了较好的试验结果；李丽等^[12]测试了自适应随机搜索（ARS）算法在资料受限条件下河网模型糙率反演的表现，得到了较理想的反演效果；Lai等^[13]在以概率方法描述系统误差的基础上，利用集合卡尔曼滤波算法对河道水动力模型的参数和状态变量进行了估计；Butler等^[14]提出了新途径求解曼宁场，但计算结果具有不确定性；陈一帆等^[15]将糙率和水位状态量同时作为河网非线性系统的变量，构建了基于扩展卡尔曼滤波算法的糙率和水位同步校正的数据同化模型，分析了糙率噪声水平、水位噪声水平等因素对模型同化效果的影响；贾有本等^[16]采用粒子群算法对淮河中游河道进行了河道糙率系数校正，指出粒子群算法具有一定随机性和不确定性。

河道糙率作为河网水力计算的重要水力参数，是时空相关变量，且具有一定的自然规律可循。作者以糙率空间分布平滑性作为先验知识约束条件，采用结合软约束的扩展卡尔曼滤波算法，构建了平原河网糙率参数动态反演模型。通过经典算例计算分析，系统检验了糙率空间分布平滑度权重、糙率初始估值和测站个数对模型动态反演效果的影响，证明了模型的实用性。

1 结合糙率空间分布平滑性的平原河网水动力参数动态反演模型 1.1 平原河网水流运动控制方程

城镇化平原区河道人工渠化明显，通常采用圣维南方程组描述其1维非恒定水流运动：

$\left\{ \begin{aligned} & \frac{{\partial Q}}{{\partial x}} + B\frac{{\partial {\textit{Z}}}}{{\partial t}} = {q_{\rm{L}}}, \\ & \frac{{\partial Q}}{{\partial t}} + \frac{\partial }{{\partial x}}(\frac{{a {Q^2}}}{A}) + gA\frac{{\partial {\textit{Z}}}}{{\partial x}} + gA\frac{{Q|Q|}}{{{K^2}}} = {q_{\rm{L}}}{v_x} \end{aligned} \right.$

(1)

式中：Q、Z分别为断面流量和水位； $A $ 、 $B $ 分别为断面过水面积和宽度；x为沿河方向位移；t为时间；q_L为河道单位长度旁侧入流；a为动量修正系数；g为重力加速度；K为流量模数，其表达式为 $K = A{R^{2/3}}/n $ ，其中，R为水力半径，n为糙率曼宁系数；v_x为旁侧入流沿河道水流方向的速度分量，通常近似为零。

圣维南方程组的时空离散采用Preissamann四点隐式差分格式（具体离散方法详见文献[17]），最终可得河道首尾两端水位与流量关系式：

$\left\{ \begin{aligned} & Q_{\rm{s}}^t + {a_1}{\textit{Z}}_{\rm{s}}^t + {b_1}{\textit{Z}}_{\rm{e}}^t = {c_1}, \\ & Q_{\rm{e}}^t + {a_2}{\textit{Z}}_{\rm{e}}^t + {b_2}{\textit{Z}}_{\rm{s}}^t = {c_2} \end{aligned} \right.$

(2)

式中： $Q_{\rm{s}}^t $ 、 $Z_{\rm{s}}^t $ 分别为河道首端断面流量和水位； $Q_{\rm{e}}^t $ 、 $Z_{\rm{e}}^t $ 分别为河道末端断面流量和水位；a₁、a₂、b₁、b₂、c₁和c₂均为系数，详见文献[17]。

1.2 糙率先验知识——糙率空间分布平滑性

通常，平原河网区河道的糙率变化在空间上表现为连续性和缓变性。基于这种自然特性，作者采用河网空间分布平滑度矩阵 $ {{S}}$ 与糙率向量n的相关公式表征糙率的空间变化量。

假设某一河网的河道糙率为N个，即糙率向量为n=[n₁，n₂，n₃， $\cdots$ ，n_N]^T。任取河道i，与其相连的河道有j₁，j₂， $\cdots$ ，j_M，则矩阵 ${{S}}$ 中对角线元素S_ii取为M，元素S_ij1，S_ij2， $\cdots$ ，S_ijM取为–1，第i行其余元素取为0。依次遍历所有河道，可得空间分布平滑度矩阵 ${{S}} $ ，该矩阵反映了河道之间的空间连接关系，进而可用(Sn)^TSn统计相邻河道间的糙率变化量。

1.3 糙率动态反演模型

当把河道糙率作为时空变量时，可采用如下通用形式描述糙率方程：

${{{n}}^t} = {{{\varPhi }}_1}({{{n}}^{t - 1}}) + {{\varGamma }}({{{n}}^{t - 1}}){{W}}$

(3)

式中：n为河道糙率变量组；t为当前时刻，t–1为上一时刻； ${{{\varPhi }}_1} $ 为已知t–1时刻条件下求t时刻的糙率函数，函数具体形式将在第1.4节中给出； ${{\varGamma }} $ 为糙率白噪声权重系数矩阵，通常可取为单位矩阵 ${{I}} $ ； ${{W}} $ 为均值等于零的糙率白噪声。

引入河网测站量测方程：

${{{Y}}^t} = {{{\varPhi }}_2}({{{n}}^t}) + {{V}}$

(4)

式中： ${{{Y}} ^t}$ 为量测变量，通常为量测断面的水位或流量； ${{{\varPhi }}_2} $ 为关于糙率n^t的非线性量测函数，即量测断面处的水流运动控制方程； $ {{V}}$ 为均值等于零的量测白噪声。

合并式（3）和（4），可建立平原河网糙率参数动态反演模型的系统方程组：

$\left\{ \begin{aligned} & {{{n}}^t} = {{{\varPhi }}_1}({{{n}}^{t - 1}}) + {{\varGamma }}({{{n}}^{t - 1}}){{W}} {\text{，}}\\ & {{{Y}}^t} = {{{\varPhi }}_2}({{{n}}^t}) + {{V}} \end{aligned} \right.$

(5)

应用扩展卡尔曼滤波递推算法，可建立平原河网糙率参数动态反演模型的求解算法如下（计算流程如图1所示）：

1）糙率预测

${{{\hat n}}^{t/t - 1}} = {{{\varPhi }}_1}({{{\hat n}}^{t - 1}})$

(6)

2）预测误差方差阵

${{{P}}^{t/t - 1}} = \frac{{\partial {{{\varPhi }}_1}}}{{\partial {{{{\hat n}}}^{t - 1}}}}{{{P}}^{t - 1}}{\left[ {\frac{{\partial {{{\varPhi }}_1}}}{{\partial {{{{\hat n}}}^{t - 1}}}}} \right]^{\rm{T}}} + {{\varGamma }}({{{\hat n}}^{t - 1}}){{{D}}^t}{{\varGamma }}{({{{\hat n}}^{t - 1}})^{\rm{T}}}$

(7)

式中，D^t为糙率白噪声协方差阵。

3）增益矩阵

${{{K}}^t} = {{{P}}^{t/t - 1}}{\left(\frac{{\partial {{{\varPhi }}_2}}}{{\partial {{{{\hat n}}}^{t/t - 1}}}}\right)^{\rm{T}}}{\left[\left(\frac{{\partial {{{\varPhi }}_2}}}{{\partial {{{{\hat n}}}^{t/t - 1}}}}\right){{{P}}^{t/t - 1}}{\left(\frac{{\partial {{{\varPhi }}_2}}}{{\partial {{{{\hat n}}}^{t/t - 1}}}}\right)^{\rm{T}}} + {{{E}}^t}\right]^{ - 1}}$

(8)

式中， ${{{E}}^t} $ 为量测白噪声协方差阵。

4）滤波校正

${{{\hat n}}^t} = {{{\hat n}}^{t/t - 1}} + {{{K}}^t}[{{{Y}}^t} - {{{\varPhi }}_2}({{{\hat n}}^{t/t - 1}})]$

(9)

5）滤波误差方差阵

${{{P}}^t} = \left({{I}} - {{{K}}^t}\frac{{\partial {{{\varPhi }}_2}}}{{\partial {{{{\hat n}}}^{t/t - 1}}}}\right){{{P}}^{t/t - 1}}$

(10)

糙率参数动态反演模型的基本思想包括预测和校正两步：在预测阶段，根据先前时刻的参数值生成当前时刻的预测值；在校正阶段，引入观测数据，采用最小方差估计方法对预测值重新分析，并修正参数。

1.4 糙率函数 ${{{\varPhi }}_1} $ 的选取

通常，前后时刻糙率变化很小，可忽略不计，即糙率函数 ${{{\varPhi }}_1} $ 可取为：

${{{n}}^t} = {{{n}}^{t - 1}}$

(11)

式（11）是基于朴素思想的糙率函数，未用到糙率先验知识，相应的糙率滤波解记为 ${{\hat n}}_{\rm{bas}}^t$ 。若要应用糙率先验知识，需将糙率先验知识以软约束条件的形式结合到模型中。当采用糙率空间分布平滑性作为软约束条件时，可建立优化目标函数：

$\varOmega = {{{(}}{{\hat n}}_{\rm{opt}}^t - {{\hat n}}_{\rm{bas}}^t{{)}}^{\rm{T}}}{{(}}{{\hat n}}_{\rm{opt}}^t - {{\hat n}}_{\rm{bas}}^t{{)}} + {{{(}}{{S\hat n}}_{\rm{opt}}^t{{)}}^{\rm{T}}}{{{\theta }}_n}{{S\hat n}}_{\rm{opt}}^t$

(12)

式中， ${{\hat n}}_{\rm{opt}}^t$ 为结合糙率先验知识的滤波解， ${{S}} $ 为河网空间分布平滑度矩阵， ${{{\theta }}_n}$ 为糙率空间分布平滑度权重系数矩阵。

将目标函数 $\varOmega $ 对 ${{\hat n}}_{\rm{opt}}^t$ 进行1阶求导得：

$\frac{{\partial \varOmega }}{{\partial {{\hat n}}_{\rm{opt}}^t}} = 2\left( {{{\hat n}}_{\rm{opt}}^t - {{\hat n}}_{\rm{bas}}^t} \right) + 2{{{S}}^{\rm{T}}}{{{\theta }}_n}{{S\hat n}}_{\rm{opt}}^t$

(13)

并令 $\displaystyle\frac{{\partial \varOmega }}{{\partial {{\hat n}}_{\rm{opt}}^t}} = 0$ ，可得 $\varOmega $ 最优时对应的糙率滤波解：

${{\hat n}}_{\rm{opt}}^t = {\left( {{{I}} + {{{S}}^{\rm{T}}}{{{\theta }}_n}{{S}}} \right)^{{{ - }}1}}{{\hat n}}_{\rm{bas}}^t$

(14)

式（14）等价于糙率函数：

${{{n}}^t} = {\left( {{{I}} + {{{S}}^{\rm{T}}}{{{\theta }}_n}{{S}}} \right)^{{{ - }}1}}{{{n}}^{t - 1}} + {{\varGamma }}({{{n}}^{t - 1}}){{W}}$

(15)

式（15）即为结合糙率空间分布平滑性的糙率函数。

2 算例分析

选用一个具有四级河道的河网进行算例分析，如图2所示。水力计算为非恒定流，河道断面资料、初始条件和边界条件等资料详见文献[18]。为了检验模型对糙率突变的空间识别能力，设置15号河道糙率为0.026，27号河道糙率为0.018，其余河道糙率均取0.023。

针对糙率空间分布平滑度权重、糙率初始估值及测站个数分别进行计算分析，结果如下：

2.1 糙率空间分布平滑度权重影响分析

测试设置为：1）在4号节点（位于一级河道）、7号节点（位于二级河道）和26号节点（位于三级河道）上布置水位监测点；2）各河道曼宁糙率初始估值均取为0.018；3）分别选取糙率空间分布平滑度权重系数θ_n=0、0.1和1.0进行测试比较，其中θ_n=0为不考虑糙率空间分布平滑性。不同糙率空间分布平滑度权重下，各河道最终的糙率反演结果见表1和图3，并选取靠近水位监测点的3号河道和远离水位监测点的16号河道展示糙率动态反演过程，见图4。

结果显示：1）在动态反演开始阶段有一个波动调整期，随后动态反演趋于稳定，最终靠近监测点的河道其糙率反演值趋于真值，远离监测点的河道其糙率反演值趋于空间分布平滑，体现了信息在空间传递上的衰减性。2）θ_n=0的糙率误差平方和的开方为0.020 94，相应的糙率空间变化量(Sn)^TSn值的开方为0.028 81；θ_n=0.1的糙率误差平方和的开方为0.013 69，相应的糙率空间变化量(Sn)^TSn值的开方为0.017 09；θ_n=1.0的糙率误差平方和的开方为0.008 98，相应的糙率空间变化量(Sn)^TSn值的开方为0.005 49；即随着糙率空间分布平滑度权重θ_n的增大，信息在空间上的影响范围增大，更多的河道糙率得到有效反演，糙率系统误差减小，糙率反演值趋于空间分布平滑，且动态反演开始阶段的波动调整幅度减小。3）随着糙率空间分布平滑度权重θ_n的增大，虽然糙率系统误差减小，但15号河道和27号河道的糙率反演误差增大，意味着模型对糙率突变的空间识别能力有所下降。因此，糙率空间分布平滑度权重θ_n不宜取太大，也不宜取太小。

2.2 糙率初始估值影响分析

测试设置为：1）在4号节点（位于一级河道）、7号节点（位于二级河道）和26号节点（位于三级河道）上布置水位监测点；2）取糙率空间分布平滑度权重θ_n=0.1；3）各河道糙率初始估值取值相同，分别取曼宁糙率n=0.018、0.025和0.032进行测试比较。

不同糙率初始估值下，各河道最终的糙率反演结果见表2和图5，并选取靠近水位监测点的3号河道和远离水位监测点的16号河道展示糙率动态反演过程，如图6所示。

结果显示：1）靠近监测点的河道其糙率反演值受初始估值影响较小，远离监测点的河道其糙率反演值受初始估值影响较大，符合信息空间传递衰减性。2）n=0.018的糙率误差平方和的开方为0.013 69，相应的糙率空间变化量(Sn)^TSn值的开方为0.017 09；n=0.025的糙率误差平方和的开方为0.007 42，相应的糙率空间变化量(Sn)^TSn值的开方为0.010 42；n=0.032的糙率误差平方和的开方为0.018 38，相应的糙率空间变化量(Sn)^TSn值的开方为0.013 94；这意味着糙率初始设置越接近真值，糙率系统误差越小，动态反演开始阶段的波动调整幅度越小，反演效果越好。

2.3 测站个数影响分析

测试设置为：1）取糙率空间分布平滑度权重θ_n=0.1；2）各河道曼宁糙率初始估值均取为0.018；3）分别选取测点数为3（布置在节点4、7和26）、4（布置在节点4、7、23和26）和5（布置在节点4、7、14、23和26）进行测试比较。

不同测站个数下，各河道最终的糙率反演结果见表3和图7，并选取靠近水位监测点的3号河道和远离水位监测点的16号河道展示糙率动态反演过程，如图8所示。

结果显示：测点数为3的糙率误差平方和的开方为0.013 69，相应的糙率空间变化量(Sn)^TSn值的开方为0.017 09；测点数为4的糙率误差平方和的开方为0.011 72，相应的糙率空间变化量(Sn)^TSn值的开方为0.016 94；测点数为5的糙率误差平方和的开方为0.009 37，相应的糙率空间变化量(Sn)^TSn值的开方为0.016 74；这意味着随着测点个数增加，可用于反演的信息量增多，更多的河道糙率得到有效反演，模型对糙率突变的空间识别能力有所提高，糙率系统误差减小，反演效果越好。2）3号河道随着附近监测点的增多，动态反演开始阶段的波动调整幅度减小，且反演值趋于真值；16号河道由于远离新增的测点，反演效果没有明显改善。

3 结　论

以河道糙率空间分布平滑性作为先验知识，结合扩展卡尔曼滤波器，构建了一个稳健的城镇化平原河网水动力参数动态反演模型，通过算例分析系统检验了糙率空间分布平滑度权重、糙率初始估值和测站个数对模型动态反演效果的影响，得到如下结论：

1）通过调节糙率空间分布平滑度权重，可以有效控制动态反演开始阶段的不稳定性，并在保障模型具有参数突变识别能力的基础上，使动态反演结果符合糙率空间分布平滑性。

2）靠近监测点的河道糙率反演值趋于真值，且受初始估值影响较小，远离监测点的河道糙率反演值趋于空间分布平滑。

3）糙率初始估值越精确，反演用的监测点越多，则模型动态反演越稳定，反演结果越合理，更多河道糙率得到有效反演。