CN105115573A - Correction method and device for flood flow forecasting - Google Patents

Abstract

The invention discloses a correction method and device for flood flow forecasting, which are used for weakening the influences on a linear Kalman filtering real-time correction method by observation errors and guaranteeing the correction precision of the flood flow forecasting. The correction method for the flood flow forecasting comprises the steps: establishing a first-order autoregression model of remnant information between observed flow and forecasted flow; establishing a state equation according to parameters of the first-order autoregression model; determining a state parameter estimated value of the current moment from the first-order autoregression model and the state equation by adopting a robust M-estimator; obtaining a remnant forecasted value of next moment according to the determined state parameter estimated value of the current moment and the remnant information of the current moment; correcting forecasted flow of the next moment according to the remnant forecasted value of the next moment to obtain corrected flow of the next moment.

Description

A kind of bearing calibration of flood discharge forecast and device

Technical field

The present invention relates to flood forecasting technical field, be specifically related to a kind of bearing calibration and device of flood discharge forecast.

Background technology

Residual information between the forecasting runoff that Real-time Flood Forecasting system utilizes the observed volume of continuous renewal and Flood Forecasting Model to obtain, carries out real time correction, improves forecast precision.A kind of linear Kalman filter (English name Kalmanfilter adopting recursive algorithm is there is in prior art, be called for short KF) method, because it is without the need to storing all conception of history measurement informations, only need store previous moment state estimation and the new value of observation, save and calculate the used time, be widely used in Real-time Flood Forecasting correction.

At present, also widespread use is there is in telemetry system in Real-time Flood Forecasting, can realize automatically testing hydrographic data by telemetry system, ensure the ageing of observation data simultaneously, but use telemetry system also to bring telemetry system observational error to hydrographic data, telemetry system observational error can produce a very large impact the monitoring of waterlevel data.The reservoir inflow of Real-time Flood Forecasting system, mostly obtains waterlevel data by on-line monitoring, then obtains observed volume data according to water level-flow corresponding relation or reservoir area water balance calculation method.Overflow stage, adjoint strong wind also can make waterlevel data there is exceptional value, and when abnormal waterlevel data is transformed into observed volume data, the normal distribution of observational error is often contaminated.Desirable assumed condition due to linear Kalman filter method requires that observational error meets normal distribution, and the freshwater monitoring situation of reality can not meet this ideal hypothesis, according to this kalman filter method, the result estimated is influenced comparatively large, even occurs the phenomenon correcting poorer and poorer.

Summary of the invention

The object of the present invention is to provide a kind of bearing calibration and device of flood discharge forecast, for weakening the impact of observational error on linear Kalman filter real-time correction method, ensureing the correction accuracy of flood discharge forecast.

In order to achieve the above object, the present invention adopts so following technical scheme:

On the one hand, the invention provides a kind of bearing calibration of flood discharge forecast, comprising:

Set up first order autoregressive model to the residual information between observed volume and forecasting runoff, as measurement equation, the parameter according to described first order autoregressive model sets up state equation;

Adopt robust M-estimator from the state parameter valuation of described measurement equation and described state equation determination current time;

The residual error predicted value of subsequent time is obtained according to the state parameter valuation of the described current time determined and the residual information of current time;

The forecasting runoff of residual error predicted value to subsequent time according to described subsequent time corrects, and obtains the calibrated flow of subsequent time.

Wherein, described first order autoregressive model is set up to the residual information between observed volume and forecasting runoff, comprising:

Set up the measurement equation of described first order autoregressive model in the following way:

e _t＝φ _te _t-1+ε _t；

Wherein, described e _tfor the residual information of t, described e _t-1for the residual information in t-1 moment, described φ _tfor the state parameter of t, described ε _tfor the observation noise of t, described t is natural number;

Described e _tcalculate in the following way:

e _t＝Q _ot-Q _ct；

Wherein, described Q _otfor the observed volume of t, described Q _ctfor the forecasting runoff of t.

Wherein, the described parameter according to described first order autoregressive model sets up state equation, comprising:

Set up the state equation of described dynamic model in the following way:

φ _t＝Iφ _t-1+ω _t；

Wherein, described φ _tfor the state parameter of t, described I is the transition matrix of state equation and described I representation unit matrix, described φ _t-1for the state parameter in t-1 moment, described ω _tfor the process noise of t, described t is natural number.

Wherein, described employing robust M-estimator, from the state parameter valuation of described measurement equation and state equation determination current time, comprising:

Robust Modified Equivalent Weight Function is calculated by the robust characteristic variable designed in advance;

Iterated adjustment calculating is carried out according to described robust Modified Equivalent Weight Function parity price weight factor, when the equivalence weight factor does not meet preset power threshold condition in an iterative process, continue parity price weight factor and carry out Iterated adjustment calculating, only have when weight factor of equal value in iterative process meets described power threshold condition, export the robust gain matrix at the end of iteration;

The state parameter valuation of current time is calculated according to the robust gain matrix at the end of iteration.

Wherein, describedly carry out Iterated adjustment calculating according to described robust Modified Equivalent Weight Function parity price weight factor, comprising:

In the iterative process of kth step, according to the equivalence weight factor of (k-1) step of described robust Modified Equivalent Weight Function and t calculate the equivalence weight factor of the kth step of t described t and k is natural number;

According to the equivalence weight factor that the kth of described t walks the residual information in t-1 moment and the robust error covariance in t-1 moment calculate the robust gain matrix of t;

According to the equivalence weight factor that the kth of described t walks the robust gain matrix of the residual information in t-1 moment and the robust error covariance in t-1 moment, described t calculates the robust error covariance of t;

The state parameter valuation of the kth step of t is calculated according to the residual information of the robust gain matrix of t, t, the residual information in t-1 moment and the state parameter valuation in t-1 moment.

Wherein, the described state parameter valuation calculating current time according to the robust gain matrix at the end of iteration, comprising:

In described iterative process, weight factor of equal value meets described power threshold condition, comprising: described μ is described power threshold condition, and the iteration of kth step is the final step that iteration terminates, and current time is t, the robust gain matrix that at the end of obtaining iteration, the equivalence weight factor pair of the kth step of t is answered;

The state parameter valuation of robust gain matrix, the residual information of current time, the residual information of the previous moment of current time and the previous moment of current time that the equivalence weight factor pair walked according to the kth of t at the end of iteration is answered calculates the state parameter valuation of current time.

Wherein, described method also comprises:

The criterion of described robust M-estimator is:

$({\widehat{\mathrm{\φ}}}_t{e}_{t-1}-e_t){\stackrel{\‾}{p}}_t^{\left(1\right)}({\widehat{\mathrm{\φ}}}_t{e}_{t-1}-e_t)+({\widehat{\mathrm{\φ}}}_t-{\widehat{\mathrm{\φ}}}_{t-1}){\left(P_{t-1}\right)}^{-1}({\widehat{\mathrm{\φ}}}_t-{\widehat{\mathrm{\φ}}}_{t-1})=min;$

Wherein, described in for the initial equivalence weight factor of t, described in as the weight of t observational error, described in for the valuation of t state parameter, described in for the valuation of t-1 moment state parameter, described e _tfor the residual information of t, described e _t-1for the residual information in t-1 moment, described P _t-1for the covariance of the state parameter sequence in t-1 moment, as the weight of t-1 etching process error.

Wherein, the state parameter valuation of the described current time that described basis is determined and the residual information of current time obtain the residual error predicted value of subsequent time, comprising:

Calculate the residual error predicted value of the subsequent time of current time in the following way;

${\stackrel{\‾}{e}}_{t+1}={\widehat{\mathrm{\φ}}}_t\·e_t;$

Wherein, described in for residual error predicted value when subsequent time is t+1, described in for the state parameter valuation of t, described e _tfor residual information when current time is t.

Wherein, the described forecasting runoff of residual error predicted value to subsequent time according to described subsequent time corrects, and obtains the calibrated flow of subsequent time, comprising:

Calculate the calibrated flow of subsequent time in the following way:

$Q_{t+1}=Q_{c(t+1)}+{\stackrel{\‾}{e}}_{t+1};$

Wherein, described Q _t+1for calibrated flow when subsequent time is t+1, described Q _{c (t+1)}for forecasting runoff when subsequent time is t+1, described in for residual error predicted value when subsequent time is t+1.

On the other hand, the invention provides a kind of means for correcting of flood discharge forecast, comprising:

MBM, for setting up first order autoregressive model to the residual information between observed volume and forecasting runoff, as measurement equation, the parameter according to described first order autoregressive model sets up state equation;

State estimation module, for adopting robust M-estimator from the state parameter valuation of described measurement equation and described state equation determination current time;

Residual error forecast module, for obtaining the residual error predicted value of subsequent time according to the state parameter valuation of the described current time determined and the residual information of current time;

Flux modification module, corrects for the forecasting runoff of residual error predicted value to subsequent time according to described subsequent time, obtains the calibrated flow of subsequent time.

Wherein, described MBM, specifically for setting up the measurement equation of described first order autoregressive model in the following way:

e _t＝φ _te _t-1+ε _t；

Wherein, described e _tfor the residual information of t, described e _t-1for the residual information in t-1 moment, described φ _tfor the state parameter of t, described ε _tfor the observation noise of t, described t is natural number;

Described e _tcalculate in the following way:

e _t＝Q _ot-Q _ct；

Wherein, described Q _otfor the observed volume of t, described Q _ctfor the forecasting runoff of t.

Wherein, described MBM, specifically for setting up described state equation in the following way:

φ _t＝Iφ _t-1+ω _t；

Wherein, described φ _tfor the state parameter of t, described I is the transition matrix of state equation and described I representation unit matrix, described φ _t-1for the state parameter in t-1 moment, described ω _tfor the process noise of t, described t is natural number.

Wherein, described state estimation module, comprising:

First computing module, calculates robust Modified Equivalent Weight Function for the robust characteristic variable by designing in advance;

Second computing module, for carrying out Iterated adjustment calculating according to described robust Modified Equivalent Weight Function parity price weight factor, when the equivalence weight factor does not meet preset power threshold condition in an iterative process, continue parity price weight factor and carry out Iterated adjustment calculating, only have when weight factor of equal value in iterative process meets described power threshold condition, export the robust gain matrix at the end of iteration;

3rd computing module, for calculating the state parameter valuation of current time according to the robust gain matrix at the end of iteration.

Wherein, described second computing module, comprising:

First computing unit, in the iterative process that walks in kth, according to the equivalence weight factor of (k-1) step of described robust Modified Equivalent Weight Function and t calculate the equivalence weight factor of the kth step of t described t and k is natural number;

Second computing unit, for the equivalence weight factor walked according to the kth of described t the residual information in t-1 moment and the robust error covariance in t-1 moment calculate the robust gain matrix of t;

3rd computing unit, for the equivalence weight factor walked according to the kth of described t the robust gain matrix of the residual information in t-1 moment and the robust error covariance in t-1 moment, described t calculates the robust error covariance of t;

4th computing unit, calculates the state parameter valuation of the kth step of t for the residual information of the robust gain matrix according to t, t, the residual information in t-1 moment and the state parameter valuation in t-1 moment.

Wherein, in described iterative process, weight factor of equal value meets described power threshold condition, comprising: described μ is described power threshold condition, and the iteration of kth step is the final step that iteration terminates, and current time is t, the robust gain matrix that at the end of obtaining iteration, the equivalence weight factor pair of the kth step of t is answered; Described 3rd computing module, the state parameter valuation of the previous moment of the robust gain matrix that the equivalence weight factor pair specifically for walking according to the kth of t at the end of iteration is answered, the residual information of previous moment of current time, the residual information of current time and current time calculates the state parameter valuation of current time.

Wherein, the criterion of described robust M-estimator is:

$({\widehat{\mathrm{\φ}}}_t{e}_{t-1}-e_t){\stackrel{\‾}{p}}_t^{\left(1\right)}({\widehat{\mathrm{\φ}}}_t{e}_{t-1}-e_t)+({\widehat{\mathrm{\φ}}}_t-{\widehat{\mathrm{\φ}}}_{t-1}){\left(P_{t-1}\right)}^{-1}({\widehat{\mathrm{\φ}}}_t-{\widehat{\mathrm{\φ}}}_{t-1})=min;$

Wherein, described in for the initial equivalence weight factor of t, described in as the weight of t observational error, described in for the valuation of t state parameter, described in for the valuation of t-1 moment state parameter, described e _tfor the residual information of t, described e _t-1for the residual information in t-1 moment, described P _t-1for the covariance of the state parameter sequence in t-1 moment, (P _t-1) ^-1as the weight of t-1 etching process error.

Wherein, described residual error forecast module, specifically for calculating the residual error predicted value of the subsequent time of current time in the following way;

${\stackrel{\‾}{e}}_{t+1}={\widehat{\mathrm{\φ}}}_t\·e_t;$

Wherein, described in for residual error predicted value when subsequent time is t+1, described in for the state parameter valuation of t, described e _tfor residual information when current time is t.

Wherein, described flux modification module, specifically for calculating the calibrated flow of subsequent time in the following way:

$Q_{t+1}=Q_{c(t+1)}+{\stackrel{\‾}{e}}_{t+1};$

Wherein, described Q _t+1for calibrated flow when subsequent time is t+1, described Q _{c (t+1)}for forecasting runoff when subsequent time is t+1, described in for residual error predicted value when subsequent time is t+1.

Adopt after technique scheme, technical scheme provided by the invention will by following advantage:

First first order autoregressive model is set up to the residual information between observed volume and forecasting runoff, as measurement equation, and set up state equation according to the parameter of first order autoregressive model, then adopt robust M-estimator from the state parameter valuation of measurement equation and state equation determination current time, the state parameter valuation of the current time that following basis is determined and the residual information of current time obtain the residual error predicted value of subsequent time, and correct according to the forecasting runoff of residual error predicted value to subsequent time of subsequent time, obtain the calibrated flow of subsequent time.Robust M-estimator is introduced in linear Kalman filter real-time correction method by the present invention, utilizes robust M-estimator to reduce the state parameter valuation of unknown quantity as far as possible.Abnormal observational error is reduced as much as possible, weakens observational error to the impact of linear Kalman filter real-time correction method, ensure the correction accuracy of flood discharge forecast.

Term " first ", " second " etc. in instructions of the present invention and claims and above-mentioned accompanying drawing are for distinguishing similar object, and need not be used for describing specific order or precedence.Should be appreciated that the term used like this can exchange in the appropriate case, this is only describe in embodiments of the invention the differentiation mode that the object of same alike result adopts when describing.In addition, term " comprises " and " having " and their any distortion, intention is to cover not exclusive comprising, to comprise the process of a series of unit, method, system, product or equipment being not necessarily limited to those unit, but can comprise clearly do not list or for intrinsic other unit of these processes, method, product or equipment.

Below be described in detail respectively.

An embodiment of the bearing calibration of flood discharge forecast of the present invention, can be applicable to, in the monitoring correction to flood, refer to shown in Fig. 1, the bearing calibration of flood discharge forecast provided by the invention, can comprise the steps:

101, first order autoregressive model is set up to the residual information between observed volume and forecasting runoff, as measurement equation, and set up state equation according to the parameter of first order autoregressive model.

In embodiments of the present invention, the observed volume of flood refers to the flood discharge measured in real time by telemetry system, and telemetry system can get observed volume based on the method for telemetering of prior art to the measurement of flood.The forecasting runoff of flood refers to the traffic forecast value of the future time instance using flood real-time predicting method to obtain.

In the present invention from initial time to current time, set up first order autoregressive model, as measurement equation according to the residual information between observed volume in the same time and forecasting runoff.Can calculate the residual information in this moment according to observed volume and forecasting runoff for synchronization, wherein residual information refers to the difference between observed volume value and forecasting runoff value, is namely the difference of actual observation value and forecasting runoff.

There is observational error in the observed volume got due to telemetry system in the present invention, this observational error of telemetry system does not meet normal distribution, concrete, observational error is made up of two class errors, one class is the random observational error meeting normal distribution, another kind of is anomalous differences, and cause primarily of reasons such as dropout, signal disturbing and mechanical faults, this kind of error disobeys normal distribution.The superposition of two class errors makes the observational error of flood flow observation flow no longer meet normal distribution, but main part is normal distribution, and secondary part is the contaminated normal distribution of anomalous differences distribution.Conventional linear Kalman filter method cannot be re-used in this case and carry out real time correction, desirable assumed condition due to existing kalman filter method requires that observational error meets normal distribution, and the freshwater monitoring situation of reality can not meet this ideal hypothesis, according to existing linear Kalman filter method, the result estimated is influenced comparatively large, even occurs the phenomenon correcting poorer and poorer.The embodiment of the present invention proposes robust Kalman filtering method for the real time correction to flood forecasting flow for this reason.First autoregression (Auto-regressive is set up according to residual information, AR) model and state equation, concrete, first order autoregressive model is set up to the residual information between observed volume and forecasting runoff, and sets up state equation according to the parameter of first order autoregressive model.Set up autoregressive model and state equation in the present invention, employing be robust Kalman filtering algorithm, be a recursive feedback algorithm in essence.Its point two parts: time update equation (autoregressive model of residual error) and measuring state renewal equation (state equation).Wherein, the recursion of time update equation deadline, the feedback of measuring state renewal equation completion status value.

Concrete, in some embodiments of the invention, in step 101, first order autoregressive model is set up to the residual information between observed volume and forecasting runoff, can comprise the steps:

A1, set up the measurement equation of first order autoregressive model in the following way:

e _t＝φ _te _t-1+ε _t；

Wherein, e _tfor the residual information of t, e _t-1for the residual information in t-1 moment, φ _tfor the state parameter of t, ε _tfor the observation noise of t, t is natural number;

E _tcalculate in the following way:

e _t＝Q _ot-Q _ct；

Wherein, Q _otfor the observed volume of t, Q _ctfor the forecasting runoff of t.

Realize in scene at above-mentioned A1, describe the measurement equation in autoregressive model for moment t, by the measurement equation of above-mentioned record, the concrete value of t is by needing to determine the flood calibrated flow in which moment and selecting the intermediate variable in which moment to decide.In addition according to above-mentioned e _tcomputing formula, the residual information of t can be calculated according to the residual information in t-1 moment.

It should be noted that, in above-mentioned measurement equation, have anomalous differences due in the observed volume that telemetry system gets, i.e. ε _tno longer Normal Distribution, but obedience main part is normal distribution, secondary part is the contaminated normal distribution of anomalous differences distribution.

In some embodiments of the invention, the state parameter according to first order autoregressive model in step 101 sets up state equation, specifically can comprise the steps:

A2, set up state equation in the following way:

φ _t＝Iφ _t-1+ω _t；

Wherein, φ _tfor the state parameter of t, I is the transition matrix of state equation and I representation unit matrix, φ _t-1for the state parameter in t-1 moment, ω _tfor the process noise of t, t is natural number.

Realize in scene at above-mentioned A2, describe the state equation in dynamic model for moment t, by the state equation of above-mentioned record, the concrete value of t is by needing to determine the flood calibrated flow in which moment and selecting the intermediate variable in which moment to decide.In addition according to above-mentioned φ _tcomputing formula, the state parameter of t can be calculated according to the state parameter in t-1 moment.

It should be noted that, in above-mentioned state equation, ω _tstill meet normal distribution, therefore to ω _testimation can adopt least square principle.

102, adopt robust M-estimator from the state parameter valuation of measurement equation and state equation determination current time.

In embodiments of the present invention, after setting up first order autoregressive model and state equation, the state parameter valuation of current time is estimated by recursive algorithm, known by Such analysis, there is observational error in the observed volume that telemetry system gets, observational error directly affects the correction to flood forecasting flow, use robust M-estimator from the state parameter valuation of measurement equation and state equation determination current time in the present invention, introduced by robust estimation theory in linear kalman filter method in the present invention, robust M-estimator is that a kind of generalized likelihood is estimated.

In an embodiment of the present invention, the criterion of robust M-estimator has done appropriate adjustment on M estimation criterion basis.Wherein M estimation criterion is:

$({\widehat{\mathrm{\φ}}}_t{e}_{t-1}-e_t){\left(Z_{t-1}\right)}^{-1}({\widehat{\mathrm{\φ}}}_t{e}_{t-1}-e_t)+({\widehat{\mathrm{\φ}}}_t-{\widehat{\mathrm{\φ}}}_{t-1}){\left(P_{t-1}\right)}^{-1}({\widehat{\mathrm{\φ}}}_t-{\widehat{\mathrm{\φ}}}_{t-1})=min;$

Wherein, for the valuation of t state parameter, for the valuation of t-1 moment state parameter, e _tfor the residual information of t, e _t-1for the residual information in t-1 moment, Z _t-1for the covariance of t observation sequence, (Z _t-1) ^-1as the weight of t observational error, P _t-1for the covariance of the state parameter sequence in t-1 moment, (P _t-1) ^-1as the weight of t-1 etching process error.

When there is exceptional value in t observed volume, e _tthere is exceptional value, use (Z _t-1) ^-1no longer applicable as weight.Robust M-estimator utilizes the equivalence weight factor replace (Z _t-1) ^-1as the weight of observed volume, criterion becomes:

$({\widehat{\mathrm{\φ}}}_t{e}_{t-1}-e_t){\stackrel{\‾}{p}}_t^{\left(1\right)}({\widehat{\mathrm{\φ}}}_t{e}_{t-1}-e_t)+({\widehat{\mathrm{\φ}}}_t-{\widehat{\mathrm{\φ}}}_{t-1}){\left(P_{t-1}\right)}^{-1}({\widehat{\mathrm{\φ}}}_t-{\widehat{\mathrm{\φ}}}_{t-1})=min;$

Wherein, for the initial equivalence weight factor of t, as the weight of t observational error, for the valuation of t state parameter, for the valuation of t-1 moment state parameter, e _tfor the residual information of t, e _t-1for the residual information in t-1 moment, P _t-1for the covariance of the state parameter sequence in t-1 moment, (P _t-1) ^-1as the weight of t-1 etching process error.

The criterion of robust M-estimator is differentiated according to state parameter and equals 0 and can obtain:

$e_{t-1}{\stackrel{\‾}{p}}_t^{\left(1\right)}({\widehat{\mathrm{\φ}}}_t{e}_{t-1}-e_t)+{\left(P_{t-1}\right)}^{-1}({\widehat{\mathrm{\φ}}}_t-{\widehat{\mathrm{\φ}}}_{t-1})=0;$

Wherein, P _t-1for the robust error covariance in t-1 moment, e _t-1for the residual information in t-1 moment, for the initial equivalence weight factor of t, for the state parameter valuation of t, for the state parameter valuation in t-1 moment.

Further, above formula is out of shape, obtains:

According to recurrence relation, can derive and obtain following expression:

P _t＝(I-K _te _t-1)P _t-1；

Now, P in above-mentioned formula _tcontaining the equivalence weight factor in the state parameter error covariance represented, be called robust error covariance.

It should be noted that, robust Kalman filtering bearing calibration is adopted in the embodiment of the present invention, by robust M-estimator can realize parity price power three sections power take respectively protect power, fall power, zero power process, retain normal data, weaken the impact of suspicious data, eliminate bad data, thus make observational error limited to valuation Influence on test result, high degree overcomes existing Kalman filtering bearing calibration affects excessive problem when processing non-normal errors by exceptional value, ensure that Real-time Flood Forecasting precision.

In some embodiments of the invention, step 102 adopts robust M-estimator from the state parameter valuation of measurement equation and state equation determination current time, specifically can comprise the steps:

C1, the robust characteristic variable passing through to design in advance calculate robust Modified Equivalent Weight Function;

C2, carry out Iterated adjustment calculating according to robust Modified Equivalent Weight Function parity price weight factor, when the equivalence weight factor does not meet preset power threshold condition in an iterative process, continue parity price weight factor and carry out Iterated adjustment calculating, only have when weight factor of equal value in iterative process meets power threshold condition, export the robust gain matrix at the end of iteration;

C3, calculate the state parameter valuation of current time according to the robust gain matrix at the end of iteration.

Wherein, in step C1, the embodiment of the present invention adopts the parameter of robust M-estimator to the current time of system state to estimate, first need to design robust characteristic variable (representing with alphabetical M) in advance, then robust Modified Equivalent Weight Function (representing with letter w) is calculated by robust characteristic variable, the robust M-estimator adopted in the present invention needs to be applicable to flood forecasting, so just needs in conjunction with flood forecasting feature to design robust Modified Equivalent Weight Function.Robust Modified Equivalent Weight Function is the core point of robust estimation theory, and the robust Modified Equivalent Weight Function for different problem is different, and the design of robust Modified Equivalent Weight Function must meet the actual features of the required flood forecasting problem solved.In Real-time Flood Forecasting process, careful attitude to be kept to the rejecting of suspected outlier, namely need to fall power to suspicious data, normal data maintains former power, and bad data is eliminated, based on such actual features, mode at the IGG syllogic Modified Equivalent Weight Function of survey field Application comparison maturation can be selected to design robust Modified Equivalent Weight Function, such as can use IGGI or IGGIII syllogic Modified Equivalent Weight Function, specifically connected applications scene can need, do not limit herein.

After designing robust Modified Equivalent Weight Function by step C1, continue to perform step C2, carry out Iterated adjustment calculating according to robust Modified Equivalent Weight Function parity price weight factor, utilize robust equivalence weight characteristic variable whether to judge in fresh information containing anomalous differences.If containing anomalous differences, reduce fresh information weight factor according to robust Modified Equivalent Weight Function, weaken it to the impact correcting result; If not containing anomalous differences, weight factor is constant, is similar to and is equal to linear Kalman filter.Robust Kalman filtering real-time correction method utilizes robust Modified Equivalent Weight Function to improve fresh information weight factor, exceptional value observational error not being met to normal distribution reduces power even zero power, make it limited to valuation Influence on test result, high degree overcomes linear Kalman filter bearing calibration affects excessive problem when processing non-normal errors by exceptional value, ensure that Real-time Flood Forecasting precision.Whether it should be noted that, the Iterated adjustment of parity price power weight factor is the derivation of a poll, can pre-set power threshold condition, can the basis for estimation that calculates of finishing iteration as.In the iterative computation of each step, all need to calculate the robust gain matrix in this step iterative process, and when satisfied power threshold condition, export the robust gain matrix at the end of iteration.

Further, in some embodiments of the invention, step C2 carries out Iterated adjustment calculating according to robust Modified Equivalent Weight Function parity price weight factor, specifically can comprise the steps:

C21, in the iterative process of kth step, according to the equivalence weight factor of (k-1) step of robust Modified Equivalent Weight Function and t calculate the equivalence weight factor of the kth step of t t and k is natural number;

C22, the equivalence weight factor walked according to the kth of t the residual information in t-1 moment and the robust error covariance in t-1 moment calculate the robust gain matrix of t;

C23, the equivalence weight factor walked according to the kth of t the robust gain matrix of the residual information in t-1 moment and the robust error covariance in t-1 moment, t calculates the robust error covariance of t;

C24, residual information according to the robust gain matrix of t, t, the residual information in t-1 moment and the state parameter valuation in t-1 moment calculate the state parameter valuation of the kth step of t.

Wherein, what describe during the Iterated adjustment of step C21 to C24 calculates be the iterative process of any one step is example, wherein, k in the iterative process of kth step can refer to any one step, adopt iterative algorithm in robust Kalman filtering bearing calibration in the present invention, calculated the equivalence weight factor that can get (k-1) step of t by last round of Iterated adjustment combine the robust Modified Equivalent Weight Function precomputed again, the equivalence weight factor of the kth step of rear t can be calculated according to the equivalence weight factor that the kth of t walks the residual information in t-1 moment and the robust error covariance in t-1 moment calculate the robust gain matrix of t, wherein, the residual information in t-1 moment is calculated by the recursion formula of aforementioned residual information, the robust error covariance in t-1 moment can be obtained, according to the equivalence weight factor that the kth of t walks by the robust error covariance formula in the last moment before t the robust gain matrix of the residual information in t-1 moment and the robust error covariance in t-1 moment, t calculates the robust error covariance of t, the state parameter valuation of the kth step of t is calculated according to the residual information of the robust gain matrix of t, t, the residual information in t-1 moment and the state parameter valuation in t-1 moment, wherein, the state parameter valuation in t-1 moment can be obtained by the state parameter valuation computing formula in the last moment before t.Be understandable that, the choosing power that what abovementioned steps C21 to step C24 described is in an iterative process calculates, sequencing between each computation process of weight factor of equal value, robust gain matrix, robust error covariance, state parameter valuation in an iterative process can not limit, and also needs according to successively relevant with the acquisition of some parameter before current time.

In step C3, obtain robust gain matrix can calculate current time state parameter valuation according to the robust gain matrix at the end of iteration, it should be noted that, when adopting robust gain matrix to calculate the state parameter valuation of current time, do not limit and need to adopt other parameter except robust gain matrix to carry out computing mode parameter estimation, in actual applications, the state parameter valuation of the previous moment using current time is also needed in the computation process of recursion, the residual information etc. of the residual information of current time and the previous moment of current time, specifically can need to obtain concrete parameter value according to actual computation.

In some embodiments of the invention, step C3 calculates the state parameter valuation of current time according to the robust gain matrix at the end of iteration, specifically can comprise the steps:

In C31, iterative process, weight factor of equal value meets power threshold condition, comprising: μ is power threshold condition, and the iteration of kth step is the final step that iteration terminates, and current time is t, the robust gain matrix that at the end of obtaining iteration, the equivalence weight factor pair of the kth step of t is answered;

The state parameter valuation of robust gain matrix, the residual information of current time, the residual information of the previous moment of current time and the previous moment of current time that C32, the equivalence weight factor pair walked according to the kth of t at the end of iteration are answered calculates the state parameter valuation of current time.

Wherein, in step C31, weight factor of equal value meets power threshold condition, between the equivalence weight factor of kth-1 step of the equivalence weight Summing Factor t that the kth referring to t walks, difference is very little, be less than power threshold condition μ, in this case, can think that walking iteration in kth terminates, the robust gain matrix that the equivalence weight factor pair exporting now kth step is answered.If current time is t, the state parameter valuation of the robust gain matrix that the equivalence weight factor pair walked according to the kth of t at the end of iteration is answered, the residual information of t, the residual information in t-1 moment and t-1 calculates the state parameter valuation of t.

It should be noted that, previous embodiment realize scene under, for selected different robust Modified Equivalent Weight Function, such as Huber function, Hampel function, can also amplify out the step of different Iterated adjustment, the mode of concrete Iterated adjustment can be determined according to specifically choosing of robust Modified Equivalent Weight Function.

103, the residual error predicted value of subsequent time is obtained according to the state parameter valuation of the current time determined and the residual information of current time.

In embodiments of the present invention, after employing robust M-estimator obtains the state parameter valuation of current time, residual information in conjunction with current time can calculate the residual error predicted value of subsequent time, wherein, residual error predicted value is the predicted value of the residual information to certain moment following, describe the account form of residual information in foregoing teachings, utilize the state parameter valuation of the current time got in the residual information of current time and step 102 can estimate the residual error predicted value of the subsequent time of current time.

In some embodiments of the invention, step 102 obtains the residual error predicted value of subsequent time according to the residual information of the state parameter valuation of the current time determined and current time, specifically can comprise the steps:

Calculate the residual error predicted value of the subsequent time of current time in the following way;

${\stackrel{\‾}{e}}_{t+1}={\widehat{\mathrm{\φ}}}_t\·e_t;$

Wherein, for residual error predicted value when subsequent time is t+1, for the state parameter valuation of t, e _tfor residual information when current time is t.

104, correct according to the forecasting runoff of residual error predicted value to subsequent time of subsequent time, obtain the calibrated flow of subsequent time.

In embodiments of the present invention, the residual error predicted value of subsequent time is estimated by robust M-estimator, then the residual error predicted value of subsequent time may be used for correcting the forecasting runoff of the subsequent time that Flood Forecasting Model dopes, the forecasting runoff of subsequent time use the residual error predicted value of subsequent time correct after exports and shows the calibrated flow of subsequent time, what the calibrated flow of the subsequent time exported in the present invention represented is obtain correction result according to robust Kalman filtering bearing calibration of the present invention.Owing to robust estimation theory being introduced in linear Kalman filter real-time correction method in the present invention, utilize robust M-estimator to have when anomalous differences is inevitable, make full use of valid data, restriction utilizes data available, get rid of the characteristic of exceptional value, make unknown quantity valuation reduce or remit the characteristic of anomalous differences impact as far as possible.Distribute according to the pollution of error reality in the present invention, select the Modified Equivalent Weight Function being applicable to flood forecasting real-time correction method, in conjunction with Kalman filtering, derive robust Kalman filtering real-time correction method.Further, can also be calculated by the Iterated adjustment of robust Modified Equivalent Weight Function, the equivalence weight factor containing anomalous differences observed reading is progressively reduced, weaken its impact on bearing calibration, ensure the correction accuracy to flood forecasting flow.

In some embodiments of the invention, step 104 corrects according to the forecasting runoff of residual error predicted value to subsequent time of subsequent time, obtains the calibrated flow of subsequent time, specifically can comprise the steps:

Calculate the calibrated flow of subsequent time in the following way:

$Q_{t+1}=Q_{c(t+1)}+{\stackrel{\‾}{e}}_{t+1};$

Wherein, Q _t+1for calibrated flow when subsequent time is t+1, Q _{c (t+1)}for forecasting runoff when subsequent time is t+1, for residual error predicted value when subsequent time is t+1.

It should be noted that, use to Q _{c (t+1)}correction, be that the state parameter valuation obtained by robust M-estimator is calculated, this robust Kalman filtering account form of the present invention weakens the impact of observational error, can ensure correction accuracy.

Known by the above embodiment description of this invention, first first order autoregressive model is set up to the residual information between observed volume and forecasting runoff, as measurement equation, and set up state equation according to the parameter of first order autoregressive model, then adopt robust M-estimator from the state parameter valuation of measurement equation and state equation determination current time, the state parameter valuation of the current time that following basis is determined and the residual information of current time obtain the residual error predicted value of subsequent time, and correct according to the forecasting runoff of residual error predicted value to subsequent time of subsequent time, obtain the calibrated flow of subsequent time.Robust M-estimator is introduced in linear Kalman filter real-time correction method by the present invention, utilize robust M-estimator to reduce the weight of exceptional value as far as possible, the valuation impact of abnormal observational error on state parameter is reduced, weaken observational error to the impact of linear Kalman filter real-time correction method, ensure the correction accuracy of forecasting runoff.

