现代洪水预报技术简介

洪水预报员当前迫切需要掌握的是两类技术：一类是现代系统预测科学的一套基本方法，促使现行传统洪水预报方法发生质的变化；另一类则属计算机技术，熟练地将图形交互技术应用于预报方法之中。

6.3.2.1　线性最小二乘估计

线性最小二乘估计是迄今应用最广泛、研究最深入的系统参数估计方法，且派生出来的新算法繁多。使用的方式又可区分为离线方式、在线方式和动态方式。

1.离线最小的二乘估计

离线最小二乘估计是最基本的参数估计方式，也是其他方式的基础。

（1）线性系统参数估计的通用方程组。洪水预报中使用的线性系统数学模型有差分方程、响应函数、状态空间等形式。用于参数估计时，也可将它们转化为通用的线性方程组的形式，即：

式中：ФN为系统输入矩阵，N×n维；zN为系统输出向量，N维；θ为系统参数向量，n维；eN为模型噪声向量，N维。

n是模型参数的个数，N是系统的输入、输出变量的观测数据的组数。也就是说，系统依据N组输入、输出数据去估计系统的n个参数。一般来说，N≫n。eN则是估计值与观测值之间不闭合的残差，也称为模型的噪声。因为估计的参数是真实系统的一种近似，故凡是系统估计（识别）模型都有这个残差存在。

对于任何一个实际的洪水系统，只要把系统的数学模型转换为式（6.19）的形式，就可以使用本节介绍的算法，实现参数的最小二乘估计（识别）。需注意，Ф阵中的元素有时只包括系统的输入变量（例如用于响应函数估计时），有时也包括部分已知的输出变量（例如用于差方程参数估计时），这随着对于问题的归纳方式不同而改变。但是，对式（6.19）中θ的估计方法则是完全相同的、通用的。

（2）最小二乘估计准则。确认估计的参数结果是“最优”的最小二乘准则是：用估计的参数还原计算系统的输出估计值与已知值的误差平方总和达到最小。相应的估计参数就称为最小二乘估计值。由此导出的参数计算方法即为最小二乘估计方法。

据此，最小二乘估计准则的目标函数为：

（3）最小二乘估计计算公式。根据最小二乘准则式（6.20），可以推导出模型式（6.19）的参数计算公式。

将式（6.19）中eN项移到等式左端，其余项移到等式右端，并代入式（6.20），得到：

依据极值理论，E{θ}达到极小值的位置在∂E/∂θ＝0处，即：

则：

式（6.23）称为式（6.19）的正规方程组或法方程组。它也是一个线性方程组，容易求出它的解向量为：

这就是模型式（6.19）参数的离线最小的二乘估计公式。如用式（6.23）法方程直接求解，容易夸大方程组的病态问题，采用正交化分解方法，算法性能较好，效果可靠。

2.在线最小二乘估计

对于一个运行中的系统，如果每增加一组新观测数据后，企图计算系统参数的变化，使用离线估计方法就需全部数据重新进行一次计算，不仅计算工作量大，还需把历史数据全部保存起来，很不方便。在线递推算法可以解决这一问题。

通用方程式（6.19）的矩阵、向量的具体表达形式：

其中

符号Ψi代表Φ阵的第i个n维行向量，记为：ψi＝[φi1，φi2，…，φin]T，也就是说，输入矩阵由N个（观测数据的组数）行向量构成。

当有N组观测数据进入程序之后，可以得到它的离线最小二乘估计值，即式（6.24），如使用符号，则有。当第N＋1组新观测数据进入程序之后，ΦN＋1在ΦN的基础上增加了新的行向量：

输出向量ZN＋1比ZN也增加了一个新元素zN＋1。这时，可以推导出根据N组观测求出的最小二乘估计值和根据N＋1组观测求出新的估值的递推关系。

3.有限记忆最小二乘估计

