用最小二乘法擬合水位流量關(guān)系線的方法

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

流量是過(guò)水?dāng)嗝嫔系钠骄魉俸瓦^(guò)水?dāng)嗝婷娣e的乘積。對(duì)于河流及任何明渠水流來(lái)說(shuō)，過(guò)水?dāng)嗝婷娣e由水位決定。水位可隨時(shí)準(zhǔn)確測(cè)記，但過(guò)水?dāng)嗝嫔细鼽c(diǎn)流速不一樣，要實(shí)時(shí)測(cè)得流量仍有一定困難，還難免有誤差。日常工作中多采用在一年中的不同時(shí)期或時(shí)段，在各級(jí)水位實(shí)測(cè)流量，點(diǎn)繪水位、流量關(guān)系圖，通過(guò)點(diǎn)群重心建立全年或多個(gè)時(shí)段的水位流量關(guān)系線，以水位推求任何時(shí)刻的流量的方法。此方法采取目估并用專門的曲線板在圖上畫出關(guān)系線，再?gòu)木€上查讀一些點(diǎn)次做成水位流量關(guān)系節(jié)點(diǎn)表，用于推求需要時(shí)刻的流量。但此項(xiàng)工作不僅量大且易受人為因素影響，因此，利用某種相關(guān)函數(shù)式表達(dá)水位流量關(guān)系線，以水位為自變量推定任何時(shí)刻的流量就顯得尤為必要。而高等數(shù)學(xué)中的最小二乘法就是求解最優(yōu)數(shù)值和相關(guān)函數(shù)式的可靠又方便的方法。

1 最小二乘法概述

1.1 最小二乘法的原理

最小二乘法(又稱最小平方法)是一種數(shù)學(xué)優(yōu)化技術(shù)，它通過(guò)最小化誤差的平方和尋找數(shù)據(jù)的最佳函數(shù)匹配。利用最小二乘法可以簡(jiǎn)便地求得未知的數(shù)據(jù)，并使得這些求得的數(shù)據(jù)與實(shí)際數(shù)據(jù)之間誤差的平方和為最小。

從幾何意義上講，就是尋求與給定點(diǎn)(xi，yi)(i=0，1，…，m)的距離平方和為最小的曲線y=p(x)。函數(shù)p(x)稱為擬合函數(shù)或最小二乘解，求擬合函數(shù)p(x)的方法稱為曲線擬合的最小二乘法[1]。

1.2 利用最小二乘法進(jìn)行曲線擬合的基本方法

利用最小二乘法進(jìn)行多項(xiàng)式擬合的一般方法可歸納為以下幾步：

(1)由已知數(shù)據(jù)繪制出函數(shù)粗略的圖形——散點(diǎn)圖，確定擬合多項(xiàng)式的次數(shù)n。建議用2次或3次多項(xiàng)式，即n=2或n=3。多數(shù)以n=3為好，個(gè)別n=2時(shí)也很好，不妨都試一下，選用最優(yōu)者。

(2)根據(jù)最小二乘法原理確定函數(shù)式各項(xiàng)系數(shù)，如當(dāng)函數(shù)為n次多項(xiàng)式，即Y=A0+A1X+A2X2+…+AnXn，則使下式

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

(1)

(3)對(duì)式(1)中的A0、A1、A2、A3…An分別求偏導(dǎo)數(shù)并令其等于0，列出方程組求解出A0、A1、A2、A3…An等系數(shù)值。

(4)寫出擬合多項(xiàng)式。

2 最小二乘法在水位流量關(guān)系式推定中的應(yīng)用

本文中所用水文資料來(lái)源于四川省岷江水文水資源局彭山水文站。

2.1 最小二乘法的曲線擬合

根據(jù)彭山站2018年1線水文測(cè)量資料(表1)繪制散點(diǎn)圖，并根據(jù)散點(diǎn)圖的最佳擬合效果(圖1)，確定擬合多項(xiàng)式的次數(shù)n=3。

表1 彭山站2018年實(shí)測(cè)流量的1線資料、最小二乘法定線及誤差

根據(jù)式(1)，當(dāng)n=3時(shí)，求解系數(shù)A0、A1、A2、A3的方程組為式(2)

