复杂水情下H-ADCP流量在线监测推流方法

刘墨阳，蒋四维，林云发，张生稳，晋 涛，黄迎春， 兰蓓蓓，周晓英，蒋飞卿，童冰星， 周文静

(1.长江水利委员会水文局汉江水文水资源勘测局，湖北 襄阳 441025； 2.河海大学水文水资源学院，江苏 南京 210098)

为提高流量自动化监测水平，国内部分水文站陆续投入定点式声学多普勒流速剖面仪(ADCP)进行比测。常见的定点式ADCP为水平式，也称水平声学多普勒流速剖面仪(H-ADCP)。该仪器安装在水中某一高度，能连续不断地发射声波信号，获取其安装高程对应水层沿横断面流速分布，具有精度高、速度快、数据传输频次高、数据量大的特点[1-2]。

H-ADCP测得的断面某一水层流速分布与断面平均流速有相关性但不等价，因此数据获取后需经过分析处理才能用于推求断面流量。断面流量推求方法主要有数值法和代表流速法，数值法也称流速剖面法，该方法基于H-ADCP获取的流速数据，根据理论垂向流速分布，推算得到全断面流量[3]；代表流速法也称指标流速法，通过分析流速网格数据关系找到合适区间得到代表流速，将代表流速与该测站常规测流方法得到的断面平均流速建立函数关系，进而推求流量。代表流速法实际应用时通过代表流速推求断面平均流速，与相应时刻水位下查算出的断面面积相乘得到断面流量。国内部分水文站如高坝洲[4]、宜昌[5]、黄陵庙[6]、长沙坝[7]、小河坝[8]等基于代表流速法建立了比较简单实用的推流公式。这些应用成果表明，实测断面平均流速与代表流速有较强的线性相关性，一般用一元一次或者二次方程就能达到满足规范和生产要求的推流精度。

近年来由于测流断面受上下游水利工程调蓄影响，部分水文站点测流断面水情复杂，代表流速与断面平均流速的关系无法采用代表流速-断面平均流速线性方程或一元二次方程描述；此外，其他影响因素如洪水涨落、回水顶托、连续洪水壅水、分洪溃口等引起水面比降改变，影响断面流速时空分布的均匀性，导致断面流速的非线性变化，增大了H-ADCP应用的难度。为解决该问题，国内外不少学者做了相关尝试，主要有两种途径：①结合先进测量仪器更好地控制断面流速监测质量。例如：韦立新等[9]将H-ADCP和垂向定点式ADCP相结合，分别得到水平平均流速和垂线平均流速，建立多元线性回归模型，得到南京水文实验站不同流速级代表流速与断面平均流速的关系；肖林[10]为解决滇池口流量在线监测问题，在新运粮河布设流量在线监测系统，通过非接触式雷达波传感器与H-ADCP分别监测表面流速和水下代表流速，采用监测数据模拟断面流速分布。②探索新方法率定基于代表流速的断面平均流速计算公式。这类方法主要是改进拟合方法，充分利用多种数据分析方案得到合理的推流公式。例如：杜兴强等[11]分不同水位级率定了高坝洲水文站的代表流速与断面平均流速的关系，建立了与低枯水水位-流量相结合的推流方案；李辉等[12]考虑断面水位变化，引入变化水位流速标定系数，建立了虚拟梯形断面平均流速与测验断面平均流速相关关系，并与直接建立代表流速与断面平均流速定线关系的传统方法进行对比，在上海苏州河温州路水文站取得了较好的应用效果；徐刚等[13-14]采用小波分析法进行滤波处理，通过BP神经网络建立断面平均流速与H-ADCP代表流速关系,实现了复杂水流河段(三峡—葛洲坝段)的流量自动测报及资料整编；邓山等[15]从粗大误差剔除、线性回归模型损失函数、代表单元格选取三方面优化代表流速法，并以小河坝站为例，证明了所提出的优化方法可显著提高模型精度；袁德忠等[16]采用支持向量机、BP神经网络、极限学习机等机器学习方法，根据清泉沟水文站H-ADCP数据模拟其断面流速分布，探究了机器学习方法与传统水文测验结合的可能；Hoitink等[17]结合代表流速法和流速剖面法，采用边界层模型基于点流速计算局部流量，通过回归模型得到断面流量并应用于感潮河段；Sassi等[18-19]在Hoitink工作的基础上，进一步考虑边壁与水位的影响推求断面流量，并在受变动回水影响的内陆河流进行了应用推广。

