Flood forecasting error correction method based on Volterra series filter

Abstract

The present invention relates to the fields of hydrology and water resource, water conservancy and hydropower engineering, and emergency disaster prevention. Disclosed is a flood forecasting error cor

Description

[0003] At present, most of the models used in flood forecast are conceptual models. The models are characterized by centralized inputs and lumped parameters, which is contradictory to the basin objects of scattered inputs and distributed parameters. That is to say, after the spatial average, hydrological elements result in the variation of hydrological laws. The simulated hydrological state has no real correspondence in the basin, which leads to the simulation error. Because the state variables of the model cannot be measured directly in the prototype, the input and output process of the basin or reach can only be usually observed, and the state evolution of the system is inferred through comparative analysis. In addition, the observation error of input data of the model also makes the model parameter deviate from the true value during identification. The presence of these factors causes the flood forecast results to have estimation errors. Therefore, it is necessary to correct the model forecast results in real time to improve the forecast accuracy.

[0004] The real-time correction technology of the model forecast results can be roughly divided into three categories. A first technology is to explore the regularity of hydrological phenomena, study the simulation theory, improve the model structure and improve the model simulation ability. A second technology is to overcome the forecast error resulted from input errors and use the relatively effective parameter identification technology. A third technology is to perform a secondary forecast on the forecast residual and add the forecast residual to the forecast result. Therefore, the present disclosure provides a flood forecast error correction method based on a Volterra series filter. The present disclosure uses a Volterra series filter to establish a model for the flood forecast residual part and forecast the flood forecast residual part, so as to correct the model forecast results in real time. SUMMARY

[0005] In order to solve the problem that it is still difficult for the flood forecast accuracy to meet the actual requirements at present, the purpose of the present disclosure is to provide a flood forecast error correction method based on a Volterra series filter, so as to correct flood forecast results in real time and thereby improve the forecast accuracy.

[0006] In order to achieve the above purpose, the present disclosure provides the following technical scheme.

[0007] A flood forecast error correction method based on a Volterra series filter, including the following steps of:

[0008] acquiring storm flood data and using a Xinanjiang model for simulation to obtain a model residual sequence;

[0009] modeling the Xinanjiang model residual sequence based on the Volterra series filter;

[0010] solving the Volterra series filter based on a self-organizing method; and

[0011] correcting flood forecast results of the Xinanjiang model in real time based on the Volterra series filter.

[0012] As a further technical scheme of the present disclosure, the step of acquiring storm flood data and using a Xinanjiang model for simulation to obtain a model residual sequence includes:

