1. 引言
水文现代化是新阶段水利高质量发展的基础性先行性工作,其主要特征之一就是水文要素监测自动化,要求广泛应用现代技术手段和仪器设备,全面推进水文要素全量程自动监测。河流流量是重要的水文要素之一,推进河流流量全量程自动监测,对水旱灾害防御、水资源管理、水环境保护和水生态修复具有重要意义。
河流流量自动监测方法有雷达法、影像法、声学多普勒法、超声波时差法、水工建筑物法、水位流量关系单值化法等[1] [2]。河流流量自动监测的研究,多为单一的流量自动监测方法应用[3]-[10];少数学者探索了组合式流量自动监测系统的应用[11] [12];多源流量监测数据的融合,多见用于洪水预报模型研究[13]-[16],鲜有用于流量自动监测研究[17]。
由于气候变化和人类活动的影响,以及水文监测技术发展限制,在河流小流速和高含沙量洪水的条件下实现流量自动监测面临较大挑战,多数水文站难以仅依靠单一方法实现河流流量全量程自动监测,通常需要结合两种或多种方法以确保监测的准确性;此外,流量自动监测设备偶尔会发生故障,导致数据中断,因此必须具备辅助方案以保障数据的连续性。对于采用多种监测方法的站点,还需妥善处理不同方法之间的数据衔接问题,为流量资料整编和报汛报旱提供可靠支持。基于此现状,本文提出了一种基于多源数据融合算法来实现河流流量全量程自动监测,并以西江梧州水文站作为典型案例验证了该算法的可行性。
2. 多源数据融合方法
2.1. 基本思路
假定某水文站有多种流量自动监测方法,不同方法的流量自动监测数据相互独立、互不影响。各种方法均通过与实测流量建立相关关系进行率定分析,满足符号检验、适线检验、偏离检验和系统误差、随机不确定度要求,确定了各自不同的适用范围,适用范围覆盖了从最低水位到最高水位的全量程水位变幅。
在此基础上,采用基于动态方差的卡尔曼滤波流量融合算法,对多种自动监测方法提供的实时数据进行融合,通过建立“范围适用、数值无缺、方差最小”等指标体系,实现对河流流量监测数据的实时评估与校正,从而选取最优的实时数据,生成更完整、稳定和准确的融合流量数据。
下一步则建立实时融合流量与实测流量的相关关系,并经过关系曲线的符号检验、适线检验、偏离检验后,计算相关关系的相关系数、系统误差、随机不确定度,若优于单一流量自动监测方法的相关指标,且满足《水文资料整编规范》(SL/T 247-2020)规定的精度要求,则证明融合后的流量自动监测数据为最优,可作为该水文站流量全量程自动监测数据。
2.2. 基于动态方差的卡尔曼滤波流量融合算法
2.2.1. 卡尔曼滤波介绍
卡尔曼滤波是一种利用线性系统状态方程,通过系统输入输出观测数据,对系统状态进行最优估计的算法。假设状态预测值为
,状态转移矩阵为
,控制输入矩阵为
,控制向量为
,则在第k时刻的状态预测公式为:
(1)
假设状态协方差矩阵为
,预测方差为
,则在第k时刻的方差更新公式为:
(2)
假设观测值为
,观测方差为
,测量矩阵为
,则卡尔曼增益
为:
(3)
第k时刻状态的最优估计为:
(4)
第k时刻状态协方差矩阵更新为:
(5)
卡尔曼滤波通过观测值与上述5个公式不断迭代更新误差项,使得融合结果误差最低,最终得到最优融合结果。
2.2.2. 卡尔曼滤波在流量计算中的应用
以不同测流方法对同一水体的测量结果作为卡尔曼滤波计算中的观测值
,测量误差作为观测方差
。由于观测值不止一个,卡尔曼滤波公式需要转变为多观测值的卡尔曼滤波公式:
(6)
(7)
(8)
其中,n为观测值数量,i为第i个观测量。
本文选用趋势预测法作为流量预测模型,因此状态转移矩阵
为单位阵,控制项
等于上一时刻与上上时刻流量值之差
,由于观测值和预测值都是流量值,测量矩阵
也是单位阵。所以卡尔曼滤波的5个公式经过简化变为:
(9)
(10)
(11)
(12)
(13)
2.2.3. 方差动态调整
原卡尔曼滤波算法中,观测误差
为一个定值,本文采用方差动态调整方法,在实时计算中根据收集到的数据情况动态调整观测方差
和预测方差
,使得算法能适应水体环境的不断变化,计算出最优融合流量。
实时计算得到各个观测值
,其均值为
、方差为
,预测值为
,则理论上,观测均值
越接近预测值
、观测方差
越小,最后融合得到的流量结果越准确、稳定。
当观测方差不变,观测均值与预测值偏差变大,我们认为是预测值不准确,更偏向于相信观测值,算法上增大了预测方差
的值;当观测均值与预测值偏差不变,观测方差变大,我们认为是观测值稳定性下降,更偏向于相信预测值,算法上增大了观测方差
的值,见图1。
图1. 方差调整原理
为了保证数据变化的连续性,防止异常数据导致方差突变,造成计算结果有较大的跳动,本文在预测方差
和观测方差
的调整中添加了渐变系数
和
,更新公式变为:
(14)
(15)
2.3. 融合步骤
算法的完整流程如下:
1、数据预处理,利用上下限过滤方法对各测流方法结果去除粗大误差。
2、计算各观测值的权重
,可以根据观测值变化曲线的抖动程度计算得到,也可以在有比测数据的前提下,根据观测值与比测值的偏差计算得到。观测值权重作为对不同观测值质量的先验知识,可以应用于后续计算过程。若不作权重分配,默认等分权重即可;
3、初始化参数,状态协方差
、预测方差
,初始化上一时刻流量值
等于所有观测值的均值(加权平均值)
、流量变化趋势
;
4、计算流量预测值:
(16)
5、计算观测均值:
(17)
6、计算预测方差:
(18)
7、计算观测方差:
(19)
8、根据计算公式(9)~(13)计算出最优融合流量
和状态方差
;
9、更新流量变化趋势,用于下一次运算:
(20)
10、重复步骤4~9,实现实时的流量融合。
2.4. 缺值处理
由于水体环境变化、通讯中断、设备故障等因素,测流结果可能会出现数据缺失。常规卡尔曼滤波算法无法处理缺值数据,因此本小节将对算法进行针对数据缺值的优化调整,使其能应对缺值情况,保证算法稳定。
首先,存在缺值数据时,各观测值的权重将要重新计算,即去掉缺值对应的权重,其他权重归一化。第i种非缺值测量结果对应的权重变为:
(21)
其中,D是非缺值测流方法集合,j是D中每一种测流方法。
然后,在卡尔曼滤波算法中,卡尔曼增益、最优估计值和状态协方差矩阵的计算都需要作相对应的调整:
(22)
(23)
(24)
其中,D是非缺值集合。
修改后,算法将能应对部分测量数据缺失的情况,并能保证最终结果受到的影响最小。但若是所有测量数据都缺失,本算法也无法得出计算结果,只能使用后处理方法进行插补处理。
2.5. 方法筛选
不同的测流方法在不同的水体环境下有不同的表现,当某种测流方法在当前环境下表现不佳,我们应该考虑过滤掉该方法得到的流量结果,不参与流量融合直到其测量效果恢复正常为止。
判断一种测流方法的测流效果是否正常,可以通过先验筛选、自检验筛选和互检验筛选三种方法来判断:
1、先验筛选是指基于以往大量的实测数据,为每种测流方法设置特定的筛选条件。例如当某种方法在特定时段内受环境或设备因素影响导致测流效果不佳时,可以暂时停用该方法的流量,直到其恢复稳定。
2、自检验筛选是指每种测流方法用当前的测量结果与历史数据进行对比,以识别是否存在显著的突变或异常值,从而对测量结果进行筛选。
3、互检验筛选是通过比较不同测流方法之间的测量结果,若某方法的结果与其他方法相比存在明显差异,或流量涨率与其它方法不一致,则可以认为该方法可能存在误差,需对结果进行过滤处理。
2.5.1. 先验筛选
先验筛选需要收集过往大量实测流量数据,通过与比测数据对比,计算出测量误差大小,再与日期时间、水位高低、流速大小等几个影响因素作相关性分析,以确定各测流方法的主要影响因素,从而设定合理的过滤范围。
在后续的实例分析章节中,本文将针对每种测流方法详细分析,并提供相应的先验筛选条件。
2.5.2. 自检验筛选
自检验筛选需要对比当前测流结果与过去结果。假设测流结果为f,误差阈值为t,该流域的流量上下限分别为
、
,则当第i时刻的测流结果与前m时刻的测流结果均值相比,误差小于t时,被判定为正常结果,否则被判定为异常结果,需要被过滤掉:
(25)
2.5.3. 互检验筛选
互检验筛选需要对比每一种测流方法与其他方法之间的差异。我们可以通过Hampel滤波算法对每种方法的测流结果进行离群点判断和过滤。假设各方法的测流结果为f,若符合以下条件,则可以认为该方法的测流结果异常,需将其过滤掉:
(26)
其中,Mdn为求中位数,s为标准差。
通过以上检验和筛选方法,我们可以实时判断出每种测流方法取得的测流结果是否应用于最终流量融合的计算中,从而确定流量融合所需的方法组合。
3. 应用实例分析
3.1. 典型水文站概况
梧州水文站是珠江流域西江水系西江干流的重要控制站,属于国家重要水文测站,集水面积327,006 km2,
表1. 梧州水文站流量自动监测方法适用范围情况表
设备或方法 |
单一方法适用范围 |
单一方法精度情况 |
水位(m) |
流速(m/s) |
流量(m3/s) |
相关系数 |
系统误差 |
随机不确定度(%) |
侧扫雷达 |
3.02~24.10 |
0.30~2.16 |
1980~41,000 |
0.988 |
5.70 |
11.4 |
水平式ADCP |
1.34~22.27 |
0.19~1.88 |
1180~34,000 |
0.996 |
0.18 |
3.46 |
落差指数法1 (梧州–桂粤) |
2.10~24.10 |
0.24~2.16 |
1440~41,000 |
0.999 |
−0.73 |
10.2 |
落差指数法2 (梧州–龙圩) |
3.00~24.10 |
0.30~2.16 |
2490~41,000 |
0.999 |
−0.03 |
9.20 |
图2. 各测流方法得到的流量过程线
图3. 各测流方法与走航式ADCP结果对比情况
占西江流域集水面积的94.6%,控制了广西境内85%的来水。梧州水文站实测最高水位为27.80 m (1985国家高程基准),最大流量为54,500 m3/s,最大平均流速3.39 m/s;实测最低水位为1.33 m,最小流量为656 m3/s,最小平均流速0.12 m/s;历年最高水位与最低水位之差为26.47 m。
梧州水文站流量自动监测方法有侧扫雷达自动测流法、水平式ADCP自动测流法、落差指数自动推流法,经率定分析和符号检验、适线检验、偏离数值检验,其适用范围均为部分水位(流速、流量)量程范围,见表1。
3.2. 各测流方法情况分析
本文选取梧州水文站2024年的侧扫雷达自动测流法、水平式ADCP自动测流法和两种落差指数自动推流法得到的流量自动监测数据,覆盖水位1.34 m~24.1 m、流速0.19 m/s~2.16 m/s、流量1180 m3/s~41,000 m3/s,约占历年最大流量变幅的73%,具有较好代表性;与走航式ADCP的实测流量进行对比分析,相关系数均达到0.997以上,见图2、图3。
图4. 侧扫雷达流量过程线
图5. 侧扫雷达流量误差分析(水位、流速、时间)
从整体过程线与水位、流速、时间三个方面因素,分别对每种测流方法进行误差分析如下。
1、侧扫雷达自动测流法:测流结果缺值比较多,在低流速、低水位时跳变较多,在高水位时也有一定测量误差。由于发生跳变时误差都较大,在融合之前通过互检验筛选过滤掉错误数据,见图4、图5。
图6. 水平式ADCP流量过程线
图7. 水平式ADCP流量误差分析(水位、流速、时间)
图8. 落差指数自动推流法1流量过程线
图9. 落差指数自动推流法1流量误差分析(水位、流速、时间)
图10. 落差指数自动推流法2流量过程线
图11. 落差指数自动推流法2流量误差分析(水位、流速、时间)
总体来看,四种方法的测流结果与实测流量均表现出较高的相关性,但在特定水位级别上仍存在相对较大的偏差,在一定程度上影响了数据的准确性,因此这些方法均难以作为流量全量程自动监测成果。
3.3. 融合结果比较分析
3.3.1. 融合过程与融合结果
通过分析各测流方法率定时的实时流量与实测流量的均方根误差,取其倒数归一化后,得到各测流方法的权重,分别为0.21、0.50、0.13和0.16。这些权重将作为融合算法中初始权重值使用。
根据文中提到的三种检验筛选方法(先验检验、自检验和互检验),对各测流方法得出的实时流量结果进行了筛选和处理,作为融合算法的源数据。随后,依据融合算法流程对实时流量数据进行融合计算,得到最终的融合流量。最后,根据计算结果动态更新算法模型参数,为下一时刻的计算做准备,直至所有数据处理完毕,算法终止。
图12. 融合流量过程线
图12是最终融合结果的流量过程线与实测流量值对比。可以看出,融合流量过程线相较于单一监测方法的过程线更加平滑、完整,消除了明显的跳动和错误数据,且与实测点高度吻合,表明融合后的提升效果显著。
融合算法在保持水平式ADCP结果的准确性前提下,利用低水位的落差指数2结果、中水位的侧扫雷达结果和高水位的落差指数1结果,进行流量值的去噪、平滑,使得最终的融合流量结果在准确性、稳定性、完整性上都有很好的表现。
3.3.2. 精度比较分析
将融合流量数据与所有单一监测方法结果进行对比,对比结果见表2。其中,我们使用了流量过程线所有点的三阶导均值作为评价过程线平滑程度的指标,该指标值越小,流量过程线越平滑。具体计算公式如下,其中f为流量值,n为样本总数。
(27)
表2. 融合流量结果与各监测方法结果对比分析表
测流方法 |
完整率(%) |
相关系数 |
系统误差(%) |
随机不确定度(%) |
平滑程度 |
侧扫雷达 |
87.7 |
0.9976 |
−2.18 |
19.17 |
8535 |
水平式ADCP |
99.5 |
0.9996 |
−0.68 |
11.86 |
986,212 |
落差指数法1 (梧州–桂粤) |
99.8 |
0.9986 |
−4.79 |
19.34 |
43,126 |
落差指数法2 (梧州–龙圩) |
99.8 |
0.9985 |
1.75 |
16.59 |
33,147 |
融合法 |
100 |
0.9994 |
0.85 |
10.92 |
93 |
将融合流量与实测流量建立相关关系,对比点分布及精度比较分析结果见图13、表3。
图13. 融合流量与走航式ADCP对比结果
表3. 2024年梧州水文站融合法流量与实测流量相关关系分析表
序号 |
时间 |
水位(m) |
实测流量(m3/s) |
融合流量(m3/s) |
1 |
2024/1/2 16:04 |
2.54 |
1910 |
2090 |
2 |
2024/1/5 10:14 |
2.12 |
1770 |
1700 |
3 |
2024/1/10 11:00 |
2.59 |
2090 |
2220 |
… |
… |
… |
… |
… |
20 |
2024/06/20 09:34 |
23.34 |
39,300 |
38,600 |
21 |
2024/06/21 09:06 |
24.08 |
41,000 |
41,000 |
22 |
2024/06/21 12:34 |
24.10 |
40,600 |
41,000 |
… |
… |
… |
… |
… |
61 |
2024/12/16 16:47 |
2.81 |
2850 |
2610 |
62 |
2024/12/20 10:01 |
2.86 |
1970 |
2000 |
63 |
2024/12/30 12:03 |
2.50 |
2170 |
2220 |
样本容量: |
N = 63 |
正号个数:33 |
符号交换次数:27 |
符号检验: |
u = 0.25 |
允许:1.15 (显著性水平a = 0.25) |
合格 |
适线检验: |
U = 0.89 |
允许:1.28 (显著性水平a = 0.10) |
合格 |
偏离数值检验: |
|t| = 1.26 |
允许:1.30 (显著性水平a = 0.20) |
合格 |
标准差: |
Se (%) = 5.46 |
随机不确定度(%):10.92 |
系统误差(%):0.85 |
相关系数: |
R = 0.9994 |
数据完整率(%):100 |
|
融合后的梧州水文站最小流量为1440 m3/s、最大流量41,700 m3/s,完整覆盖了实测流量范围1440 m3/s~41,000 m3/s,这一改进实现了流量全量程自动监测,确保了数据的连续性,通过了符号检验、适线检验和偏离数值检验。融合方法不仅解决了单一方法的缺值问题,还在保证相关系数、系统误差、随机不确定度等误差指标基本优于所有单一监测方法的前提下,显著提升了数据的平滑程度,有效减少了跳动误差对结果的影响。融合后的流量自动监测数据质量明显优于单一自动监测方法。
3.3.3. 特征值比较分析
利用融合法得到的梧州水文站各月平均流量与资料整编成果特征值进行比较分析,各月的平均流量相对误差位于−2.0%~2.1%之间,均在允许相对误差±3%以内,符合要求,结果见表4。
表4. 2024年梧州水文站融合法流量与整编成果特征值比较分析表
比较项目 |
一月 |
二月 |
三月 |
四月 |
五月 |
六月 |
七月 |
八月 |
九月 |
十月 |
十一月 |
十二月 |
融合法 |
2260 |
2500 |
2860 |
7290 |
10,400 |
24,900 |
12,900 |
9930 |
7500 |
3410 |
2520 |
2180 |
整编成果 |
2260 |
2530 |
2840 |
7140 |
10,200 |
24,400 |
12,800 |
9780 |
7640 |
3450 |
2570 |
2210 |
相对误差(%) |
0.0 |
−1.2 |
0.7 |
2.1 |
2.0 |
2.1 |
0.8 |
1.5 |
−1.8 |
−1.2 |
−2.0 |
−1.4 |
允许误差(%) |
±3.0 |
±3.0 |
±3.0 |
±3.0 |
±3.0 |
±3.0 |
±3.0 |
±3.0 |
±3.0 |
±3.0 |
±3.0 |
±3.0 |
是否合格 |
合格 |
合格 |
合格 |
合格 |
合格 |
合格 |
合格 |
合格 |
合格 |
合格 |
合格 |
合格 |
3.3.4. 径流总量比较分析
利用融合法得到的梧州水文站年径流总量、汛期径流总量、次洪水径流总量与资料整编成果进行比较分析,年径流总量相对误差在±3%以内、汛期径流总量相对误差在±3.5%以内、次洪水径流总量相对误差在±6%以内,符合要求,结果见表5。
表5. 2024年梧州水文站融合法流量与整编成果径流总量比较分析表
比较项目 |
融合法 |
整编成果 |
相对误差(%) |
允许相对误差(%) |
是否合格 |
年径流总量(108 m3) |
2340 |
2313 |
1.17 |
±3.0 |
合格 |
汛期径流总量(108 m3) |
1920 |
1894 |
1.37 |
±3.5 |
合格 |
次洪水径流总量(108 m3) |
303.7 |
297.2 |
2.18 |
±6.0 |
合格 |
注:次洪水径流总量为6月15日至25日的径流总量。
3.3.5. 总体评价
经比较分析,利用融合法得到的梧州水文站流量自动监测数据与实测流量相关系数达到0.9994,系统误差为0.85%,随机不确定度10.92%,符号检验、适线检验和偏离数值检验合格,测验精度优于单一的自动监测方法;各月平均流量、年径流总量、汛期径流总量、次洪水径流总量与资料整编成果相比,相对误差均符合《水文资料整编规范》(SL/T 247-2020)规定的精度要求。
3.3.6. 适用性分析
利用郁江南宁水文站、红水河都安水文站、黔江大藤峡水文站的2024年数据进行适用性测试,结果见表6。通过符号检验、适线检验、偏离数值检验和完整率、相关系数、系统误差、随机不确定度、平滑程度等指标分析,结果表明,融合法结果均通过检验,总体上优于单一方法,达到了流量全量程自动监测和提升结果完整性、稳定性、准确性的目的,说明该方法可适用于其他河流流域和不同水文条件的水文站。
表6. 融合法在不同水文条件下的适用性分析表
水文站名称 |
测流 方法 |
符号检验
(<1.15) |
适线检验
(<1.28) |
偏离数值检验
(<1.30) |
完整率(%) |
相关 系数 |
系统误差(%) |
随机不确定度(%) |
平滑程度 |
南宁水文站 |
声层析 时差法 |
通过(0.14) |
不通过(1.43) |
通过(0.54) |
57.09 |
0.9993 |
−0.72 |
13.04 |
2371 |
水平式 ADCP |
通过(0.47) |
通过(−0.71) |
通过(0.70) |
98.54 |
0.9985 |
1.21 |
13.32 |
52,036 |
融合法 |
通过(0.00) |
通过(0.81) |
通过(0.73) |
99.16 |
0.9992 |
1.02 |
12.37 |
35 |
都安水文站 |
水平式 ADCP |
通过(0.00) |
通过(−2.21) |
通过(0.20) |
94.09 |
0.9995 |
0.29 |
3.73 |
24,817 |
单一线法 |
通过(−0.27) |
通过(−0.55) |
通过(1.01) |
98.62 |
0.9985 |
0.69 |
6.45 |
402 |
融合法 |
通过(0.80) |
通过(0.55) |
通过(0.90) |
99.57 |
0.9996 |
0.81 |
3.69 |
11.18 |
大藤峡水文站 |
落差 指数法 |
通过(0.75) |
通过(−1.52) |
通过(0.55) |
89.74 |
0.9974 |
−0.62 |
15.06 |
173,642 |
定点 雷达波 |
通过(1.06) |
不通过(2.74) |
不通过(2.07) |
90.02 |
0.9960 |
5.06 |
34.44 |
460,809 |
融合法 |
通过(0.00) |
通过(0.44) |
不通过(1.38) |
97.09 |
0.9989 |
1.29 |
13.20 |
318 |
4. 结语
1、本文提出一种基于多源数据融合的河流流量全量程自动监测方法,该方法以动态方差的卡尔曼滤波算法为核心,辅以缺值处理和方法筛选等技术。以率定后的多种流量自动监测数据为分析对象,利用监测数据与实测流量的均方根误差建立融合算法的各种方法初始权重值,实时分析不同时间、不同水位、不同流速下各方法监测数据的优劣,选出不同阶段下的最优监测数据组合,得到完整性、稳定性、准确性最佳的流量结果作为该水文站全量程在线自动监测数据。
2、利用西江梧州水文站2024年数据的4种自动监测方法,对基于多源数据融合的河流流量全量程自动监测方法进行了验证,结果表明,相较于单一方法,融合方法在数据完整率、相关系数、系统误差、随机不确定度和平滑程度等多个关键指标上表现更为优异;以整编资料的计算结果为对比对象,融合方法在各月平均流量、年径流总量、汛期径流总量、次洪水径流总量的相对误差均符合《水文资料整编规范》(SL/T 247-2020)规定的精度要求;此外,融合方法确保了水文站在不同时间、水位和流速条件下,均能准确且稳定地获取流量数据,成功实现了河流流量全量程自动监测。
3、提出的基于多源数据融合的河流流量全量程自动监测方法,相较于单一监测方法或简单融合方法,优势在于:不同监测方法之间可以互相监督、检验,并且能识别出误差大、跳动多的方法,大大降低奇异值对融合结果的误差影响;系统处于动态调整状态,而非固定参数,能有效应对环境突变、设备故障等突发情况,最大限度地减少误差,增强系统稳定性;若后续增加实测数据,可以持续优化更新模型参数,使模型适应不断变化的河流环境;当水文站接入新的流量监测方法,通过调整系统参数即可融合新方法,无需重构系统。
4、提出的基于多源数据融合的河流流量全量程自动监测方法还存在一定局限性:融合算法结果受各测流方法结果优劣的影响,若所有测流方法得到的结果都不太准确,融合算法也无法得到一个比较可靠的结果;融合算法结果的精度符合性只是事后分析验证,还未能引进实测流量对自动监测数据进行实时验证。在下阶段研究中,可以确保各测流方法结果相对较好的基础上,引进实测流量对自动监测数据进行实时验证,从而保障融合算法结果的精度符合性。
基金项目
感谢广西壮族自治区水利厅水利科技推广及“工程带科研”项目SK2021-3-15对本研究的支持。
NOTES
作者简介:文宏展(1969-),男,广西玉林人,研究生,高级工程师,主要从事水文水资源监测、评价与管理工作,Email: 609579263@qq.com