在线最小二乘估计与离线估计完全等价，仅在实现方式上运用了递推的技巧。每一组新的观测数据ΨN＋1、zN＋1最后都进入到记忆阵PN＋1和的计算结果中去，其作用在后续递推中永远存在，不会消失。故这种方式被称为具有“无限记忆”的递推方法，当N不大时，每一组新数据都对具有一定影响和修正作用。但是，当N到达数百、千乃至万的数量级之后，ΨN＋1与PN记忆的历史数据的作用相比就微小了，这时系统的增益、新息的修正作用也愈来愈小，与的差异微乎其微。称估计达到了“数据饱和”，系统进入了“稳态”。

如果一个客观系统本身的参数就是接近于定常的，系统的稳态参数已趋近于系统的真实参数。这种系统可称为“稳态系统”。但是，如果系统本身是一个真正的动态系统，这种数据饱和现象就会使递推最小二乘估计丧失继续跟踪系统参数变化的能力。所以，实现系统动态跟踪还需要对普通的递推估计进行改进。

从“无限记忆”这一特征出发，对在线最小二乘估计方法的第一种改进就是实现“有限记忆”的递推。它的思路是：设定一个系统记忆长度M（即记忆观测数据的组数，如果是等距采样的对象，它也代表着一个时间长度），当有一组新的观测数据进入识别程序后，把原先记忆M组数据中的最早进入的一组数据同时舍弃，保持用M组数据进行最小二乘估计。早于M组的观测数据对于新的参数估计就不再起作用，这种方式估计所得参数永远是反映当前最新的M组观测数据的结果，故称为“有限记忆最小二乘递推估计”。它所得到的参数估计值就有不断地跟踪系统新变化的性能。

4.衰减记忆最小二乘估计

将最小二乘目标函数的等权结构改为加权结构，可以实现加权递推最小二乘法。

将等权目标函数式（6.20）改写为：

式中：ρ为加权因子，0＜ρ≤1。若ρ＝1，就是等权的情况。所以，加权最小二乘法在更高一级层次上概括了普通等权方法。

从式（6.27）可见，当k＝N，也就是资料序列当中（k＝1～N）最新的一组数据所求出的残差eN的权重ρN－N＝1为最大，k逐个向前减少时，ek的权重则从ρ，ρ2，ρ3，…，ρN－1逐个依次减小，由于减小呈指数方次关系，ρ的大小对于权重的衰减速度影响很大。

6.3.2.2　最大似然估计

最大似然估计方法是根据数理统计理论中似然函数极大化原则构建起来的。所谓似然函数就是条件概率密度函数。若对式（6.19）所示系统给出一组输入的观测值{uN}，如果系统参数θ采用不相同的一组值，则上式可产生出一组不相同的观测计算结果，由此可以建立一个条件概率密度函数P（{zN}|{uN}；θ），使估计值随着θ而变化时，其出现概率不同，按统计推断的原则，概率大者出现可能性大，可以认为实测{zN}是该系统的条件概率达到极大时的结果，这时相应的θ值作为似然函数极大化的选择结果就是参数θ的一种估计结果，称为最大似然估计值。最大似然估计方法也分为离线和在线两种形式。

1.离线估计方法

由于似然函数的建立与模型的具体表达形式有关，现使用带有相关残差ek的差分方程表达形式进行推导。当模型形式不同时，方法原理仍同。当ek是m阶相关残差时，用自回归模型表示，即：

记：θ＝[a1…an，b0…bn，c1…cm]T，{zk}＝{zk，zk－1，…，z0}，{uk}＝{uk，uk－1，…，u0}，依照前述最大似然估计的意义，对参数θ的最大似然估计分两步进行，首先依照式（6.28）给出似然函数的表达式，然后对似然函数进行极大的优化计算。

利用多事件概率乘法公式可得：，其中P（z0）是系统输入开始之前，输出的概率分布，由于系统未受激发，可认为处于稳定状态，作已知处理，而：

若为zk的估计值，据式（6.28）有：

显然，ξk是θ确定后zk的估计误差，因为它是θ的函数，又记为ξk（θ）。因为当{zk－1}、{uk}、θ确定后，亦确定，与ξk无关。故ξk（θ）可认为是服从正态分布N（0，σξ）的白噪声序列。

为方便，取条件概率密度函数的对数作为似然函数：

其中，是不依赖于σξ与θ的常数。

