用最小二乘法拟合水位流量关系线的方法

(四川省岷江水文水资源局彭山水文站，四川 彭山，620800)

流量是过水断面上的平均流速和过水断面面积的乘积。对于河流及任何明渠水流来说，过水断面面积由水位决定。水位可随时准确测记，但过水断面上各点流速不一样，要实时测得流量仍有一定困难，还难免有误差。日常工作中多采用在一年中的不同时期或时段，在各级水位实测流量，点绘水位、流量关系图，通过点群重心建立全年或多个时段的水位流量关系线，以水位推求任何时刻的流量的方法。此方法采取目估并用专门的曲线板在图上画出关系线，再从线上查读一些点次做成水位流量关系节点表，用于推求需要时刻的流量。但此项工作不仅量大且易受人为因素影响，因此，利用某种相关函数式表达水位流量关系线，以水位为自变量推定任何时刻的流量就显得尤为必要。而高等数学中的最小二乘法就是求解最优数值和相关函数式的可靠又方便的方法。

1 最小二乘法概述

1.1 最小二乘法的原理

最小二乘法(又称最小平方法)是一种数学优化技术，它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。

从几何意义上讲，就是寻求与给定点(xi，yi)(i=0，1，…，m)的距离平方和为最小的曲线y=p(x)。函数p(x)称为拟合函数或最小二乘解，求拟合函数p(x)的方法称为曲线拟合的最小二乘法[1]。

1.2 利用最小二乘法进行曲线拟合的基本方法

利用最小二乘法进行多项式拟合的一般方法可归纳为以下几步：

(1)由已知数据绘制出函数粗略的图形——散点图，确定拟合多项式的次数n。建议用2次或3次多项式，即n=2或n=3。多数以n=3为好，个别n=2时也很好，不妨都试一下，选用最优者。

(2)根据最小二乘法原理确定函数式各项系数，如当函数为n次多项式，即Y=A0+A1X+A2X2+…+AnXn，则使下式