For ease of better understanding and implement the such scheme of the embodiment of the present invention, corresponding application scenarios of illustrating below is specifically described.In prior art, Kalman filtering real-time correction method needs the restriction of Normal Distribution to error, the kalman filter method namely described in background technology, also can be referred to as linear Kalman filter, also extending outside linear kalman filter method has EKF and Ensemble Kalman Filter etc.This kind of Kalman filtering real time correction algorithm extended all is supposed under a kind of ideal conditions, namely requires observational error and state error Normal Distribution, if hypothesis is set up, then correction accuracy is higher, and supposing is false, and correction accuracy will be greatly affected even distortion.Refer to as shown in Fig. 2-a, the bearing calibration of the flood discharge forecast provided for the embodiment of the present invention realize principle schematic, for overcoming Kalman filtering real-time correction method in prior art error needed to the restriction of Normal Distribution, the present invention proposes a kind of robust Kalman filtering real-time correction method (RobustKalmanfilter is called for short RKF).Robust estimation theory is introduced in linear Kalman filter real-time correction method, utilize Robust filter to have when anomalous differences is inevitable, make full use of valid data, restriction utilizes data available, get rid of the method for exceptional value, make unknown quantity valuation reduce or remit the characteristic of anomalous differences impact as far as possible.Distribute according to the pollution of error reality, select the Modified Equivalent Weight Function being applicable to flood forecasting real-time correction method, in conjunction with Kalman filtering, derive robust Kalman filtering real-time correction method.Calculated by the Iterated adjustment of Modified Equivalent Weight Function, the equivalence weight containing anomalous differences observed reading is progressively reduced, weaken its impact on bearing calibration, ensure correction accuracy.

The bearing calibration of flood discharge forecast provided by the invention realize principle, specifically can comprise the steps:

First, what the present invention adopted is robust Kalman filtering real time correction, needs the iterative computation by the equivalence weight factor in the present invention, in conjunction with the recursive algorithm of Kalman filtering, obtains real time correction valuation.

A (), current time (t) (use Q from telemetry system Real-time Obtaining to observed volume _otrepresent), then calculate and obtain residual information e _t=Q _ot-Q _ct, Q _ctbe the forecasting runoff that t Flood Forecasting Model gets, and compose Q _ctthe initial equivalence weight factor i.e. e _tthe initial equivalence weight factor

(b), set up first order autoregressive model according to residual information, following measurement equation can be obtained:

e _t＝φ _te _t-1+ε _t(1)

In formula (1): φ _tfor the state parameter of t, ε _tfor the observation noise of t.

To the state parameter φ in first order autoregressive model _tset up dynamic model, following state equation can be obtained:

φ _t＝Iφ _t-1+ω _t(2)

In formula (2): ω _tbe the process noise of t, the transition matrix I of state equation is unit battle array.

If (c) ω _t, ε _tnormal Distribution, can meet standard Kalman filtering design conditions.But there is anomalous differences in remote measurement observed reading, i.e. ε _tno longer Normal Distribution, but obedience main part is normal distribution, secondary part is the contaminated normal distribution of anomalous differences distribution.If now ε _tcontinue the least square principle adopting standard Kalman filtering, exceptional value will produce extreme influence to valuation result.Robust estimation theory is introduced, to ε in the present invention _temploying generalized likelihood is estimated, i.e. robust M-estimator.And the ω in state equation _tstill meet normal distribution, therefore it can continue to adopt least square principle.Now, the extremum conditions of robust Kalman filtering becomes:

$({\widehat{\mathrm{\φ}}}_t{e}_{t-1}-e_t){\stackrel{\‾}{p}}_t^{\left(1\right)}({\widehat{\mathrm{\φ}}}_t{e}_{t-1}-e_t)+({\widehat{\mathrm{\φ}}}_t-{\widehat{\mathrm{\φ}}}_{t-1}){\left(P_{t-1}\right)}^{-1}({\widehat{\mathrm{\φ}}}_t-{\widehat{\mathrm{\φ}}}_{t-1})=min---\left(3\right)$

In formula (3) formula, for the state parameter valuation of t, for the state parameter valuation in t-1 moment, P _t-1for the robust error covariance in t-1 moment.

To state parameter, extreme value is asked for formula (3):

$e_{t-1}{{\stackrel{\‾}{p}}_t}^{\left(1\right)}({\widehat{\mathrm{\φ}}}_t{e}_{t-1}-e_t)+{\left(P_{t-1}\right)}^{-1}({\widehat{\mathrm{\φ}}}_t-{\widehat{\mathrm{\φ}}}_{t-1})=0---\left(4\right)$

D (), iterative computation start, the value of iterative steps k=1, k increases progressively along with the carrying out of iteration.

(e), the equivalence weight factor walked according to the kth of t the residual information in t-1 moment and the robust error covariance in t-1 moment calculate the robust gain matrix of t, and robust gain matrix can be following formula:

$K_t=\frac{e_{t-1}{\stackrel{\‾}{p}}_t^{\left(k\right)}{P}_{t-1}}{1+e_{t-1}{\stackrel{\‾}{p}}_t^{\left(k\right)}{P}_{t-1}{e}_{t-1}}---\left(5\right)$

In formula (5), e _t-1for the residual information in t-1 moment, P _t-1for the robust error covariance in t-1 moment.

(f), the equivalence weight factor walked according to the kth of t the robust gain matrix of the residual information in t-1 moment and the robust error covariance in t-1 moment, t calculates the robust error covariance of t, and robust error covariance can be following formula:

P _t＝(I-K _te _t-1)P _t-1(6)

In formula (6), e _t-1for the residual information in t-1 moment, I is unit battle array, P _t-1for the robust error covariance in t-1 moment, K _tfor the robust gain matrix of t.

G (), residual information according to the robust gain matrix of t, t, the residual information in t-1 moment and the state parameter valuation in t-1 moment calculate the state parameter valuation of the kth step of t, state parameter valuation can be following formula:

${\widehat{\mathrm{\φ}}}_t^{\left(k\right)}={\widehat{\mathrm{\φ}}}_{t-1}+K_t(e_t-e_{t-1}{\widehat{\mathrm{\φ}}}_{t-1})---\left(7\right)$

In formula (7), for the state parameter valuation that the kth of t walks, e _tfor the residual information of t, e _t-1for the residual information in t-1 moment, K _tfor the robust gain matrix of t, for the state parameter valuation in t-1 moment.

(h), calculate robust characteristic variable M, then the equivalence weight factor of (k-1) step according to robust Modified Equivalent Weight Function W and t calculate the equivalence weight factor of kth step, the calculating of M and W is described in detail below.

${{\stackrel{\‾}{p}}_t}^{\left(k\right)}={{\stackrel{\‾}{p}}_t}^{(k-1)}\·W---\left(8\right)$

In above-mentioned formula, by the iteration of k=k+1, the equivalence weight factor of each step can be calculated.

If be the final step that iteration terminates in the iteration of kth step, then the state parameter valuation of the kth step of t be the state parameter valuation of t therefore there is following relation:

(i), repetition above-mentioned steps (d), (e), (f), g () carries out Iterated adjustment, judge whether the equivalence weight factor meets power threshold condition, concrete, in order to ensure computational accuracy, power threshold condition μ can get 0.0001, if meet iteration terminates.μ can also get other numerical value, herein just citing.

(j), calculate the residual error predicted value of the subsequent time of current time in the following way:

${\stackrel{\‾}{e}}_{t+1}={\widehat{\mathrm{\φ}}}_t{e}_t---\left(9\right)$

The calibrated flow of (k), in the following way calculating subsequent time t+1:

$Q_{t+1}=Q_{C(t+1)}+{\stackrel{\‾}{e}}_{t+1}---\left(10\right)$

According to final weight factor judge whether the observed volume of this period exists exceptional value.If for valid data; for useful data; for exceptional value.

Following design is applicable to the robust Modified Equivalent Weight Function W of flood forecasting feature.Robust Modified Equivalent Weight Function is the core of robust theory, and the Modified Equivalent Weight Function of different problem is different, the actual features of the necessary compliance problem of design of Modified Equivalent Weight Function.In Real-time Flood Forecasting process, careful attitude to be kept to the rejecting of suspected outlier.Such as, IGG I syllogic Modified Equivalent Weight Function form is selected, such as formula (11).

$W=\left\{\begin{array}{cc}1.0& M\&le;1.5\\ 1.5/M& 1.5<M\&le;2.5\\ 0& 2.5<M\end{array}\right.---\left(11\right)$

In above-mentioned formula, M is robust characteristic variable, and the specific implementation of M is as follows:

$M=\frac{(e_t-{\widehat{\mathrm{\&phi;}}}_t^{\left(k\right)}{e}_{t-1})}{{\mathrm{\&sigma;}}_v}=\frac{(Q_{ot}-Q_{ct})-(Q_t-Q_{ct})}{{\mathrm{\&sigma;}}_v}=\frac{(Q_{ot}-Q_t)}{{\mathrm{\&sigma;}}_v}---\left(12\right)$

${\mathrm{\&sigma;}}_v=\frac{\sqrt{\underset{i=1}{\overset{t}{\&Sigma;}}(e_t-{\widehat{\mathrm{\&phi;}}}_t^{\left(t\right)}{e}_{t-1}){\stackrel{\&OverBar;}{p}}_i^{\left(k\right)}(e_t-{\widehat{\mathrm{\&phi;}}}_t^{\left(k\right)}{e}_{t-1})}}{t-m-1}=\sqrt{\underset{i=1}{\overset{t}{\&Sigma;}}{{\stackrel{\&OverBar;}{p}}_i}^{\left(k\right)}\&CenterDot;{(Q_i-Q_{Oi})}^2)}/(t-m-1)---\left(13\right)$