[0013] acquiring storm flood data, where the storm flood data includes a rainfall value Pij(t) and a basin outlet flow (^(t);

[0014] inputting the rainfall value Pij(t) into the Xinanjiang model f() to obtain a forecast basin outlet flow Qt(t); and

[0015] obtaining the model residual sequence according to a difference between the forecast basin outlet flow and the basin outlet flow.

[0016] In some embodiments, the storm flood data is [P, Q] , P = [^(1),^(2),...,^(0 Pi,jW]T , Q = [^(1),^(2),...,^(0 Q^nyT , P^ denotes a rainfall value of an / -th flood at aj-th rainfall station at a time / , and Qi(t) denotes a basin outlet flow of the / -th flood at the time t, the forecast basin outlet flow is = f^Pij^y and the model residual sequence is x(0 = Qi(t) — Qt(t).

[0017] In some embodiments, the step of modeling the Xinanjiang model residual sequence based on the Volterra series filter includes:

[0018] reconstructing a phase space of a model residual sequence x(t) and forecasting xi+(m_1yr+1 = f(xi,Xi+T, , where r is a time delay, and m is an embedding dimension, Y = [y^Yz, - ,yN]T, Yi = >xi+(m_1)T]fi = 1,2,, ...,1V, and N denotes the length of the residual sequence;

[0019] selecting the time delay t;

[0020] selecting the embedding dimension m, and

[0021] using a second-order Volterra adaptive filter for fitting by f( ) according to = f^i> ^i + T> ^i+(7n-l)r+l a0 + Xi=o - i) + S™o1 a2x(t - l)x(t - D-

[0022] In some embodiments, the time delay t is obtained according to a mutual information method. The graph of the mutual information I(t) with respect to the delay time r is drawn in the actual calculation. When the mutual information decreases to a minimum value for the first time, the corresponding time Tmin is taken as the time delay t of reconstructing the phase space.

[0023] In some embodiments, the embedding dimension m is selected by using a trial calculation method, that is, a value of m gradually increases from 1 according to m = ml. a corresponding maximum Lyapunov exponent is calculated until the maximum Lyapunov exponent stops changing, to obtain a value of m as the embedding dimension, and the maximum Lyapunov exponent is obtained by using a Rosenstein method.

[0024] In some embodiments, the step of solving the Volterra series filter based on a selforganizing method includes:

[0025] selecting a combination of any two variables from m independent variables determined by the time delay r to construct an intermediate variable,

[0026] 2 2 Z, = Cn +C,X. + C^X. + CX + C.X- +C.X.X,, K U 1 / — J Z> I -r J J 1 1 '

[0027] selecting preparation data, and determining co, ci, ..., according to a least square mmY principle: 1 , where I denotes a number of data pairs, where the preparation data is a part of data selected from a complete set of phase points in a reconstructed phase space;

[0028] obtaining a least square error between the calculated value of zt and the actual check data according to the obtained coefficient of c, using the check data, ranking the errors in an increasing order, extracting a first m %as intermediate variables, and discarding other intermediate variables, where the check data is remaining phase point data;

[0029] deeming these intermediate variables as xh repeating a check until the least square error of the check data is no longer improved, and stopping calculation; 3

[0030] substituting coefficients back layer by layer, and obtaining a forecast expression.

[0031] In some embodiments, the step of correcting flood forecast results of the Xinanjiang model in real time based on the Volterra series filter includes:

[0032] inputting rainfall data into the Xinanjiang model in real time, and outputting a basin outlet flow process Qxi(ty

[0033] forecasting a residual sequence x^t);

[0034] correcting flood forecast results in real time through the basin outlet flow process Qxt(t) and the forecast residual sequence Xi(t), Qi(t) = + xt(ty

[0035] The present disclosure has the following beneficial effects. The flood forecast error correction method based on the Volterra series filter further supplements and improves the forecast results of the Xinanjiang model, improves the accuracy of flood forecast, and has strong operability and a wide application prospect. BRIEF DESCRIPTION OF THE DRAWINGS

[0036] FIG. lisa flowchart of a flood forecast error correction method based on a Volterra series filter.

[0037] FIG. 2 is a curve graph of a / (t) versus t in a flood forecast error correction method based on a Volterra series filter. DETAILED DESCRIPTION OF THE EMBODIMENTS

[0038] As shown in FIG. 1, a flood forecast error correction method based on a Volterra series filter includes the following steps.

[0039] In Step S100, storm flood data is acquired, and a Xinanjiang model is used for simulation to obtain a model residual sequence.

[0040] In Step S200, a Xinanjiang model residual sequence is modelled based on the Volterra series filter.

[0041] In Step S300, the Volterra series filter is solved based on a self-organizing method.

[0042] In Step S400, flood forecast results are corrected in real time.

[0043] The step of acquiring storm flood data and using a Xinanjiang model for simulation to obtain a model residual sequence includes: acquiring the storm flood data, where the storm flood data includes a rainfall value Pij(t) and a basin outlet flow Qyty inputting the rainfall value Ptj(t) into the Xinanjiang model f() to obtain a forecast basin outlet flow Qt(ty, and obtaining the model residual sequence according to a difference between the forecast basin outlet flow and the basin outlet flow.

[0044] In some embodiments of the present disclosure, the storm flood data is [P, Q], P = [^-(1),^-(2),...,P^.....P^n^ , Q = [¢ / 1),^(2),...,^(1)....., P^ denotes a rainfall value of an z-th flood at a j-th rainfall station at a time t, and ¢,(1) denotes a basin outlet flow ofthez-th flood at the time / , the forecast basin outlet flow is Qt(t) = ffPijCt)), and the model residual sequence is %(t) = ¢ / 1) — Qi(t).

[0045] The step of modeling the Xinanjiang model residual sequence based on the Volterra series filter includes:

[0046] reconstructing a phase space of the model residual sequence x(t), Y = [yi,y2, — >Yn]T^ Yi [-^i» ^t+T> >x i+(ni-i)r], i 1,2,,..., N. ^1+(777-1)7+1 - / (^1,3(+7,...,^1+(777-1)7), where N denotes a length of the residual sequence, r is a time delay, and m is an embedding dimension;

[0047] selecting the time delay t;

[0048] selecting the embedding dimension m;

[0049] using a second-order Volterra adaptive filter for fitting by f( ) according to xl+^m_^T+1 = f (xit xi+T,..., xi+(m_1)T), that is, xi+(m_1)r+1 = a0 + a^t - i') + 2 / / ' a2x(t - 0x(t -j).

[0050] The time delay t is obtained according to a mutual information method. The graph of the mutual information I(t) with respect to the delay time t is drawn in the actual calculation. When the mutual information decreases to a minimum value for the first time, the corresponding time Tmin is taken as the time delay r for reconstructing the phase space. Now, the specific steps are described as follows.

[0051] Assuming that the system is S, the time series consists of a series of values containing the information of the system: ” 2’ and the probability corresponding to each information P (5 ) 7 = 19 n value is A , 1.----, / 0 j|ie average mutual information obtained from each measured value Si is the entropy H of the system: i

[0052] Considering the dependence of the quantity sj+r, obtained after a time delay t from st on the quantity .v , q, = sj+T, there is a series of \x,q.} \s ,q and the formed system is referred to as a total coupling system [5(0]. When the measured value of s is known as sj, the uncertainty of q can be measured by conditional entropy H(Q | sj, which is expressed as:

[0053] ^21 ■’,) = -IX(?„■’, )l°g|A(¾)] =H(S,Q)-H(S) J

[0054] where , the measurement of the decrease of the uncertainty of due to S can be expressed as mutual information ),

[0055] 8) = H(Q) -H(Q\S) =H(Q) + H(S)-H(S, Q)

[0056] It is necessary to estimate the probability / f7 in the calculation of mutual information, in which a Box-counting method is usually used. If the point (s, q) on the plane of s, q is in a box N with the size of , P =---—As^q, where Nsq denotes a number of phase points in the Ntotal box, and Ntotal denotes a total number of phase points.

[0057] From the curve / (t) versus t in FIG. 2, when the time delay r = 4, the mutual information 1(4) = 0.9884, when r = 3,1(3) = 1.1044, and when t = 5,1(5) = 1.0703. Therefore, a first local minimum solution is r = 4, and the time delay is r= 4.

[0058] The embedding dimension m is selected by using a trial calculation method, that is, the value of m gradually increases from 1 according to m = m+1, a corresponding maximum Lyapunov exponent is calculated until the maximum Lyapunov exponent stops changing, at this time, the value of m is the embedding dimension, and the maximum Lyapunov exponent is obtained by using a Rosenstein method. The specific steps are described as follows.

[0059] (1) The time delay t is selected by using time series, and an m-dimensional phase space is constructed according to a total number n of observed data samples, where a number of the phase points isN: N= ) r;

[0060] (2) An initial phase point y0 is taken as a base point, and the point y} closest to y0 is selected as an endpoint from the remaining phase points in a point set {y,! to form an initial vector, in which the Euclidean distance between y0 and y. can be represented as £(^);

[0061] (3) The time step is taken as s, tx=t0+s, in which an initial vector evolves forward along a trajectory to obtain a new vector, and the Euclidean distance between the corresponding base point and the endpoint can be represented as A( / ); and in the corresponding time period, the exponential growth rate of the system linearity is represented as ; 1. Ai — — In 5

[0062] (4) The process is continued in this way until all phase points are traversed, and then the average value of various exponential growth rates is taken as a maximum Lyapunov exponent , , 1A1, £( / ) estimated value A, (N) = — >—In----—. n^s

[0063] (5) The embedding dimension m is increased gradually, and the processes (2), (3) and (4) are repeated until the Lyapunov exponent estimated value becomes relatively stable with the change of m, and the calculated result at this time is the maximum Lyapunov exponent value. When the embedding dimension m gradually increases from 3 to 7, the maximum Lyapunov exponent does not change, and the embedding dimension m is 7.

[0064] Table 1: Relational table between the embedding dimension m and the maximum Lyapunov exponent value embedding maximum Lyapunov embedding maximum Lyapunov dimension exponent value dimension exponent value 3 0.026 6 0.016 4 0.018 7 0.012 5 0.017 8 0.012

[0065] The step of solving the Vol terra series filter based on a sell -organizing method includes: selecting a combination of any two variables from m independent variables determined by the time delay t to construct an intermediate variable, zk =c0 + clxj +c2x. + c3x^ +c4Xj k = 1,2,---, / / / ( / / / -1) / 2 selecting preparation data, and determining co, ci, ..., a according to a least square principle, i.e., 1 , where I denotes a number of data pairs, where the preparation data is a part of data selected from a complete set of phase points in a reconstructed phase space; obtaining a least square error between the calculated value of zk and the actual check data according to the obtained coefficient of c, using the check data, ranking the errors in an increasing order, taking out a first m Zk^ intermediate variables, and discarding other intermediate variables, where the check data is the remaining phase point data; deeming these intermediate variables as x;, repeating above process until the least square error of the check data is no longer improved, and stopping calculation; substituting coefficients back layer by layer, and obtaining a forecast expression.

[0066] Specifically, through programming calculation, the least square error will not be improved after seven-layer iterative calculation. Starting from x(n-(m-1 )r), assuming that the residual values of respective time delay points [x(n-(m-l)z), x(n-(m-2)r), ..., x(n)] are [xo, xi..., xs], the expressions of respective layers are as follows:

[0067] A first layer: z^ = 0.125 + 0.134x4 - 0.337x6 - 4.6 x 10-5x42 -8.7 x 10-5x62 -1.1 x 10-5 x4x6 zW = 0.526 -0.084x2 -0.364x6 - 6.1 x l 0 \d - 8.1 +10 s x2 - 1.4 x l 0 4x2x6 z4 = 0.214-0.029^ -0.371x6 + 4.6 + 10 ’.x' A z 10 \v3 +1.710 ;x:x.

[0068] zF = 0.197 + 0.015.+ -0.371+,.+^ zj1’ = 0.522 -0.1 08x5 -0.395x6 -5.3 x 10-5x52 -7.3 x 10 5xg + 7.9 x 10-5x5x6 z2 ’ = 0.264 + O.O59xo - 0.382x6 + 7 x 10-6 x02 -1.1 x 10-4 xf +1.2 x 10-4x0x6 z^ = -0.522-0.016xj + 0.202x4+9xl0-6x2+7.9 + 10-5x4+3.04xl0-4x1x4

[0069] A second layer: z^ = -0.116 + 0.18z^ + 0.640zW + 5.4x 10^41^ -8.1 x 10 4z, f +1.4xlO^z^z^ 42) = -0.142 + 0.928z; + 0.692Z;1’ + 7.8x 10^1)2 -6.7x 10^2 -1.4xlO^'M” z^ =-0.426 + 0.141z2U0.88z2M.2xl0-3z(22+6.4 + 10^^-2.0 + 10-^2¼0

[0070] z;: 0.46 0.036z;: • 1.04z;:) I 3.0^10 '!z:;: + 1.9 + 10-¼°5 4+ 10¼1¼0 z4 = -0.133 + 0.9294° + 0.6554} +7.3x10 4Z|I)2 -4.96xlQ-4z2)2 -1.77xlQ-J41¼° 42) = -0.184 + 0.979z2) +0.036z2} -3.5x10-3z^ -3.2+10^2^+7.6 + 10¼1¼ z^ = -0.335 + O.873z2’ + 0.157z°} -8.1 x 10-5z2)2 + 2.13 x 10-¼°2 -1.33 x 10-3z^

[0071] A third layer: z2° = 0.158 + 0.742¼ +0.317zI2) -1.4 + 10^^-2.9 + 10¼°2 + 2.9 + 10¼2¼° 4° 0.453 I 0.207z;3 • 0.84z? 5.040¼3.5 + 10 + 7.9 + 10¼2¼° z^ = 0.147 + 0.733z4 +O.325z2° +2.5^.10-4^2 -1.8 + 10¼°2 +1.3xW^z^z^