从此可直观地看出，L（θ，σξ）取极大，也就是和取得极小值（是正态分布的ξk（θ）的方差，就是zk的预测误差的总和），这时就达到了对zk的最佳模拟。

要使L极大，必有∂L/∂σξ＝0，对式（6.30）求偏导数可得：，因此，L（θ，σξ）极大值求解问题转化为用式（6.29）定义的估计误差所构造的目标函数：

这是一个非线性最小二乘问题，求解方法甚多。以下介绍一种最速下降法。

对式（6.29）求偏导数，得：

使用迭代公式：

ρ（k）的选择应使在ρ＝ρ（k）达到极小。记目标函数的梯度向量为：，则为E在θ＝θ（k）处的梯度向量。

根据式（6.31）可得：

其中ξk（θ）由式6.29计算，而：

则可用式（6.31）求出，故∂E/∂θ也相应求出。

理论上已证明，只要系统满足稳定性条件且是完全可控的，输入信号{uk}是2n阶持续激发的，最大似然识别方法有良好的渐近收敛性。本法求出的θ是无偏的强相容的估计值，且由于方法处理了相关残差，故排除噪声干扰的能力强，对系统信噪比很小（噪声强）的情况，仍能得到较好的结果。缺点是需要反复迭代，计算较繁，计算的精度受最优化方法性能的影响。

2.在线估计方法

由于似然函数极大是从概率意义上获得的，因此在大量数据基础上增加一组新数据对似然函数的影响是很难估计的。所以，最大似然法没有严格意义上的在线算法。但是，如果是新增加了一组M个数据{uk，zk}，k＝N＋1，…，N＋M，则可以导出对参数估计的修正，得到一种递推估计算法。

根据这批M个数据，使：

在满足式（6.29）下极小的θ记为，相应于数据1～N组估计的参数记为，合起来最终修正的参数估计为，它应当更逼近真值θ。

对式（6.29）中ξk（θ）在处作泰勒级数展开，取前两项近似：

这时　

式（6.36）的梯度向量为：

令，解出：

其中偏转矩阵为：

梯度向量在处的值为：

则可以采用：

作为在满足式（6.29）条件时，整个1～N＋M区间的参数最大似然估计值。即是说，以新加入的M组数据的目标函数的偏转矩阵与梯度向量的乘积作为参数估计的修正值，可以直接完成参数的递推修正，不再重复已做计算，不存储历史数据，也无需进行优化迭代。

仿照离线方法，式（6.41）可略作改进。令：

这里ρ（N＋M）是使：

达到极小的ρ，仍用二次曲线近似，则：

6.3.2.3　线性无偏差最小方差估计

理论意义上最佳的参数估计方法是线性无偏最小方差估计。在线性估计的诸方法中，由于它保证能达到无偏和方差最小两项目标，故被称为“最优线性估计”。

1.离线估计方法

将式（6.19）表达的一般线性系统简记为：

作如下假定后讨论线性无偏最小方差估计的方法。

（1）模型的噪声e是一个零均值，协方差阵为R的噪声，即：

（2）模型噪声e与参数互不相关，即：

对于式（6.45），把对参数的所有线性估计结果均记为：

显然，若L＝（ΦTΦ）－1ΦT，式（6.48）就是最小二乘估计。

因无偏估计的条件是E＝{}＝E{θ}＝θ，即：

于是，只有在LΦ＝I（单位阵）时，估计无偏。

从式（6.49）得到，在无偏估计条件满足时，＝Lz＝LΦθ＋Le＝θ＋Le，有：

故按照参数估计协方差的定义表达形式为：

这样，推求无偏最小方差估计的问题就归纳为在线性无偏约束条件下：

寻求能使：　

的最优化问题的解（找出L的表达形式）。

这个数学问题的直接求解推导涉及一些复杂的数学问题。从应用的目的出发，只需了解其解的结果，并进行验证证明其正确即可以。

满足式（6.52）、式（6.53）条件的解L*为：

或者说，式（6.45）的无偏最小方差解为：

下面，验证式（6.54）给出的L*所得到L*RL*T确是所有LRLT中最小者，即：

