- Email: [email protected]
COMPARISON OF ROBUST ALGORITHMS FOR AIRBORNE GRAVIMETRY
Preprints of the 5th IFAC Symposium on Robust Control Design
ROCOND'06, Toulouse, France, July 5-7, 2006
COMPARISON OF ROBUST ALGORITHMS FOR AIRBORNE GRAVIMETRY V.I. Kulakova,* A.V. Nebylov *, O.A. Stepanov** * State University of Aerospace Instrumentation 67, Bolshaya Morskaya, St. Petersburg, 190000, Russia email: [email protected] , [email protected]
**
Central Scientific & Research Institute Elektropribor 30, Malaya Posadskaya, St. Petersburg, 197046, Russia email: [email protected]
Abstract: The problem of a gravity anomaly (GA) estimation aboard an aircraft solved by using the data from a gravimeter and phase measurements of altitude from the differential satellite navigation system is discussed. The problem is considered from the viewpoints of different robust estimation algorithms: an H∞ estimation algorithm, a mixed H2 /H∞ estimation algorithm, a suboptimal algorithm based on unsteady models represented as integrals of a white noise and a cost-guaranteed algorithm based on the information about the GA variance and the variance of its derivative. A comprehensive comparison of robust estimators is presented. In addition, robust estimators are compared with the optimal estimator, which is obtained for the widely used GA spectral density. Keywords: robust and optimal estimation, accuracy, sensitivity, airborne gravimetry.
1. INTRODUCTION The survey of a gravity anomaly (GA) from an aircraft (airborne gravimetry) is widely used now. Airborne gravimetry is a highly productive and inexpensive method of survey, which can be employed even in hard-to-reach areas (Blaznov, et al., 2002; Stepanov, et al., 2002). However, GA measuring aboard a moving vehicle involves two main problems. The first is the problem of compensation for the vertical acceleration caused by aircraft dynamics. This problem is of vital importance for airborne gravimetry as, due to high speed of the aircraft, the spectrum of the signal to be estimated is displaced to a high frequency domain. Hence, the spectrum of GA and the spectrum of the vertical accelerations overlap and their separation by filtering is impossible. Consequently, external high-
precision information about the aircraft dynamics is needed to exclude an unknown component of the aircraft vertical acceleration. Compensation for the vertical accelerations is based on the data about the aircraft altitude from the differential satellite navigation system (DSNS). Subcentimeter altitude accuracy is achieved by using phase measurements (Stepanov, et al, 2002). The second problem consists in specifying a stochastic model for the GA. A lot of various GA spectral densities are used to describe it (Jordan, 1972). At the same time the adequacy of the model is particularly relevant for airborne gravimetry, as the GA changes with a change of altitude (Afanasieva, et al., 2002). Note that only two parameters are usually preset for the GA: the variance of GA and the GA gradient (or the variance of the GA derivative).
Preprints of the 5th IFAC Symposium on Robust Control Design
ROCOND'06, Toulouse, France, July 5-7, 2006
The aim of this paper is to investigate the efficiency of applying different robust approaches to the problem of airborne gravimetry. In order to take a better account of the properties of robust estimators the latter will be compared with the optimal estimator (Kalman filter) (Van Trees, 1968) obtained for the widely used GA spectral density. The following robust filters will be examined: an H∞ filter (Nagpal and Khargonekar, 1991; Grimble and Elsayed, 1990), which minimizes the maximum peak of the error spectral density and thus possesses an inherent robustness to a possible change of the properties of input signals; a mixed H 2 /H∞ filter (Khargonekar, et al., 1996; Rotstein and Sznaier, 1996), which is expected to perform better than either a Kalman filter or an H∞ filter. Then the problem will be simplified by straightening in the logarithmic scale of the GA spectral density in the vicinity of its cross point with measurement errors and thus the suboptimal algorithm will be derived. In this case the input signal for the filter is represented in the form of integrals of the white noise (Chelpanov, et.al., 1978). The solutions of applied problems show that the application of models in the form of integrals of the white noise to the description of signals is very efficient because it allows the development of reliable and simple estimation algorithms (Tupysev, 2004). And at last, instead of the spectral density to describe the GA, a priori information about the GA will be represented as two parameters: the variance of GA and the variance of its derivative. In this case only these numerical characteristics can be used in designing an estimator (Nebylov, 2004; Kulakova, et al., 2005). In the subsequent discussion such estimator will be called a cost-guaranteed estimator (the reason for choosing this term will be substantiated a bit later).
hSNS+ δh g~ + δg + h&&
s
The GA is often described by a process with the spectral density (Jordan, 1972):
where
σ 2g~
g~ gr g~ + δg SNS − h = − δh . s2 s2
(1)
The estimation problem consists in obtaining GA using differential measurement (1), i.e. in the filter transfer function H(s) designing. For this purpose it is necessary to specify models for g~ , δh, and δg .
(ω 2 + α 2 ) 3
is the variance;
,
(2)
α = V / µ; V − the
vehicle speed; µ - the correlation distance. If the GA gradient ∇g~ = σ ∂g~ / ∂ l (rms value of the gravity increment at a distance of l , which is assumed to be 1 km) is given, α can be found from the equation (∇g~ )2 = 2α 2 σ 2~ . GA spectral density (2) is shown in g
Fig. 2. It is not difficult to show that this process can be described by using a third-order shaping filter: g~& 1 (t ) = −αg~1 (t ) + g~2 (t ), g~& (t ) = −αg~ (t ) + g~ (t ), 2
3
g~& 3 (t ) = − αg~3 (t ) + 10α 3σ 2g~ ν(t ), g~(t ) = −αζg~ (t ) + g~ (t ), 1
2
with ζ = ( 5−1) 5 , where ν(t ) is the is the zeromean white noise with the unit intensity. 10-2
S
(1 ) g~ r
(ω )
S
g~ s
(ω )
10-4 10-6
(m/s2)2s
The difference between the second integral of the gravimeter readings and the altitude from the DSNS is formed in order to eliminate unknown altitude h (see Fig. 1). The differential measurement can be represented as (Kulakova, et al., 2004):
5 ⋅ ω2 + α 2
S g~ (ω) = 2α 3 ⋅ σ 2g~ ⋅
The gravimeter readings can be represented as a sum of three components: g~ gr = g~ + δg~ + h&& . Here g~ is
h SNS = h + δh , here h is unknown altitude and δh is the errors in DSNS measurements.
g~ˆ
H(s)
The survey data for the investigations was obtained with the use of the gravimeter developed in the CSRI Elektropribor (Blazhnov, et al., 2002; Stepanov, et al, 2002) and dual-frequency geodetic Novatel receivers. It has been shown that the errors in phase measurements are mainly of a white-noise character with the intensity Rh =(0.005m)2s and the gravimeter errors can be described as a white noise with the intensity R gr = (5mGal)2s.
2
gravity anomaly; δg is the gravimeter error and h&& is the vertical acceleration caused by the vertical motion of the aircraft. The data from the DSNS are
zh
2
Fig.1. A block diagram for the GA estimation.
2. STATEMENT OF THE PROBLEM
zh =
-
1
S g~ ( ω )
10-8
S
( 0 ,1 ) g~ r
(ω )
Rh ω4
10-10 10-12 10-14 10-6
10-5
10-4
10-3
.
10-2
10-1
1
Rad/s Fig.2. GA spectral densities and measurement errors for ∇g~ =10 mGal/km, σ g~ =30 mGal, V =50 m/s.
Preprints of the 5th IFAC Symposium on Robust Control Design
ROCOND'06, Toulouse, France, July 5-7, 2006
[
3. ROBUST SOLUTIONS OF THE PROBLEM Under the assumptions made the problem of the GA estimation reduces to estimation of g~ using equations as follows: h&(t ) = V z (t ), V&z (t ) = g~(t ) + R gr n gr (t ),
where E is the mathematical expectation, e(t ) = x(t ) − xˆ (t ) is the estimation error of x(t ) . The Kalman filter provides minimum error variance estimation for every component of the vector x . The GA value itself can be estimated as zˆ(t ) = C1 xˆ (t ) .
z − zˆ
y (t ) = h(t ) + Rh n h (t ),
[
2 2
]T ,
where VZ (t ) is the vertical velocity, n gr (t ) is the
where n = n gr
noise that describes the gravimeter error, y (t ) is the
noise attenuation and
differential measurement, n h (t ) is a measurement noise caused by the errors in DSNS measurements. For the GA model (2) the estimation problem is formulated in this way: x&1 (t ) = x 2 (t ),
nh
< γ2 ⎡ n ⎢⎣
(3)
x& 5 (t ) = −αx 5 (t ) + 10α 3σ g~2 ν (t ), y (t ) = x1 (t ) + Rh n h (t ),
(7)
γ > 0 is a prescribed level of
⋅ 2 denotes the norm in L2 [0, ∞ ) . It is essential that for the H∞-filter the noises in (4) are modeled as signals with a bounded energy and an arbitrary spectrum. Note that the H∞ filter provides optimal estimation in sense of (7) only for the signal z (t ) = C1 x(t ) .
T = [T~z ν
T~z n ] ,
(8)
T~z ν ( s ) = C1 (sI − A1) B2 , T~z n ( s ) = C1 (sI − A1) Bn −1
−1
2
where x(t ) is the state, ν(t ) is the scalar process noise, z (t ) is the signal to be estimated from the measurement y (t ) . The system (3) can be written in the following form: + B2 ν(t ),
(4)
y (t ) = C 2 x (t ) + Dn h (t ) ,
where B 2T = (0,0,0,0, 10α 3 σ 2g~ ) ,
C1 = (0,0,−αξ ,1,0) , C 2 = (1,0,0,0,0), D = Rh .
Both the Kalman and H∞ filters are designed for the GA model (2). These filters are of the form: xˆ& (t ) = Axˆ (t ) + K [y (t ) − C 2 xˆ (t )], zˆ (t ) = C1 xˆ (t ),
2 + ν 2⎤ , ⎥⎦
with A1 = A − KC 2 and Bn = [B1 − KD ] . Then the performance measurer (7) can be rewritten as
z (t ) = −αζ x3 (t ) + x 4 (t ),
B1T = (0, R gr ,0,0,0) ,
2 2
Let T denote the transfer matrix from noises to estimation error ~ z = z − zˆ :
x& 2 (t ) = −αζ x3 (t ) + x 4 (t ) + R gr n gr (t ),
x& (t ) = Ax (t ) + B1 n gr (t ) z (t ) = C1 x (t ),
(6)
The H∞-filter is designed to provide
g~& (t ) → model,
x& 3 (t ) = −αx3 (t ) + x 4 (t ), x& 4 (t ) = −αx 4 (t ) + x5 (t ),
]
De = E e T (t )e(t ) ,
(5)
where xˆ (t ) is the estimate of x(t ) , zˆ (t ) is the estimate of z (t ) and K is the filter gain. For the Kalman filter the noises n gr (t ) , ν(t ) , n h (t ) are assumed to be zero-mean unit intensity white noise signals. The gain K is chosen to achieve the minimum for the following criterion:
T ∞ < γ2 ,
where
T
∞
(9)
denotes the H∞-norm of the transfer
function T . Thus the H∞-filter minimizes the H ∞ norm of the transfer function from noises to the estimation error. Note that the Kalman filter minimizes the H 2 -norm of this transfer function. In the case when the noises in (4) are standard white noise signals the performance measurer (9) has such a meaning that the maximum peak of the error spectral density will not exceed γ 2 . It is important for many applications as the Kalman filter, which minimizes the integral of this spectral density, often leads to a large error spectrum peak that causes bad robustness. It should be noted that choosing the value of γ is a trade-off between the H 2 and H ∞ performances for the H ∞ -filter: near γ 0 , which is a minimal value of γ for which (9) can be met, the error spectrum peak will be reduced to the lowest possible level, at γ → ∞ the H ∞ filter will reduce to Kalman filter. The rms (root-mean-square) estimation errors for the Kalman filter for different values of GA gradient – 3, 5 and 10 mGal/km are presented in Table 1. In calculations it was assumed that σ g~ =30mGal. The aircraft speed was assumed to be equal to V =50m/s, which is typical for an aerogravimetric survey.
Preprints of the 5th IFAC Symposium on Robust Control Design
ROCOND'06, Toulouse, France, July 5-7, 2006
Table 1. Rms estimation errors for GA, mGal
6 . 10
8
3
5
10
5 . 10
8
Kalman filter, σ e
2
3.1
5.5
4 . 10
8
H ∞ -filter for γ ∗ , σ e
3.2
4.4
8.4
3 . 10
8
H 2 / H ∞ filter, σ e
2.9
4.0
6.0
Suboptimal filter, σ es
2.2
3.6
6.9
Suboptimal filter, σ 'es
2.1
3.4
6.4
Cost-guaranteed filter, σ e
2.8
4.1
6.7
Cost-guaranteed filter, σ er
2.8
3.9
6.4
During the investigation of the H ∞ -filter it was found that for γ near γ 0 the H ∞ -filter provides an inadmissible large error (about 25 mGal) when the system (4) is driven by the white noise signals (this means the estimation error that provides the H ∞ filter for the GA described by model (2)). That’s why the value γ was increased in order to achieve a reasonable estimation error. In Fig. 3 the error spectral densities are shown for the H ∞ -filters designed for γ = γ 0 and γ = γ ∗ , where γ ∗ is a chosen level of noise attenuation. Fig. 3 corresponds to the value of the GA gradient 10 mGal/km, in this case γ 0 ≅ 2 ⋅ 10
−8
(m / s ) s 2 2
(
(the error is about 25
)
2
mGal ) and γ * ≅ 2.5 ⋅ 10 −8 m / s 2 s (the error is about 8 mGal). In addition, Fig. 3 shows the error spectral density for the Kalman filter. It is clearly seen in the picture that decrease of the value γ results in decrease of the error spectrum peak on the one hand, but in broadening of the error spectra on the other hand. The estimation errors obtained for the H ∞ -filter for γ ∗ are presented in Table 1. Note that for γ > γ 0 there exist many filters that achieve (9), and a specific filter may be chosen to improve the H 2 performance of the filter. The problem of H 2 performance optimization subject to H ∞ constraint is referred to as a mixed H 2 / H ∞ problem (Khargonekar, et al., 1996). At the moment no analytical solution to this problem is known. In the present paper the mixed H 2 / H ∞ filtering problem will be considered in this way: find a gain K of the filter (5) that solves:
{
}
inf T 2 : subject to T ∞ < γ ,
(10)
where γ is given. The value γ for the H 2 / H ∞ filter synthesis is chosen approximately in the middle between γ 0 and γ K , where γ K is the H ∞ constraint for the Kalman filter (according to Fig. 3
(m/s2)2s
∇g~ , mGal/km
∇g~ σ g~
4
= 10 mGal/km = 30 mGal
3 2
1
8 2 . 10
1 . 10
8
0
0.2 .
0.1
0.3
0.4
Rad/s
Fig.3. Error spectral densities: 1 – the H∞-filter for γ0; 2 - the H∞-filter for γ*; 3 – the H2/H∞ filter; 4 – the Kalman filter. 2.5
∇ g~
4 2
3
2
σ g~
= 10 mGal/km = 30 mGal
1.5
1
5
0.5
0
1
1 3
0
0.1
0.2
. Rad/s
0.3
0.4
Fig.4. The AFC for the vertical acceleration: 1 – the H∞-filter; 2 – the H2/H∞ filter; 3 – the Kalman filter; 4 – the suboptimal filter; 5 – the cost-guaranteed filter.
(
)
2
γ K ≅ 5 ⋅ 10 −8 m / s 2 s ). Fig. 3 shows the error spectral density corresponding to the H 2 / H ∞ filter. The rms errors for the H 2 / H ∞ filter are presented in Table 1. The amplitude-frequency characteristics (AFC) ~ H ( s ) = H ( s ) s 2 for these three filters are shown in Fig. 4. Note that all these filters are designed for the GA model (2) and they have the fifth order. It is clear that the accuracy obtained as well as the structure of the filter depend on the type of the models used for the measurement errors and GA. The models describing the measurement errors were derived from the analysis of the real data. Model (2) is widely used in the problems that require a stochastic description of the GA. At the same time many examples show that there is a strong reason for adjusting filters not to real models of signals to be estimated but to unsteady models represented as integrals of a white noise (Tupysev, 2004). The fact is that it allows improving filter robustness and achieving guaranteed accuracy. One of the possible ways for specifying a model of a signal as integrals
Preprints of the 5th IFAC Symposium on Robust Control Design
ROCOND'06, Toulouse, France, July 5-7, 2006
of a white noise is to rectify (to replace by straight line) the spectral density of the signal in the vicinity of its cross point with error spectral density (Chelpanov, et.al., 1978). The filter which is obtained for the rectified spectral density of the real signal is referred to as a suboptimal filter. The GA spectral density (2) can be approximated, in the logarithmic scale, in the vicinity of its cross point with ω 4 Rh (see Fig. 2) by the second integral of the white noise with the intensity q g2~ = 10α 3 σ 2g~ : S g~s (ω) ≈ q g~ / ω 4 .
(11)
Thus, the suboptimal filter is a Kalman filter that uses (11) as the GA model. It is easy to show that the suboptimal filter transfer function can be represented as H s ( s ) = s 2W B ( s ) , where WB (s ) corresponds to the Butterworth 4th-order filter with the cutoff frequency
(q g~
)
Rh 1 6 .
The
AFC
of
described by the model (2) σ 'es are also presented. The final step will be rejection of the GA description by any spectral density and representation of a priori information about the GA as two parameters: the variances of GA and its derivative: (12)
It is suggested that only these numerical characteristics are known (Kulakova, et al., 2005). In the case that the signal has bounded variances of the signal itself and its derivative the following expression holds for its spectral density S (ω) : 1 2π
∫ [c0 S (ω) + c1ω
∞
2
]
S (ω) dω = c 0 D0 + c1 D1 ,
−∞
where c 0 , c1 are the real, positive, arbitrary parameters. Then S (ω) can be written as: S (ω) = S r(0,1) S w (ω) ,
(13)
with S r(0,1) =
c 2 D0 + D1 c 2 + ω2
and
1 2π
∞
∫ S w (ω)dω = 1 ,
−∞
where c 2 = c 0 c1 . It is true because the spectral S w (ω) can be chosen as density S w (ω) = S (ω) / S r(0,1) .
x& 2 (t ) = x3 (t ) + R gr ngr (t ), x& 3 (t ) = − cx 3 (t ) + c 2 D0 + D1 w(t ),
Thus the signal with the bounded variances (12) can be obtained from the signal with unit variance and unknown spectrum. It is possible to show that the problem of the GA estimation when only numerical characteristics (12) are known can be formulated as follows:
(14)
z (t ) = x3 (t ), y (t ) = x1 (t ) + Rh n h (t ),
where n gr (t ) , n h (t ) are assumed to be standard white noise signals, w(t ) is the zero mean signal with unit variance and unknown spectrum. The system (14) can be written in the form similar to (4). In the case that the spectrum of the process noise w is unknown it is impossible to find error variance. However, since w is of bounded power, it is possible to find an upper bound on the error variance and thus to design a filter that minimizes this upper bound. This means to find a filter that provides:
{
inf T~z w
this
Butterworth filter is shown in Fig. 4. The estimation errors for the suboptimal filter σ es are presented in Table 1. It should be noted that the suboptimal filter is adjusted to model (11) that cannot represent the real dynamics of the GA signal. That’s why the errors that provides suboptimal filter for the signal
D0 = σ 2g~ , D1 = (∇g~ )2 .
x&1 (t ) = x 2 (t ),
2 ∞
+ T~z n
2 2
}.
(15)
The simplest admissible filter in this case has a structure like the one in (5) (Khargonekar, et al., 1996) and it is filter (5) that is used. Thus, the function under optimization J ( K , c ) = T~z w
2 ∞
+ T~z n
2 2
depends on the filter gain
matrix K and the parameter c in the plant model. Their optimal values depend on the specific values D0 , D1 , Rh , R gr and can only be obtained during the function J ( K , c ) minimization. As the upper bound for the error variance is minimized for such filter, it will provide the guaranteed estimation of the signals which are represented by the numerical characteristics (12). That’s why this filter is referred to as a costguaranteed filter. It should be noted that problem (15) can also be regarded as a mixed H 2 / H ∞ problem but formulated in a way different from (10). It is important to understand that problem (15) was derived from other assumptions and for other purposes. In Table 1 the cost-guaranteed filter is characterized by two parameters: the upper bound of the estimation error σ e (the error value which will not be exceeded for the entire class of input signals with the preset values D0 , D1 accepted in designing a filter) and the error σ er that provides the cost-guaranteed filter for the GA described by model (2). It is interesting to note that the parameter D0 had little influence on the accuracy of the GA estimation. This means that if only the parameter D1 is preset the upper bound σ e will vary slightly in comparison with the case when two parameters are preset. In the case that only D1 is used the GA spectral density is
Preprints of the 5th IFAC Symposium on Robust Control Design
represented
as
S (ω) = S r(1) S w (ω) ,
ROCOND'06, Toulouse, France, July 5-7, 2006
where
S r(1) = D1 / ω 2 ( c = 0 ). This means that in this case the cost-guaranteed filter is adjusted to random walk.
Fig. 2 shows the spectral densities S r(1) and S r(0,1) , with the optimal parameter c for S r(0,1) (the value of c for which the upper bound on the error variance is minimal). In order to estimate the weight of the variances (13), special consideration must be given to how close the curves S r(1) and S r(0,1) are in the vicinity of their cross point with the measurement errors. It is clear that the curve
S r(1) represents the
rectified curve S r(0,1) in the logarithmic scale at its cross point with the spectral density of the measurement errors. Thus, the cost-guaranteed filter is mainly determined by the variance of the GA derivative. Fig. 2 also shows that the cost-guaranteed filter is adjusted to the spectral density that has a wider frequency bound. Therefore, it will be lowsensitive to the change of the acceptable model that causes extension of the spectral density. Therewith, the cost-guaranteed filter will provide a higher error in comparison with the filter adjusted to model (2) when the GA spectral density is narrowed. The cost-guaranteed filter AFC is shown in Fig. 4 from which it follows that the filter turns out to be similar to the classical low-pass filter. Its AFC has flat vertex and an inclination of 60 dB/dec in the high-frequency domain. The existence of peaks for the AFC of the other filters makes them too sensitive to an increased noise power at the frequencies where the AFC is high. It is also sensitive to the change of the accepted model that causes decrease of the valid signal energy in this frequency range. Therefore, from the intuitive point of view, the cost-guaranteed estimator is low-sensitive. Although no special measures were taken with this aim. CONCLUSIONS Five different filtering approaches to GA estimation aboard an aircraft have been analyzed and compared with respect to the appropriate a priori information and their numerical performance. Several tests have been carried in order to investigate their efficiency for various operating conditions. The results obtained suggest that the cost-guaranteed filter will provide reliable and efficient estimation of a GA. The costguaranteed filter makes use of two parameters (the GA variance and the GA gradient) which can be easily obtained from the available experimental data; it has the features similar to those of the mixed H 2 / H ∞ filter; it is adjusted to the model represented as the integral of the white noise; the cost-guaranteed filter loses to the optimal filter in accuracy to a small extent. In conclusion, it should be noted that the paper presents a novel approach to the cost-guaranteed filter designing, the one when only numerical characteristics are given for signals.
REFERENCES Afanasjeva, L.V., E.A. Bojarskij, N.V. Drobyshev and V.N. Konushev (2002). About the errors of gravimetry anomaly images obtained from the aircraft readings. Geophizicheskij vestnik. No 4, pp. 13-15 (in Russian). Blazhnov, B.A., L.P. Nesenjuk, V.G. Peshekhonov, A.V. Sokolov, L.S. Elinson, L.K. Zhelesnyak (2002). An integrated mobile gravimetric system. Development and test results. Proceedings of the 9th Saint-Peterburg International Conference on Integrated Navigation Systems. Chelpanov, I.B., L.P. Nesenjuk and М.V. Braginsky (1978). Computation of navigational gyro devices characteristics. Sudostroyeniye, Leningrad (in Russian). Grimble, M.J. and A. Elsayed (1990). Solution of the H∞ optimal linear filtering problem for discretetime systems. IEEE Trans. Acoustics, Speech, Signal Processing, Vol. 38, pp.1092-1104. Jordan, S.K. (1972). Self-consistent statistical models for gravity anomaly and undulation of the geoid. Geophys. Res. Vol. 77, No 20, pp. 3660-3670. Khargonekar, P.P., M.A. Rotea and E. Baeyens (1996). Mixed H2/H∞ filtering. Internat. J. Robust Nonlinear Control, Vol. 6, pp. 313–330. Kulakova, V.I., A.V. Nebylov and O.A. Stepanov (2005). Robust estimation versus optimal estimation and their application for airborne gravimetry. Proceedings of the 16th IFAC World Congress. Prague. Kulakova, V.I., A.V. Nebylov and O.A. Stepanov (2004). Application of the robust approach to the problem of airborne gravimetry. Proceedings of the 16th IFAC Symposium on Automatic Control in Aerospace. St. Petersburg, Russia. Pp.354-359. Nagpal, K.M. and P.P. Khargonekar (1991). Filtering and Smoothing in an H∞ Setting. IEEE Trans. Automatic Control, Vol. 36, pp.152-166. Nebylov, А.V. (2004). Ensuring control accuracy. Lecture Notes in Control and Information Sciences, Springer-Verlag, Heidelberg, Germany. Rotstein, H. and M. Sznaier (1996). H2/H∞ filtering theory and an aerospace application. Internat. J. Robust Nonlinear Control, Vol. 6, pp. 347–366. Stepanov, O.A., B.A. Blazhnov and D.A. Koshaev (2002). The efficiency of using velocity and coordinate satellite measurements in determining gravity aboard an aircraft. Proceedings of the 9-th Saint Petersburg International Conference on Integrated Navigation Systems. St. Petersburg, Russia. Pp.255-264. Tupysev, V.A. (2004). Guaranteed estimation of the dynamic system state in the case of uncertainty in the description of disturbance and measurement errors. Proceedings of the 11th Saint Petersburg International Conference on Integrated Navigation Systems. St. Petersburg, Russia. Pp.86-89. Van Trees, H. (1968). Detection, estimation and modulation theory. Part I. Detection, estimation and linear modulation theory. New-York.
ROCOND'06, Toulouse, France, July 5-7, 2006
COMPARISON OF ROBUST ALGORITHMS FOR AIRBORNE GRAVIMETRY V.I. Kulakova,* A.V. Nebylov *, O.A. Stepanov** * State University of Aerospace Instrumentation 67, Bolshaya Morskaya, St. Petersburg, 190000, Russia email: [email protected] , [email protected]
**
Central Scientific & Research Institute Elektropribor 30, Malaya Posadskaya, St. Petersburg, 197046, Russia email: [email protected]
Abstract: The problem of a gravity anomaly (GA) estimation aboard an aircraft solved by using the data from a gravimeter and phase measurements of altitude from the differential satellite navigation system is discussed. The problem is considered from the viewpoints of different robust estimation algorithms: an H∞ estimation algorithm, a mixed H2 /H∞ estimation algorithm, a suboptimal algorithm based on unsteady models represented as integrals of a white noise and a cost-guaranteed algorithm based on the information about the GA variance and the variance of its derivative. A comprehensive comparison of robust estimators is presented. In addition, robust estimators are compared with the optimal estimator, which is obtained for the widely used GA spectral density. Keywords: robust and optimal estimation, accuracy, sensitivity, airborne gravimetry.
1. INTRODUCTION The survey of a gravity anomaly (GA) from an aircraft (airborne gravimetry) is widely used now. Airborne gravimetry is a highly productive and inexpensive method of survey, which can be employed even in hard-to-reach areas (Blaznov, et al., 2002; Stepanov, et al., 2002). However, GA measuring aboard a moving vehicle involves two main problems. The first is the problem of compensation for the vertical acceleration caused by aircraft dynamics. This problem is of vital importance for airborne gravimetry as, due to high speed of the aircraft, the spectrum of the signal to be estimated is displaced to a high frequency domain. Hence, the spectrum of GA and the spectrum of the vertical accelerations overlap and their separation by filtering is impossible. Consequently, external high-
precision information about the aircraft dynamics is needed to exclude an unknown component of the aircraft vertical acceleration. Compensation for the vertical accelerations is based on the data about the aircraft altitude from the differential satellite navigation system (DSNS). Subcentimeter altitude accuracy is achieved by using phase measurements (Stepanov, et al, 2002). The second problem consists in specifying a stochastic model for the GA. A lot of various GA spectral densities are used to describe it (Jordan, 1972). At the same time the adequacy of the model is particularly relevant for airborne gravimetry, as the GA changes with a change of altitude (Afanasieva, et al., 2002). Note that only two parameters are usually preset for the GA: the variance of GA and the GA gradient (or the variance of the GA derivative).
Preprints of the 5th IFAC Symposium on Robust Control Design
ROCOND'06, Toulouse, France, July 5-7, 2006
The aim of this paper is to investigate the efficiency of applying different robust approaches to the problem of airborne gravimetry. In order to take a better account of the properties of robust estimators the latter will be compared with the optimal estimator (Kalman filter) (Van Trees, 1968) obtained for the widely used GA spectral density. The following robust filters will be examined: an H∞ filter (Nagpal and Khargonekar, 1991; Grimble and Elsayed, 1990), which minimizes the maximum peak of the error spectral density and thus possesses an inherent robustness to a possible change of the properties of input signals; a mixed H 2 /H∞ filter (Khargonekar, et al., 1996; Rotstein and Sznaier, 1996), which is expected to perform better than either a Kalman filter or an H∞ filter. Then the problem will be simplified by straightening in the logarithmic scale of the GA spectral density in the vicinity of its cross point with measurement errors and thus the suboptimal algorithm will be derived. In this case the input signal for the filter is represented in the form of integrals of the white noise (Chelpanov, et.al., 1978). The solutions of applied problems show that the application of models in the form of integrals of the white noise to the description of signals is very efficient because it allows the development of reliable and simple estimation algorithms (Tupysev, 2004). And at last, instead of the spectral density to describe the GA, a priori information about the GA will be represented as two parameters: the variance of GA and the variance of its derivative. In this case only these numerical characteristics can be used in designing an estimator (Nebylov, 2004; Kulakova, et al., 2005). In the subsequent discussion such estimator will be called a cost-guaranteed estimator (the reason for choosing this term will be substantiated a bit later).
hSNS+ δh g~ + δg + h&&
s
The GA is often described by a process with the spectral density (Jordan, 1972):
where
σ 2g~
g~ gr g~ + δg SNS − h = − δh . s2 s2
(1)
The estimation problem consists in obtaining GA using differential measurement (1), i.e. in the filter transfer function H(s) designing. For this purpose it is necessary to specify models for g~ , δh, and δg .
(ω 2 + α 2 ) 3
is the variance;
,
(2)
α = V / µ; V − the
vehicle speed; µ - the correlation distance. If the GA gradient ∇g~ = σ ∂g~ / ∂ l (rms value of the gravity increment at a distance of l , which is assumed to be 1 km) is given, α can be found from the equation (∇g~ )2 = 2α 2 σ 2~ . GA spectral density (2) is shown in g
Fig. 2. It is not difficult to show that this process can be described by using a third-order shaping filter: g~& 1 (t ) = −αg~1 (t ) + g~2 (t ), g~& (t ) = −αg~ (t ) + g~ (t ), 2
3
g~& 3 (t ) = − αg~3 (t ) + 10α 3σ 2g~ ν(t ), g~(t ) = −αζg~ (t ) + g~ (t ), 1
2
with ζ = ( 5−1) 5 , where ν(t ) is the is the zeromean white noise with the unit intensity. 10-2
S
(1 ) g~ r
(ω )
S
g~ s
(ω )
10-4 10-6
(m/s2)2s
The difference between the second integral of the gravimeter readings and the altitude from the DSNS is formed in order to eliminate unknown altitude h (see Fig. 1). The differential measurement can be represented as (Kulakova, et al., 2004):
5 ⋅ ω2 + α 2
S g~ (ω) = 2α 3 ⋅ σ 2g~ ⋅
The gravimeter readings can be represented as a sum of three components: g~ gr = g~ + δg~ + h&& . Here g~ is
h SNS = h + δh , here h is unknown altitude and δh is the errors in DSNS measurements.
g~ˆ
H(s)
The survey data for the investigations was obtained with the use of the gravimeter developed in the CSRI Elektropribor (Blazhnov, et al., 2002; Stepanov, et al, 2002) and dual-frequency geodetic Novatel receivers. It has been shown that the errors in phase measurements are mainly of a white-noise character with the intensity Rh =(0.005m)2s and the gravimeter errors can be described as a white noise with the intensity R gr = (5mGal)2s.
2
gravity anomaly; δg is the gravimeter error and h&& is the vertical acceleration caused by the vertical motion of the aircraft. The data from the DSNS are
zh
2
Fig.1. A block diagram for the GA estimation.
2. STATEMENT OF THE PROBLEM
zh =
-
1
S g~ ( ω )
10-8
S
( 0 ,1 ) g~ r
(ω )
Rh ω4
10-10 10-12 10-14 10-6
10-5
10-4
10-3
.
10-2
10-1
1
Rad/s Fig.2. GA spectral densities and measurement errors for ∇g~ =10 mGal/km, σ g~ =30 mGal, V =50 m/s.
Preprints of the 5th IFAC Symposium on Robust Control Design
ROCOND'06, Toulouse, France, July 5-7, 2006
[
3. ROBUST SOLUTIONS OF THE PROBLEM Under the assumptions made the problem of the GA estimation reduces to estimation of g~ using equations as follows: h&(t ) = V z (t ), V&z (t ) = g~(t ) + R gr n gr (t ),
where E is the mathematical expectation, e(t ) = x(t ) − xˆ (t ) is the estimation error of x(t ) . The Kalman filter provides minimum error variance estimation for every component of the vector x . The GA value itself can be estimated as zˆ(t ) = C1 xˆ (t ) .
z − zˆ
y (t ) = h(t ) + Rh n h (t ),
[
2 2
]T ,
where VZ (t ) is the vertical velocity, n gr (t ) is the
where n = n gr
noise that describes the gravimeter error, y (t ) is the
noise attenuation and
differential measurement, n h (t ) is a measurement noise caused by the errors in DSNS measurements. For the GA model (2) the estimation problem is formulated in this way: x&1 (t ) = x 2 (t ),
nh
< γ2 ⎡ n ⎢⎣
(3)
x& 5 (t ) = −αx 5 (t ) + 10α 3σ g~2 ν (t ), y (t ) = x1 (t ) + Rh n h (t ),
(7)
γ > 0 is a prescribed level of
⋅ 2 denotes the norm in L2 [0, ∞ ) . It is essential that for the H∞-filter the noises in (4) are modeled as signals with a bounded energy and an arbitrary spectrum. Note that the H∞ filter provides optimal estimation in sense of (7) only for the signal z (t ) = C1 x(t ) .
T = [T~z ν
T~z n ] ,
(8)
T~z ν ( s ) = C1 (sI − A1) B2 , T~z n ( s ) = C1 (sI − A1) Bn −1
−1
2
where x(t ) is the state, ν(t ) is the scalar process noise, z (t ) is the signal to be estimated from the measurement y (t ) . The system (3) can be written in the following form: + B2 ν(t ),
(4)
y (t ) = C 2 x (t ) + Dn h (t ) ,
where B 2T = (0,0,0,0, 10α 3 σ 2g~ ) ,
C1 = (0,0,−αξ ,1,0) , C 2 = (1,0,0,0,0), D = Rh .
Both the Kalman and H∞ filters are designed for the GA model (2). These filters are of the form: xˆ& (t ) = Axˆ (t ) + K [y (t ) − C 2 xˆ (t )], zˆ (t ) = C1 xˆ (t ),
2 + ν 2⎤ , ⎥⎦
with A1 = A − KC 2 and Bn = [B1 − KD ] . Then the performance measurer (7) can be rewritten as
z (t ) = −αζ x3 (t ) + x 4 (t ),
B1T = (0, R gr ,0,0,0) ,
2 2
Let T denote the transfer matrix from noises to estimation error ~ z = z − zˆ :
x& 2 (t ) = −αζ x3 (t ) + x 4 (t ) + R gr n gr (t ),
x& (t ) = Ax (t ) + B1 n gr (t ) z (t ) = C1 x (t ),
(6)
The H∞-filter is designed to provide
g~& (t ) → model,
x& 3 (t ) = −αx3 (t ) + x 4 (t ), x& 4 (t ) = −αx 4 (t ) + x5 (t ),
]
De = E e T (t )e(t ) ,
(5)
where xˆ (t ) is the estimate of x(t ) , zˆ (t ) is the estimate of z (t ) and K is the filter gain. For the Kalman filter the noises n gr (t ) , ν(t ) , n h (t ) are assumed to be zero-mean unit intensity white noise signals. The gain K is chosen to achieve the minimum for the following criterion:
T ∞ < γ2 ,
where
T
∞
(9)
denotes the H∞-norm of the transfer
function T . Thus the H∞-filter minimizes the H ∞ norm of the transfer function from noises to the estimation error. Note that the Kalman filter minimizes the H 2 -norm of this transfer function. In the case when the noises in (4) are standard white noise signals the performance measurer (9) has such a meaning that the maximum peak of the error spectral density will not exceed γ 2 . It is important for many applications as the Kalman filter, which minimizes the integral of this spectral density, often leads to a large error spectrum peak that causes bad robustness. It should be noted that choosing the value of γ is a trade-off between the H 2 and H ∞ performances for the H ∞ -filter: near γ 0 , which is a minimal value of γ for which (9) can be met, the error spectrum peak will be reduced to the lowest possible level, at γ → ∞ the H ∞ filter will reduce to Kalman filter. The rms (root-mean-square) estimation errors for the Kalman filter for different values of GA gradient – 3, 5 and 10 mGal/km are presented in Table 1. In calculations it was assumed that σ g~ =30mGal. The aircraft speed was assumed to be equal to V =50m/s, which is typical for an aerogravimetric survey.
Preprints of the 5th IFAC Symposium on Robust Control Design
ROCOND'06, Toulouse, France, July 5-7, 2006
Table 1. Rms estimation errors for GA, mGal
6 . 10
8
3
5
10
5 . 10
8
Kalman filter, σ e
2
3.1
5.5
4 . 10
8
H ∞ -filter for γ ∗ , σ e
3.2
4.4
8.4
3 . 10
8
H 2 / H ∞ filter, σ e
2.9
4.0
6.0
Suboptimal filter, σ es
2.2
3.6
6.9
Suboptimal filter, σ 'es
2.1
3.4
6.4
Cost-guaranteed filter, σ e
2.8
4.1
6.7
Cost-guaranteed filter, σ er
2.8
3.9
6.4
During the investigation of the H ∞ -filter it was found that for γ near γ 0 the H ∞ -filter provides an inadmissible large error (about 25 mGal) when the system (4) is driven by the white noise signals (this means the estimation error that provides the H ∞ filter for the GA described by model (2)). That’s why the value γ was increased in order to achieve a reasonable estimation error. In Fig. 3 the error spectral densities are shown for the H ∞ -filters designed for γ = γ 0 and γ = γ ∗ , where γ ∗ is a chosen level of noise attenuation. Fig. 3 corresponds to the value of the GA gradient 10 mGal/km, in this case γ 0 ≅ 2 ⋅ 10
−8
(m / s ) s 2 2
(
(the error is about 25
)
2
mGal ) and γ * ≅ 2.5 ⋅ 10 −8 m / s 2 s (the error is about 8 mGal). In addition, Fig. 3 shows the error spectral density for the Kalman filter. It is clearly seen in the picture that decrease of the value γ results in decrease of the error spectrum peak on the one hand, but in broadening of the error spectra on the other hand. The estimation errors obtained for the H ∞ -filter for γ ∗ are presented in Table 1. Note that for γ > γ 0 there exist many filters that achieve (9), and a specific filter may be chosen to improve the H 2 performance of the filter. The problem of H 2 performance optimization subject to H ∞ constraint is referred to as a mixed H 2 / H ∞ problem (Khargonekar, et al., 1996). At the moment no analytical solution to this problem is known. In the present paper the mixed H 2 / H ∞ filtering problem will be considered in this way: find a gain K of the filter (5) that solves:
{
}
inf T 2 : subject to T ∞ < γ ,
(10)
where γ is given. The value γ for the H 2 / H ∞ filter synthesis is chosen approximately in the middle between γ 0 and γ K , where γ K is the H ∞ constraint for the Kalman filter (according to Fig. 3
(m/s2)2s
∇g~ , mGal/km
∇g~ σ g~
4
= 10 mGal/km = 30 mGal
3 2
1
8 2 . 10
1 . 10
8
0
0.2 .
0.1
0.3
0.4
Rad/s
Fig.3. Error spectral densities: 1 – the H∞-filter for γ0; 2 - the H∞-filter for γ*; 3 – the H2/H∞ filter; 4 – the Kalman filter. 2.5
∇ g~
4 2
3
2
σ g~
= 10 mGal/km = 30 mGal
1.5
1
5
0.5
0
1
1 3
0
0.1
0.2
. Rad/s
0.3
0.4
Fig.4. The AFC for the vertical acceleration: 1 – the H∞-filter; 2 – the H2/H∞ filter; 3 – the Kalman filter; 4 – the suboptimal filter; 5 – the cost-guaranteed filter.
(
)
2
γ K ≅ 5 ⋅ 10 −8 m / s 2 s ). Fig. 3 shows the error spectral density corresponding to the H 2 / H ∞ filter. The rms errors for the H 2 / H ∞ filter are presented in Table 1. The amplitude-frequency characteristics (AFC) ~ H ( s ) = H ( s ) s 2 for these three filters are shown in Fig. 4. Note that all these filters are designed for the GA model (2) and they have the fifth order. It is clear that the accuracy obtained as well as the structure of the filter depend on the type of the models used for the measurement errors and GA. The models describing the measurement errors were derived from the analysis of the real data. Model (2) is widely used in the problems that require a stochastic description of the GA. At the same time many examples show that there is a strong reason for adjusting filters not to real models of signals to be estimated but to unsteady models represented as integrals of a white noise (Tupysev, 2004). The fact is that it allows improving filter robustness and achieving guaranteed accuracy. One of the possible ways for specifying a model of a signal as integrals
Preprints of the 5th IFAC Symposium on Robust Control Design
ROCOND'06, Toulouse, France, July 5-7, 2006
of a white noise is to rectify (to replace by straight line) the spectral density of the signal in the vicinity of its cross point with error spectral density (Chelpanov, et.al., 1978). The filter which is obtained for the rectified spectral density of the real signal is referred to as a suboptimal filter. The GA spectral density (2) can be approximated, in the logarithmic scale, in the vicinity of its cross point with ω 4 Rh (see Fig. 2) by the second integral of the white noise with the intensity q g2~ = 10α 3 σ 2g~ : S g~s (ω) ≈ q g~ / ω 4 .
(11)
Thus, the suboptimal filter is a Kalman filter that uses (11) as the GA model. It is easy to show that the suboptimal filter transfer function can be represented as H s ( s ) = s 2W B ( s ) , where WB (s ) corresponds to the Butterworth 4th-order filter with the cutoff frequency
(q g~
)
Rh 1 6 .
The
AFC
of
described by the model (2) σ 'es are also presented. The final step will be rejection of the GA description by any spectral density and representation of a priori information about the GA as two parameters: the variances of GA and its derivative: (12)
It is suggested that only these numerical characteristics are known (Kulakova, et al., 2005). In the case that the signal has bounded variances of the signal itself and its derivative the following expression holds for its spectral density S (ω) : 1 2π
∫ [c0 S (ω) + c1ω
∞
2
]
S (ω) dω = c 0 D0 + c1 D1 ,
−∞
where c 0 , c1 are the real, positive, arbitrary parameters. Then S (ω) can be written as: S (ω) = S r(0,1) S w (ω) ,
(13)
with S r(0,1) =
c 2 D0 + D1 c 2 + ω2
and
1 2π
∞
∫ S w (ω)dω = 1 ,
−∞
where c 2 = c 0 c1 . It is true because the spectral S w (ω) can be chosen as density S w (ω) = S (ω) / S r(0,1) .
x& 2 (t ) = x3 (t ) + R gr ngr (t ), x& 3 (t ) = − cx 3 (t ) + c 2 D0 + D1 w(t ),
Thus the signal with the bounded variances (12) can be obtained from the signal with unit variance and unknown spectrum. It is possible to show that the problem of the GA estimation when only numerical characteristics (12) are known can be formulated as follows:
(14)
z (t ) = x3 (t ), y (t ) = x1 (t ) + Rh n h (t ),
where n gr (t ) , n h (t ) are assumed to be standard white noise signals, w(t ) is the zero mean signal with unit variance and unknown spectrum. The system (14) can be written in the form similar to (4). In the case that the spectrum of the process noise w is unknown it is impossible to find error variance. However, since w is of bounded power, it is possible to find an upper bound on the error variance and thus to design a filter that minimizes this upper bound. This means to find a filter that provides:
{
inf T~z w
this
Butterworth filter is shown in Fig. 4. The estimation errors for the suboptimal filter σ es are presented in Table 1. It should be noted that the suboptimal filter is adjusted to model (11) that cannot represent the real dynamics of the GA signal. That’s why the errors that provides suboptimal filter for the signal
D0 = σ 2g~ , D1 = (∇g~ )2 .
x&1 (t ) = x 2 (t ),
2 ∞
+ T~z n
2 2
}.
(15)
The simplest admissible filter in this case has a structure like the one in (5) (Khargonekar, et al., 1996) and it is filter (5) that is used. Thus, the function under optimization J ( K , c ) = T~z w
2 ∞
+ T~z n
2 2
depends on the filter gain
matrix K and the parameter c in the plant model. Their optimal values depend on the specific values D0 , D1 , Rh , R gr and can only be obtained during the function J ( K , c ) minimization. As the upper bound for the error variance is minimized for such filter, it will provide the guaranteed estimation of the signals which are represented by the numerical characteristics (12). That’s why this filter is referred to as a costguaranteed filter. It should be noted that problem (15) can also be regarded as a mixed H 2 / H ∞ problem but formulated in a way different from (10). It is important to understand that problem (15) was derived from other assumptions and for other purposes. In Table 1 the cost-guaranteed filter is characterized by two parameters: the upper bound of the estimation error σ e (the error value which will not be exceeded for the entire class of input signals with the preset values D0 , D1 accepted in designing a filter) and the error σ er that provides the cost-guaranteed filter for the GA described by model (2). It is interesting to note that the parameter D0 had little influence on the accuracy of the GA estimation. This means that if only the parameter D1 is preset the upper bound σ e will vary slightly in comparison with the case when two parameters are preset. In the case that only D1 is used the GA spectral density is
Preprints of the 5th IFAC Symposium on Robust Control Design
represented
as
S (ω) = S r(1) S w (ω) ,
ROCOND'06, Toulouse, France, July 5-7, 2006
where
S r(1) = D1 / ω 2 ( c = 0 ). This means that in this case the cost-guaranteed filter is adjusted to random walk.
Fig. 2 shows the spectral densities S r(1) and S r(0,1) , with the optimal parameter c for S r(0,1) (the value of c for which the upper bound on the error variance is minimal). In order to estimate the weight of the variances (13), special consideration must be given to how close the curves S r(1) and S r(0,1) are in the vicinity of their cross point with the measurement errors. It is clear that the curve
S r(1) represents the
rectified curve S r(0,1) in the logarithmic scale at its cross point with the spectral density of the measurement errors. Thus, the cost-guaranteed filter is mainly determined by the variance of the GA derivative. Fig. 2 also shows that the cost-guaranteed filter is adjusted to the spectral density that has a wider frequency bound. Therefore, it will be lowsensitive to the change of the acceptable model that causes extension of the spectral density. Therewith, the cost-guaranteed filter will provide a higher error in comparison with the filter adjusted to model (2) when the GA spectral density is narrowed. The cost-guaranteed filter AFC is shown in Fig. 4 from which it follows that the filter turns out to be similar to the classical low-pass filter. Its AFC has flat vertex and an inclination of 60 dB/dec in the high-frequency domain. The existence of peaks for the AFC of the other filters makes them too sensitive to an increased noise power at the frequencies where the AFC is high. It is also sensitive to the change of the accepted model that causes decrease of the valid signal energy in this frequency range. Therefore, from the intuitive point of view, the cost-guaranteed estimator is low-sensitive. Although no special measures were taken with this aim. CONCLUSIONS Five different filtering approaches to GA estimation aboard an aircraft have been analyzed and compared with respect to the appropriate a priori information and their numerical performance. Several tests have been carried in order to investigate their efficiency for various operating conditions. The results obtained suggest that the cost-guaranteed filter will provide reliable and efficient estimation of a GA. The costguaranteed filter makes use of two parameters (the GA variance and the GA gradient) which can be easily obtained from the available experimental data; it has the features similar to those of the mixed H 2 / H ∞ filter; it is adjusted to the model represented as the integral of the white noise; the cost-guaranteed filter loses to the optimal filter in accuracy to a small extent. In conclusion, it should be noted that the paper presents a novel approach to the cost-guaranteed filter designing, the one when only numerical characteristics are given for signals.
REFERENCES Afanasjeva, L.V., E.A. Bojarskij, N.V. Drobyshev and V.N. Konushev (2002). About the errors of gravimetry anomaly images obtained from the aircraft readings. Geophizicheskij vestnik. No 4, pp. 13-15 (in Russian). Blazhnov, B.A., L.P. Nesenjuk, V.G. Peshekhonov, A.V. Sokolov, L.S. Elinson, L.K. Zhelesnyak (2002). An integrated mobile gravimetric system. Development and test results. Proceedings of the 9th Saint-Peterburg International Conference on Integrated Navigation Systems. Chelpanov, I.B., L.P. Nesenjuk and М.V. Braginsky (1978). Computation of navigational gyro devices characteristics. Sudostroyeniye, Leningrad (in Russian). Grimble, M.J. and A. Elsayed (1990). Solution of the H∞ optimal linear filtering problem for discretetime systems. IEEE Trans. Acoustics, Speech, Signal Processing, Vol. 38, pp.1092-1104. Jordan, S.K. (1972). Self-consistent statistical models for gravity anomaly and undulation of the geoid. Geophys. Res. Vol. 77, No 20, pp. 3660-3670. Khargonekar, P.P., M.A. Rotea and E. Baeyens (1996). Mixed H2/H∞ filtering. Internat. J. Robust Nonlinear Control, Vol. 6, pp. 313–330. Kulakova, V.I., A.V. Nebylov and O.A. Stepanov (2005). Robust estimation versus optimal estimation and their application for airborne gravimetry. Proceedings of the 16th IFAC World Congress. Prague. Kulakova, V.I., A.V. Nebylov and O.A. Stepanov (2004). Application of the robust approach to the problem of airborne gravimetry. Proceedings of the 16th IFAC Symposium on Automatic Control in Aerospace. St. Petersburg, Russia. Pp.354-359. Nagpal, K.M. and P.P. Khargonekar (1991). Filtering and Smoothing in an H∞ Setting. IEEE Trans. Automatic Control, Vol. 36, pp.152-166. Nebylov, А.V. (2004). Ensuring control accuracy. Lecture Notes in Control and Information Sciences, Springer-Verlag, Heidelberg, Germany. Rotstein, H. and M. Sznaier (1996). H2/H∞ filtering theory and an aerospace application. Internat. J. Robust Nonlinear Control, Vol. 6, pp. 347–366. Stepanov, O.A., B.A. Blazhnov and D.A. Koshaev (2002). The efficiency of using velocity and coordinate satellite measurements in determining gravity aboard an aircraft. Proceedings of the 9-th Saint Petersburg International Conference on Integrated Navigation Systems. St. Petersburg, Russia. Pp.255-264. Tupysev, V.A. (2004). Guaranteed estimation of the dynamic system state in the case of uncertainty in the description of disturbance and measurement errors. Proceedings of the 11th Saint Petersburg International Conference on Integrated Navigation Systems. St. Petersburg, Russia. Pp.86-89. Van Trees, H. (1968). Detection, estimation and modulation theory. Part I. Detection, estimation and linear modulation theory. New-York.