[0072] 4° = 0.267 + 0.7664° + 0.29¼° + 5.1 + 10^^-1.0 + 10¼°2 -1.0 + 10¼2¼° zi3) = 0.209 + 0.36z2° + 0.694° +3.0 + 10^^+1.88 + 10¼°2 -2.7 + 10¼¼° 4° = 0.268 + 0.7764° +Q.2Sz^ + 1.4 + 10¼°2 -1.5 + 10¼°2 + 7.6 + 10¼2¼° z^ =0.375+ 0.245z4 +0.803zF -2.8 + 10¼°2 -1.7 + 10¼°2 + 3.9xl0-3z442)

[0073] A fourth layer: z^ = -0.347+ 0.897z4 +0.0944° d,5d0443-2.3xl0~3^2 +3.4 + 10-3z^z® z^ = -0.337+ 0.985z4 +0.005z4 +8.75 + 10^2 -TAxlQ-4^2 +4.8xl0-4zfzf> z^ =-0.29 + 1.O5z2)-0.06z4 +1AxW-3z(^2 +1.6 + 10-3 z^-2.5xlQ-3 z^z^

[0074] z^ = -0.29 + 0.0824^ +0.916z4 -3.0x l0 3Z|3'2 -2.8z lO ’zf'2 +6.4xl0-343)z^ z4 = -0.306 +1.25 lz4 - 0.25z4 + 6.3 x 10-3z^ + 8A x 1Q-3zff -1A x 1Q-2 z^zf 44) = -0.272 +1 ,2z4 -0.197z^ +7.4x10-3z^ +1.2 x 10-2-1.9 x 10-2z^z^ z4 = 0.415 I 0.29z : + 0.68z4 -8.1 x 10:z3' -1.7 x 10^2 + 2.6x W2^z^