据式（6.54）有：

根据著名的矩阵许瓦兹不等式：

令

则

这就证明了L*确实满足式（6.54），是LRLT中最小者。

同时，检验式（6.52）的条件为：L*Φ＝（ΦTR－1Φ）－1ΦTR－1Φ＝I，因此，式（6.56）确是系统（6.45）参数θ的无偏最小方差估计值。

比较式（6.56）与最小二乘估计值式（6.22），发现二者相差一个R－1。

R是模型噪声{e}的协方差矩阵。如果假定噪声e具有相关性，最小二乘估计也就有偏，无偏最小方差估计则利用R阵实现了对最小二乘估计“纠偏”和达到方差最小。这种情况下，无偏最小方差估计将优于最小二乘估计。但是，这种改善的条件是模型噪声的统计特性—协方差矩阵R需已知。

如果噪声e具有完全的独立性，即：

其中，I是单位阵。将式（6.51）代入式（6.56），即得到式（6.22）。这就是说，在噪声完全不相关时，最小二乘估计与无偏最小方差估计等价。

2.在线估计方法

线性无偏最小方差估计也可递推实现。但是，由于R阵的存在，这种递推必须假定：由N组观测获得的噪声向量eN和N＋1时段eN＋1之间必须保持独立（不相关），也就是在矩阵：

各元素eiej中，凡是i≠j且i＝N＋1或j＝N＋1的元素全部为零，于是式（6.62）可记为分块矩阵：

在从N向N＋1递推时，新进入的观测值为ΨN＋1、zN＋1，递推时，其他矩阵、向量记号定义为：

令

则对应N组观测数据的无偏最小方差估计为：

当第N＋1组观测值进入后，由式（6.64）可得到：

对此式使用矩阵求逆公式，可得：

其中，KN＋1称为增益阵：(www.xing528.com)

按的离线计算公式，有：

式（6.66）、式（6.67）和式（6.68）就构成了无偏最小方差递推估计公式。式（6.68）中就是新息，它由最新输出观测zN＋1与用老参数和新输入观测ψN＋1计算的预测与观测值之差求出。新参数的修正量由老参数与新息加权和确定，权重大小由系统的增益阵大小而确定，这种递推算法的形式与最小二乘法递推完全相同。

但是，必须要注意，无偏最小方差的递推估计与离线方法不是等价关系，新噪声与老噪声之间需保持独立性，RN和rN＋1需作为先验知识已知，这增加了方法在实用时的处理难度。

6.3.2.4　卡尔曼滤波

1.卡尔曼滤波的假定

卡尔曼滤波突破了经典控制理论频域分析的局限性，是一种可用于平稳非平稳、线性非线性、集总和分布，多输入、出系统等相当广泛领域内的现代动态系统分析技术。卡尔曼滤波以系统状态空间模型为分析对象，根据受噪声干扰的系统模型和包含噪声干扰的系统观测量，运用现代随机估计理论给出了系统状态的无偏最小方差的递推估计值。如果估计的是系统状态现状刻的值，这称为“滤波”，如果是系统状态过去时间的值，则称为“内插”（或“平滑”），如果是系统状态将来的值，则称为“预报”。

估计理论已经证明，无偏最小方差估计是线性估计中性能最佳的一种估计。卡尔曼滤波即是一种最佳线性估计方法，故又称卡尔曼滤波为最佳线性滤波。

卡尔曼滤波方法的基本假定是：模型噪声和观测噪声都是高斯白噪声，其误差协方差矩阵Q、R为已知，模型噪声与观测噪声互相独立而且也与初始状态无关。以此为条件建立了“正规卡尔曼滤波”方法。但实际问题自然复杂得多，为了解决实际使用中存在的各种问题，各种滤波处理技术相应产生：如果噪声不是白噪声，则建立了有色噪声改正处理技术；如果Q、R阵未知，就可使用在滤波器运行中自动进行Q、R估计校正的自适应滤波技术；如果滤波出现模型误差积累而真实发散，则可使用误差补偿的各种技术；用于非线性模型时则用推广的卡尔曼波滤器技术。此处仅初步介绍正规卡尔曼滤波。