从现有研究可以看出引入机器学习以及深度学习等新方法不失为提升H-ADCP在生产实践中应用精度的可行途径，然而已有研究方法需要大量数据支撑，若参与训练样本量较小，模型容易出现过拟合现象，且构建推流方案时甚至会出现震荡现象。为解决该问题，本文采用较为简单实用的公式结构，在不同水情下选择合适的回归模型确定模型参数，探索了复杂水情下受水利工程影响测站的H-ADCP推流方案，可为实现流量在线监测提供参考。

1 研究方法

本文通过构建以下两种模型来研究复杂水情下受水利工程影响测站的H-ADCP推流方案：①除代表流速外，引入落差、H-ADCP仪器入水深等相关因素，建立断面平均流速与代表流速、落差、H-ADCP仪器入水深等相关因素之间的多元线性回归模型，根据最小二乘法得到模型系数；②考虑所有相关性较强的流速单元格，即每个参与率定的流速单元格在公式中具有不同的权重，应用LASSO回归模型求解参数，从而拟合出多因素影响下的断面平均流速回归方程。此外，率定模型前，分析H-ADCP流速与断面平均流速相关性，为挑选流速区间提供依据，根据实际水情选择相应模型建立推流公式，得到研究时段的推流结果并进行精度评估。

1.1 最小二乘回归模型

最小二乘回归模型遵循以残差平方和最小确定直线位置的原则，除了计算比较方便外，得到的估计量还具有无偏性、最小方差性等优良特性，但这种方法对异常值非常敏感。最小二乘法是求解线性回归模型应用最广泛的方法，其优化的目标函数为

(1)

式中：y1为输出，代表断面平均流速；x1为输入，代表流速、落差、水深等；β1为线性回归方程的系数。该方法试图找到一条直线，使所有样本到直线上的欧氏距离之和最小，可求得平方损失函数的极值点，进而得到待估参数[20]。在变量数目小于样本数目时，应用效果较好。

H-ADCP代表流速大多只能反映安装高度对应水层的流速情况。当河流流态稳定且水情单一时，可通过H-ADCP单元格流速和断面平均流速建立相关关系推流。然而水情复杂时，仅采用单一固定区间代表流速作为自变量无法满足推流精度，也很难反映测验断面受顶托影响的程度和上下游水利工程状况等，因此有必要综合考虑多种水力因素的影响。本文参考落差指数法原理，除选取代表流速外，加入总落差、上游落差、下游落差等因素反映蓄放水及回水顶托影响，引入H-ADCP入水深反映本断面流速分布代表性，建立多元线性回归模型，采用最小二乘法求解模型参数。

1.2 LASSO回归模型

LASSO回归模型最初由Tibshirani[21]于1996年提出，目的是为了消除最小二乘回归中的冗余变量，适用于参数数目缩减和参数选择。在水文气象领域，Hammami等[22-23]最早将其应用于模型变量选择和模型简化，而在国内的应用较为少见。

LASSO回归模型损失函数为L1范数，其优化的目标函数为

(2)

式中：y2为实测输出序列矩阵；x2为实测输入序列矩阵；λ为控制LASSO回归模型复杂度的正则参数，一般通过迭代选择最佳的λ，本文采用网格搜索优化该参数；β2为模型回归系数。

LASSO回归模型的L1正则化项为向量中各个元素绝对值之和，导致损失函数在零点位置不可导，其常见求解方法为坐标轴下降法和最小角回归法，本文采用坐标轴下降法求解回归系数β2[24]。坐标轴下降法沿某一坐标方向进行搜索，固定其他坐标方向，找到函数的局部极小值，不需要对目标函数求导，计算效率较高，最大迭代次数设置为2 000。

在流量较小时，引入多个水力因素后仍然无法较好拟合实测流量样本。小流量时流速较小且流场紊乱，如果将固定区间流速取均值后可能无法反映实际情况，故将所有相关性较强的单个流速单元格纳入考虑，采用LASSO回归模型求解。

1.3 评价指标

采用相对偏差、系统误差、随机不确定度等评价指标评价流量模拟精度，并进行符号检验、适线检验和偏离数值检验[25]。除评估实测点拟合效果外，还收集上下游水文站点及水利工程资料进行逐日平均流量对比，并对照上下游水文站点的年径流量，综合评估推流成果合理性。

2 应用实例