∑[Yi-(A0+A1Xi+A2Xi2+……+AnXin〗2=Min

(1)

(3)对式(1)中的A0、A1、A2、A3…An分别求偏导数并令其等于0，列出方程组求解出A0、A1、A2、A3…An等系数值。

(4)写出拟合多项式。

2 最小二乘法在水位流量关系式推定中的应用

本文中所用水文资料来源于四川省岷江水文水资源局彭山水文站。

2.1 最小二乘法的曲线拟合

根据彭山站2018年1线水文测量资料(表1)绘制散点图，并根据散点图的最佳拟合效果(图1)，确定拟合多项式的次数n=3。

表1 彭山站2018年实测流量的1线资料、最小二乘法定线及误差

根据式(1)，当n=3时，求解系数A0、A1、A2、A3的方程组为式(2)

(2)

图1 彭山站2018年实测流量散点及曲线拟合

上述拟合计算可用Excel(电子表格)软件完成。经Excel拟合后图1的关系式为：

Y=404.34X3-323.06X2+2536X-2023(R2=0.9987)

(3)

根据式(3)计算出各施测号数基本水尺水位对应的流量(表1中的机定流量)，经误差(表1中的机线误差)分析可知，其最大相对误差绝对值为4.2，而人工定线误差绝对值为4.9，机线误差总体上小于人工定线误差，即最小二乘法定线优于人工定线。

用Excel拟合水位流量关系式的具体方法如下：

(1)在表格中输入用于定线的实测点的施测号数、基本水尺水位和流量(任何一次基本水尺水位的整米数值都不能省略)；

(2)选出基本水尺水位最低的水位数值，将其整米数值记为Z0，再将基本水尺水位中各水位值减去Z0，这一列的数值后面用于计算流量，叫它计算水位(见表1)。因基本水尺水位规定取至0.01m，一般有4～6位数字。如直接用此数值，当n=3时，就会计算到此数值的6次方，数值的位数太多，不利于计算，因此，须减一个我们称它为Z0的数值；

(3)将计算水位数值和流量数值作为做图数据，在“图表”中选“散点图”，将出现以计算水位和实测流量为数据源的点群关系图。此图的横坐标是水位，纵坐标是流量，与水文部门点绘的水位流量关系图相反。这是由于数学上一般采用横轴为自变量，纵轴为函数值，其不影响最后结果；

(4)右击图上任一关系点，点击“増添趋势线”，再选“多项式”及“次数”后，点击“选项”栏，确定显示“公式”和“R平方值”，最后点“确定”。这时图上方将显示水位流量关系式及R平方值。式中Y为Q，X是计算流量用的水位数值，称为计算水位，由实时水位值减去Z0所得；

(5)将水位流量关系式代入各实测流量计算水位值即可计算相应水位在线上的流量值，并可算出实测流量与线上流量的误差；

(6)如认为所得函数式符合要求，即可将推流时段内水位值减去Z0后，代入式中计算其相应流量。亦可用函数式制作水位流量关系节点表，供现行资料整编的推流用。

2.2 最小二乘法的改进应用

一般情况下，用Excel软件拟合的水位流量关系线可完全满足推流要求。但若水位变幅大，时段内最大最小流量相差几十倍甚至上百倍，如2011年彭山站当年最大实测流量为7850m3/s，最小流量仅为353m3/s，两者相差达22倍(表2)。而实际工作中一年只定一条线，这时用Excel定出的关系线就不太理想。

以彭山站2011年3线水文测量资料为例，绘制彭山站2011年3线水文测量资料散点图，并根据散点图的最佳拟合效果，确定拟合多项式的次数n=3，见图2。经Excel拟合后图2的关系式为：

表2 彭山站2011年站实测流量的3线资料、最小二乘法定线及误差

Y=-20.231X3+340.96X2+116.16X+72.861(R2=0.9989)

(4)

根据上式(4)计算出各施测号数基本水尺水位对应的流量(表2中的机定流量)，经误差(表2中的机线误差)分析可知，其绝对误差均<5%，总体上小于人工定线误差。

图2 彭山站2011年3线水位流量散点及曲线拟合

但由式(1)可知，求得的函数式的各系数(A0、A1、A2、A3…An)是使实测值与函数值之绝对误差的平方和达最小。在水位流量关系线的定线实践中，低水部分关系点距线很近但相对误差却很大甚至不合格；高水部分的关系点有些看似离线远但相对误差并不大，显然由式(1)确定的函数式对高水部分的关系点考虑多些。如果以下式

(5)

就是使实测值与函数值的相对误差之平方和为最小来求得函数式各系数。但n=3时，上式中的A0、A1、A2、A3等由解下列方程组求得：

(6)

以彭山站2011年3线水文测量数据为例，其方程组为：

(7)

注：为方便计算，式(7)方程组中等式两边均乘以106。

解得函数式为：

Y=48.46+154.09X+324.40X2-18.32X3(R2=0.9932)

(8)

据上式(8)求得各点线上流量(表2中的相对误差定线)及与实测流量对比，其误差(表2中的相对误差定线误差)总体上小于人工定线误差和用式(1)拟合的定线误差(表2中的机定误差)，相对误差均<5%，优于人工定线和用式(1)拟合的定线。高原、山区的水文站多数年水位变幅不太大，过水断面稳定，河道坡降也大，往往一年只定一条水位流量关系线，此方法更适用这类测站。

3 结论与讨论

(1)利用最小二乘法原理的水文基础测量的水位流量推定，流量测验精度高，计算可利用Excel软件，简便易行，可在实际工作中推广应用；

(2)通常以多项式函数表达水位流量关系式。日常工作中利用最小二乘法计算得出的水位流量关系函数式一般取用3次多项式。个别关系线也有2次多项式比3次多项式还稍好，建议用Excel软件时可分别用2次多项式和3次多项式求出结果进行比较后选定函数式。3次以上的多项式不宜采用；

(3)在使实测值与函数值的相对误差之平方和为最小来求得函数式各系数的最小二乘法方法中，解得的A0、A1、A2、A3等系数建议取至3位小数。此方法更适用于年水位变幅不太大，过水断面稳定，河道坡降也大的高原、山区水文站；

(4)若在拟定水位流量关系线时测点数较少且水位变幅不大，可试选幂函数或指数函数拟合曲线；

(5)为了避免计算水位出现负值，前述的Z0建议选定为推流时段中最低水位的整米数值。