若用马尔可夫型状态空间表达线性系统，如下式所示：

2.卡尔曼滤波问题所需的已知条件

（1）模型系数矩阵：状态转移矩阵Φ，模型噪声分配阵Г，观测阵H。

（2）噪声统计特征：模型误差协方阵Q、观测误差协方阵R。

（3）初始状态向量x^0，初始状态向量估计的误差协方阵

（4）模型噪声wk，观测噪声vk均为高斯白噪声，即：E{wk}＝0，E{vk}＝0，对一切k成立，，其中δk，j为克罗内克算符，

模型噪声协方差阵Qk，j和观测噪声协方差阵Rk为半正定阵。

（5）模型噪声与观测噪声不相关，也与初始状态无关，即：

，对一切k，j成立，

，对一切k成立。

3.卡尔曼滤波解决的问题

依据带噪的系统观测量z1，z2，…，zk推求以下各项：

（1）xk＋i的线性无偏最小方差的估计值（i＜0内插，i＝0滤波，i＞0预报）。

（2）同时给出此估计的误差大小。用状态向量预测误差协方阵Pk|k－1，状态滤波误差协方阵用Pk|k表示，有：

其中所有下标中的“|”划线前表示该值所在的时刻，其后代表推求该值所依据的时间。卡尔曼对上述问题导出了以下5个递推估计公式，合称为“正规卡尔曼滤波器”。

（1）状态预测：

（2）状态预测误差协方阵：

（3）增益矩阵：

（4）状态滤波：

（5）状态滤波误差协方阵：

根据已知的（即x0|0），P0（即P0|0），以k＝1，2，…，顺序代入式（6.71）～式（6.75），就实现了正规卡尔曼滤波。

4.使用正规卡尔曼滤波器的注意点

（1）作为线性无偏最小方差的递推估计，卡尔曼滤波器性能良好，但导出的依据是已知条件（4）、（5）所作的基本假定，使用时不可忽视。

（2）如果系统稳定，系统初值、P0的假定对滤波器影响不大，它将会很快收敛到与、P0无关的系统估值上。

（3）Q、R阵是滤波器设计的关键参数，但对实际系统它的正确估计又是一个比较困难的问题，以下将进一步对此进行讨论。

（4）式（6.75）在实用上还有一种改进形式（证明从略）：

由于它的表达Pk|k为一个对称二次型矩阵，可以使数值计算的稳定性有所改进。

（5）如果系统模型使用非马尔可夫型的算式，即有系统的输入控制项uk，这时滤波器作两处修改，式（6.71）和式（6.74）分别改为：

6.3.2.5　实用最优化方法

最优化技术是在高速计算机出现以后，由于科学实验、社会生产等大量实践的客观需要推动而迅速发展起来的。最优化技术范围甚大，前面所述最小二乘等方法，实际上也是最优化方法，只是它们实现了完全的解析分析，最优解可以用数学表达式直接进行计算。运筹学、决策论讨论的课题也属于最优化技术。最优化技术在数学规划中使用最多。数学规划包括线性规划、非线性规划、几何规划和动态规划。习惯上，将上述学科使用的在函数极值理论基础上建立和发展起来计算方法称为最优化方法。

最优化技术是系统分析的基本工具。数学模型的建立，归根到底就是参数的“优选”，它运用的就是最优化技术。在洪水预报应用中经反复检验，证明其可靠的几种实用性强的最优化方法有：网格逐步缩小法，模式搜索法，LM（麦夸特）法，惩罚函数法。

1.网格逐步缩小法

传统的网格法就是一种直观的优劣比较法。如果系统数学模型的参数用向量记为：

系统的目标函数记为：

只需要将θ的每个分量按一定的步长划分，即相当于设定计算采用的不同的θi方案m个，直接进行模型试算，求出θi对应的m个目标函数，再直接比较其大小，目标函数最小的一组参数就是最优的结果。

网格法的主要优点是可靠，因为无需设立参数初值，也不涉及搜索方向等任何判断，不会在算法的组织环节上发生问题。缺点在于如果步长选择较大，可能会漏掉最优的参数，但步长选择较小，则计算量剧增，效率很低。在实际使用时，一般而言，每逢其他算法出现问题之后，反过来用本法都可顺利完成优化任务。

