Inverse Methods - Optimal Estimation

Estimating Surface/Internal Sources and Sinks Measurement Equation In Lagrangian framework: (0,0) Observing Station s',t', j Time v(s',t') s, t, l, k Back Trajectory Velocity Position Change in mole fraction y ( s, t ) from its initial condition y * ( 0, 0 ) is given by y ( s, t ) = y * ( s, t ) − y * ( 0, 0 ) t = ∫ xdt ' 0 x ds ' ⎞ ⎛ = ∫ ds ' ⎜v = ⎟ v dt ' ⎠ ⎝ 0 where x is the net chemical production (ignoring molecular diffusion and assuming perfect definition of back trajectory) s Real measurements of y have errors ε s xt 0 y ( s, t ) = ∫ ds '+ ε ( s, t ) v 0 Real measurement True value if perfect model Measurement error In discrete form: yi0 = ∑ h ij x tj + εi j In vector-matrix form for multiple observing sites i y0 = Hx t + ε measurement equation where H = ⎡⎣ h ij ⎤⎦ is the partial derivative (or observation) matrix, and ∂yi0 ∂yi ⎫⎪ = ⎬ compute in model ∂x tj ∂x j ⎪⎭ where ( yi , x j ) denotes an estimate of ( yi0 , x tj ) . h ij = Eulerian models: Here x refers to grid points in the Eulerian model rather than to points along a back trajectory. To utilize the above equations we must define: x = x ( grid point ) − x ( ref , grid point ) y = y ( grid point ) − y ( ref , grid point ) where y(ref) is the mole fraction computed in a reference run using best available estimates x ( ref ) of the state vector prior to estimating x t . The measurement equation expresses an apparent linear relation between the observation vector y0 and the unknowns contained in the state vector x t . It is an expression of the forward problem. The theory of the Linear Kalman Filter allows us to perform the inverse problem in which we optimally estimate x t given a model and a time series of observations y0 at one or more sites. Optimal estimation – the “cost” function (J) ∂J =0) ∂x ) = probability distribution function of ( ) ) Seek the estimate x of x t that minimizes J (i.e. (Notation: p ( (a) No knowledge of p ( x t ) , p ( y0 ) J = sum of squares of differences between predictions ( y = Hx ) and observations ( y 0 ) = ( y0 − Hx ) T (y 0 − Hx ) (called “least squares” minimization) (b) Know p ( y0 ) and hence R = expectation ⎡⎣εεT ⎤⎦ ( ) Then seek estimate x of x t that maximizes the conditional probability p y0 x t : J = ( y 0 − Hx ) R −1 ( y 0 − Hx ) T ( R −1 weights each squared difference by inverse of variance of each relevant observation) (called “weighted least squares” minimization or “maximum likelihood”) (c) Know p ( y0 ) and p ( x t ) ( ) Then seek x that maximizes p x t y0 which is defined from Bayes Theorem: p(y o | x t )p(x t ) p(y o ) (called the “minimum variance Bayes estimate” – see below for the definition of J in this case for the Kalman filter) p(x t | y o ) = KALMAN FILTER Suppose we have an optimal estimate of the state vector available prior to consideration of the kth measurement yko in a data series and we wish to obtain a new optimal estimate xka and its error νka using this measurement. We can thus propose in general xka = Κ k′ xkf + Κ k yko (13) νka = Κ k′ νkf + Κ k εk (14) and seek to specify the matrices Κ k and Κ k′ . Using the measurement equation (5) yo = Hxt + ε to define ykf = Ηk xkf , the definition xk = xkt + νk , the demand that the random measurement errors have zero mean (Ε[εk ] = 0) , and finally, the demand that estimations are unbiased (Ε[νka ] = 0) , we can show that Κ k′ = Ι − Κ k Ηk . Hence the new estimate is xka = xkf + Κ k [yko − Ηk xkf ] (15) with an error νka = (Ι − Κ k Ηk )νkf + Κ k εk (16) and estimation error covariance matrix Ρka = Ε[νka (νka )Τ ] (17) Substituting (16) into (17), using the definition of the measurement error covariance matrix Rk = Ε[εk εkT ] , and demanding that measurement errors and state errors are uncorrelated (so that Ε[νkf εkT ] = Ε[εk (νkf )Τ ] = 0 we obtain Ρka = (Ι − Κ k Ηk )Ρkf (Ι − Κ k Ηk )T + Κ kRk Κ kT (18) We now use the criterion of optimality to determine Κ k . Since we will assume we know p(yo ) and p(xt ) , we will choose a value for Κ k which minimizes the cost function J (equation 12) for the minimum variance Bayes estimate. Specifically Jk = Ε[(νka )T νka ] = trace [Ρka ] (19) Evaluating ∂Jk / ∂Κ k = 0 and solving for the so-called “Kalman Gain” matrix Κ k we have Κ k = Ρkf ΗkT [Ηk Ρkf ΗkT + Rk ]−1 (20) Substituting (20) into (18) then yields Ρka = [Ι − Κ k Ηk ]Ρkf (21) Finally, using the state space equation (7) x(t) = M(t, t o )x(t o ) + η(t, t o ) we then obtain the estimates of xkf needed in (15) and Ρkf needed in (21) xkf = Mk −1xka−1 (22) Ρkf = Μk −1Ρka−1ΜkT−1 + Qk −1 (23) where Qk −1 = Ε[ηk −1ηkT−1 ] , and xka −1 and Ρka −1 are the optimal outputs from the previous iteration of the filter. From our earlier discussion (Section 3), Q could represent random forcing in the system model due to transport model errors. To use the filter we must provide initial (a priori) estimates for x and P. Then from any prior output estimates (xka −1 , Ρka −1 ) , we use measurement k information ( yko ,Rk ) and model information ( Hk , Qk ) together with equations (22), (23), (20), (15), and (21) to provide outputs xka and Ρka for inputs to the next step. The filter equations are summarized in Table 1. Some intuitive concepts regarding the DKF are useful in understanding its operation. First, from equation (20), the gain matrix Κ k → Ηk−1 (its “maximum” value) as the measurement error covariance (noise) matrix R k → 0, and Κ k → Ρkf ΗkTRk−1 (its “minimum” value) as R k → ∞ . Since the update in the state vector xka − xkf varies linearly with Κ k , it is clear that measurements noisy enough so that R k much exceeds Ηk Ρkf ΗkT , will contribute much less to improvement of the state vector estimation. In this respect we can usefully consider Ηk Ρkf ΗkT as the error covariance matrix for the measurement estimates yk . This emphasizes the importance of the weighting of the data inherent in R k and the distortions created if erroneous R k are used. Note that R k can include model error, mismatch error, and instrumental error as noted earlier. Second, using (21), and recognizing that the maximum value of Κ kHk = Ι , we a k see Ρ ≤ Ρkf with equality occurring for infinitely noisy measurements. Hence, the error covariance matrix Ρk (whose diagonal elements are the variances of the state vector element estimates) decreases by amounts sensitively dependent on the measurement errors. Third, we note from (23), that random forcings η in the system (state-space) model [equation (7)], which are represented here by Q, will increase the extrapolated error covariance matrix Ρkf by amounts depending on the relative values of Qk −1 and the extrapolation matrix Μk −1Ρka−1ΜkT−1 in the absence of system (statespace) model noise. The inclusion of Q lessens the influence (or memory) of previous iterations in the filter operation. In the extreme, sufficiently large values of Q will prevent the capability of even non-noisy measurements to decrease Ρk and hence increase the confidence in the state vector estimate. In other words excellent (nonnoisy) measurements are of little use if the system (state-space) model is very noisy (e.g., through random variations η introduced by random transport errors). Table 1: Kalman Filter Equations* Definition Equation Measurement equation (model) yko = Hk xkt + εk ; System (state) equation (model) xk = Mk −1xk −1 + ηk −1 State update xka − xkf = Κ k (yko − yk ) Error Update Ρka = (1 − Κ k Ηk )Ρkf Kalman gain update Κ k = Ρkf ΗkT (Ηk Ρkf ΗkT + Rk )−1 State time extrapolation xkf = Μk −1xka −1 Error time extrapolation Ρkf = Μk −1Ρka −1ΜkT−1 + Qk −1 System random forcing covariance Qk = Ε(ηk ηkT ) Measurement error covariance Rk = Ε(εk εk T ) Estimation error covariance Ρk = Ε(νk νk T ) Input measurement matrix = Ηk = ∂yk / ∂xk Input system random forcing covariance = Qk Input state extrapolation = Μk Input measurement yko Input measurement error covariance = Rk Filter iteration − − − → (k − 1)f , → estimate yk = Hk xkf → (k − 1)a , → extrapolate → (k)f , → − − − ___________________________________________________________________ *A superscript a or superscript f denotes respectively the value before (f) or after (a) an update of an estimate using measurements, and k denotes the measurement number. In general, errors are assumed random with zero mean and measurement and estimation errors are uncorrelated.