2.1 研究区概况

白河水文站是汉江流域的重要水文站点，以往水位流量关系呈一条或多条有规律可循的临时曲线。近年来，上游受白河电站蓄放水调蓄影响，且上游200 m处新建白郧大桥影响测验断面流速分布，下游受新建孤山电站试验性蓄水顶托影响，水情复杂，白河水文站水位流量关系点散乱(图1)，特别是小流量时期，同水位级对应多个流量，无法应用临时曲线法推流。

图1 白河水文站7—12月实测水位流量

白河水文站于2020年7月3日安装H-ADCP，选用600 kHz换能器收集数据，设置128个单元格，单元格宽度为1 m，安装高程为172.50 m。8月左右白河进入主汛期，出现流量峰值4 000 m3/s左右的洪水过程，伴随着含沙量增大，600 kHz换能器声波能量衰减较快，其剖面有效测量范围缩短至原有的1/3，所收集到的流速数据代表性减弱，无法满足数据采集需要，改为300 kHz换能器，设置90个单元格，单元格宽度为2m。随着下游孤山电站蓄水水位逐渐升高，H-ADCP安装相对位置不断下降，经分析发现该水层流速数据逐渐失去代表性，10月14日将仪器安装高程提高至173.77 m，单元格设置为120个。安装位置及仪器测量范围如图2所示。

图2 白河水文站H-ADCP安装位置及测量范围

2.2 数据选择及预处理

选取白河水文站及其上游大王沟水位站、下游将军河水位站、下游孤山桥水位站2020年8—12月水文资料作为研究数据。数据预处理步骤如下：

a.比测数据选择。选取走航式ADCP流量测次，数据选取遵循以下原则：①回放测验时的走航式ADCP原始数据，采用在H-ADCP正常工作时段内的测次；②当走航式ADCP在GGA模式下流量小于BTM(底跟踪)模式下流量时，选用BTM模式计算流量；③对于在水情平稳时测得的数据，根据SL337—2006《声学多普勒流量测验规范》要求，施测 2 个测回，任一半测回BTM或GGA模式下测量值与平均值的相对误差不应大于5%，否则补测同向的一个测次流量。当第一测回任一次BTM或GGA模式下测量值与平均值的相对误差小于2%时，可不施测第二测回，流量取施测第一测回的平均值。满足以上条件的视为闭合流量，否则做不闭合流量处理。将走航式 ADCP 测流流量除以相应水位对应过水断面面积，得到断面平均流速作为实测数据，为下一步公式率定做准备。

b.流速单元格选取。采用走航式ADCP流量对应起始时间段内H-ADCP流速网格数据(剔除伪值)计算每个流速单元格在该时间段的流速均值；分析每个单元格流速与实测断面平均流速的相关关系，计算两者相关系数，选取与实测流量断面平均流速相关性较好的单元格区间，准备好模型率定时所需H-ADCP流速数据。

c.计算相关水力因素。引入上下游相关辅助水位站点同时刻的摘录水位(若时间不对应，对摘录水位做线性插值处理)计算落差，用于反映上下游水利工程蓄放水、回水顶托对断面流速的影响；计算H-ADCP仪器入水深和断面平均水深，用于反映变化水位下相对水深对断面流速分布的影响。

2.3 推流方案

推流方案流程如图3所示。当H-ADCP流速数据沿横断面具有较好的分布规律时，采用最小二乘回归模型推流较为简便；当H-ADCP流速数据沿横断面分布正负不定时，采用LASSO回归模型推流。该方法引入单个流速单元格，待率定参数个数较多。

图3 推流方案流程

a.最小二乘回归法。先采用互相关系数确定代表流速预选区间[20]，选择相关系数较高的网格区间计算区间代表流速均值；然后分析模型自变量与因变量相关性，为确定公式项提供依据，并以提高拟合精度为目标，逐步增加自变量个数，采用最小二乘回归法率定回归模型。

b.LASSO回归法。白河水文站流量较小时测流断面流速分布较紊乱，单个网格数值在±0.2 m/s之间波动，H-ADCP流速区间均值代表性不佳，使用最小二乘回归法很难得到满足要求的推流精度，因此改变思路，将断面相关性较强的单元格流速数据放入LASSO回归模型中进行训练。该方法在模型计算时并未采用流速区间均值，而是考虑所有相关性较强的单个流速单元格，即每个参与率定的流速单元格在公式中具有不同的权重。此外，将相关性较好的单个流速网格数据与仪器入水深、平均水深、白河孤山桥落差等相关因素组合，进一步提高拟合精度。