Wherein, m, for eliminating sample number, is namely defined as the number of superseded data.

According to Error Theory, during error Normal Distribution, the probability of error outside ± 1.5 σ is 13%, and the probability outside ± 2.5 σ is only 0.7%.What robust characteristic variable M reflected is the ratio relation of error and mean square deviation (σ).All observation datas are divided three classes by the calculated value according to robust characteristic variable M.Effective observation data (M≤1.5), the equivalence weight factor is constant; Available observation data (1.5 < M≤2.5), reduces the equivalence weight factor, weakens its impact on valuation; Eliminate observation data (M > 2.5), compose zero power, do not participate in valuation and calculate, on valuation without impact.

It should be noted that, in the present invention in aforementioned robust Kalman filtering real time correction, the iteration initial value of the equivalence weight factor of robust Kalman filtering directly affects the speed of convergence of Iterated adjustment.And the flood that same basin occurs has certain similarity, real-time flood is the fluctuation within the specific limits of historical flood average, therefore the historical flood autoregressive model do not affected by exceptional value can be adopted to pass through the initial value of parameter estimation as the equivalence weight factor of least square acquisition by the gross, the initial value of such imparting can realize convergence fast, and counting yield is high.

By aforementioned to of the present invention illustrate known, the robust Kalman filtering real-time correction method that the present invention adopts, the residual information between observed volume and forecasting runoff is utilized to set up first order autoregressive model as measurement equation, state parameter according to autoregressive model sets up state equation, obtain parameter estimation according to robust Kalman filtering real-time correction method, and then obtain calibrated flow.Whether the initial weight factor of fresh information (new residual error) is 1.0, utilize the robust characteristic variable of Modified Equivalent Weight Function to judge in fresh information containing anomalous differences.If containing anomalous differences, reduce fresh information weight factor according to Modified Equivalent Weight Function, weaken it to the impact correcting result; If not containing anomalous differences, weight factor is constant, is similar to and is equal to standard Kalman filtering.Robust Kalman filtering real-time correction method utilizes robust Modified Equivalent Weight Function to improve fresh information weight factor, exceptional value observational error not being met to normal distribution reduces power even zero power, make it limited to valuation Influence on test result, high degree overcomes existing Kalman filtering real-time correction method affects excessive problem when processing non-normal errors by exceptional value, ensure that Real-time Flood Forecasting precision.

In order to further illustrate the present invention's attainable accurate forecast in flood forecasting flow, refer to the concrete example explanation as figure embodiment, in standard Kalman filtering of the prior art and the present invention, the real-time correction method of robust Kalman filtering all applies in the flood forecasting in Qi Lijie basin, the Min River, the Xinanjiang model flood forecasting flow specifically choosing 7 floods of 2007 ~ 2008 carries out real time correction, leading time is 3 hours, is next compared the calibration result of two kinds of methods by data.

(1) for the calibration result being as good as constant value situation

According to the forecasting runoff that Xinanjiang model obtains, use robust Kalman filtering (being called for short RKF as follows) real time correction in standard Kalman filtering (being called for short KF as follows) and the present invention respectively, as provided match value in non-correction accuracy, KF correction accuracy, RKF correction accuracy three kinds of situations and error total amount in following table 1 respectively.

The effect data table of each bearing calibration when table 1 is for being as good as constant value

Flood	Dc	Dc(KF)	Dc(RKF)	Er	Er(KF)	Er(RKF)
							080808	0.870	0.966	0.962	-0.029	-0.002	-0.001
080609	0.910	0.958	0.942	0.083	0.006	0.041
							080401	0.890	0.976	0.978	0.025	0.004	0.003
080315	0.813	0.965	0.959	-0.021	0.001	0.002
							070901	0.876	0.976	0.976	-0.096	-0.002	-0.003
070718	0.89	0.969	0.978	0.083	0.006	0
							070505	0.721	0.889	0.905	-0.097	-0.002	0.024
Mean value	0.853	0.957	0.957

Wherein, the 1st be classified as and count on 7 flood samples in table 1.Dc is the match value in non-correction accuracy situation, Dc (KF) is the match value in KF correction accuracy situation, Dc (RKF) is the match value in RKF correction accuracy situation, Er is the flood volume relative error in non-correction accuracy situation, Er (KF) is the flood volume relative error after KF correction, and Er (RKF) is the flood volume relative error after RKF corrects.

for measured discharge; Q is forecasting runoff; for measured discharge average,

q ₀for measured discharge; Q _kFfor KF calibrated flow; for measured discharge average;

q ₀for measured discharge; Q _rKFfor RKF calibrated flow; for measured discharge average;

$Er=\frac{\mathrm{\&Sigma;}(Q_{0i}-Q_i)}{{\mathrm{\&Sigma;Q}}_{0i}},$

$Er\left(KF\right)=\frac{\mathrm{\&Sigma;}(Q_{0i}-Q_{KFi})}{{\mathrm{\&Sigma;Q}}_{0i}},$

$Er\left(RKF\right)=\frac{\mathrm{\&Sigma;}(Q_{0i}-Q_{RKFi})}{{\mathrm{\&Sigma;Q}}_{0i}}.$

Known by the data analysis of table 1, when being as good as constant value, existing KF and RKF of the present invention two kinds of methods have good calibration result, and effect is substantially identical.

The effect data table of each bearing calibration when table 2 is for there being an exceptional value

Flood	Dc	Dc(KF)	Dc(RKF)	Er	Er(KF)	Er(RKF)
							080808	0.870	0.876	0.967	-0.029	-0.031	-0.002
080609	0.910	0.940	0.964	0.083	0.007	0.021
							080401	0.890	0.908	0.976	0.025	0.017	0.006
080315	0.813	0.822	0.961	-0.021	0.015	0.002
							070901	0.876	0.886	0.976	-0.096	-0.060	-0.003
070718	0.89	0.890	0.927	0.083	0.017	0.006
							070505	0.721	0.721	0.839	-0.097	0.003	0.009
Mean value	0.853	0.863	0.944

In table 2, the 1st is classified as and counts on 7 flood samples.Dc is the match value in non-correction accuracy situation, Dc (KF) is the match value in KF correction accuracy situation, Dc (RKF) is the match value in RKF correction accuracy situation, Er is the error total amount in non-correction accuracy situation, Er (KF) is the error total amount in KF correction accuracy situation, and Er (RKF) is the error total amount in RKF correction accuracy situation.

Table 2 shows the calibration result of two kinds of methods in exceptional value situation, and the average fit effect that robust Kalman filtering real-time correction method obtains is 0.944, and the average fit effect of standard Kalman filtering is 0.863.From the relative error of flood volume, result each flood corresponding of robust Kalman filtering real time correction is all little compared with standard Kalman filtering.Such as flood 080808, for being-0.029 through the relative error of overcorrect, the relative error after standard Kalman filtering and calibration is-0.031, and more correction error is larger.And the relative error after robust Kalman filtering of the present invention corrects is-0.002, after correcting, precision improves.Therefore, the error total amount of RKF provided by the invention is very little, and the calibration result of RKF is better than KF.Main cause is that RKF passes through to improve the equivalence weight factor, limits the impact of exceptional value on valuation, ensure that calibration result.

Next be illustrated the calibration result of robust Kalman filtering bearing calibration in the present invention, observed volume when being error without exception for big vast number 080808, Qo, manually add anomalous differences to Qo, the mode of interpolation is

