吉林大学精品课>>专门水文地质学>>教材>>专门水文地质学  §8.5相关分析法 在地下水资源预报中，经常碰到两种类型的变量关系。一类是确定性的函数关系，如一口井的抽水量和水位降这两个变量关系，就满足一一对应的函数关系。另一类是非确定性的依赖关系，如一个开采区的开采量和降深的关系，虽然是相互依赖的，降深大开采量也多，但变量关系不是确定的，对应同样降深的开采量不一定完全相等，找不到任何函数可以表示这类变量关系。此外，开采量和降雨量的关系，开采量和地表水侧渗量或灌水回渗量的关系等，也都有类似的性质，这类既有依赖又不确定的变量关系，统称为相关关系。 相关分析法就是研究上述变量的相关关系，建立回归方程，并用来预报开采量的方法。严格说，研究随机变量和确定变量之间的相关关系称为回归分析，研究随机变量和随机变量之间的相关关系才称为相关分析。为了简单，以后不加以区分，统称为相关分析。 相关分析法的适用范围很广，只要在具备多年观测资料地区，无论水文地质条件多么复杂，条件是否查清，都可建立回归方程进行预报，用回归方程预报的效果，一般来说预报系列较短精度较高，系列较长则精度下降，但对较长系列，可以采用不断校正方程，不断延长系列的逐步预报法。 根据研究的变量数目和相关性质不同，可分为一元线性相关分析，多元线性相关分析，多元非线性相关分析，最后简介常用的逐步回归分析法。 一、一元线性相关分析 设有两个变量Y和X，存在前述的相关关系，对这两个变量进行n次独立观测，得一组观测值： （Xt,Yt）,t=1、2、3……n 如用这组观测值在直角坐标纸上做散点图，点的分布呈直线趋势时，则表示变量Y和X呈线性相关，现用Y的平均状态配一直线方程：  (8-17) 这就是Y对X的回归方程，b0，b称为回归系数，是待定值。 若使上式成为最佳的配合直线，可用最小二乘法确定b0和b值。为此，取全部观测值Yt与的偏差平方和  ① 为最小。这样确定b0、b后配出的直线。同全部观测值（Xt,Yt）的偏差将是一切直线中最小的一条，或者说，它是一切直线中同全部观测值最靠近的一条，所以最能代表X和Y之间的变化规律，其精度可由式①算出的Δ的最小值来代表，在后面检验回归方程的效果时将用到这一点。 因为Δ是b0和b的二次函数，又是非负值，故其最小值总是存在的。按数学分析中的极值原理，取Δ关于b0和b的偏导数，并令导数等于零：   表示对从1到求和，上式也可改写成：   解方程组②，得回归系数：   式中：和，为样本的算术平均值，反映样本的平均水平；SXY 和SXX为变量的协方差， 。 把④式代入（8-17）式，即得到要求的一元线性回归方程，也可简写成下列形式：  （8-18） 这说明，对于一组样本观测值来说，回归方程是一条通过散点图几何重心的直线，b是直线的斜率，明确这一点，对做散点图有帮助。 但是，建立的回归方程能否正确反映Y和X的变化规律，不能肯定。因为用最小二乘法建立回归方程时并没有用到Y和X必然存在线性相关的假定，即使观测值在散点图上呈现完全散乱的点子，没有线性相关，同样可以建立一个回归方程，只是这种方程毫无价值而已。所以还要解决两个问题：①变量之间的线性相关是否存在，即变量的相关程度问题；②用回归方程预报时，能有多大的随机误差，即预报的精度问题。这是检验回归方程有无实际价值的两个重要问题。 为了建立检验标准，我们从方差分析着手。把观测值Yt对其平均值的总偏差平方和进行分解：  最后的等式，是由于交积差为零。由③式可知： 于是得分解式及其自由度为 S总=S回+S余 ⑤  ，自由度f总=n-1,  ，自由度f回=1(自变量个数)  ，自由度f余= n-2 (取样组数减去两个变量个数) 式⑤说明，总偏差平方和是由两部分组成的。一是S回，称为回归平方和，表示由变量X取值不同引起的Y值的差异，其大小反映了变量之间的相关程度；二是S余，称剩余平方和，表示由观测因素和未加控制因素等随机误差引起的Y值的差异，其大小反映了随机误差对预报精度的影响。据此，可分别确定相关程度和预报精度的检验标准： ⑴ 变量之间的相关程度 取S回在S总中所占的比例大小，作为衡量变量之间相关程度的数量指标：  r称为相关系数，把代入S回中,同时利用④式，可得r的计算式：  (8-19) 式中： ,时，表示Y和X没有线性相关，观测值在散点图上呈散乱点子；也可能存在非线性相关，散点图呈曲线趋势。r=1时，表示Y和Z呈线性函数关系，所有观测值都应落在同一直线上。,表示Y和X存在线性相关。越大线性相关越密切，观测值在散点图上越靠近配合直线的两侧分布；越小则线性相关越微弱，观测值在散点图上越远离配合直线。 根据(8-19)式计算的值，只要在0～1区间就表示变量之间存在线性相关，但只要值大到一定程度时才有明显的相关关系。那么究竟要多大时Y和X之间才有明显的线性关系？这就必须对相关系数进行显著性检验。为此，构造统计量：  把代回S回中，再利用④和⑤式，得：  S余=S总-S回=SYY(1-r2) ⑥ 代入F式中得计算式：  （8-20） 上式服从自由度为1和n-2的F分布。在给定显著性水平α下，把计算的F值同附表I中达到显著时的临界值做比较，若，则说明在α水平上两变量之间的线性关系是显著的，否则说明线性关系不显著。 由（8-20）式知，F和r存在一一对应关系，按F分布表中的临界值可以求出相应r的临界值，如附表II。因此，在实际检验时，也可用（8-19）式计算的r值同附表II的临界值直接做比较，若r超过表II的临界值可认为变量之间的线性关系是显著的，反之，低于临界值则认为线性关系不显著。 所谓显著性水平α0是指做出“显著”这个结论时发生判断错误的可能性。如α=0.05时，表示判断错误的可能性不超过5%；α=0.01时，这种可能性不超过1%。检验时选取的α值越小，检验越严格，相应的要求r值就越大。例如，n=12时，如果计算的，则说明r值在α=0.05的水平上（概率或可信度为95%）显著。如把可信度提高到99%，计算的相关系数有，则说明在α=0.01水平上显著。α越小，r的显著程度越高，变量之间的线性关系越密切，所建的回归方程才越有价值。 (2) 回归方程的预报精度 已知S余表示随机误差，即观测值并不完全落在配合直线上，而是散布在它的两侧所存在的误差。所以S余越小，用回归方程预报的效果越好，这里的S余，实际上就是前面最小二乘法建立回归方程时所给出的最小的那个Δ值。因此，可用均方差  作为衡量回归方程预报效果好坏的指标。也称为剩余标准差，它表示观测值偏离配合直线的平均距离，把⑥式代入上式得的计算式：  （8-21） 这样计算的越小预报的效果越好，但究竟多小效果最好？仍可用显著性检验来说明。 所谓预报问题，就是在一定的显著性水平α下寻找一个偏差δ，使得按给定的X预报Y时的观测值，以（1-α）的概率落在（Y-δ, Y+δ）的区间内，即：  而δ和有如下的一般关系：  这说明，在一定的样本观测值和显著性水平下，给出的X越接近，δ越小。如在（X，Y）平面上划成曲线，则形成一个包围回归直线的带域，两头呈喇叭形，在处带域最窄。因此，如果给出X预报Y时，则预报区间越短（X越靠近），带域宽度越窄，说明预报饿精度越高。 在实际计算中，n值往往过大，要求的预报区间较短，这时，结果得预报区间的近似式： =68.3% 这说明，在给出X预报Y时的观测值中，有68.3%可能落在配合直线两侧各为一个均方差的范围内，或者说观测值Y可能落在之间的概率为68.3%，同样，观测值Y落在之间的概率为95.4%，落在之间的概率为99.7%，即： =95.4% =99.7% 称为可信区间。用回归方程预报Y时，越小，可信区间越窄，预报的精度越高。当r=1时，Y和X呈函数关系，这时=0，表示所有观测值都落在配合直线上。 控制是预报的反问题。即在一定概率P下，要求Y在区间（Y1，Y2）内取值时，X应控制的范围或求出相应的X1和X2，例如，概率为68.3%，如取   则给出了Y1和Y2，就可解出X1和X2作为控制X的上下限。 〈实例〉，某水源地已有多年的开采历史。经分析，取其中十年的开采资料进行相关分析。为了扩大开采，要求预报降深S=26m时的开采量Q。原始数据和计算数据均列入表（8-10）中，表中采用的符号Q～Y，S～X。 根据表中的数据，按下列顺序预报： 首先，计算相关系数，按（8-19）求得：  按n=10，查附表II知，在显著性水平α=0.01时相关系数的临界值为0.765。显然，计算值0.997〉0.765，所以Q和S之间的线性关系是显著的。 其次，建立预报方程，进行预报。按④式计算回归系数。  把B值和，代入（8-18）式，整理后得回归方程：  代入S=26m,，得预报的开采量Q=97.43×104m3/d。 最后，检验预报精度。按（8-21）式计算剩余均方差  表（8-10） 年份 开采资料 计算数据   Qt (104m3/d) St (m)       1959 60 16.5 -13 -3.3 169 10.89 42.9  1960 67 18 -6 -1.8 36 3.24 10.8  1961 60 16.5 -13 -3.3 169 10.89 42.9  1962 63 14.7 -10 -2.3 100 5.29 23  1970 80 21.5 7 1.7 49 2.89 11.9  1971 80 21.9 7 2.1 49 4.41 14.7  1972 78 21 5 1.2 25 1.44 6  1973 76 20.5 3 0.7 9 0.49 2.1  1974 82 22 9 2.2 81 4.84 19.8  1975 84 22.5 11 2.7 121 7.29 29.7  总和 730 197.9   808 51.67 203.8  平均 73 19.8       如取概率为68.3%，则有  这说明，预报S=26m的开采量时，Q的观测值将落在96.65-98.2区间的概率为68.3%，也就是实测开采量将在宽约15.6×104m3/d的带形区内波动，所以回归方程的预报精度是足够的。 二、多元线性相关分析 地下水资源总是在多种因素影响下形成的，所以在资源评价中采用多元线性相关分析的 效果总比一元相关分析为好。参与多元相关的各种变量可按当地具体条件确定。如在潜水区，可以考虑开采量和潜水位、开采面积、历年的降雨量、蒸发量等各变量的相关关系；在承压水区，可以建立开采量和水位、侧渗补给量、回灌量、开采时间等的相关关系；矿井疏干区，可以建立涌水量和地下水位、降雨量、巷道长度等的相关关系；在泉群地区，可以建立泉群流量和前一年的流量、前几年的降雨量、蒸发量等的相关关系。下面简介多元线性相关的原理和应用。 设有因变量Y和自变量Xi（ i=1,2,3……m）。对Y和Xi进行n次观测，得一组观测值为(X1t , X2t , ……Xmt , Yt) ， t=1,2,3……n。 如果Y和Xi之间存在线性相关，则可配一线性回归方程：  (8-22) 式中：b0, bi, i=1,2,3……m , 为回归系数,是待定值。 为使上式成为最佳的配合直线，同一元相关分析法一样，可用最小二乘法确定m+1个回归系数。因此，取全部观测值和的偏差平方和  ⑦ 达到最小。按数学分析中的极值原理，首先取  由此解得  ⑧ 式中：，； 把代入⑦式，整理后得：  再取上式对的导数。并令导数为零：   ⑨ 如果    则由⑨式得：   （8-23） 这是一个含m个待定值和m个方程的联立方程组，也称回归方程（8-22）的正规方程组。如果它的系数（）行列式不为零，则可用消元法解出。再由⑧式求出，把求出的和代回（8-22）式，即得要求的多元线性回归方程。 在实际计算中，由于变量Y和X i可能有不同的量纲和数值上的很大差异，而正规方程组（8-23）的系数矩阵（），又都是X的二次项，所以之间在数值上的差异可能非常之大，求解时难免会给结果带来较大的舍入误差。为了避免这一点，常将正规方程组（8-23）标准化。取    或  ，  ⑩ 把⑩式代入（8-23）中，得：  在等式两边同除以，得  已知：   结果得正规方程组的标准形式：   （8-24） 式中：为标准回归系数；为变量和之间的单相关系数；为变量和Y之间的单相关系数。它们都是无因次量。且恒有，方程组（8-24）也可写成矩阵形式：    相关系数矩阵是满秩的，的逆矩阵存在。所以，可以用伴随矩阵求逆矩阵的方法解出  (8-25) 式中：为伴随矩阵，为矩阵的行列式。 把所有的B式代入⑩式求出，由⑧式求出和再代回（8-22）式即得要求的回归方程。 所求的回归方程是否反映Y和Xi之间的变化规律，还要进一步做方差分析才能清楚。为此，仍取对其平均值的偏差平方和：  同一元相关的分析方法一样，上式最后一个等号是因为交叉积。由（8-23）式和⑧式知：  = = = = 于是得分解式和自由度为 =+ （11） =，自由度； =，自由度（自变量个数） =，自由度； Y和Xi整体之间的相关程度，不能再用单相关系数r来衡量，必须用复相关系数R表示。R的定义和单相关系数一样，仍取和中所占的比值：  利用（8-22）式和⑧式，可表示为： == = = （12） 把结果代入R式中，得复相关系数的计算式：  （8-26） 用上式计算R很方便，可以利用解方程组的已有结果直接计算。R的变化范围为0-1。R=0，说明Y和Xi整体之间没有线性关系；R=1，表示Y和Xi整体之间为函数关系。当R在0～1之间变化时，R越大，表示Y和Xi整体之间的相关越密切。 同一元相关一样，复相关系数的显著性，可用下列统计量来检验：  为了计算方便，利用R的定义和（11）式，上式可改用R表示： （8-27） 把值同附表Ⅱ中的临界值作比较。若计算的＜，则说明Y和Xi整体之间的线性关系不显著；若＞则表明线性关系显著，所建的回归方程基本上反映了Y和Xi整体之间的变化规律。 同一元相关一样，仍用剩余的均方差  做为衡量回归方程预报精度的指标。由（11）式可改用R表示  = （8-28） 当n值较大，预报区间较短（给出的Xi分别靠近）时，可用下列近似检验预报或控制的显著性： =68.3% =95.4% =99.7% 实例1：泉群流量的相关分析。该泉主要接受灰岩山区的降雨入渗补给和砂页岩山区河流进入灰岩地区的河水渗漏补给。经分析，不同年数的降雨量平均值和相关系数的关系证明，五年的相关系数最高，六年的相关系数下降，这说明泉流量同前五年的降雨量有关，由于储存量的调节性，泉流量还受到前几年的流量影响。为了说明计算方法，现建立受前五年平均降雨量和前一年泉流量影响的回归方程：   —一年平均泉流量；—前一年年平均泉流量；—前五年的年平均降雨量。观测数据和计算数据，均列入表（8—11）。 首先，建立多元线性回归方程。接表中数据计算协方差： Syy== — ()2=29.46—=0.529 S11==— =29.879—=0.516 S 22==— =611.764—=6.708 S 1y==— =29.595—=0.449 年度 原始数据 计算数据   Y 10m3/s X1 10m3/s X2 102mm Y2 X12 X22 yX1 yX2 X1X2  1960 1.23 1.32 5.852 1.513 1.742 34.246 1.624 7.1980 7.7246  1961 1.27 1.23 5.685 1.613 1.513 32.013 1.562 7.1857 6.9593  1962 1.21 1.27 6.136 1.464 1.613 37.650 1.573 7.4246 7.7927  1963 1.36 1.21 6.598 1.850 1.464 43.534 1.646 8.9733 7.9836  1964 1.58 1.36 6.758 2.496 1.850 45.671 2.149 10.6776 9.1909  1965 1.52 1.58 6.422 2.310 2.496 41.242 2.402 9.7614 10.1468  1966 1.45 1.52 6.626 2.103 2.310 43.904 2.204 9.6077 10.0715  1967 1.48 1.45 6.572 2.190 2.103 43.191 2.146 9.7266 9.5294  1968 1.43 1.48 5.860 2.045 2.190 34.340 2.116 8.3798 8.6728  1969 1.33 1.43 5.516 1.769 2.045 30.426 1.902 7.3363 7.8879  1970 1.28 1.33 5.796 1.638 1.769 33.594 1.702 7.4189 7.7087  1971 1.24 1.28 5.596 1.538 1.638 31.315 1.587 6.9390 7.1629  1972 1.13 1.24 4.898 1.277 1.538 23.990 1.401 5.5347 6.0735  1973 1.08 1.13 5.414 1.166 1.277 29.311 1.220 5.8471 6.1178  1974 1.07 1.08 4.938 1.145 1.166 24.384 1.156 5.2837 5.3330  1975 1.01 1.07 4.930 1.020 1.145 24.305 1.081 4.9793 5.2751  1976 1.00 1.01 4.928 1.000 1.020 24.285 1.010 4.9280 4.9773  1977 1.15 1.00 5.862 1.323 1.000 34.363 1.150 6.7413 5.8620  总和∑ 22.82 22.99 104.360 29.460 29.879 611.764 29.595 133.943 134.4698  平均值 1.268 1.277 5.7978   表(8-11) S 2y==—  =133.943—=1.6377 S 12= S 21==—  =134.47—=1.1789 于是，按（2～23）式得正规方 程组0.516b1+1.1789b2=0.449 1.17896 b1+6.708 b2=1.6377 按(2～24)式变换为标准方程组： b/1+0.6337 b/2=0.8593 0.6337 b/1+ b/2=0. 8694 联解得：b/1=0.5154, b/2=0.5429 按(10)式和(8)式求得：b1= b/1=0.5291，b2= b/2=0.1525 b0=－b1－b1=1.268－0.5219×1.277－0.1525×5.7978=－0.2827 由此，得线性回归方程：=－0.2827+0.5219X1+0.1525X2 从上式可做近期预报(从略) 其次，检验。按(2～26)式得复相关系数： R===0.9568 按(2～27)式检验其显著性： F=(n－m－1)=(18－3)=80.63 由附表Ⅰ中查得α=0.01时,临界值F0.01(218－)=６.３６，显然F>Fα0， σ余===0.0548 如取可信度为68.3%，则预报时的观测值将有68.3%落在宽约2σ余=0.1096的带域内，可见预报方程有足够精度。 实例2：※济南地区地下水开采量的相关分析。 该区地下水流域面积约为1502平方公里，其中灰岩面积占60%，变质岩面积占7%，第四系覆盖面积占33%。因此大气降水是该区地下水的主要补给源。根据长期观测资料的分析和筛选证明，开采量和四个因子存在明显相关：地区灰岩地下水年平均水位Ht，前一年降水量Pt-1、前二年降水量Pt-2和前三年降水量Pt-3。n=10年的观测资料列入表(8-12)中。 表(8-12) 序号  项目 年份 Qt(y) （104m3/s） Ht(X1) (m) Pt-1 (X1) Pt-2(X2) Pt-3(X3)  1 1968 36.28 28.32 552 508 444  2 1969 42.58 28.19 385 552 508  3 1970 48.66 28.17 712 385 552  4 1971 55.86 27.94 534 712 385  5 1972 62.33 27.83 765 534 712  6 1973 62.76 27.71 569 765 534  7 1974 62.20 28.51 854 569 765  8 1975 75.27 27.85 737 854 569  9 1976 79.79 27.18 578 737 854  10 1977 84.56 26.72 703 578 737   ∑ 610.47 278.42 6389 6194 6060   平均 610.47 27.842 638.9 619.4 606.0   Syy 2253.92 26204 172860.9 182084.4 209340.0  ※ 摘自山东省地质局第三水文地质队 根据所给资料，可建立下列多元线性回归方程：= b0+ b1Ht+ b2 Pt-1＋b3 Pt-2+ b4 Pt-3 计算协方差后按(2～23)式求得正规方程组：Si1 b1+ Si2 b2+ Si3 b3 +Si4 b4= SiQ 2.6204 b1+4.382 b2-225.898 b3 -347.92 b4=-61.11 4.3820 b1+172860.9 b2-18643.6 b3 +100460 b4=851076 -225.898 b1-18643.6 b2+182084.4 b3 +4357 b4=11097.27 -347.92 b1+100460 b2+4357 b3 +209340 b4=15080.09 按(2～24)式变换为标准方程组 +++= b1/+0.00651 b2/-0.32703 b3/-0.46975 b4/=-0.79511 0.00651 b1/+ b2/-0.10509 b3/+0.5281 b4/=0.43163 -0.32703 b1/-0.10509 b2/+ b3/+0.02232 b4/=0.54779 -0.46975 b1/+0.52810 b2/+0.02232 b3/ +b4/=0.69424 由这个方程组可解得：b1/ =-0.55465，b2/ =0.35044，b3/ =0.39790，b4/=0.23975 因而求得：b1= b1/=-0.55465=-16.25604 b2= b2/=0.35044=0.04002 以及b3=0.04427，b4=0.02488，b0=445.57965 于是得回归方程：=445.57965－16.25604Ht+0.04002 Pt－1+0.04427 Pt－2+0.02488 Pt－3 复相关系数为： R== =0.98827 F=(n－m－1)=48.75 由附表Ⅰ中查得α=0.01时的临界值F0.01(4，10－5)=11.39，F> F0.01说明是显著的。 剩余均分查：α余===3.24226 说明预报时的偏差较大，预报时间不宜过长。 用方程计算的1968～1977年地下水开采量曲线和实测曲线吻合尚好，证明利用所建的回归方程，对济南地区地下水开采量可以做近期预报。 三、多元非线性相关分析 在地下水资源评价中，变量间的非线性相关比线性相关具有更大的适应性。因为后者只是前者的特例。但是，非线性相关比较复杂，直到现在还有许多值得探讨的问题。 有些非线性相关问题，可表示为多项式或幂函数形式： ＝b0＋b1Ｘ１＋b２＋b３　　（１３） ＝b0 （１４） 式中b0 ，b1 ，b２，b３是待定系数。 这类非线性问题，通过简单交换，可以化为线性相关问题来处理。 做下列变换：对（１３）式，令：u１＝Ｘ１，u２＝，u３＝ 对（１４）式，令：＝lgy， ＝lgＸ１，＝lgＸ２ 式（１３）和（１４）变为：＝b0＋b1 u１＋b２u２＋b３u３　　（１５） ＝lg b0＋b1＋b２ （１６） 这说明，只要把观测值（），（i＝１，２，３；t＝１，２，……n）分别换成新的观测值（）和（），则非线性的多项式和幂函数便可转化为一般的多元线性回归方程。因此，前两段介绍的线性相关分析法，便可全部用来解决多项式或幂函数这类非线性问题。 数学分析告诉我们，多项式可以分段逼近任何函数，幂函数也有满足多数曲线（或曲面）的性质，所以，利用（１３）和（１４）式可以解决相当数量的非线性相关问题，这两类函数在相关分析中占重要地位。但是，它们只是非线性相关中的特例，还有一定的局限性。 这种方法在地下水动态预报中已显示了计算简单、预报准确的优点，具体方法如下： 设有因变量X0和自变量X1，X1，……Xl，统一记为Xi(I=0,1,2……l)。对全部变量进行n次观测，得一组观测值： (Xit), t=1,2,……n 如果变量间存在单调的（无极值）非线性相关关系，则回归方程可表示为：  (8-29) 这是由非线性单值函数组成的线性形式的复回归方程。其中：ui(xi)为非线性单值函数，正态变量，b1,b2,……bl为回归系数，待定值。 为了客观地确定l+1个非线性单值函数，可对全部变量的观测值(Xit,Yt)按由小到大的顺序排队，标上序号，再按每一观测值的序号求出经验频率  (17) m(Xit)为的Xit序号，n是观测次数或序列长度。 同所求的概率相当的正态变量ui(Xit)，可按正态积分分布函数确定  (18) 其反函数为：  (19) 这样，变量的每一观测值Xit，都对应一个概率Pm，而每一概率Pm又和一个正态变量um相对应，即把所有变量的观测值Xit都转化为相应的正态变量ui(Xit)。这样转化的结果有一特点，因为um只同n和m序号有关，所以和各变量每一观测值对应的um可以不同，但各变量的完全相等。根据这个特点，附表III中列出了n=15～25时，对应不同序号的Pm和um值，它对统一观测系列中的全部变量都适用，可供计算中查用。 以所求的正态变量ui(Xit)作为观测值，再用最小二乘法确定(8-29)式的回归系数b1,b2,……bl。因为该式是线性形式，可以采用多元线性相关分析的方法，取：  达到最小。令，整理后得  (8-30) 其中的单相关系数为 (20) i=1, 2, ……l, j=0, 1, 2, ……l uit 、ujt－－系列中序号为的观测值的正态变量值。i=j时。 式（8-30）是一含个方程和个未知数的方程组，只要它的系数行列式不为零，就可解出回归系数。把已知代回（8-29）式，即得要求的多元非线性回归方程。 因为分析非线性过程中引用的是正态分布函数，所以线性相关计算中有关回归方程可信度的检验理论，均可在这里应用。故复相关系数可按下式计算；  （8-31） 均方差为  （8-32） 根据（8-29）式预报时，可先选定变量，绘制-之间的关系曲线，如图（8-17）。然后用预报年份的已知变量、……值，从关系曲线上查得对应的（）、（）……（）值，代入（8-29）式求得（）值，再从-关系曲线上反查出值，即为预报值。 实例：对某水源地孔，预报开采条件下年最低水位。根据当地水文地质条件和多年地下水水位动态过程线的分析认为，影响历年最低水位的因子主要有三：一是历年最高水位，反映了前一年降雨量的大小；二是历年农灌水开采系数，衡量当年开采量的相对大小；三是历年1-5月降雨量，反映当年雨季来到的时间。这样就归结为一个三元非线性相关问题。 有关的观测资料和Xit转化为等的计算过程，列入表（8-12）和表（8-13）中。 首先，计算单相关系数，按表8-13中的总和值得： = = = = = = 其次，组成方程组：    解之得回归系数  ，  ，  于是得回归方程：  最后，预报和检验。已知1977年的m，，㎜。由关系Xi～ui曲线查得：    把U1，U2，U3代入回归方程求得 -1.1032-0.2042+0.1860=-1.12 从关系曲线上反查，当=-1.12时，对应的X0=20.15m，这就是预报的1977年的最低水位。当年实测最低水位为19.9m，预报误差为0.25m。 复相关系数为R=（）1/2=0.995 预报的标准误差δ余= （1-R2）1/2 =0.1m 表（8-12） Xit-Pit-Ui表 年度 变量观测值 观测值序号 正态变量值   X0t X1t X2t X3t m m m m U0（X0t） U1（X1t） U2（X2t） U3（X3t）       （X0t） （X1t） （X2t） （X3t）      1958 29.8 30.26 1.11 74.6 18 17 3 15 -1.343 -1.076 1.076 -0.693  1959 29.46 29.75 0.85 55.4 15 14 1 10 -0.693 -0.539 1.77 0  1960 29.67 30.22 0.95 19.4 16 16 2 3 -0.871 -0.871 1.343 1.076  1961 29.2 29.77 1.52 34.7 14 15 6 8 -0.539 -0.693 0.539 0.265  1962 28.87 29.56 1.46 32.6 13 12 5 7 -0.396 -0.261 0.693 0.396  1963 28.57 29.5 1.41 183.9 12 10 4 19 -0.261 0 0.871 -1.77  1964 30.18 30.75 1.64 177.2 19 19 8 18 -0.771 -1.77 0.261 -1.343  1965 29.73 30.64 1.59 60.2 17 18 7 11 -1.076 -1.343 0.396 -0.128  1966 27.69 29.36 3.2 65.8 11 9 10 14 -0.128 0.128 0 -0.539  1967 27.51 29.54 6.84 63.5 10 11 12 13 0 -0.128 -0.261 -0.396  1968 26.51 29.74 3.05 30.8 9 13 9 6 0.128 -0.396 0.128 0.539  1969 25.75 27.65 6.44 93.9 8 7 11 17 0.261 0.396 -0.128 -1.076  1970 25.69 28.12 8.18 61.8 7 8 13 12 0.396 0.216 -0.369 -0.261  1971 24.57 26.96 14.2 24.9 6 6 16 5 0.539 0.539 -0.871 0.693  1972 21.93 25.07 13.62 38.4 5 5 15 9 0.693 0.693 -0.693 0.128  1973 20.66 22.76 12.6 78 4 2 14 16 0.871 1.342 -0.539 -0.871  1974 20.25 24.04 17.03 17.3 3 3 18 2 1.076 1.076 -1.343 1.343  1975 19.81 24.41 14.67 17.1 2 4 17 1 1.343 0.871 -1.076 1.77  1976 18.6 22.02 21.08 24.4 1 1 19 4 1.77 1.77 -1.77 0.871  1977  22.66 16.26 126.5          ※ ——历年最低水位，——历年最高水位，——历年农灌开采系数，——历年1～5月的降雨量。 ※根据观测值的m在n=19行中由附表Ⅱ中查得。 表（8-13） 相关系数计算表 年度 正态变量的平方 正态变量的交差   U02 U12 U22 U32 U0U1 U0U2 U0U3 U1U2 U1U3 U2U3  1958 1.804 1.158 1.185 0.48 1.4451 -1.4451 0.9307 -1.1578 0.7457 -0.7457  1959 0.48 0.291 3.133 0 0.3735 -1.2266 0 -0.954 0 0  1960 0.759 0.759 1.804 1.158 0.7586 -1.1698 -0.9372 -1.1698 -0.9372 1.4451  1961 0.291 0.48 0.291 0.068 0.3735 -0.2985 -0.1407 -0.3735 -0.1809 0.1407  1962 0.157 0.068 0.48 0.157 0.1034 -0.2744 -0.1568 -0.1809 -0.1034 0.2744  1963 0.068 0 0.759 3.133 0 -0.2273 0.462 0 0 -1.5417  1964 3.133 3.133 0.008 1.804 3.1329 -0.462 2.3771 -0.462 2.3771 -0.3505  1965 1.158 1.804 0.157 0.016 1.4451 -0.4621 0.1377 -0.5318 0.1919 -0.0507  1966 0.016 0.016 0 0.291 -0.0164 0 0.069 0 -0.069 0  1967 0 0.016 0.068 0.157 0 0 0 0.0334 0.0577 0.01034  1968 0.016 0.157 0.016 0.291 -0.0507 0.0164 0.069 -0.0507 -0.2134 0.069  1969 0.068 0.157 0.016 0.158 0.1034 -0.0334 -0.2808 -0.0507 -0.4261 0.1377  1970 0.157 0.068 0.157 0.068 0.1034 -0.1568 -0.1034 -0.1034 -0.0681 0.1034  1971 0.291 0.291 0.759 0.48 0.2905 -0.4695 0.3735 -0.4695 0.3735 -0.6036  1972 0.48 0.48 0.48 0.016 0.4802 -0.4802 0.0887 -0.4802 0.0887 -0.0887  1973 0.759 1.804 0.291 0.759 1.1698 -0.4695 -0.7586 -0.7239 -1.1698 0.4695  1974 1.158 1.158 1.804 1.804 1.1578 -1.4451 1.4451 -1.4451 1.4451 -1.8036  1975 1.804 0.759 1.158 3.133 1.1698 -1.4451 2.3771 -0.9372 1.5417 -1.9045  1976 3.133 3.133 3.133 0.759 3.1329 -3.1329 1.5417 -3.1329 1.5417 -1.5417  总和 15.732 15.732 15.732 15.732 15.1728 -13.1379 7.4941 -12.19 5.1682 -5.8875  四、逐步回归分析 在前两段介绍的多元相关分析中，都是根据水文地质条件和相关动态资料，尽可能选取较多的自变量建立多元回归方程。一般而言，考虑的自变量越多，S回越大，S余越小，建立的回归方程越有价值。但是，如果在回归式中包含了某些对因变量没有作用或作用不大的自变量，则S余不会因它们的加入而减少多少，反会因S余的自由度减少，引起F值变小，使回归效果在一定程度上变坏。其次，为了计算方便，也要求回归式中尽量包含较少的自变量为佳，以便减少观测项目，减少无效计算量，增大方程的稳定性。如果两个因子完全重叠，则正规方程组就无法求解。由此可见，在回归方程中包含过多而作用不大的自变量，不仅无益，反而有害。因此，在地下水资源评价中，总希望建立一个“最优”的回归方程，以便进行有效的预报。这种“最优”的回归方程，应当包含一切对因变量作用显著的因子，而不包括任何作用不显著的自变量，它的均方差又必须相对最小。 现以数值实例说明“最优”回归方程的直接选优法。初选四个自变量，对这四个自变量的可能组合，可建立15个偏回归方程，（包含一个自变量的四个，两个自变量的六个，三个自变量的四个，四个自变量的一个）。然后，对每一方程中包含的自变量做显著性检验，标**者表示在α=0.01水平上显著，标*者表示在α=0.1水平上显著。全部计算结果列入表（8-14）中。 从表中数据可见，按S余衡量，以四元全回归方程（15）为最小（（47.86），而从均方差来看，则是三元偏回归方程（12）为最小（23）。但这两个方程中都含有作用不显著的自变量，所以都不是“最优”的回归方程。综合观察15个方程，其中包含的自变量都起显著作用，同时均方差又相对最小者是方程（5），即 y=52.5773+1.4683X1+0.6623X2 这就是要求的“最优”回归方程。 这种直观的优选法，虽然可以从初选的许多自变量中找到一个最优的回归方程。但计算量实在太大。如果初步选10个自变量，则要建立210-1=1023个方程。这样大的计算量显然很不实用。所以，在相关分析中发展了逐步回归法。 逐步回归法，是在初选若干自变量的基础上，把其中作用显著的自变量逐个引入回归方程，同时对误选入而作用变得不显著的自变量随时剔除回归方程的一种“优选”法。优选法的基本过程是：先从初选的若干自变量中找一偏回归平方和为最大的自变量，通过F检验，作用显著时可引入回归防横；再从余下的自变量中找一偏回归平方和为最大的自变量。经F检验，证明显著时引入，同时如果因这一自变量的引入而使早引入自变量由显著变为不显著（经F检验）。则应随时将它从回归方程中剔除。这样反复进行，可把一切作用显著的自变量逐个引入方程，而把任何作用不显著的自变量随时剔除方程，直到没有一个作用显著的自变量的可以引入，也没有一个作用不显著的自变量可以剔除时为止，最后的偏回归方程就是要求的读某些批数据讲为最优的回归方程。 实现上述的逐步回归过程，利用相关矩阵的消元变换法比较有效。因此，先复习一下这方面的有关知识。 有方程组（4～24），得标准方程  （21） 式中：rij为自变量Xi和Xj之间的单相关系数；riy为自变量Xi和因变量y之间的单相关系数；是标准回归系数，待求量，j=1、2、……m。 表（8-14） 系数 序号 b0 b1 b2 b3 b4 S余 F余 σ余  1 81.4793 1.8687    126569 11 10.73  2 57.4236  0.7891   906.34 11 9.08  3 110.2026   -12.558  1939.4 11 13.28  4 117.5679    -0.7382 883.87 11 8.96  5 52.5773 1.4683 0.6623   57.9 11 7.61  6 72.349 2.3125  0.4945  1227.07 10 11.08  7 103.0973 1.44   -0.614 74.76 10 2.74  8 72.0747  0.7313 -1.0084  415.44 10 6.45  9 94.16  0.3109  -0.4569 868.88 10 9.32  10 131.2824   -1.1999 -0.7246 175.74 10 4.19  11 48.1396 1.695 0.6569 0.25  48.11 9 2.31  12 71.6482 1.4519 0.4617  -0.2365 47.97 9 2.3  13 111.6844 1.0519  -0.41 -0.6428 50.97 9 2.38  14 203.6418  -0.2934 -1.448 -1.557 73.82 9 2.86  15 62.4052 1.5511 0.5101 0.1019 -0.1441 47.86 8 2.45  由（21）式中的单相关系数可以组成增广矩阵。为了检验的方便，在增广矩阵的下面再填一行（ry1、ry2……rym、ryy），形成一个对称的初始方阵R（0），以后简称相关矩阵，记为：  （22） 由线性代数知，如在R（0）中对待求量进行下列消元变换：  （23） 则依次可得新的相关矩阵：  （24） 和，，…………。这种变换过程有下列重要性质： 1.每一步消元变换的结果可得到一个子方程组的解及其相关矩阵的逆矩阵。例如，对待求量做消元变换得，则得子方程组的解为及相关矩阵的逆矩阵-1=；对待求量、做消元变换得，则得子方程  的解=，=及相关矩阵的递矩阵：  如此类推，对待求量、……做消元变换后，在所得矩阵R（m）中的最后一列（最后一行除外）就是m个方程组的解，而前m列就是相关矩阵的递矩阵。 在逐步回归分析中，利用这种每一步变换正好提供偏回归系数的性质，可从初选的m个自变量中逐个引入作用显著的自变量，建立相应的偏回归方程。 2.对同一待求量实行重复的消元变换，其效果相当没有变换。例如，对做了的消元变换得，如果对再做的消元变换得，则。一般来说，已对再做的消元变换为，如对再做变换得，同样有这一性质可从式（23）直接得到证明： =     可见对一切ij都有，所以。 在逐步回归分析中，利用这一复原性质，可把误差引入而作用变得不显著的自变量从回归方程中剔出来，把相关矩阵复原到对未做消元变换前的状态，再从这个状态开始引入新变量。 3.根据每一步变换的矩阵元素，可以直接计算自变量的偏回归平方和并可检验其显著性。设在矩阵中，已把P个自变量引入回归方程，则容易证明（从略），P个自变量的剩余平方和可用的右下角元素表示：  (25)  根据这个关系，P个自变量的回归平方和可表示为：  先取对（对应自变量Xk）的消元变换，得。设这一变换是引入自变量Xk，则得P+1个自变量的回归平方和  显然＞。取二者之差，得自变量Xk的偏回归平方和，记为： == 或取下列比值：  （8-33） 这里的，表示引入或剔出一个自变量Xk时引起回归平方和的增减量，称为Xk的偏回归平方和。其值越大，说明Xk在回归中起的作用越大，而且可用变换前的元素计算。 自变量Xk在回归中作用的显著性，可按已知的统计量（4～27）式做检验。如为剔除变量，利用（25）和（26）式，则得  （8-34） 这时的=1。如为引入变量，考虑到自变量增加，剩余平方和将减少，故得  （8-35） 根据（25）式，也可以直接计算p个自变量的复相关系数和剩余均方差：  （8-36）  （8-37） 4.矩阵具有如下的对称性和反对称性：即  这说明，计算的值有正有负，标志着自变量Xk不在或已在偏回归方程中。即当＞0时，表示Xk尚未引入回归方程，如果作用显著则是应当引入的变量；当＜0时，表示Xk已经引入回归方程，如果作用不显著，则是应该剔除的变量。 具备上述知识后，可按下列步骤进行逐步回归分析： 第一阶段，准备工作。 首先，按问题性质和所处的水文地质条件，考虑一切可能的相关因素。初选因变量和自变量的个数。根据这些变量的长观资料和试验资料，计算均值、均方差和各变量之间的单相关系数，组成初始康德相关矩阵。通常是一对称阵，只要计算它的上三角部或下三角部就可以了。 其次，根据可能选入回归方程的自变量的个数p，规定一个临界值Fα（1，n-p-1），作为检验自变量在回归中的作用是否显著的标准。Fα和要求的检验水平α有关，为使最终的回归方程包含较多的自变量，α不宜取得过高。Fα还随自由度f余=n-p-1发生变化，随着引入或剔除自变量个数的增减，f余将不断变化。但从Fα的分布表中可以看出，在n＞p的条件下，对于一定的α来讲，Fα只随f余的变化略有增减。因此，从较保险的角度出发，可以近似规定一个常值Fα作为标准。例如，在a=0.05条件下，如果n＞15，在初选的全部自变量中可能有2-3个引入回归方程时，则可取f回=1，f余=13，查附表（Ⅰ）得Fα=4.67。 第二阶段，逐步回归。 第一步:考虑引入第一个自变量。为此，根据R（0）中的元素，按（8-33）式计算所有自变量的偏回归平方和（k=1，2，…………m）。 要引入的自变量应是值最大的自变量，它能使剩余平方和减少得最多，因而作用最大。例如max=，对应自变量X2。对它按（8-35）式做显著性检验得： F1= （p=0）， 如果F1＞Fα，证明X2作用显著，则可以引入回归方程。一般来讲，这第一个自变量总是可以引入的。 取对做消元变换，得。由中得解，再按（8）式和（10）式求出和，结果得关于变量的一元偏回归方程：。 回归方程的相关系数和剩余均方差，可按（8-36）式和（8-37）式计算，式中的，P=1。 第二步:考虑引入第二个自变量。根据中的元素按（8-33）式计算所有自变量的偏回归方程和。 根据反对称性，中只有一个为负值，它对应第一步引入的自变量，不必考虑剔除变量问题。 从的正值中选一最大值。例如，对应自变量。对它做F检验，按（8-35）式得：  P=1 如果＞，证明变量可以引入回归方程。 取对做消元变换，得。由中得解，，再计算对应的、、，结果得二元偏回归方程：  这个方程的复相关系数和剩余均方差，仍按（8-36）式和（8-37）式计算，式中的，P=2。 第三步:考虑剔除变量和引入第三个自变量。根据中的元素按（8-36）式计算所有自变量的偏回归平方和。 根据反对称性，中的和都是负值，同已经引入的两个自变量和相对应。单独检验和时对因变量的作用是显著的，但合起来一个变量的作用就可能被另一个变量所代替，而变的不显著，这时应当随时剔除它。 剔除变量时，无须对已经引入回归方程的所有变量进行检验。可找出这些变量中偏回归平方和的绝对值最小者做检验，如果作用显著，则证明所有引入的自变量都不用剔除。否则，应首先剔除偏回归平方和绝对值最小的那个变量。因此，从的负值中找出绝对值最小者，例如，对应自变量。按（8-34）式做F检验：  P=2 如果＞,则引入。取取对做消元变换，得。由此得解、、，进而求得、、、，结果得三元偏回归方程：  复相关系数和剩余均方差，可按（8-36）式（8-37）式计算，式中的。P=3。 重复上述步骤，可以继续剔除旧变量引入新变量。一般讲，在R（0）经过1次消元变换得R（0），已经引入p个自变量，建立了p元回归方程。 根据R’元素，按（8-38）式计算偏回归。其中有p个负值，同已经引入的自变量相对应。其余为正值，属于尚为云如的自变量。从的负值中可以找出绝对值最小者，按（8-34）式做检验。如果作用显著，则所有已经引入的自变量不必剔除。 从的正值中找出最大值，按（8-35）式做检验。如果作用显著，应引入新变量。做消元变换LP+1，R（l）=R(l+1),，解出新的偏回归系数，建立P+1元回归方程。如果作用不显著，则逐步回归就此结束，已经建立的P元回归方程  就是要求的最优的回归方程。它包含了一切作用显著的自变量，不显著者已经随时剔除。为了进一步检验方程的回归效果，还要对整个回归方程进行F检验。先按（8-36）式和（8-37）式求得复相关系数R和剩余均方差σ余。再按下式计算  如果F＞Fα，则说明因变量和p个自变量之间的相关都是显著的，因而预报效果一定教好。 实例：娘子关泉流量的主动部回归分析。该泉的泉域面积约3601Km2，其中可溶岩的分布面积为1675 Km2，主要含水层厚度约470～630m。在这种水文地质条件下，泉水直接或间接都是接受降水入渗补给。而且，从泉水流量的继承性可以想象泉流量不仅受当年的降水量影响，也受前一年、前两年……直至前n-1年的降水量影响。泉域的调节能力越强，这个前往逆推的年数就越多，逆推年数称为调节期或滞后期。为了分析泉流量和调节期内各年降水量的相关关系，列入表（8-15）中，预报的调节期为9年（包括当年在内）。 首先，确定泉流量和可能选入回归方程的自变量，按前面分析，预计泉流量Qt可能受当年的和前8年的降水量影响。设当年降水量Pt的影响比为a0；前一年降水量Pt-1的影响为a1；前二年降水量Pt-2的影响比为a2，如此类推，前8年降水量Pt-8的影响比为a8。如取a0+ a1+……+ a8=1，则得 a0Pt+ a1Pt-1+ a2Pt-2+……+ a8Pt-8 这里的为引入降水量，是指各年降水量对泉流量影响的总合。如果近似认为和 表（8-15）  年份 1955 1956 1957 1958 1959 1960 1961 1962  年平均降水量（mm） 514 690 392 570 707 568 593 632  年平均泉流量（m3/s） — — — — — — — —  年份 1963 1964 1965 1966 1967 1968 1969 1970  年平均降水量（mm） 802 787 399 696 605 445 615 539  年平均泉流量（m3/s） 13.6 15.8 15.2 14.5 14.8 14.4 13.3 12.3  年份 1971 1972 1973 1974 1975 1976 1977 1978  年平均降水量（mm） 595 255 703 378 535 594 722 481  年平均泉流量（m3/s） 12.4 11.3 10.8 10.7 10.1 10 11.5 11.8  注:根据山西省地质局水文地质队资料   表（8-16） 序号 年份 项目 Qt Pt Pt-1 Pt—2 Pt—3 Pt—4 Pt—5 Pt—6 Pt—7 Pt—8  1 1963 12.6 802 632 593 568 707 570 392 690 514  2 1964 15.8 787 802 632 593 568 707 570 392 690  3 1965 15.2 399 787 802 632 593 568 707 570 392  4 1966 14.5 696 399 787 802 632 593 568 707 570  5 1967 14.8 605 696 399 787 802 632 593 568 707  6 1968 14.4 445 605 696 399 787 802 632 593 568  7 1969 13.3 615 445 605 696 399 787 802 632 593  8 1970 12.8 539 615 445 605 696 399 787 802 632  9 1971 12.4 595 539 615 445 605 696 399 787 802  10 1972 11.3 255 595 539 615 445 605 696 399 787  11 1973 10.8 703 255 595 539 615 445 605 696 399  12 1974 10.7 378 703 255 595 539 615 445 605 696  13 1975 10.1 535 378 703 255 595 539 615 445 605  14 1976 10 594 535 378 703 255 595 539 615 445  15 1977 11.5 722 594 535 378 703 255 595 539 615  16 1978 11.8 481 722 594 535 378 703 255 595 539   ∑ 203 9151 9302 9173 9147 9319 9511 9200 9635 9554   平均值 12.688 57.198 581.38 573.313 571.688 582.438 594.438 575 620.188 597.12    7.317 594.167 589.023 517.198 510.811 583.901 544.973 565.814 468.129 475.878  Qt呈线性相关，则得～Qt的关系式是： Qt=b0+bPt 将前式带入上式，整理后得 Qt=b0+b1Pt+ b2Pt-1+ b3Pt-2+……+ b9Pt-8 结果，问题归结为一个9元线性回归方程，其中b0、b1、……b9为10个待定的回归系数。 其次，按已给资料计算待定系数，构成正规方程组。按上式需要进行观测时对应Q1的降水量应为P1、P1-1、P1-2、P1-3……P1-8；对应Q2的降水量为P2、P2-1、P2-2、……P2-8；如此类推，对应Q16的降水量为P16、P16-1、P16-2……P16-8等。按这个规律可把表（8-15）中的资料重新组合成适合计算用的数据，如表（8-16）。根据表（8-16）的数据计算各变量的总和及平均值，进一步求出包含未知量bj（j=1、2、……9）的正规方程组，结果列入表（8-17）中。 最后，实行逐步回归，初选的九个自变量对泉流量的作用是否都显著？这要通过逐步回归分析来确定。因此，按（8-24）式的方法。 把正规方程组标准化，然后组成相关矩阵R（0），结果如表（8-18）。根据泉域条件，检验用的临界值取为Fα=3 第一步:考虑引入第一自变量。为此，按R（0）元素计算初选自变量的偏归平方和VK得 V1=rQ1（0）r1Q（0）/r11（0）=0.05886， V4=0.06829 V7=0.03485 V2=0.21536 V5=0.21793 V8=0.00001 V3=0.19660 V6=0.12965 V9=0.00093 所有VK均为正值，表明全部自变量尚为引入回归方程。VK中的最大值为V5，对应自变量为Pt—4。做F检验： F1=V5/（νQ（1）Q-V5）（n-2）=3.90124＞Fα=3 证明Pt-4的作用显著，可以引入回归方程。取R（0）对做消元变换，得L5R（0）=R（1）（见表8-19），因此得解： =r5Q（1）=0.46683，再按（10）式和（8）式计算：   结果得一元偏回归方程为  这时的复相关系数和剩余均方差分别为   表（8-17） 正规方程组  Si1 b1 + Si2 b2 + Si3 b3 + Si4 b4 + Si5 b5 + Si6 b6 + Si7 b7 + Si8 b8 + Si9 b9 = Si Q  353034.59 b1 81343.63 b2 48048.31 b3 18584.69 b4 39853.44 b5 -70724.56 b6 -54156.00 b7 65977.19 b8 -57312.88 b9 1054.79  -81343.63 b1 346947.8 b2 -65883.88 b3 49855.88 b4 21752.38 b5 64143.38 b6 -70076.00 b7 -80454.13 b8 60791.25 b9 2000.08  48048.31 b1 -65883.88 b2 326267.4 b3 -69820.44 b4 59150.82 b5 51904.81 b6 44016.0 b7 -20231.94 b8 -81121.63 b9 1853.16  18584.69 b1 49855.88 b2 -65883.88 b3 325825.4 b4 -69580.81 b5 62582.19 b6 50229.00 b7 50731.94 b8 -19746.38 b9 1091.40  39853.44 b1 21752.38 b2 59160.82 b3 -69580.81 b4 340939.9 b5 -80748.06 b6 42516.00 b7 49373.69 b8 39406.13 b9 1994.49  -70724.56 b1 64143.38 b2 51904.81 b3 62582.19 b4 -80748.06 b5 296996.9 b6 -52304.00 b7 -29554.31 b8 49973.13 b9 1435.79  -54156.00 b1 -70076.00 b2 44016.0 b3 50229.00 b4 42516.00 b5 -52304.00 b6 320146.0 b7 -32389.00 b8 -12496.00 b9 772.90  65977.19 b1 -80454.13 b2 -20231.94 b3 50731.94 b4 49373.69 b5 -29554.31 b6 -32389.00 b7 21914434 b8 -42628.38 b9 -12.36  -57312.88 b1 60791.25 b2 -81121.63 b3 -19746.38 b4 39406.13 b5 49973.13 b6 -12496.00 b7 -42628.38 b8 226459.8 b9 106.23   表（8-18） 初始相关矩阵R（0）  ri1（0） ri2（0） ri3（0） ri4（0） ri5（0） ri6（0） ri7（0） ri8（0） ri9（0） rQ（0）  1.00000 -0.23243 0.14157 0.54800 0.11437 -0.21842 -0.16109 0.23720 -0.20270 0.24262  -0.23243 1.00000 -0.19582 0.14823 0.06325 0.19982 -0.21026 -0.29178 0.21688 0.46407  0.14157 -0.19582 1.00000 -0.21414 0.17738 0.16674 0.13619 -0.07566 -0.29844 0.44340  0.54800 0.14828 -0.21414 1.00000 -0.20877 0.20118 0.15552 0.18986 -0.07269 0.26132  0.11487 0.06325 0.17738 -0.20877 1.00000 -0.25376 0.12869 0.18063 0.14182 0.46683  -0.21842 0.19982 0.16674 0.20118 -0.25376 1.00000 -0.16962 -0.11585 0.19269 0.36007  0.10109 -0.21026 0.13619 0.15552 0.12869 -0.16962 1.00000 -0.12228 -0.04644 0.18669  0.23720 -0.29178 -0.07566 0.18986 0.18063 -0.11585 -0.12228 1.00000 0.19135 0.00361  -0.20270 0.21688 -0.29844 -0.07269 0.14182 0.19269 -0.04644 0.19135 1.00000 0.30510  0.24262 0.46407 0.44340 0.26132 0.46683 0.36007 0.18669 0.00361 0.03051 1.00000   表（8-19） L5R（0）=R（1）  ri1（1） ri2（1） ri3（1） ri4（1） ri5（1） ri6（1） ri7（1） ri8（1） ri9（1） riQ（1）  0.98681 -0.23970 0.12119 0.07878 0.11487 -0.18928 -0.17587 0.21645 -0.21899 0.18899  -0.23970 0.99600 -0.20704 0.16149 0.06325 0.21587 -0.21840 -0.30321 0.20791 0.43454  0.12119 -0.20704 0.96854 -0.17711 0.17738 0.21175 0.11336 -0.10770 -0.32060 0.36059  0.07878 0.16149 -0.17711 0.95642 0.20877 0.14820 0.18239 0.22757 -0.04308 0.35878  0.11487 0.06325 0.17738 -0.20877 1.00000 -0.25376 0.12869 0.18063 0.14182 0.46688  -0.18927 0.21587 0.21175 0.14820 0.25376 0.93561 -0.13696 -0.07001 0.22868 0.47853  -0.17587 -0.21840 0.11336 0.18239 -0.12869 -0.13696 0.98344 -0.14553 -0.06469 0.12661  0.21645 -0.30321 -0.10770 0.22757 -0.18063 -0.07001 -0.14553 0.96737 0.16573 -0.08793  -0.21899 0.20791 -0.32060 -0.04308 -0.14182 0.22868 -0.06469 0.16573 0.97989 -0.03570  0.18899 0.43454 0.36059 0.35878 -0.46683 0.47853 0.12661 -0.08793 -0.03570 0.78207   表（8-20） L6R（1）=R（2）  ri1（2） ri2（2） ri3（2） ri4（2） ri5（2） ri6（2） ri7（2） ri8（2） ri9（2） riQ（2）  0.94852 -0.19603 0.16403 0.10876 -0.06353 0.20231 -0.20358 0.20229 -0.17273 0.28580  -0.19602 0.91319 -0.25590 0.12729 -0.12180 -0.23073 -0.18680 -0.28706 0.15515 0.32413  0.16403 -0.25590 0.92062 -0.21065 -0.23481 -0.22632 0.14436 -0.09086 -0.37536 0.25229  0.10876 0.32730 -0.21065 0.93295 0.16858 -0.15840 0.20408 0.23866 -0.07930 0.28298  0.06354 0.12180 0.23481 -0.16858 1.06883 0.27122 0.09154 0.16164 0.20384 0.59662  -0.20230 0.23073 0.22632 0.85840 0.27122 1.06882 -0.14639 -0.07483 0.24443 0.51146  -0.20358 -0.18680 0.14436 0.20408 -0.09154 0.14639 0.96339 -0.15578 -0.03122 0.19600  0.20229 -0.28706 -0.09086 0.23866 -0.16164 0.07483 -0.15578 -0.03122 0.18284 -0.05212  -0.17273 0.15515 -0.37536 -0.07930 -0.20384 0.24442 -0.03121 0.18284 0.92400 -0.15266  0.28580 0.32413 0.25229 0.28298 -0.59662 -0.51146 0.19666 -0.05212 -0.15266 0.53732   第二步：考虑引入第二个自变量。根据R（1）元素计算VK值，得 V1=0.03619 V4=0.13459 V7=0.01630 V2=0.18958 V5=-0.21793 V8=0.00799 V3=0.13425 V6=0.24475 V9=0.00130 VK值中只有一个负值，对应刚才引入的Pt-4，故不必剔除变量。值的VK最大者为V6，对应自变量Pt-5,作F检验。  决定引入Pt-5，取R（1）对b6ˊ做消元变换，L6R（1）=R（2），见表（8-20）得解：，，继而求得，，  于是得二元偏回归方程 这时的复相关系数和剩余均方差分别为  第三步：考虑剔除变量和引入第三个自变量。按R（2）中的元素计算VK值得： V1=0.08611 V4=0.08583 V7=0.04015 V2=0.11104=max V5=-0.23303 V8=0.000282 V3=0.06914 V6=-0.24475=min V9=0.02522 先考虑剔除变量的可能性。从负值的两个VK中找出绝对值最小者V6对应变量Pt—5做F检验： 证明不必剔除变量Pt-5。再考虑可否引入新变量，正值的六个VK中最大者为V2=0.11104，对应变量Pt-1、检验： 应该引入变量Pt-1。在R（2）中对b2ˊ做消元变换，L2R（2）=R（3），见表（8-21）。 由此得解 ，，  ， 结果得三元偏回归方程  复相关系数和剩余均方差为  第四步：继续剔除变量和引入第四个自变量。按R（3）的元素计算VK值得 V1=0.13735 V4=0.06257 V7=0.07333 V2=-0.11103=min V5=-0.28392 V8=0.00244 V3=0.13573=max V6=-0.16620 V9=0.04717 先考虑剔除变量的可能性。在负值的三个VK中以为最小，对应变量Pt-1，做F检验  证明不必剔除。再考虑引入新变量的可能性。在正值的VK中，最大值为V3=0.12573，对应变量Pt—2。做F检验：  可以引入变量Pt-2。因此在R（3）中对b2’做消元变换，L2R（3）=R（4），见表（8-22）。 由此得解 ，， ，继而求得：，，， ，  结果得四元偏回归方程： 表（8-21） L2R（2）=R（3）  ri1（3） ri2（3） ri3（3） ri4（3） ri5（3） ri6（3） ri7（3） ri8（3） ri9（3） riQ（3）  0.90701 0.20517 0.10102 0.13313 -0.08877 0.65450 -0.24224 0.14282 -0.14059 0.35295  -0.20717 1.05687 -0.27046 0.13454 -0.12873 -0.24385 -0.19843 -0.30339 0.16397 0.34256  0.11102 0.27946 0.33144 -0.17622 -0.26775 -0.28372 0.09840 -0.26950 -0.33340 0.33995  0.13513 -0.13454 -0.17622 0.91582 0.16497 -0.12736 0.22921 0.27728 -0.10017 0.23937  0.08877 -0.12873 0.26775 -0.18497 1.08457 0.30092 0.11559 0.19859 0.18387 0.59490  -0.15450 -0.24385 0.28872 0.12736 0.31192 1.12508 -0.10084 -0.00483 0.20659 0.43242  -0.24228 0.19742 0.09384 0.22921 -0.11559 0.10084 0.92651 -0.21245 -0.00059 0.26065  0.14282 0.30339 -0.16950 0.27728 -0.19859 0.00483 -0.21245 0.37505 0.22991 0.04622  -0.14059 -0.16397 -0.33340 -0.10017 -0.18387 -0.20659 -0.00059 0.22991 0.89856 -0.20581  0.35295 -0.34256 0.33995 0.23937 -0.55440 -0.43242 0.26065 0.04622 -0.20581 0.42629   表（8-22） L3R（3）=R（4）  ri1（4） ri2（4） ri3（4） ri4（4） ri5（4） ri6（4） ri7（4） ri8（4） ri9（4） riQ（4）  0.89253 0.17191 -0.13039 0.15811 -0.05386 0.19215 -0.25452 0.16492 -0.09712 0.30862  -0.17191 1.14278 0.31765 0.07856 -0.21378 -0.33556 -0.16761 -0.35723 0.05807 0.45055  0.13039 0.31765 1.17448 -0.20697 -0.31447 -0.33910 0.11021 -0.19908 -0.39157 0.39927  0.15811 -0.07856 0.20697 0.87935 0.12956 -0.18712 0.24863 0.24220 -0.16917 0.30973  0.05386 -0.21378 -0.31447 -0.12956 1.16871 0.39171 0.08608 0.25190 0.28871 0.44800  -0.19215 0.33556 -0.33910 -0.18712 0.39171 1.22298 -0.13266 0.05265 0.31965 0.31714  -0.25452 0.16761 -0.11021 0.24863 -0.08608 0.13266 0.91617 -0.19377 0.03616 0.22318  0.16492 0.35723 0.19908 0.24220 -0.25189 -0.05265 -0.19377 0.14131 0.16354 0.11390  -0.09712 -0.05807 0.39157 -0.16917 -0.28871 -0.31965 0.03616 0.16354 0.76801 -0.07270  0.30862 -0.45055 -0.39927 0.30973 -0.44800 -0.31714 0.22318 0.11390 -0.07270 0.29056   复相关系数和剩余均方差为， 第五步：继续考虑剔除变量的可能性，从负值的VK中取绝对值最小者为，对应变量Pt-5。做F检验： 证明不必剔除。再考虑引入新变量问题。在正值的VK中，最大值V4=0.10910，对应变量Pt-3。做F检验： 可以引入变量Pt-3。在R（4）中对做消元变换，L4R（4）=R（5），见表（8-23）。得解： ，， ， 继而求得：，， ， ， 结果得五元偏回归方程： 复相关系数和剩余均方差为， 第六步:继续考虑剔除变量和引入第六个自变量。按R（5）的元素计算VK值得 V1=0.07404=max V4=-0.10910 V7=0.02174 V2=-0.15553 V5=-0.20515 V8=0.00106 V3=-0.18227 V6=-0.05826=min V9=0.00023 考虑剔除变量，取负值VK中的绝对值最小者为，对应变量为Pt—5。做F检验：不必剔除变量。再考虑引入变量，在正值的VK中，最大值为V1=0.07404，对应变量Pt。做F检验： 可以引入Pt0在R（5）中对b1’做消元变换，L1R（5）=R（6），见表8-24。于是得解： b1’=0.29272 b2’=0.47734 b3’=0.42311 b4’=0.29960 b5’=0.47104 b6’=0.33733. 于是求得 ，    表（8-23） L4R（4）=R（5）  ri1（5） ri2（5） ri3（5） ri4（5） ri5（5） ri6（5） ri7（5） ri8（5） ri9（5） riQ（5）  0.86408 0.18604 -0.16760 -0.17980 -0.07716 0.22580 -0.29922 0.12137 -0.06670 0.25293  -0.18604 1.14980 0.29916 -0.08934 -0.22536 -0.31884 -0.18982 -0.37887 0.07318 0.42288  0.16760 0.29916 1.22319 0.23537 -0.28398 -0.38314 0.16870 -0.14207 -0.43139 0.47217  0.17980 -0.08934 0.23537 1.13720 0.14734 -0.21279 0.28274 0.27543 -0.19238 0.35223  0.07716 -0.22536 -0.28398 0.14734 1.18780 0.36414 0.12271 0.28759 0.26379 0.49363  -0.22586 -0.31884 -0.38310 -0.21279 0.36414 1.26280 -0.18557 0.00111 0.35565 0.27123  -0.29922 0.18982 -0.16873 -0.28274 -0.12271 0.18557 0.84587 -0.26225 0.08399 0.13561  0.12137 0.37887 0.14207 -0.27543 -0.28759 -0.00111 -0.26225 0.77460 0.21014 0.02895  -0.06670 -0.07318 0.43139 0.19238 -0.26379 -0.35565 0.08399 0.21014 0.73547 0.01311  0.25293 -0.42288 -0.47217 -0.35223 -0.49363 -0.27123 0.13561 0.02859 -0.01311 0.18146   表（8-24） L1R（5）=R（6）  ri1（6） ri2（6） ri3（6） ri4（6） ri5（6） ri6（6） ri7（6） ri8（6） ri9（6） riQ（6）  1.15730 0.21530 -0.19396 -0.20808 -0.08930 0.26132 -0.34630 0.14046 -0.07719 0.29272  0.21530 1.18986 0.26308 -0.12805 -0.24197 -0.27022 -0.25424 -0.35274 0.05882 0.47734  -0.19396 0.26308 1.25570 0.27025 -0.26901 -0.42694 0.22677 -0.16561 -0.41845 0.42311  -0.20808 -0.12805 0.27025 1.17461 0.16340 0.25978 0.34500 0.25018 -0.17850 0.29960  -0.08930 -0.24197 -0.26901 0.16340 0.19469 0.34398 0.14943 0.27675 0.26975 0.47104  0.26132 -0.27022 -0.42694 -0.25978 0.34398 1.32181 -0.26376 0.03283 0.33822 0.33733  0.34630 0.25424 -0.42677 -0.34500 -0.14943 0.26376 0.74225 -0.22022 0.06089 0.22320  -0.14046 0.35274 0.16561 -0.25018 -0.27675 -0.03283 -0.22022 0.75756 0.21951 -0.00694  0.07719 -0.05882 0.41845 0.17850 -0.26975 -0.33822 0.06089 0.21951 0.73032 0.00641  -0.29272 -0.04773 -0.42311 -0.29960 -0.47104 -0.33733 0.22320 -0.00694 0.00641 0.10742     结果得六元偏回归方程：    第七步:继续剔除和引入变量。按R（6）的元素计VK值。得 V1=-0.07404=min V4=-0.07642 V7=0.06712=max V2=-0.19150 V5=-0.13572 V8=0.000064 V3=-0.14488 V6=-0.08609 V9=0.000056 考虑剔除变量。负值VK中最小者为，对应变量Pt，检查：  证明不必剔除。再考虑引入变量。在正值VK中的最大值为V7=0.06712，对应变量Pt-6检验：  可以引入变量Pt-6。在R（6）中对做消元变换，L7R（6）=R（7），见表（8-25）。得解 bˊ1=0.39686, bˊ2=0.55379 , bˊ3=0.35492, bˊ4=0.19586 ,bˊ5=0.42611 , bˊ6=0.41665 ,bˊ7=0.30071, 于是求得         最后得七元偏回归方程：  R=0.97964 σ余≈0.51933 第八步:继续剔除和引入变量。按R（7）的元素计算VK值，得： V1=-0.11942 V4=-0.02874=min V7=-0.06712 V2=-0.24017 V5=-0.14825 V8=0.00508=max V3=-0.09507 V6=-0.12264 V9=0.000195 在负值VK中绝对值最小者为， 对应变量Pt-3。检验得无剔除。在正值中最大值为V8=0.00508，对应变量Pt-7。检验得 证明初选的变量中已经没有再能引入的变量，逐步回归就此结束。方程（5）就是要求的最优回归方程。实际调节期为7年（包括当年在内）。整个方程的显著性为  证明因变量和引入自变量之间的相关都是显著的，预报效果不会差。 值得指出的是，如果不经这样筛选，而直接对Q求关于9个变量的全回归，可得回归方程  这时的复相关系数和剩余均方差为 R=0.98305 σ余=0.54764 结果说明，自变量增加两个，R增大量很小，而σ余不仅没有减少，反而由于自由度变小而增大，预报效果变差。这就是逐步回归比普通回归优越的地方。