需要注意的是，两种推流方法均属于多元线性回归模型求解方法，主要区别是对代表流速处理方法不同、公式形式不同(LASSO回归法引入单个流速网格，其公式项数较多)和模型参数率定方法不同。

2.4 结果分析

由于仪器型号、仪器安装位置改变以及水情变化的影响，H-ADCP获取的流速资料不具有一致性，故根据仪器型号、水情和安装位置，将资料划分为7月3日到8月9日(600 kHz换能器，安装高程172.50 m)、8月10日至9月24日(300 kHz换能器，安装高程172.50 m)、9月25日至10月13日(300 kHz换能器，安装高程172.50 m)和10月14日至12月31日(300 kHz，安装高程173.77 m)4个时段。其中，8月10日至10月13日仪器型号与安装高程未改变，但考虑到9月22—24日受孤山电站调节影响，白河水文站水位下降，退至H-ADCP安装高度以下，9月24日后电站恢复蓄水，断面水位上升至正常蓄水位，H-ADCP开始正常运行，水情有较大改变，故以9月24日为分界点再次进行时段划分。

600 kHz换能器正常运行时间较短，因此主要对8月10日至12月31日收集的后3个时段300 kHz流速数据进行分析推算。

2.4.1时段1(8月10日至9月24日)

8月10日至9月24日共测流49次，测次均闭合，实测断面平均流速范围为0.547～2.19 m/s。去除测流期间H-ADCP流速数据异常或中断的测次，共选用34次实测流量进行相关关系拟合。

拟合前对H-ADCP单元格流速数据进行相关性分析，计算每个单元格与实测断面平均流速相关性(图4)，以此作为单元格挑选依据。

图4 时段1单元格流速与实测断面平均流速相关系数

由图4可知，从单元格20开始，相关性逐渐增强，经过主泓段后，相关性逐渐减弱，到单元格67以后相关性锐减，流速相关性系数分布与断面形状较为一致。选取相关系数在0.8以上的单元格42～57区间流速均值作为代表流速，该区间基本位于主泓段，代表性较好，同时考虑上下游落差、断面形状、换能器入水深等相关因素，通过分析这些因素与断面平均流速的相关性，以提高拟合精度为原则逐步引入相关变量，与断面平均流速进行多元回归分析。拟合公式为

(3)

表1 时段1单元格系数

使用拟合公式推算时段内逐时流量，并进行平滑摘录，得到该时段流量过程线。以实测流量作为真值，推算流量过程线相应时段流量均值作为推求流量(下同)，进行误差检验。时段1使用拟合关系推算流量的系统误差为0.41%，随机不确定度为8.68%，符号检验、适线检验、偏离数值检验均通过，满足规范要求。

2.4.2时段2(9月25日至10月13日)

9月25日至10月13日共测流15次，其中闭合流量14次，实测断面平均流速范围为0.011～0.808 m/s。去除测流期间H-ADCP流速数据缺失或异常的测次后，共有闭合流量9次，不闭合流量1次。由于该时段样本数量较少，无法拟合出精度较好的相关关系，因此采用单测回流量进行拟合。同样，首先进行单元格流速相关性分析，如图5所示。

图5 时段2单元格流速与实测断面平均流速相关系数

从图5可以看出，靠近左岸的主泓段相关性较好，中间部分受上游在建白郧大桥桥墩影响相关系数较为散乱，因此选取左岸相关性较好的14～33单元格进行分析计算。因该时段总体流量较小，上下游落差变幅较小，使用时段1的分析方法拟合时精度不佳。为提高拟合精度，将14～33单元格流速值乘以H-ADCP相对水深，增强其代表性，并采用LASSO回归法进行优化求解。拟合公式为

(4)

式中：bj为H-ADCP流速14～33单元格系数，取值见表2；vind,j为H-ADCP单元格j的流速。由表2可知，参与率定的20个单元格系数不为0的共有17个。

表2 时段2单元格系数

对时段2拟合的相关关系进行误差分析和检验，对于不闭合流量测次，检验时取测次均值进行误差统计(下同)。推算流量的系统误差为-0.55%，随机不确定度为9.02%，符号检验、适线检验、偏离数值检验均通过，满足规范要求。

2.4.3时段3(10月14日至12月31日)

