CN102693215A - 幂函数型水位流量关系的一种拟合方法 - Google Patents

Abstract

幂函数型水位流量关系的一种拟合方法，属于水文水资源技术领域。根据某水文站待拟合的水位流量关系值(Z_el，Q_l)，i＝1～n，首先利用n组关系值中较大、较小的两组关系值，确定幂函数型水位流量关系式中待估参数的初值θ₍₀₎＝(β₀₍₀₎，β₍₀₎)′，或将化非线性为线性回归的传统方法求得的非线性回归系数作为待估参数β₀，β的初值；然后采用高斯-牛顿法确定回归系数β₀，β。水位流量关系是水文资料整编的重要内容之一；是水利水电工程规划设计、水文预报和水库调度的重要依据之一。

Description

幂函数型水位流量关系的一种拟合方法

一、技术领域

本发明属于水文水资源技术领域。水位流量关系是水文资料整编的重要内容之一；是水利水电工程规划设计、水文预报和水库调度的重要依据之一。 

二、背景技术

水位流量关系常选配成幂函数型^[1]： 

$Q={\mathrm{\β}}_0{(Z-Z_0)}^{\mathrm{\β}}={\mathrm{\β}}_0{Z}_e^{\mathrm{\β}}---\left(1\right)$

式中Q——流量，m³/s；Z——水位，m；Z₀——断流水位，即Q＝0时的水位，m；Ze——有效水位，m；β₀，β——待估参数。 

为叙述方便，将式(1)表示为一般形式 

y＝ax^b            (2) 

目前常用的方法^[2，3]是首先通过变量代换将式(2)化为线性模型，即对式(2)两边取对数，得 

lgy＝lga+blgx    (3) 

并令Y＝lgy，X＝lgx，a′＝lg a，b′＝b，则有线性模型： 

Y＝a′+b′X      (4) 

其次，利用线性回归方法推求式(4)的线性回归系数a′，b′，然后再由线性回归系数反求式(2)的非线性回归系数：a＝10^a′，b＝b′。这是数理统计教科书以及生产实际中的经典方法，称其为传统方法。上述方法看起来是合理的，其实不然，说明如下。 

设幂函数型因变量的残差平方和为 

线性化后的因变量Y的残差平方和为 

文献[4](张子贤.化非线性为线性的加权回归方法及其应用[J].水文，1991，(2)：46-48.)导出二者关系为： 

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(Y_i-{\hat{Y}}_i)}^2{\left(\frac{\mathrm{dy}}{\mathrm{dY}}\right)}_i^2\≈\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-{\hat{y}}_i)}^2---\left(5\right)$

式(5)表明，传统方法所求回归系数a′，b′满足 

为最小，而由此求得的非线性回归系数a，b不满足幂函数型因变量的残差平方和 

为最小。 

文献[5](周世良，尚明芳，李怡，等.拉格朗日乘子法在水位流量关系拟合中的应用[J].水文，2011，31(5)：15-17.)、文献[6](杨晓华，陆桂华，郦建强.自适应加速遗传算法及其在水位流量关系拟合中的应用[J].水文，2002，22(2)：14-18.)分别采用拉格朗日乘子法、自适应加速遗传算法拟合水位流量关系。 

本项发明利用高斯-牛顿法，拟合幂函数型水位流量关系。 

高斯-牛顿法^[7](方开泰.实用回归分析[M].北京：科学出版社，1988.168-172.)： 

设一般的非线性方程为 

y＝f(X，θ)+ε 

式中f为一般函数；X可以是单个自变量，也可以是r个自变量X＝(x₁，x₂，…，x_r)；θ为P维参数向量，即θ＝(θ₁，θ₂，…，θ_p)′；ε为随机误差项，且∈N(0，σ²)。设y和X具有n组观察值( 

…，x_rl；y_l)，i＝1～n。求非线性模型参数θ的“最小二乘”估计，就是求θ的估计值使残差平方和 

$S\left(\mathrm{\θ}\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{\ϵ}}_i^2=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{[y_i-f(X_i,\mathrm{\θ})]}^2---\left(6\right)$

