1  引言

流域水文系统因其输入输出和系统输送既受确定性因素的作用又受随机性因素的影响，因此流域水文过程的变化十分复杂，其精确描述显然是十分困难的。在当前一种既考虑水文过程形成的物理背景又不苛求预报量与预报因子间确切的运行机制的系统建模思想比较符合人们对客观水文规律认识的现实，因而得到了众多学者和专家的认可。多维混合回归系统模型是将回归与自回归有机结合起来的混合回归模型，它既考虑了预报量自身的演变规律又充分利用了预报量对预报因子的依赖关系，预报模型具有一定的物理基础又不过分强调每一个参数的物理意义^[1~3]。其建模方式灵活，使用方便，精度较高。

2       多维混合回归系统水文模型的结构和算法

2．1 结构^[1]

设流域水文系统的输出为{Z}，输入为{y₁}，{y₂}，…，{y_s}，则经过流域输送作用后，其系统描述可表示为

     Z_t+1=b₀⁽⁰⁾+ b₁⁽⁰⁾Z_t + b₂⁽⁰⁾Z_t-1+ … + b_p0⁽⁰⁾Z_t-p0

              + b₁⁽¹⁾y_1,t+ b₂⁽¹⁾y_1,t-1 + … + b_p1⁽¹⁾y_1,t-p1

            + b₁⁽²⁾y_2,t + b₂⁽²⁾y_2,t-1+ … + b_p2⁽²⁾y_2,t-p2

                     +… + b₁^(s)y_s,t + b₂^(s)y_s,t-1+ … + b_ps^(s)y_s,t-ps+ε_t+1            (1）

上式即为多维混合回归系统水文的结构。式中y_1,t ，y_2,t ，…，y_s,t分别为系统输入{y₁}，{y₂}，…，{y_s}在t时刻的数值，Z_t为系统输入{Z_t}在t 时刻的数值，{b⁽⁰⁾}，{b⁽¹⁾}，{b⁽²⁾}，…，{b^(s)}为系统响应函数。

2．2 算法

式（1）中包含了S+1个变量，是多维的。参数共有（p₀+p₁+p₂+…+p_s+1）个。实际上p₀为系统自回归阶数，p₁，p₂，…，p_s为因子{y₁}，{y₂}，…，{y_s}的阶数。模型阶数和参数识别算法如下。

设已知实测数据{Z_t}，{y_1,t}，{y_2,t}，…，{y_s,t}，其中 t=1，2，…，m。令p=max{p₀，p₁，…，p_s}，在式（1）中分别取t=p+1，p+2，…，m，由此可得m-p 个等式，其缩写形式为

     Z=AB+ε                                              ( 2 )

式中

         Z=[Z_p+1 Z_p+2  …  Z_m]^T

        B = [b₀ b₁  b₂  …  b_s]^T

              b₀ = [b₀⁽⁰⁾  b₁⁽⁰⁾  b₂⁽⁰⁾  …  b_p0⁽⁰⁾]^T

        b_I = [b₁⁽ⁱ⁾  b₂⁽ⁱ⁾   b₃⁽ⁱ⁾  …  b_pi⁽ⁱ⁾]^T   (i=1，2，…，s)

        ε= [ε_p+1  ε_p+2  ε_p+3  …ε_m]^T

        A = [A₀  A₁  A₂  …  A_s]^T

                       Z_p   Z_p-1   …  Z_p-p0+1

              Z_p+1   Z_p    …   Z_p-p0+2

        A₀=   …   …    …   …

              Z_m-1  Z_m-2  …    Z_m-p0

              y_i,p    y_i,p-1   …  y_i,p-pi+1

              y_i,p+1   y_i,p    …  y_i,p-pi+2

        A_i =   …    …    …   …      ( i=1，2，…，s)

.                                 

              y_i,m-1  y_i,m-2   …   y_i,m-pi

    根据最小二乘原理，B的最小二乘估计为

         B=(A^TA)^-1A^TZ                                          (3 )

相应的残差平方和与以上的 有关，也有各变量的阶数有关，其表达式为

        S（B，p₀，p₁，…，p_s）= Z^TZ – Z^TA(A^TA)^-1A^TZ            ( 4 )

在p₀，p₁，…，p_s 未知的情况下，可按赤池信息准则（AIC准则）确定各自的阶数值。AIC准则函数为

        AIC(p₀，p₁，…，p_s)= log S（B，p₀，p₁，…，p_s）

                            + 2(p₀，p₁，…，p_s）/m              (5 )

在给定某个备选的最大延迟量 时，在所有可能的符合以下不等式的值中找出使式（5）中的AIC最小的p₀，p₁，…，p_s。

          0 ≤ p₀，p₁，…，p_s≤ p

如此选取出的p₀，p₁，…，p_s即为式（1）中各变量的阶数。

具体做法是先定最大延迟量p，然后在式（5）范围内选取一组（p₀，p₁，…，p_s），以这些阶数为基础，用最小二乘法估计出参数后用式（4）计算出残差平方和，并由式（5）计算AIC值。然后再选一组阶数，重复同样的计算步骤，比较AIC的大小，直至找出使AIC最小的阶数。

参数估算和阶数确定是相互联系，互为前提的，可编程上机一次实现。当S较大，参数个数很多时，为减少计算工作量提高模型的统计稳定性，在实际应用时可采用疏系数方法求解。 

3        朗梨站洪峰水位多维混合回归预报模型

3．1       基本情况

朗梨站系湘江一级支流浏渭河的总控制站，控制流域面积3815km²，干流长度110km，平均坡降0.27‰。其上游设有双江囗水文站，控制面积2067km²，区间有江背报汛雨量站。浏渭河在朗梨站下游约25km处汇入湘江。朗梨站洪峰主要由上游流域的暴雨形成，同时也受湘江洪水顶托。由于洪水波运动的非线性特征加之下游严重的回水顶托影响，使得朗梨站的水位流量关系异常复杂，相似的上游来水和区间降水在朗梨站形成的洪峰水位相距甚远。用传统的上下游水位相关法建立洪峰水位预报方案合格率不高，不能正常投产，作业预报的难度很大。

   为此，特考虑建立该站洪峰水位多维混合回归预报模型。

3．2       模型结构的改进和系统输入的非线性化

根据朗梨站洪峰水位预报的实际需要，系统以朗梨站洪峰水位为输出，以双江囗站出峰前相关站的水位和区间降水为输入。其混合回归预报模型的结构在式（1）的基础上略加改进，结构形式为

    ^朗Z_{m , t+τ}= b₀+ b₁^{(0) 朗}Z_t+ b₂^{(0) 朗}Z_t-1+… +b_p0^{(0) 朗}Z_t-p0

                  + b₁^{(1) 双}Z_m,t+ b₂^{(1) 双}Z_t-1+… +b_p1^{(1) 双}Z_t-p1                     

                  + b₁^{(2) 长}Z_t+ b₂^{(2) 长}Z_t-1+… +b_p2^{(2) 长}Z_t-p2

                  + b₁^{(3) 双}p’_t+b₂^{(3) 双}p’_t-1+… +b_p3^{(3) 双}p’_t-p3

                  + b₁^{(4) 江}p’_t+b₂^{(4) 江}p’_t-1+… +b_p4^{(4) 江}p’_t-p4

                  + a₁s₁ + a₂s₂ + a₃s₃+ … + a_ns_n +ε_t+τ          ( 7 )

式中：^朗Z_t为t时刻朗梨站的水位；τ为洪峰从双江囗至朗梨站的传播时间，也是本模型的预见期；^双Z_m,t 为 t 时刻双江囗站的水位，下标加注m 表示其值同时为洪峰水位；^长Z_t 为t 时刻长沙站的水位；^双p’_t为(t-1)~t时刻双江囗站的降水量非线性化因子；^江p’_t 为(t-1)~t时刻江背站的降水量非线性化因子；s_i 为第i个河段的水力特性因子；b_j^(k),a_i 为系统响应。

由于降水径流和河段水力特性的非线性影响，对系统的降水和河段水力特性输入因子均作非线性处理。即

^l          对降水输入取 p’=p^α(其中α分别取 0.1,0.2,0.3,…,2.0)，按AIC准则确定采用值；

^l          对河段水力特性因子主要考虑双江囗站出峰时双江囗~朗梨，朗梨~长沙的水面落差，由曼宁公式不难推知河道流量与水面比降的平方根成正比，故模型固定取

           s₁= (^双Z_m,t  - ^朗Z_t +α )^(1/2)

              s₂=(^朗Z_t   - ^长Z_t  +β)^(1/2)

其中α和β为基面换算系数，以保证s₁和 s₂大于零。

3．3       模型的建立

因式（7）较之式（1）有较大改进，在资料处理上已将连续水文系统截断为以单次洪峰为单位的多个序列资料，以之分别代入式（7）根据最小二乘法并结合AIC 准则即可得到最优预报模型。

采用朗梨站1965，1973~1995年共34 次洪峰资料，并给定最大延迟量p=5，得到洪峰水位最优预报模型为

    ^朗Z_{m , t+τ}=  45.7604– 0.6304^朗Z_t – 0.6913 ^朗Z_t-1+ 0.2188^朗Z_t-2

                 +1.8674^双Z_m,t +0.2970^双Z_t-1 - 0.5683^双Z_t-2 +0.0007 ^双Z_t-3

                + 0.0409^长Z_t  +0.3240 ^长Z_t-1+0.1854^长Z_t-2

                +0.0588 ^双p’_t +0.0458 ^双p’_t-1+ 0.0713^江p’_t

                 - 7.8302 s₁+0.5204 s₂

式中：      ^双p’_t =(^双p_t )^0.5 ;

           ^双p’_t-1=(^双p_t-1 )^0.5 ;

           ^江p’_t=(^江p_t)^0.5 ;

           s₁= (^双Z_m,t  - ^朗Z_t + 16.0)^(1/2);

           s₂= (^朗Z_t   - ^长Z_t  + 3.4 )^(1/2);

           △t=6h

   峰现时间的预报采用同期资料建立双江囗站洪峰水位与洪峰传播历时相关图，即（^双Z_m,t~τ）经验相关图查得。限于篇幅（^双Z_m,t~τ）关系图从略。

3．4模型的检验与评价

根据《水文情报预报规范》，模型的评定指标主要有合格率和有效性两项^[4]。

            合格率 = 合格场次/总场次×100%

模型率定期内34 场洪水洪峰水位变幅6.37m，模型拟合均方σ_△为1.61m，评定时允许误差最大值取1.00m，最小值取0.30m，单次洪水洪峰水位允许预报误差据预见期水位变化幅度（△H）在0.30~1.00m间内插确定，允许误差为

△_允许= 0.30 +(0.70△H)/6.37

有效性评定采用效率系数

       d_△ = 1 - S_△²/σ_△²

其中

           S_△² = ∑(△_i-△）²/n       σ_△² =∑(△_I -△）²/n

并且以1966~1971年共10 场洪水对模型作验证预报。率定期 34 场洪水和验证期 10 场洪水的回顾与检验预报精度的综合评价详见表1。

   表1             模型预报精度的综合评价

          单次预报误差幅度   确定性系数  单次预报优等率  单次预报良好率 单次预报合格率

资料分段    |△|_min~|△|_max          d_△              (%)         (%)         (%)

率定期      0~0.67           0.97          41.1            64.7            94.1

验证期    0.02~0.59          0.93          50.0           70.0            90.0

模型可综合评定为优级，可用于正式预报作业。

2．5 实际预报作业

模型用于1997年和1998年三次洪水预报，预报结果列于表2。

     表2              模型作业预报精度评定

  洪  号    实测洪峰水位（m）   预报洪峰水位(m)   预报误差(m)    允许误差(m)     评定

 970608         38.65              38.75            -0.10          0.76          优秀    

 980614         36.58              36.44            0.14           0.45          良好    

 980627         40.23              40.18            0.05           0.42          优秀

  三次预报二优一良，可见预报效果是令人满意的。

参考文献

1．李纪人，刘德平。水文时间序列模型及预报方法。南京：河海大学出版社，1991：40~81

2．李正最，吴雅琴。混合门限自回归模型在长期水文预报中的应用。江西水利科技，1993，19（4）：309~314

3．徐向阳。具有物理基础的太湖水位多元回归模型。河海大学学报，1990，18（6）：48~54

4．中华人民共和国水利电力部标准。水文情报预报规范。北京：水利电力出版社，1985：1~8