2.模式搜索法

本算法属于直接寻优类的方法。通过对参数的反复迭代计算，对目标函数大小进行比较而搜索到最优参数。方法的优点是勿需对目标函数求导，也不限目标函数有无数学表达式，很实用。缺点则是当参数向量维数较高时计算量偏大，收敛的速度较慢，效率不够高。

方法的基本思路是：寻找目标函数的“极小值点”，犹如一个人站在山坡上“下山”（找极大值则为“爬山”），从站立点（参数的初值）出发，向四周探索，找到一个比站立点低得多的方向，就向下走出一步。反复这种步骤，一定可以下到一个再也找不到比他站立处更低的位置。所以，模式搜索法包括两个步骤，首先进行“探测性移动”（寻找目标函数在多大步长上可能减值的方向），然后进行“模式性移动”，按找到的步长（它由前一步探测得知）让探测基点（出发比较点）移动到新位置。然后，反复进行这一过程，直至找到最优参数为止。

3.LM（麦夸特）法

在最优化方法中，如果碰到目标函数是平方和函数的形式，就可以利用它构造更有效（收敛快）的特殊算法。即当模型的目标函数为：

如果ri（θ）是θ的线性函数，则属线性最小二乘问题，可采用前面的方法处理，否则属非线性最小二乘问题。这里的讨论针对非线性问题。

采用向量函数记法：

则目标函数式（6.80）可记为：

求式（6.82）的最优解θ的方法非常多，常用而著名的有高斯—牛顿法（简称G-N法），最速下降法等。以后，Levenberg-Marguarbt改进了G-N法，构造了性能更好的LM（麦夸特）法，并使G-N法和最速下降法成为LM法的特例。

4.惩罚函数法

前述各种方法，参数优选都没有对参数取值的约束条件，属无约束最优化方法。在实际水文模型中有一部分要求参数取值有一定的范围约束，在指定的约束条件下来推求参数最优解，属于有约束最优化方法。

有约束最优化方法种类繁多。“惩罚函数法”是迄今诸多算法中最有实效的算法。惩罚函数法最大的优点是可靠，适用条件宽松，及少出现求解失败的情形，因而对于有约束的水文模型的最优化计算很适用。

惩罚函数法的思路如下：改造有约束条件的目标函数的结构形式，在等价的条件下，将各个约束条件变为“罚函数”纳入新的增广目标函数之中，从而把有约束条件的优化求解转换为等价的无约束条件的优化求解。

对于n维自变量（可视作模型参数）：

模型的目标函数为：

求解x要求服从m个等式约束条件：

以及（p－m）个不等式约束条件：

建立罚函数：

其中ρi≥0称为罚因子，H、G表示用hi（x）或gi（x）构造的某种形式的函数（泛函）。

建立一个搜索逼近的过程，令增广的目标函数为：

如果逼近的极限过程满足条件：

于是就有：

式（6.91）表明，当k→∞，无约束条件的式（6.88）的解就逼近了有约束条件的式（6.84）的解，实现了惩罚函数法的目标。

6.3.2.6　逐步回归方法

由于逐步回归方法在函数拟合时可以自动挑选最佳的因子来组合拟合函数，具有很大的弹性和较高的拟合精度，故在洪水预报上可成为一种很有用的系统显式数学模型的离线建模工具。

1.n元线性回归方法及其推广

若函数y是由n个自变量x1，x2，…，xn的线性组合所决定，即：

式中：β0、β1、…、βn为方程自变量的系数；ε为方程的随机误差，一般假定随机变量ε为正态分布N（0，σ2）。

要求用m组互相独立的观测值（m≫n）

（x1t，x2t，…，xnt，yt） （t＝1，2，…，m）

在最小二乘准则下，求式（6.92）的n＋1个未知数的估计值b0，b1，…，bn，得出回归方程：

这种问题称为n元线性回归。式（6.94）代表n维空间的一个平面方程，称为回归平面。