为最小。 

对于非线性模型，无法直接求参数的“最小二乘”解。高斯-牛顿法将f(X，θ)在参数初值θ₍₀₎＝(θ₁₍₀₎，θ₂₍₀₎，…，θ_p(0))′处展开成只包含一次项的泰勒级数，从而使非线性模型线性化。为方便起见，用记号f_l(θ)代替f(x_l，θ)，可导出满足式(6)S(θ)为最小的参数递推公式，并写成矩阵形式为^[7]

θ_(k+1)＝θ_(k)+[J′(θ_(k))J(θ_(k))]^-1J′(θ_(k))[(y-f(θ_(k))]    (7) 

式中k为递推次数。 

$J\left({\mathrm{\θ}}_{\left(k\right)}\right)={\left[\begin{array}{cccccc}\frac{\∂f_1\left(\mathrm{\θ}\right)}{\∂{\mathrm{\θ}}_1}& \frac{\∂f_1\left(\mathrm{\θ}\right)}{\∂{\mathrm{\θ}}_2}& \·& \·& \·& \frac{\∂f_1\left(\mathrm{\θ}\right)}{\∂{\mathrm{\θ}}_p}\\ \frac{\∂f_2\left(\mathrm{\θ}\right)}{\∂{\mathrm{\θ}}_1}& \frac{\∂f_2\left(\mathrm{\θ}\right)}{\∂{\mathrm{\θ}}_2}& \·& \·& \·& \frac{\∂f_2\left(\mathrm{\θ}\right)}{\∂{\mathrm{\θ}}_p}\\ \·& \·& \·& \·& \·& \·\\ \·& \·& \·& \·& \·& \·\\ \frac{\∂f_n\left(\mathrm{\θ}\right)}{\∂{\mathrm{\θ}}_1}& \frac{\∂f_n\left(\mathrm{\θ}\right)}{\∂{\mathrm{\θ}}_2}& \·& \·& \·& \frac{\∂f_n\left(\mathrm{\θ}\right)}{\∂{\mathrm{\θ}}_p}\end{array}\right]}_{\mathrm{\θ}={\mathrm{\θ}}_{\left(k\right)}}$

y＝(y₁，y₂，…，y_n)′，f(θ_(k))＝[f₁(θ_(k))，f₂(θ_(k))，…，f_n(θ_(k))]′。 

利用式(7)从参数初值θ₍₀₎开始，一步步递推下去，直到θ_(k)收敛稳定，即|θ_(k+1)-θ_(k)|小于或等于预先指定的小正数δ，从而得到θ的估计值。 

三、发明内容

为了克服传统方法拟合幂函数型水位流量关系，所求拟合参数不满足因变量即流量的残差平方和为最小，所得回归线并非最佳拟合曲线的缺陷，本发明采用高斯-牛顿法，拟合幂函数型水位流量关系。 

本发明采用的技术方案是：根据水位流量关系值(Z_el，Q_l)，i＝1～n，首先确定式(1)中待估参数的初值，并记θ₍₀₎＝(β₀₍₀₎，β₍₀₎)′，然后采用高斯-牛顿法确定式(1)中的参数β₀，β。 

本项发明的有益效果是，所求幂函数型水位流量关系的拟合参数满足流量的残差平方和为最小，所得水位流量关系曲线为最佳拟合曲线，且计算方法易于实现。本项发明已利用实测水位流量资料进行了验证，效果很好。 

四、具体实施方式

已知某水文站待拟合的水位流量关系值(Z_el，Q_l)，i＝1～n。 

1.确定式(1)中待估参数β₀，β的初值。在(Z_el，Q_l)，i＝1～n中，选数值较大、较小的两组关系值，建立方程组： 

$\left\{\begin{array}{c}{Q}_1={\mathrm{\β}}_0{Z}_{e1}^{\mathrm{\β}}\\ Q_2={\mathrm{\β}}_0{Z}_{e2}^{\mathrm{\β}}\end{array}\right.---\left(8\right)$

求解方程组(8)中的两个未知参数，并将其作为式(1)参数β₀，β的初值θ₍₀₎＝(β₀₍₀₎，β₍₀₎)′。 

或将传统方法求得式(1)的回归系数作为式(1)参数β₀，β的初值。 

2.对式(1)中β₀，β分别求偏导数，即 

$\∂Q/\∂{\mathrm{\β}}_0=Z_e^{\mathrm{\β}}---\left(9\right)$

$\∂Q/\∂\mathrm{\β}={\mathrm{\β}}_0{Z}_e^{\mathrm{\β}}\ln Z---\left(10\right)$

3.应用式(9)、式(10)以及n组实测值(Z_el，Q_l)，i＝1～n，可计算偏导数矩阵J(θ₍₀₎)以及f(θ₍₀₎)。 

4.根据式(7)可计算θ₍₁₎。 

5.再以θ₍₁₎作为初始值θ₍₀₎，重复步骤3、步骤4，一步步递推下去，直到θ_(k)收敛稳定， 即|θ_(k+1)-θ_(k)|小于或等于预先指定的小正数δ(例如，δ＝0.0005)，从而得到非线性回归系数β₀，β的估计值。 

应用实例1：某水文站水位流量关系值如表1^[5]。根据散点图，可选配幂函数。分别采用化非线性为线性的传统方法、高斯-牛顿法推求β₀，β的估计值、残差平方和 相关指数R²，并将其以及文献[5]、文献[6]中本例的计算结果进行比较。 

对于高斯-牛顿法，在实测值中利用(9.48，258)、(14.55，574)两组关系值确定待估参数的初值θ₍₀₎＝(β₀₍₀₎，β₍₀₎)′＝(3.87533，1.86660)′，迭代结果见表2，可见收敛较快。不同方法的计算结果见表3。 

表1某水文站水位流量关系值 

表2高斯-牛顿法参数的递推结果 

表3不同回归方法的计算结果 

实例1表明，高斯-牛顿法显著优于传统方法，拟合精度明显提高；高斯-牛顿法与文献[5]方法所得结果相同、优于文献[6]的结果，而高斯-牛顿法比文献[5]、文献[6]方法简便。 

应用实例2：某水文站流量与相应水位数据如表4。根据散点图，可选配幂函数。分别采用化非线性为线性的传统方法、高斯-牛顿法推求非线性回归系数的β₀，β的估值、残差平方和 

相关指数R²。 

高斯-牛顿法所求待估参数的迭代结果见表5，其中待估参数的初值利用实测值中的两组关系值(1.05，28.2)、(1.5，93.5)求得。不同方法的计算结果见表6。 

表4某水文站水位流量数据 

表5高斯-牛顿法参数的递推结果 

表6不同回归方法的计算结果 

实例2表明，高斯-牛顿法显著优于传统方法，拟合精度显著提高。 

有关文献 

[1]SL247-1999.水文资料整编规范[s]. 

[2]袁志发，周静芋.多元统计分析[M].北京：科学出版社，2002.98. 

[3]金炳陶.概率论与数理统计[M].北京：高等教育出版社，2000.159-160. 

[4]张子贤.化非线性为线性的加权回归方法及其应用[J].水文，1991，(2)：46-48. 

[5]周世良，尚明芳，李怡，等.拉格朗日乘子法在水位流量关系拟合中的应用[J].水文，2011，31(5)：15-17. 

[6]杨晓华，陆桂华，郦建强.自适应加速遗传算法及其在水位流量关系拟合中的应用[J].水文，2002，22(2)：14-18. 

[7]方开泰.实用回归分析[M].北京：科学出版社，1988.168-172 。 

Claims (1)

1.幂函数型水位流量关系的一种拟合方法，根据某水文站待拟合的水位流量关系值(Z_el，Q_l)，i＝1～n，拟合幂函数型水位流量关系：

$Q={\mathrm{\β}}_0{Z}_e^{\mathrm{\β}}$

其特征是：首先确定非线性回归系数β₀，β的初值θ₍₀₎＝(β₀₍₀₎，β₍₀₎)′，然后采用高斯-牛顿法确定非线性回归系数β₀，β。