Complementary filter for wind and position

Background

Wind measurements combine a measurement of relative wind with a measurement of aircraft motion to determine the air motion relative to the ground. The aircraft motion has long been measured by an IRS, and recently also by a GPS. These have complementary strengths: The IRS provides very good information on short-term motion but drifts with a characteristic period of more than an hour, while the GPS provides good absolute accuracy but sometimes is unable to receive the GPS signals and (except in differential-GPS mode) can have short-term errors that make short segments of the track look jagged. To take advantage of the strengths of each, a complementary-filter calculation was developed and implemented in the 1980s, but it was never documented publicly. This memo is partly an attempt to remedy that and partly a suggestion to make some minor changes to how it is implemented. This note is complementary to the information in the document on RAF Processing Algorithms, section 3.4, and contains additional detail as well as notes regarding implementation of changes.

Aircraft velocity

To accomplish this, a low-pass filter,
$F_{L} ({G V N S, G V E W})$

, is applied to the GPS measurements of groundspeed, {GVNS,GVEW}, which are assumed to be valid for frequencies at or lower than the cutoff frequency
$f_{c}$

of the filter. Then the complementary high-pass filter, denoted (
$1 - F_{L}$

)(
${V N S, V E W}$

), is applied to the IRS measurements of groundspeed, {VNS,VEW}, which are assumed valid for frequencies at or higher than
$f_{c}$

. Ideally, the transition frequency would be selected where the GPS errors (increasing with frequency) are equal to the IRS errors (decreasing with frequency). The filter used is a three-pole Butterworth lowpass filter, coded following the algorithm described in Bosic, S. M., 1980: Digital and Kalman filtering : An Introduction to Discrete-Time Filtering and Optimum Linear Estimation, p. 49. The digital filter used is recursive, not centered, to permit calculation during a single pass through the data. If the cutoff frequency lies where both the GPS and INS measurements are valid and are almost the same, then the detailed characteristics of the filter in the transition region (e.g., phase shift) do not matter because the complementary filters have cancelling effects when applied to the same signal. The transition frequency
$f_{c}$

was chosen to be (1/600) Hz (but this value can be overridden via the “defaults” file). The Butterworth filter was chosen because it provides flat response away from the transition. The resulting variables for aircraft motion, {VNSC,VEWC}, are then each the sum of two filtered signals, calculated as described in the following box:

VEW = IRS-measured east component of the aircraft ground speed
VNS = IRS-measured north component of the aircraft ground speed
GVEW = GPS-measured east component of the aircraft ground speed
GVNS = GPS-measured north component of the aircraft ground speed

F_{L} ()

= three-pole Butterworth lowpass recursive digital filter

{V N S C} = {V N S} + F_{L} ({G V N S} - {V N S})

{V E W C} = {V E W} + F_{L} ({G V E W} - {V E W})

This is straightforward and effective when both sets of measurements (IRS and GPS) are available. The approach in use becomes more complicated when the GPS signals are lost, as sometimes happens in sharp turns. Then some means is needed to avoid sudden discontinuities in velocity (and hence windspeed), which would introduce spurious effects into variance spectra and other properties dependent on a continuously valid measurement of wind. To extrapolate measurements through periods when the GPS measurements are not available, a fit is determined to the difference between the best-estimate variables {VNSC,VEWC} and the IRS variables {VNS,VEW} for the period before GPS reception was lost, and that fit is used to extrapolate through periods when GPS reception is not available. The procedure is described in section 3.4 of ProcessingAlgorithms.pdf.

The changes proposed in that section are as follows:

Change FCTRF to 0.995 from the present value of 0.997. Because in this case the exponential filter is applied constantly to the correction factors like dvy, it is best to choose the time constant for the transition filter to have small impact on any frequencies below the cutoff frequency of the low-pass filter, which in the standard implementation is 600 s. The constant used in standard processing, 0.997, has about 17% residual at 600 s (or 1.8 time constants), so it is better to choose a smaller value. For 0.995, the residual is 5%, so this would be a better choice.
Correct the filter constant
$b_{1}$
. The factor
$\sqrt{\frac{3}{2}}$
was erroneously
$\frac{\sqrt{3}}{2}$
. See the Processing Algorithms discussion.
To be consistent with standard values of constants used elsewhere in processing, the Schuler oscillation period should be adjusted from 5040 s to 5067 s.
The present code (gpsc.c; see /scr/raf/cjw/builds/tikal/raf/nimbus/src/amlib/std/gpsc.c) has, at line 345, decay factors to use when GPS is lost and there is no accumulated Schuler-oscillation fit to use. The code adjusts the correction factors coming from the complementary filter toward the IRU solution. There does not seem to be any benefit from doing this, instead of just continuing to apply the offset determined from the complementary filter. Implementation could be by commenting lines 345–348.
At present, the code makes a smooth transition to the Schuler-oscillation fit when GPS measurements become invalid to avoid introducing discontinuity into the new variables. However, when the GPS reception again becomes valid, a discontinuous transition is made back to the complementary-filter solution. This should be avoided by a smooth transition. Code similar to that used to update dvx (for example) for the transition to the Schuler-oscillation solution can be used for this transition. However, to avoid introducing a new branch in the code and applying another exponential filter for this transition, it seems preferable to apply this filter to all good measurements but choose a filter constant that does not have significant influence at the low-pass-filter cutoff frequency. A suggested way to implement this would be to replace lines 359–367 by these:

Note that it is not necessary to filter gvns and vns separately because only the difference between the filtered results is used.

The following provides more documentation of the fit procedure used to determine the Schuler oscillation. The errors are assumed to result primarily from this oscillation, so the three-term fit is of the form
$Δ = c_{1} + c_{2} sin (Ω_{S c h} t) + c_{3} cos (Ω_{S c h} t)$

, where
$Ω_{S c h}$

is the angular frequency of the Schuler oscillation (taken to be
$2 π / (5067 s))$

and
$t$

is the time since the start of the flight. A separate fit is used for each component of the velocity and each component of the position (discussed below under LATC and LONC). The fit matrix used to determine these coefficients is updated each time step but the accumulated fit factors decay exponentially with about 30-min decay constant, so the terms used to determine the fit are exponentially weighted over the period of valid data with a time constant that decays exponentially into the past with a characteristic time of 30 min. This is long enough to determine a significant portion of the Schuler oscillation but short enough to emphasize recent measurements of the correction. The procedures for accumulating the matrices for the fit are as follows:

Define

u_{G}

as the aircraft eastward velocity measured by the GPS and
$u_{I}$

the corresponding velocity measured by the IRS, so that the difference is

δ u = u_{G} - u_{I}

If
$Ω_{S}$

is the Schuler oscillation period, with
$Ω_{S} = 2 π t / T_{s}$

where
$T_{s} = 5067$

s,
$τ_{u}$

is the time constant for the update (1800 s),
$t$

is the time from the start of the flight, and the measurement matrix is
$A_{i j}$

, then updated terms of the measurement matrix each sample period (
$A_{i, j}^{'}$

) are (for
$δ u$

A_{0, 1}^{'} = A_{0.1} (1 - \frac{1}{τ_{u}}) + δ u

A_{1, 1}^{'} = A_{1, 1} (1 - \frac{1}{τ_{u}}) + δ u sin (Ω_{S} t)

A_{2, 1}^{'} = A_{2, 1} (1 - \frac{1}{τ_{u}}) + δ u cos (Ω_{S} t)

The matrix components

A_{j, 0}

apply to the northward velocity component and so are represented by the same equations with
$δ u$

replaced by
$δ v$

. Similar matrices are calculated for latitude
$θ$

and longitude
$φ$

, based on the differences
$δ θ$

and
$δ φ$

between GPS and IRS measurements. The information matrix
$H_{i j}$

is calculated via

H_{i, j} = H_{i, j} (1 - \frac{1}{τ_{u}}) + V_{i, j}

where

V_{0, 0} = 1

,
$V_{0, 1} = V_{1, 0} = sin (Ω_{S} t)$

,
$V_{0, 2} = V_{2, 0} = cos (Ω_{S} t)$

,
$V_{1, 1} = {sin}^{2} (Ω_{S} t)$

,
$V_{1, 2} = V_{2, 1} = sin (Ω_{S} t) cos (Ω_{S} t)$

, and
$V_{2, 2} = {cos}^{2} (Ω_{S} t)$

. When the fit is needed, the matrix
$H_{i, j}$

is inverted and the result multiplied by the measurement matrix
$A_{i, j}$

to get the fit coefficients
$C_{i j}$

to use for predicting the results for
$δ u$

,
$δ v$

,
$δ θ$

, and
$δ φ$

via equations like
$δ u = C_{0, 1} + C_{1, 1} sin (Ω_{S} t) + C_{2, 1} cos (Ω t)$

Aircraft position

A simpler approach using an exponential filter is possible with aircraft position, as follows:

LAT = latitude measured by the IRS
LON = longitude measured by the IRS
GLAT = latitude measured by the GPS
GLON = longitude measured by the GPS
VNSC = aircraft ground speed, north component, corrected
VEWC = aircraft ground speed, east component, corrected

Initialize the corrected position at the IRS position at the start of the flight or after any large change (>5
$^{○}$
) in the IRS position.
Integrate forward from that position using the aircraft groundspeed with components {VNSC,VEWC}. Note that in the absence of GPS information this will introduce long-term errors because it does not account for the Earth's spherical geometry. It provides good short-term accuracy, but the GPS updating in the next step is needed to compensate for the difference between a rectilinear frame and the Earth's spherical coordinate frame and provides a smooth yet accurate track.
Use an exponential adjustment to the GPS position, with time constant that is typically about 100 s.

1

specifically, LATC +=
$η$
(GLAT-LATC) with
$η = 2 π / (600 s)$
To handle periods when the GPS becomes invalid, use an approach analogous to that for groundspeed, whereby a Schuler-oscillation fit to the difference between the GPS and IRS measurements is developed and used to extrapolate through periods when the GPS is invalid.

Earth Observing Laboratory

Background

Aircraft velocity

Aircraft position

NCAR

UCAR