(2)

圖1 彭山站2018年實(shí)測(cè)流量散點(diǎn)及曲線擬合

上述擬合計(jì)算可用Excel(電子表格)軟件完成。經(jīng)Excel擬合后圖1的關(guān)系式為：

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

(3)

根據(jù)式(3)計(jì)算出各施測(cè)號(hào)數(shù)基本水尺水位對(duì)應(yīng)的流量(表1中的機(jī)定流量)，經(jīng)誤差(表1中的機(jī)線誤差)分析可知，其最大相對(duì)誤差絕對(duì)值為4.2，而人工定線誤差絕對(duì)值為4.9，機(jī)線誤差總體上小于人工定線誤差，即最小二乘法定線優(yōu)于人工定線。

用Excel擬合水位流量關(guān)系式的具體方法如下：

(1)在表格中輸入用于定線的實(shí)測(cè)點(diǎn)的施測(cè)號(hào)數(shù)、基本水尺水位和流量(任何一次基本水尺水位的整米數(shù)值都不能省略)；

(2)選出基本水尺水位最低的水位數(shù)值，將其整米數(shù)值記為Z0，再將基本水尺水位中各水位值減去Z0，這一列的數(shù)值后面用于計(jì)算流量，叫它計(jì)算水位(見(jiàn)表1)。因基本水尺水位規(guī)定取至0.01m，一般有4～6位數(shù)字。如直接用此數(shù)值，當(dāng)n=3時(shí)，就會(huì)計(jì)算到此數(shù)值的6次方，數(shù)值的位數(shù)太多，不利于計(jì)算，因此，須減一個(gè)我們稱它為Z0的數(shù)值；

(3)將計(jì)算水位數(shù)值和流量數(shù)值作為做圖數(shù)據(jù)，在“圖表”中選“散點(diǎn)圖”，將出現(xiàn)以計(jì)算水位和實(shí)測(cè)流量為數(shù)據(jù)源的點(diǎn)群關(guān)系圖。此圖的橫坐標(biāo)是水位，縱坐標(biāo)是流量，與水文部門點(diǎn)繪的水位流量關(guān)系圖相反。這是由于數(shù)學(xué)上一般采用橫軸為自變量，縱軸為函數(shù)值，其不影響最后結(jié)果；

(4)右擊圖上任一關(guān)系點(diǎn)，點(diǎn)擊“増?zhí)碲厔?shì)線”，再選“多項(xiàng)式”及“次數(shù)”后，點(diǎn)擊“選項(xiàng)”欄，確定顯示“公式”和“R平方值”，最后點(diǎn)“確定”。這時(shí)圖上方將顯示水位流量關(guān)系式及R平方值。式中Y為Q，X是計(jì)算流量用的水位數(shù)值，稱為計(jì)算水位，由實(shí)時(shí)水位值減去Z0所得；

(5)將水位流量關(guān)系式代入各實(shí)測(cè)流量計(jì)算水位值即可計(jì)算相應(yīng)水位在線上的流量值，并可算出實(shí)測(cè)流量與線上流量的誤差；

(6)如認(rèn)為所得函數(shù)式符合要求，即可將推流時(shí)段內(nèi)水位值減去Z0后，代入式中計(jì)算其相應(yīng)流量。亦可用函數(shù)式制作水位流量關(guān)系節(jié)點(diǎn)表，供現(xiàn)行資料整編的推流用。

2.2 最小二乘法的改進(jìn)應(yīng)用

一般情況下，用Excel軟件擬合的水位流量關(guān)系線可完全滿足推流要求。但若水位變幅大，時(shí)段內(nèi)最大最小流量相差幾十倍甚至上百倍，如2011年彭山站當(dāng)年最大實(shí)測(cè)流量為7850m3/s，最小流量?jī)H為353m3/s，兩者相差達(dá)22倍(表2)。而實(shí)際工作中一年只定一條線，這時(shí)用Excel定出的關(guān)系線就不太理想。

以彭山站2011年3線水文測(cè)量資料為例，繪制彭山站2011年3線水文測(cè)量資料散點(diǎn)圖，并根據(jù)散點(diǎn)圖的最佳擬合效果，確定擬合多項(xiàng)式的次數(shù)n=3，見(jiàn)圖2。經(jīng)Excel擬合后圖2的關(guān)系式為：