$e_i=\left\{\begin{array}{cc}(r-0.5)*{\stackrel{\&OverBar;}{Q}}_O*3& i=\left(\mathrm{int}\right(i/40))*40\\ 0& i\&NotEqual;\left(\mathrm{int}\right(i/40))*40\end{array}\right.$

In above formula, r is the random number that satisfied (0,1) distributes, for the mean value of Qo.

The abnormal flow of t is Q _oCt=Q _ot+ e _t.

Through existing standard Kalman filtering and calibration, and through the flood discharge process that robust Kalman filtering of the present invention corrects, as shown in Fig. 2-b, be the calibration result analogous diagram of the flood discharge forecast that the embodiment of the present invention provides.In Fig. 2-b, the dotted line of long line segment composition is the real curve (dotted line that namely in figure, Qo represents) of flood discharge, the curve of solid line composition is the calibration curve (namely in figure with the solid line that Qkf represents) of the flood discharge after correcting according to existing standard Kalman filtering, the dotted line of short line segment composition is the calibration curve (namely in figure with the dotted line that Qrkf represents) of flood discharge after correcting according to robust Kalman filtering of the present invention, can determine accurately according to the discharge curve after two kinds of modes correct, the solid line that Qkf represents there will be the fluctuation of exceptional value at peak value place, and the influence of fluctuations of exceptional value has filtered by Qrkf, from flood peak situation, the flood peak that robust Kalman filtering corrects is closer to real traffic, forecast precision is higher.

Illustrate known by above content to of the present invention, Kalman filtering widely uses in real-time flood rectification, but when observational data error disobeys normal distribution, calibration result is subject to anomalous differences impact, not steadily and surely.Robust estimation theory is introduced in linear Kalman filter real time correction by the present invention, adopts robust M-estimator, improves the equivalence weight factor of anomalous differences, and limit its impact on calibration result, forecast precision is sane.

It should be noted that, for aforesaid each embodiment of the method, in order to simple description, therefore it is all expressed as a series of combination of actions, but those skilled in the art should know, the present invention is not by the restriction of described sequence of movement, because according to the present invention, some step can adopt other orders or carry out simultaneously.Secondly, those skilled in the art also should know, the embodiment described in instructions all belongs to preferred embodiment, and involved action and module might not be that the present invention is necessary.

For ease of better implementing the such scheme of the embodiment of the present invention, be also provided for the relevant apparatus implementing such scheme below.

Refer to shown in Fig. 3-a, the means for correcting 300 of a kind of flood discharge forecast that the embodiment of the present invention provides, can comprise: MBM 301, state estimation module 302, residual error forecast module 303, flux modification module 304, wherein,

MBM 301, for setting up first order autoregressive model to the residual information between observed volume and forecasting runoff, as measurement equation, the parameter according to described first order autoregressive model sets up state equation;

State estimation module 302, for adopting robust M-estimator from the state parameter valuation of described measurement equation and described state equation determination current time;

Residual error forecast module 303, for obtaining the residual error predicted value of subsequent time according to the state parameter valuation of the described current time determined and the residual information of current time;

Flux modification module 304, corrects for the forecasting runoff of residual error predicted value to subsequent time according to described subsequent time, obtains the calibrated flow of subsequent time.

In some embodiments of the invention, described MBM 301, specifically for setting up the measurement equation of described first order autoregressive model in the following way:

e _t＝φ _te _t-1+ε _t；

Wherein, described e _tfor the residual information of t, described e _t-1for the residual information in t-1 moment, described φ _tfor the state parameter of t, described ε _tfor the observation noise of t, described t is natural number;

Described e _tcalculate in the following way:

e _t＝Q _ot-Q _ct；

Wherein, described Q _otfor the observed volume of t, described Q _ctfor the forecasting runoff of t.

In some embodiments of the invention, described MBM 301, specifically for setting up described state equation in the following way:

φ _t＝Iφ _t-1+ω _t；

Wherein, described φ _tfor the state parameter of t, described I is the transition matrix of state equation and described I representation unit matrix, described φ _t-1for the state parameter in t-1 moment, described ω _tfor the process noise of t, described t is natural number.

In some embodiments of the invention, refer to as shown in Fig. 3-b, described state estimation module 302, comprising:

First computing module 3021, calculates robust Modified Equivalent Weight Function for the robust characteristic variable by designing in advance;

Second computing module 3022, for carrying out Iterated adjustment calculating according to described robust Modified Equivalent Weight Function parity price weight factor, when the equivalence weight factor does not meet preset power threshold condition in an iterative process, continue parity price weight factor and carry out Iterated adjustment calculating, only have when weight factor of equal value in iterative process meets described power threshold condition, export the robust gain matrix at the end of iteration;

3rd computing module 3023, for calculating the state parameter valuation of current time according to the robust gain matrix at the end of iteration.

Further, in some embodiments of the invention, refer to as shown in Fig. 3-c, described second computing module 3022, comprising:

First computing unit 30221, in the iterative process that walks in kth, according to the equivalence weight factor of (k-1) step of described robust Modified Equivalent Weight Function and t calculate the equivalence weight factor of the kth step of t described t and k is natural number;

Second computing unit 30222, for the equivalence weight factor walked according to the kth of described t the residual information in t-1 moment and the robust error covariance in t-1 moment calculate the robust gain matrix of t;

3rd computing unit 30223, for the equivalence weight factor walked according to the kth of described t the robust gain matrix of the residual information in t-1 moment and the robust error covariance in t-1 moment, described t calculates the robust error covariance of t;

4th computing unit 30224, calculates the state parameter valuation of the kth step of t for the residual information of the robust gain matrix according to t, t, the residual information in t-1 moment and the state parameter valuation in t-1 moment.

In some embodiments of the invention, in described iterative process, weight factor of equal value meets described power threshold condition, comprising: described μ is described power threshold condition, and the iteration of kth step is the final step that iteration terminates, and current time is t, the robust gain matrix that at the end of obtaining iteration, the equivalence weight factor pair of the kth step of t is answered; Described 3rd computing module 3023, the state parameter valuation of robust gain matrix, the residual information of current time, the residual information of the previous moment of current time and the previous moment of current time that the equivalence weight factor pair specifically for walking according to the kth of t at the end of iteration is answered calculates the state parameter valuation of current time.

In some embodiments of the invention, the criterion of described robust M-estimator is:

$({\widehat{\mathrm{\&phi;}}}_t{e}_{t-1}-e_t){\stackrel{\&OverBar;}{p}}_t^{\left(1\right)}({\widehat{\mathrm{\&phi;}}}_t{e}_{t-1}-e_t)+({\widehat{\mathrm{\&phi;}}}_t-{\widehat{\mathrm{\&phi;}}}_{t-1}){\left(P_{t-1}\right)}^{-1}({\widehat{\mathrm{\&phi;}}}_t-{\widehat{\mathrm{\&phi;}}}_{t-1})=min;$

Wherein, described in for the initial equivalence weight factor of t, described in as the weight of t observational error, described in for the valuation of t state parameter, described in for the valuation of t-1 moment state parameter, described e _tfor the residual information of t, described e _t-1for the residual information in t-1 moment, described P _t-1for the covariance of the state parameter sequence in t-1 moment, (P _t-1) ^-1as the weight of t-1 etching process error.

In some embodiments of the invention, described residual error forecast module 303, specifically for calculating the residual error predicted value of the subsequent time of current time in the following way;

${\stackrel{\&OverBar;}{e}}_{t+1}={\widehat{\mathrm{\&phi;}}}_t\&CenterDot;e_t;$

Wherein, described in for residual error predicted value when subsequent time is t+1, described in for the state parameter valuation of t, described e _tfor residual information when current time is t.

In some embodiments of the invention, described flux modification module 304, specifically for calculating the calibrated flow of subsequent time in the following way:

$Q_{t+1}=Q_{c(t+1)}+{\stackrel{\&OverBar;}{e}}_{t+1};$

Wherein, described Q _t+1for calibrated flow when subsequent time is t+1, described Q _{c (t+1)}for forecasting runoff when subsequent time is t+1, described in for residual error predicted value when subsequent time is t+1.

Known by the above embodiment description of this invention, first first order autoregressive model is set up to the residual information between observed volume and forecasting runoff, as measurement equation, and set up state equation according to the parameter of first order autoregressive model, then adopt robust M-estimator from the state parameter valuation of measurement equation and state equation determination current time, the state parameter valuation of the current time that following basis is determined and the residual information of current time obtain the residual error predicted value of subsequent time, and correct according to the forecasting runoff of residual error predicted value to subsequent time of subsequent time, obtain the calibrated flow of subsequent time.Robust M-estimator is introduced in linear Kalman filter real-time correction method by the present invention, utilizes robust M-estimator to reduce the state parameter valuation of unknown quantity as far as possible.Abnormal observational error is reduced as much as possible, weakens observational error to the impact of linear Kalman filter real-time correction method, ensure the correction accuracy of flood discharge forecast.

It should be noted that in addition, device embodiment described above is only schematic, the wherein said unit illustrated as separating component or can may not be and physically separates, parts as unit display can be or may not be physical location, namely can be positioned at a place, or also can be distributed in multiple network element.Some or all of module wherein can be selected according to the actual needs to realize the object of the present embodiment scheme.In addition, in device embodiment accompanying drawing provided by the invention, the annexation between module represents to have communication connection between them, specifically can be implemented as one or more communication bus or signal wire.Those of ordinary skill in the art, when not paying creative work, are namely appreciated that and implement.

Through the above description of the embodiments, those skilled in the art can be well understood to the mode that the present invention can add required common hardware by software and realize, and can certainly comprise special IC, dedicated cpu, private memory, special components and parts etc. realize by specialized hardware.Generally, all functions completed by computer program can realize with corresponding hardware easily, and the particular hardware structure being used for realizing same function also can be diversified, such as mimic channel, digital circuit or special circuit etc.But under more susceptible for the purpose of the present invention condition, software program realizes is better embodiment.Based on such understanding, technical scheme of the present invention can embody with the form of software product the part that prior art contributes in essence in other words, this computer software product is stored in the storage medium that can read, as the floppy disk of computing machine, USB flash disk, portable hard drive, ROM (read-only memory) (ROM, Read-OnlyMemory), random access memory (RAM, RandomAccessMemory), magnetic disc or CD etc., comprising some instructions in order to make a computer equipment (can be personal computer, server, or the network equipment etc.) perform method described in the present invention each embodiment.

In sum, above embodiment only in order to technical scheme of the present invention to be described, is not intended to limit; Although with reference to above-described embodiment to invention has been detailed description, those of ordinary skill in the art is to be understood that: it still can be modified to the technical scheme described in the various embodiments described above, or carries out equivalent replacement to wherein portion of techniques feature; And these amendments or replacement, do not make the essence of appropriate technical solution depart from the spirit and scope of various embodiments of the present invention technical scheme.

Accompanying drawing explanation

Fig. 1 provides a kind of process blocks schematic diagram of bearing calibration of flood discharge forecast for the embodiment of the present invention;

The bearing calibration of the flood discharge forecast that Fig. 2-a provides for the embodiment of the present invention realize principle schematic;

The calibration result analogous diagram of the flood discharge forecast that Fig. 2-b provides for the embodiment of the present invention;

The composition structural representation of the means for correcting of the flood discharge forecast that Fig. 3-a provides for the embodiment of the present invention;

The composition structural representation of the state estimation module that Fig. 3-b provides for the embodiment of the present invention;

The composition structural representation of the second computing module that Fig. 3-c provides for the embodiment of the present invention.

Embodiment

Embodiments providing a kind of bearing calibration to flood discharge forecast and device, for weakening the impact of Outliers error on linear Kalman filter real-time correction method, ensureing the correction accuracy of flood discharge forecast.

For making goal of the invention of the present invention, feature, advantage can be more obvious and understandable, below in conjunction with the accompanying drawing in the embodiment of the present invention, technical scheme in the embodiment of the present invention is clearly and completely described, obviously, the embodiments described below are only the present invention's part embodiments, and not all embodiments.Based on the embodiment in the present invention, the every other embodiment that those skilled in the art obtains, all belongs to the scope of protection of the invention.

Claims (10)

1. a bearing calibration for flood discharge forecast, is characterized in that, comprising:

Set up first order autoregressive model to the residual information between observed volume and forecasting runoff, as measurement equation, the parameter according to described first order autoregressive model sets up state equation;

Adopt robust M-estimator from the state parameter valuation of described measurement equation and described state equation determination current time;

The residual error predicted value of subsequent time is obtained according to the state parameter valuation of the described current time determined and the residual information of current time;

The forecasting runoff of residual error predicted value to subsequent time according to described subsequent time corrects, and obtains the calibrated flow of subsequent time.

2. method according to claim 1, is characterized in that, described employing robust M-estimator, from the state parameter valuation of described measurement equation and described state equation determination current time, comprising:

Robust Modified Equivalent Weight Function is calculated by the robust characteristic variable designed in advance;

Iterated adjustment calculating is carried out according to described robust Modified Equivalent Weight Function parity price weight factor, when the equivalence weight factor does not meet preset power threshold condition in an iterative process, continue parity price weight factor and carry out Iterated adjustment calculating, only have when weight factor of equal value in iterative process meets described power threshold condition, export the robust gain matrix at the end of iteration;

The state parameter valuation of current time is calculated according to the robust gain matrix at the end of iteration.

3. method according to claim 2, is characterized in that, describedly carries out Iterated adjustment calculating according to described robust Modified Equivalent Weight Function parity price weight factor, comprising:

In the iterative process of kth step, according to the equivalence weight factor of (k-1) step of described robust Modified Equivalent Weight Function and t calculate the equivalence weight factor of the kth step of t described t and k is natural number;

According to the equivalence weight factor that the kth of described t walks the residual information in t-1 moment and the robust error covariance in t-1 moment calculate the robust gain matrix of t;

According to the equivalence weight factor that the kth of described t walks the robust gain matrix of the residual information in t-1 moment and the robust error covariance in t-1 moment, described t calculates the robust error covariance of t;

The state parameter valuation of the kth step of t is calculated according to the residual information of the robust gain matrix of t, t, the residual information in t-1 moment and the state parameter valuation in t-1 moment.

4. according to the method in claim 2 or 3, it is characterized in that, the described state parameter valuation calculating current time according to the robust gain matrix at the end of iteration, comprising:

In described iterative process, weight factor of equal value meets described power threshold condition, comprising: described μ is described power threshold condition, and the iteration of kth step is the final step that iteration terminates, and current time is t, the robust gain matrix that at the end of obtaining iteration, the equivalence weight factor pair of the kth step of t is answered;

The state parameter valuation of robust gain matrix, the residual information of current time, the residual information of the previous moment of current time and the previous moment of current time that the equivalence weight factor pair walked according to the kth of t at the end of iteration is answered calculates the state parameter valuation of current time.

5. method according to claim 1, is characterized in that, the state parameter valuation of the described current time that described basis is determined and the residual information of current time obtain the residual error predicted value of subsequent time, comprising:

Calculate the residual error predicted value of the subsequent time of current time in the following way;

${{\stackrel{\&OverBar;}{e}}_t}_{+1}={\widehat{\mathrm{\&phi;}}}_t\&CenterDot;e_t;$

Wherein, described in for residual error predicted value when subsequent time is t+1, described in for the state parameter valuation of t, described e _tfor residual information when current time is t;

The described forecasting runoff of residual error predicted value to subsequent time according to described subsequent time corrects, and obtains the calibrated flow of subsequent time, comprising:

Calculate the calibrated flow of subsequent time in the following way:

$Q_{t+1}=Q_{c(t+1)}+{\stackrel{\&OverBar;}{e}}_{t+1};$

Wherein, described Q _t+1for calibrated flow when subsequent time is t+1, described Q _{c (t+1)}for forecasting runoff when subsequent time is t+1, described in for residual error predicted value when subsequent time is t+1.

6. a means for correcting for flood discharge forecast, is characterized in that, comprising:

MBM, for setting up first order autoregressive model to the residual information between observed volume and forecasting runoff, as measurement equation, the parameter according to described first order autoregressive model sets up state equation;

State estimation module, for adopting robust M-estimator from the state parameter valuation of described measurement equation and described state equation determination current time;

Residual error forecast module, for obtaining the residual error predicted value of subsequent time according to the state parameter valuation of the described current time determined and the residual information of current time;

Flux modification module, corrects for the forecasting runoff of residual error predicted value to subsequent time according to described subsequent time, obtains the calibrated flow of subsequent time.

7. device according to claim 6, is characterized in that, described state estimation module, comprising:

First computing module, calculates robust Modified Equivalent Weight Function for the robust characteristic variable by designing in advance;

Second computing module, for carrying out Iterated adjustment calculating according to described robust Modified Equivalent Weight Function parity price weight factor, when the equivalence weight factor does not meet preset power threshold condition in an iterative process, continue parity price weight factor and carry out Iterated adjustment calculating, only have when weight factor of equal value in iterative process meets described power threshold condition, export the robust gain matrix at the end of iteration;

3rd computing module, for calculating the state parameter valuation of current time according to the robust gain matrix at the end of iteration.

8. device according to claim 7, is characterized in that, described second computing module, comprising:

First computing unit, in the iterative process that walks in kth, according to the equivalence weight factor of (k-1) step of described robust Modified Equivalent Weight Function and t calculate the equivalence weight factor of the kth step of t described t and k is natural number;

Second computing unit, for the equivalence weight factor walked according to the kth of described t the residual information in t-1 moment and the robust error covariance in t-1 moment calculate the robust gain matrix of t;

3rd computing unit, for the equivalence weight factor walked according to the kth of described t the robust gain matrix of the residual information in t-1 moment and the robust error covariance in t-1 moment, described t calculates the robust error covariance of t;

4th computing unit, calculates the state parameter valuation of the kth step of t for the residual information of the robust gain matrix according to t, t, the residual information in t-1 moment and the state parameter valuation in t-1 moment.

9. the device according to claim 7 or 8, is characterized in that, in described iterative process, weight factor of equal value meets described power threshold condition, comprising: described μ is described power threshold condition, and the iteration of kth step is the final step that iteration terminates, and current time is t, the robust gain matrix that at the end of obtaining iteration, the equivalence weight factor pair of the kth step of t is answered;

Described 3rd computing module, the state parameter valuation of robust gain matrix, the residual information of current time, the residual information of the previous moment of current time and the previous moment of current time that the equivalence weight factor pair specifically for walking according to the kth of t at the end of iteration is answered calculates the state parameter valuation of current time.

10. device according to claim 6, is characterized in that, described residual error forecast module, specifically for calculating the residual error predicted value of the subsequent time of current time in the following way;

${\stackrel{\&OverBar;}{e}}_{t+1}={\widehat{\mathrm{\&phi;}}}_t\&CenterDot;e_t;$

Wherein, described in for residual error predicted value when subsequent time is t+1, described in for the state parameter valuation of t, described e _tfor residual information when current time is t;

Described flux modification module, specifically for calculating the calibrated flow of subsequent time in the following way:

$Q_{t+1}=Q_{c(t+1)}+{\stackrel{\&OverBar;}{e}}_{t+1};$

Wherein, described Q _t+1for calibrated flow when subsequent time is t+1, described Q _{c (t+1)}for forecasting runoff when subsequent time is t+1, described in for residual error predicted value when subsequent time is t+1.