
2.2 系统微分响应修正方法
2.2.1 系统微分响应误差修正方法的实现步骤
在微分关系中,如果把模型看作一个复杂的函数,自变量的变化量就是模型自变量或影响因素误差量,而因变量的变化量即为模型计算结果与实测值之间的误差量。在一般的水文预报误差修正问题中,没有误差修正的原始模型计算系列与实测系列之间的误差系列是可以获得的。假设知道引起模型计算误差的因素或变量,那么模型或系统对于引起误差的因素或变量的相应导数或导函数是可以计算的,由微分响应关系式可知,三要素已知了两要素,则自变量的变化量就可以通过系统函数的微分响应关系反求获得。以图2-3所示的系统概化图为例,假设其系统响应函数为y=f(x),则通过系统响应理论反求自变量误差估计量的步骤如下:
(1)计算模型模拟误差y-f(x),y表示模型实测值,f(x)表示模型计算值。
(2)分析误差特点,确定引起误差的因素或自变量x。
(3)确定微分响应函数
(4)估计自变量误差Δx。
(5)修正模型预报结果f(x+Δx)。
以上步骤中,鉴于所选用模型结构的确定性,第(1)、第(5)步比较简单、直接,对于确定的模型而言没有太多需要研究的内容,而对于修正量的确定和修正过程而言,第(2)、第(3)和第(4)步是有许多需要研究的,因此函数或系统微分响应误差修正研究的重点在于确定引起误差的因素或自变量及其确定相应的微分响应函数。
2.2.2 单变量系统微分响应误差修正方法有效性证明
本书中为了分析简单且又不失一般性,这里先对单自变量时间函数不考虑导数差分误差的情况进行证明。
对于单自变量随时间变化函数一般可以表达为

假设现有样本观测系列(y1,t1,x0),(y2,t2,x0),…,(ym,tm,x0),如果自变量存在未知误差量ex,则有误差的自变量可以表示为

那么模型的计算值与实测值之间存在的误差系列为

其没有自变量误差修正的误差平方和为

以yi为已知的目标值,x0为已知的自变量值,根据yi估计自变量误差e'x,则修正后的计算结果误差平方和为

误差修正有效性证明就是要证明修正后误差系列的目标函数值小于修正前的目标函数值。即

根据微分响应误差修正方法步骤为
(1)计算模型模拟误差ei=yi-f(x0,ti)(i=1,2,…,m)。
(2)确定引起误差的自变量x。
(3)确定微分响应f'(x0,ti)(i=1,2,…,m),得关系:

(4)估计自变量误差。式(2-32)是矛盾方程组,通常采用最小二乘法获得自变量误差的估计值为

(5)修正模型计算结果为

要证明式(2-31),则将式(2-34)代入式(2-31)左边可得

再将式(2-33)代入上式右边得

问题得证。
2.2.3 一般系统微分响应误差修正方法有效性证明
研究系统响应规律最常用的工具是模型。对于一般多函数组合都可以表达为

式中:Q(t)为系统输出;X(t)为参变量,如时间、空间位置、输入变量等;Ω为所有影响系统输出结果的特征因素向量,包括模型状态变量St、中间变量Md、特征量Cr、参数Par等。
对于模型式(2-37),如果模型相对于要素向量Ω在Rn空间处处连续可导,则有微分表达式:

式中:·为向量点积运算;dQ(t)为系统因素向量改变引起的响应。
当研究系统中一个特征量Cr(其他系统因素量都为已知)变化的微分关系时,式(2-38)变为

特别当特征量增量dCr为1或一个单位时,上式简化为

式中:dQ(t)为系统特征量增量为1引起的单位响应函数。
由式(2-40)可知模型对特征量的偏导数等于特征量单位改变引起的系统单位响应函数,式(2-39)就表达了特征量增量、系统输出增量与系统单位响应函数之间的关系。类似地可以推广到一般的微分关系式(2-38),其中的偏导数函数向量表达的是因素向量单位改变引起的系统单位响应函数向量,所以说微分式(2-38)表达了因素向量增量、系统输出增量向量与系统单位响应函数向量之间的关系。对于离散系统,可表达为向量矩阵形式:

式中:ΔΩ为系统特征因素增量;ΔQ为系统响应增量;U为系统单位响应矩阵。其具体向量矩阵为

式中:Ω(j)为Ω的初值;tm为资料的系列长度。
则反演计算基本方程为

系统规律研究分正演与反演两个方向。已知模型输入、模型结构、模型参数计算模型输出为正演计算。已知模型输出、模型结构确定模型参数、模型特征量或模型中间变量等为反演计算。
研究系统反演确定特征因素量Ω的方法,如果描述系统的模型是线性的,则其反演已有成熟的理论与方法。对于非线性系统是否可以先把非线性模型线性化,然后通过逐步迫近获得其结果?这就是系统微分响应反演方法的构建思路。其实施步骤如下:
(1)给定特征因素量初值Ω(0)。
(2)据已知的特征因素量向量计算导数矩阵U和函数向量f(j)。
(3)确定新的特征因素量向量。

(4)判断Ω(j+1)是否最优值,如果是则寻优结束,否则转步骤(2)继续循环。计算流程见图2-4。
以非线性系统线性化为基础构建的系统微分响应逐步迫近反演方法成败的关键是逐步迫近过程是否收敛。
上述系统响应计算方法,关键是要证明第(3)步获得的新因素量估计Ω(j+1)代入模型计算结果f(Ω(j+1),X(t))比原结果f(Ω(j),X(t))更接近于目标值Q(t)。如果依次循环,每步迫近,最终趋于系统因素量的目标值,就说明方法是收敛的。系统响应线性迫近定理可证明其方法的正确性。

图2-4 逐步反演计算流程
系统响应线性迫近定理:对任一X(t)为参变量的函数Q(t)=f(Ω,X(t)),如果在自变量区域[Ω(j),Ω(obj)]内X(t)取任意值时相对于自变量连续、可导,则对于两组函数值:

若偏导数为

组成的导数矩阵U满秩,则变量增量矩阵为

则f(Ω,X(t))的微分表达式为

使得其函数偏差平方和:


满足关系:

以函数在Ω(obj)的值代入式(2-46)得方程组

以上表达中,Ω(j)为Ω修正前的值,Ω(obj)表示Ω的目标值。上式写成矩阵形式:

其中:

因为U满秩,据式(2-51)偏差平方和最小可得自变量增量解:

式(2-46)代入式(2-48)得

式(2-53)写成矩阵形式有

再将式(2-52)代入式(2-54)得

式(2-55)中(UTU)-1是满秩实对称矩阵,如果m=n,则

对于m>n,据满秩实对称矩阵的乔列斯基(Cholesky)分解定理,可分解为

其中,L为非奇异下三角矩阵,则有

即

证毕。
以上证明说明这种修正方法通过对于自变量误差的估计之后,重新使用模型得到的计算系列比修正前更加接近实测系列,从而也说明了此方法在理论上是有效的。