10月14日调整H-ADCP安装高度后共测流39次，其中闭合流量34次，不闭合流量5次，实测断面平均流速范围为0.037～0.799 m/s，其中11月13日测流期间H-ADCP异常(数据中断)，不参与后续分析。各单元格流速序列与实测断面平均流速的相关系数如图6所示。

图6 时段3单元格流速与实测断面平均流速相关系数

由图6可以看出，提高H-ADCP安装高度后，断面中间部分受上游在建白郧大桥影响，流速分布依然散乱。选取相关系数0.6以上单元格区间6～40、46～55、76～93进行分析，分析方法与时段2相同，拟合公式如下：

(5)

式中：bj、bm、bn分别为H-ADCP单元格6～40、46～55、76～93对应的系数，取值见表3；vind,j、vind,m、vind,n分别为H-ADCP单元格j、m、n的流速。

由表3可知，参与率定的63个单元格系数不为0的一共为37个，可见LASSO回归模型在保证计算精度的同时有效地减少了变量数目，并且主泓段对应的流速单元格参数权重较大，与实际情况相符。

表3 时段3单元格系数

对相关关系进行误差分析和检验，推算流量的系统误差为0.73%，随机不确定度为9.37%，符号检验、适线检验、偏离数值检验均通过，满足规范要求。

2.4.4推流成果

为消除流速脉动影响[26]，应用H-ADCP流速数据推流时，需对H-ADCP原始流速数据进行平滑处理后代入率定公式推流。平滑处理公式为

(6)

8月10日12:10之前，采用连时序法绘制水位-流量关系线进行推流，之后(H-ADCP正常运行期间)，由构建的相应时段公式采用连实测流量过程线法推流，H-ADCP中断期间采用连时序法绘制水位-流量关系线推流。

将白河水文站8月10日后逐日平均流量与上游来水逐日平均流量(白河电站+长沙坝)进行对照(图7)，两者趋势相应性较好。

图7 上游来水逐日平均流量与白河水文站逐日平均流量对比

白河水文站洪水水文要素摘录流量和实测流量如图8所示,实测流量基本位于逐日平均流量过程线之上。从上下游水量来看，2020年白河水文站年径流量233.7亿m3，白河水文站下游黄家港水文站以上各站合计来水238.1亿m3，白河水文站下游黄家港(二)站年径流量271.6亿m3。黄家港水文站以上区间各站集水面积合计82 708 km2，占黄家港(二)站集水面积的87%，上游来水量238.1亿m3，占黄家港(二)站年径流量的88%，上下游水量基本平衡，说明推流成果合理。

图8 白河水文站洪水水文要素摘录流量与实测流量对比

推流结果主要受到以下两个不确定性因素的影响：①下游电站试验性蓄水的影响。白河水文站水情变化较大，H-ADCP安装高程并非一直保持在较为理想的位置。仪器使用时需要根据流速剖面分布[27]，合理选择安装高程，以更好地监测代表流速。②应用LASSO回归模型时，需要确定算法本身的相关参数，参数优选范围以及优化方法对模型计算结果可能有一定的影响。

3 结 语

本文分析了不同水情下断面平均流速与H-ADCP流速数据的相关关系，考虑落差、仪器入水深等因素建立最小二乘回归模型，解决受水利工程影响测站的推流问题；引入代表性较好的单个流速网格，利用LASSO回归模型进行参数估计，解决小流量下推流问题，构建了白河水文站2020年8月10日至12月31日H-ADCP推流方案。基于本文方案推算的流量成果系统误差、随机不确定度均满足SL247—2012《水文资料整编规范》精度要求，并通过符号检验、适线检验和偏离数值检验；从日水量和年径流量来看，推算流量满足上下游水量平衡，有较好应用效果。本文所构建的推流方案可为白河水文站及类似受水利工程影响测站的H-ADCP推流方案构建提供参考。

H-ADCP流速代表性以及对应的实测流量成果精度直接影响推流关系率定，决定推流成果质量。在水利工程影响下，测流断面选择、仪器参数设置以及安装位置等均会影响数据收集成果，进而影响分析成果。未来需规范比测流程，依据定点式ADCP相关测量规范，确定H-ADCP安装高程和有效工作范围，保证数据采集质量，收集不同水情下数据进行分析对照，更好地控制H-ADCP流速代表性。同时尝试采用新方法应对水利工程等影响下的复杂水情，得到推流精度高且稳健的方案。