[0075] A fifth layer: z^ = 0.168 + 0.948z4 + 0.056z4 +0.0227^2 +0.0134°2 -0.036z44} 4’ =0.227 + 0.31474 +0.708^4 +1.6xl0~34)2 +4.0xl0~34°2 -6.2x10-3z44} 4’ = 0.272- 0.414z4 +1 -44z4 + 1.1+102+1.4 x IO"2 z4 ~ 2.6 x l 0 2z^zf

[0076] z4 = 0.129- 0.772z4 +1.784^4 -2.2x 10 2z4? -3.5 xlQ-2z'f + 5.7 x 10-2z^z^ 4’ = 0.167 - 0.01 9z4 +1 -03z4 + 0.014^ -1.1 x 10 3 z4" - 9.2 X10-3 z^ z4 = 0.013-0.374 + 1.38z4 -0.024)2 + 0.04z44’ z4 = 0.24 - 0.67z4 + 1.69z4 -0.02z42 - 0.03z42 +0.04z4z4

[0077] A sixth layer: 4’ = 0.166 + 1.453z4 -0.4z4 +0.05z4'2 +0.04z42 -0.09z44^ z4 = -0.01 + 0.549z4 + 0.467z4 +3^10^2 +2.8x104z^ -6.7xl0-3z44} z4 = -0.133 + 0.614z4 +0.395z4 + 0.012z42 -0.016z4 ' -0.03z445)