表2 彭山站2011年站實(shí)測(cè)流量的3線資料、最小二乘法定線及誤差

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

(4)

根據(jù)上式(4)計(jì)算出各施測(cè)號(hào)數(shù)基本水尺水位對(duì)應(yīng)的流量(表2中的機(jī)定流量)，經(jīng)誤差(表2中的機(jī)線誤差)分析可知，其絕對(duì)誤差均<5%，總體上小于人工定線誤差。

圖2 彭山站2011年3線水位流量散點(diǎn)及曲線擬合

但由式(1)可知，求得的函數(shù)式的各系數(shù)(A0、A1、A2、A3…An)是使實(shí)測(cè)值與函數(shù)值之絕對(duì)誤差的平方和達(dá)最小。在水位流量關(guān)系線的定線實(shí)踐中，低水部分關(guān)系點(diǎn)距線很近但相對(duì)誤差卻很大甚至不合格；高水部分的關(guān)系點(diǎn)有些看似離線遠(yuǎn)但相對(duì)誤差并不大，顯然由式(1)確定的函數(shù)式對(duì)高水部分的關(guān)系點(diǎn)考慮多些。如果以下式

(5)

就是使實(shí)測(cè)值與函數(shù)值的相對(duì)誤差之平方和為最小來(lái)求得函數(shù)式各系數(shù)。但n=3時(shí)，上式中的A0、A1、A2、A3等由解下列方程組求得：

(6)

以彭山站2011年3線水文測(cè)量數(shù)據(jù)為例，其方程組為：

(7)

注：為方便計(jì)算，式(7)方程組中等式兩邊均乘以106。

解得函數(shù)式為：

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

(8)

據(jù)上式(8)求得各點(diǎn)線上流量(表2中的相對(duì)誤差定線)及與實(shí)測(cè)流量對(duì)比，其誤差(表2中的相對(duì)誤差定線誤差)總體上小于人工定線誤差和用式(1)擬合的定線誤差(表2中的機(jī)定誤差)，相對(duì)誤差均<5%，優(yōu)于人工定線和用式(1)擬合的定線。高原、山區(qū)的水文站多數(shù)年水位變幅不太大，過(guò)水?dāng)嗝娣€(wěn)定，河道坡降也大，往往一年只定一條水位流量關(guān)系線，此方法更適用這類測(cè)站。

3 結(jié)論與討論

(1)利用最小二乘法原理的水文基礎(chǔ)測(cè)量的水位流量推定，流量測(cè)驗(yàn)精度高，計(jì)算可利用Excel軟件，簡(jiǎn)便易行，可在實(shí)際工作中推廣應(yīng)用；

(2)通常以多項(xiàng)式函數(shù)表達(dá)水位流量關(guān)系式。日常工作中利用最小二乘法計(jì)算得出的水位流量關(guān)系函數(shù)式一般取用3次多項(xiàng)式。個(gè)別關(guān)系線也有2次多項(xiàng)式比3次多項(xiàng)式還稍好，建議用Excel軟件時(shí)可分別用2次多項(xiàng)式和3次多項(xiàng)式求出結(jié)果進(jìn)行比較后選定函數(shù)式。3次以上的多項(xiàng)式不宜采用；

(3)在使實(shí)測(cè)值與函數(shù)值的相對(duì)誤差之平方和為最小來(lái)求得函數(shù)式各系數(shù)的最小二乘法方法中，解得的A0、A1、A2、A3等系數(shù)建議取至3位小數(shù)。此方法更適用于年水位變幅不太大，過(guò)水?dāng)嗝娣€(wěn)定，河道坡降也大的高原、山區(qū)水文站；

(4)若在擬定水位流量關(guān)系線時(shí)測(cè)點(diǎn)數(shù)較少且水位變幅不大，可試選冪函數(shù)或指數(shù)函數(shù)擬合曲線；

(5)為了避免計(jì)算水位出現(xiàn)負(fù)值，前述的Z0建議選定為推流時(shí)段中最低水位的整米數(shù)值。