运用变量代换的方法，可以将n元线性回归方法进行推广。例如令：x1＝x，x2＝x2，…，xn＝xn或，于是就将n元线性回归方法推广到非线性系统显式回归型模型建模应用领域。

对式（6.93）的参数估计一般是在m≫n情况下进行，这属于方程的线性最小二乘问题。由于系统矩阵阶数高，不仅计算量较大，且容易出现病态问题等麻烦，因此，在回归计算中，常常将回归方程式（6.93）进行“正规化”处理。

设式（6.93）代表的回归平面与观测值之间的残差记为et，有：

根据最小二乘准则：

按极值条件，即可求出bi值，对于：

可得　

其中　

将式（6.96）代入式（6.95），并令

即得：

令

记：

由式（6.99）可得：

其中i＝1，2，…，n，此式称为回归方程式（6.93）的“正规方程”，式（6.101）实行了“标准化”。实际上，它的系数矩阵是n个变元之间的自、互相关系数，称为相关矩阵。它的右端向量是n个自变元与应变量y的互相关系数。此式是普通的n元线性方程组，只要相关矩阵非奇异，按任何线性代数方法均可求出其解向量b＝[b1，b2，…，bn]，再由式（6.96）求出的b0，回归方程式（6.93）便可建立。

2.逐步回归的基本思想

在实际问题中，n个自变元（因子）x1，x2，…，xn与y的关系切程度并不一致，而且有时单个因子的作用与多个因子的综合效果也不一样。计算经验表明，那些与y关系不甚密切的因子不放进回归方程，常常可以提高回归方程的精度。也就是说，希望建立一个偏回归方程：

它包括一切对y作用显著的因子（共p个），而不包括任何对y作用不显著的因子，且满足它是在所有可以建立的n个因子的回归方程式（6.93）当中回归精度最高的一个。实现这一目标的n元线性回归方法，由于要通过对因子的逐步筛选来定的，故称为“逐步回归法”。

3.方差贡献及显著性检验

要建立上述目标的逐步回归计算方法，首先要解决回归方程精度的衡量标准以及如何判断各个因子对y的作用的显著性。

在数理统计方法中，是用“方差”来衡量回归精度的。对于式（6.93）来说，有：

其中　

Q总称为回归总方差，它是每个观测值对其平均值总偏离的平方和。Q回称为回归平面的方差，它是每个观测点的回归计算值与观测均值总偏离的平方和。Q余称为剩余方差，它是每个观测点观测值与回归计算值总偏离的平方和。显然，Q余代表回归平面对测观点的吻合程度，就是回归方程的一种精度指标。

因为Q余是对n个因子、m个观测点据求出，当因子的选择个数有所变化，将影响Q余的结果。为了消除自由度不同的影响，便于因子数不同的回归方程的比较，剩余方差可以改用：

或

称为回归方程的剩余方差，称为剩余标准差。

此外，Q余是Q回和Q总的代数和。因此用Q回、Q总二者也可构造出反映回归方程精度的指标。

R称为y对x1、x2，…，xn的复相关系数。当Q余＝0时，R＝1说明观测点全部落在回归平面上，属于理想的回归结果。当Q回＝0，即Q余＝Q总则R＝0，回归平面就是所有测点的均值面，说明没有回归效果。一般情况下，0≤R≤1，R愈大说明回归效果愈好。

在逐步回归方法中，某个因子对y作用显著与否是用该因子对回归方程的“方差贡献”大小来衡量的。

在回归方程式（6.93）的诸因子中，如果包括因子xj在内时，可计算方程的Q余值，再将因子xj从回归方程中剔除，又可以计算回归方程及其剩余方差，这时称：

为因子xj的“方差贡献”。实际上，它表示由于减少因子xj后回归剩余方差的增加量。显然，方差贡献为正值，表示该因子对提高回归方程精度有作用。反之，方差贡献为负值，该因子加入反而降低回归方程精度。

正方差贡献值的大小代表因子对y关系密切程度的大小。希望选入回归方程的因子都是对y关系最密切的。因此，还需要有衡量这种密切程度的方法。这在逐步回归计算中是通过“显著性”的检验来实现的。