[0078] z4 = -0.018 + 0.409z4 +0.648z4 +0.011z42 + 0.01z4^ -0.022z44’ z4 =-0.139 + 0.268z4 +0.727z4 -1.5x10-3z42 -2^10-4^+4^10^44 z4 = -0.037 + 0.485z4 +0.554z4 + 0.Qnz^ + 0.01^ -0.022z44) z4 =-0.148 + 0.745z4 + 0.253z4 -4x10-3z42 -3.6x10-3 z4" +8.2x10-34’4)

[0079] A seventh layer:

[0080] z4 =0.168 + 1.553z4 -0.5z4 + 0.07z4 +0.05z4)2 -0.08z44)

[0081] The complete expression (the forecast formula) is y=zo(7).

[0082] The step of correcting flood forecast results in real time includes:

[0083] inputting rainfall data into the Xinanjiang model in real time, and outputting the basin outl et fl ow process Qxi {t) \

[0084] forecasting the residual sequence x^t)-,

[0085] correcting flood forecast results in real time, Q((0 = +

[0086] The flood forecast results are corrected in real time. According to the comparison of forecast results (see Table 2), the peak flood flow after real-time correction is further closer to the measured value.

[0087] Table 2: Real-time correction results of flood forecast Flood frequency measured value No real-time correction of flood forecast Real-time correction of flood forecast peak flood flow (m3 / s) Peak moment peak flood flow (m3 / s)) Peak moment peak flood flow (m3 / s) Peak moment 20060519 617 2006-5-18 713 2006-5-18 664 2006-5- 20 o’clock 21 o’clock 18 21 o’clock 20070614 499 2007-6-14 16 o’clock 620 2007-6-14 16 o’clock 591 2007-6- 14 16 o’clock 20071009 1218 2007-10-8 14 o’clock 1334 2007-10-8 16 o’clock 1286 2007-10- 8 16 o’clock 20080610 641 2008-6-9 20 o’clock 697 2008-6-9 21 o’clock 668 2008-6-9 20 o’clock

Claims

1. A flood forecast error correction method based on a Volterra series filter, comprising the steps of:acquiring storm flood data and using a Xinanjiang model for simulation to obtain a model residual sequence;modeling a Xinanjiang model residual sequence based on the Volterra series filter;solving the Volterra series filter based on a self-organizing method; andcorrecting flood forecast results of the Xinanjiang model in real time based on the Volterra series filter.

2. The method according to claim 1, wherein the step of acquiring storm flood data and using a Xinanjiang model for simulation to obtain a model residual sequence comprises:acquiring storm flood data, wherein the storm flood data comprises a rainfall value and a basin outlet flow;inputting the rainfall value into the Xinanjiang model f( ) to obtain a forecast basin outlet flow; andobtaining the model residual sequence according to a difference between the forecast basin outlet flow and the basin outlet flow.

3. The method according to claim 1, wherein the step of modeling a Xinanjiang model residual sequence based on the Vol terra series filter comprises:reconstructing a phase space of a model residual sequence x(Z) and forecasting = f(Xi,Xi+T, where r is a time delay, and m is an embeddingdimension;selecting the time delay r;selecting the embedding dimension m; andsimulating and forecasting the model residual sequence using a second-order Volterra adaptive filter.

4. The method according to claim 3, wherein the time delay ris obtained according to a mutual information method.

5. The method according to claim 3, wherein the embedding dimension m is selected by using a trial calculation method, that is, a value of m gradually increases from 1 according to m = m+1, a corresponding maximum Lyapunov exponent is calculated until the maximum Lyapunov exponent stops changing, to obtain a value of m as the embedding dimension, and the maximum Lyapunov exponent is obtained by using a Rosenstein method.

6. The method according to claim 3, wherein the step of correcting flood forecast results of the Xinanjiang model in real time based on the Volterra series filter comprises:inputting rainfall data into the Xinanjiang model in real time, and outputting a basin outlet flow process Qxi(ty,forecasting a residual sequence x, (t);correcting flood forecast results in real time through the basin outlet flow process Qxi(t) and the forecast residual sequence %j(t).

7. The method according to claim 2, wherein the storm flood data is [P, Q], P = [^-(1),^(2),...,^(0.....PijWY , ¢ = ^(1),^(2),...,^(0.....¢,000 , P^t) denotes a rainfall value of an z-th flood at a / -th rainfall station at a time t, and ¢00 denotes a basin outlet flow of the z-th flood at the time t, the forecast basin outlet flow is ¢00 = f(.Pi,j(f)\ and the model residual sequence is x(0 = ¢00 — ¢00-8. The method according to claim 3, wherein the step of solving the Volterra series filter based on a self-organizing method comprises:selecting a combination of any two variables from m independent variables determined by the time delay rto construct an intermediate variable,zk = co +c\xi ^ckxi + c3xi + +c5xtxp k = 1,2,---, / / 7( / 7 / -1) / 2selecting preparation data, and determining co, ci, ..., C5 according to a least square principle: mmS(T / ~zki)21 , where I denotes a number of data pairs, wherein the preparation data is a partof data selected from a complete set of phase points in a reconstructed phase space;obtaining a least square error between the calculated value of Zk and the actual check data according to the obtained coefficient of c, using the check data, ranking the errors in an increasing order, extracting a first m z^as intermediate variables, and discarding other intermediate variables, wherein the check data is remaining phase point data;deeming these intermediate variables as Xi, repeating a check until the least square error of the check data is no longer improved, and stopping calculation;substituting coefficients back layer by layer, and obtaining a forecast expression.

Images