Extended Gauß-Markov model

Based on the above introduced principles of parameter estimation within a GMM by means of one observation type, an extended GMM now is set up in order to derive regional gravity models from the combination of heterogeneous data. The different observation types are distinguished and collected in groups, following the criteria of

• measurement system (satellite, absolute/relative terrestrial, ship-/airborne gravimetry,. . .),

• observed functional (∆V,V_ab,T,∆g, δg, . . .),

• observation height (orbit height of satellite, measurements at the Earth’s surface, at a platform,. . .),

• spectral resolution (depending on functional, measurement system,. . .),

• accuracy (e. g. externally influenced by the climate, humans,. . .),

• observation epoch (hours, months, consideration of seasonal variations,. . .),

• sampling rate (of one or several measurement campaigns,. . .),

• pre-processing (application of reductions, corrections, elimination of outliers,. . .),

• reference (height, gravity,. . .),

• spatial distribution,

• . . .

For the combination of the different groups within an appropriate estimation model, especially the first three aspects are challenging, i. e. handling a variety of functionals and observation heights from diverse measurement systems. For the different functionalsY[V] ∈ {∆V,Vab,T,∆g, δg, . . .}, different observation equations (5.15) have to be formulated. The design matrix in Eq. (5.16) then contains the correspondingly adapted basis functions ˜b(x^p,xq)from Tab. 4.7 withx= x^p. The varying observation heights are taken into account by the downward continuation term(^R/r)^l+1, cf. Tab. 2.2, within the observation equations. Bouman et al.(2013) study signal and degree RMS at different observation heights;Kern(2003), pp. 105, studies the effect of attenuation and amplification of gravitational functionals and their errors at different heights. The remaining criteria which classify the diverse observation groups, depend on the specific data sets, cf. Sec. 3.3.

Examples are given in the next chapter. In the following, components of the single-level combination are studied, referring to the two boxes of analysis and synthesis in Fig. 5.9.

Analysis: Estimation of coefficients

Within the analysis step (light blue box in Fig 5.9), the unknown scaling coefficients are estimated. They describe the amplitudes of basis functions in the spatial domain. In Fig. 5.1, it would mean to locate Shannon SBFs at the Reuter grid points (red crosses) in∂ΩC, and to multiply them with the estimated coefficients.

Choice of basis functions

The basis functions ˜bqare defined as Shannon functions, cf. Eq. (4.6), i. e. insertingBl =1 in the expressions of Tab. 4.7. For the simplest functional,T, the basis function yields

b˜L(x^p,x_q)= XL

l=0

2l+1 4πR²

R r

!l+1

Pl

(r^p)^Tr_q

, (5.21)

based on Legendre polynomialsP_l w.r.t. the spherical distance angle between observed point P^p = P(x^p) and computation point Pq = P(x_q), i. e. r^p andr_q as introduced in Eq. (4.2). The bandlimiting maximum degreeLcorresponds to the chosen maximum degree of the modeling resolution, cf. Eq. (5.12). As discussed before, the Shannon kernel ensures an exact band limitation at degree Lwithout any smoothing behavior.

Observation equations

For the combination of different observation techniques, the parameter estimation model (5.16) is extended from one observation group toKgroups with observation vectors y= y_k,k =1, . . . ,K. Next, an appropriate background model is selected and subtracted according to Eq. (5.13) from the elements of all observation vectors y_k. It yields the differential ¯n_k ×1 observation vector y_k of the k^th measurement technique with differential elements y(x^p) = y_k(x^p) according to Eq. (5.14). e_k and A_k are the corresponding ¯nk ×1 error vector and the ¯nk ×QL design matrix. For the combination of the altogetherKobservation techniques including the additional model (5.19) for the prior information, theextended GMMreads







 y₁ y₂ ... y_K µ_d







 +







 e₁ e2

... eK

e_d









=







 A₁ A2

... AK

I









d_L with D

*...

....

,







 y₁ y₂ ... y_K µ_d







 +///

////

-=









σ²₁P⁻¹₁ 0 . . . 0 0

0 σ²₂P⁻¹₂ . . . 0 0

... ... . . . ... ...

0 0 . . . σ²_KP⁻_K¹ 0

0 0 . . . 0 σ²_dP⁻_d¹









(5.22)

(Schmidt et al., 2007). It is displayed within the red box in Fig. 5.9. The design matrices A_k contain the adapted basis functions from Eq. (5.21). Each covariance block matrix D(y_k) = σ²_kP⁻_k¹ consists of an unknown Variance Component (VC)σ_k² and the inverse of a given ¯nk ×n¯k positive definite weight matrix P_k. Those block matrices, and the covariance matrixD(µ_d)of the prior information, are located on the main diagonal of the covariance matrix D([y^T₁, . . . ,y^T_K,µ^T_d]), cf. Eq. (5.22). The off-diagonal block elements are zero-matrices, assuming that the observation groupsy_k and the prior informationµ_dare uncorrelated.

Further, as mentioned in the context of Eq. (5.16), in all numerical investigations in Chapter 6, the measurement errors are assumed to be uncorrelated and the observations are assumed to have similar accuracies within one observation groupk. Thus, identity matrices can be introduced for the weight matricesP_k in Eq. (5.22).

Least-squares adjustment

Following Lieb et al. (2016), the K observation techniques with each ¯nk (k = 1, . . . ,K) observations are combined at normal equation level. A combination at observation equation level instantly delivered residuals of the observations, which are necessary for VCE. However, the size and, thus, the computation time of the linear equation systems restrict a practical realization within the scope of this work. According to the normal equation system (5.17) for one observation group, the observation equation system (5.22) forK groups and the prior information yields

* ,

XK

k=1

* ,

1

σ²_kA^T_kP_kA_k+

-+ 1

σ_d²P_d+ -Dd_L =

XK

k=1

* ,

1

σ²_kA^T_kP_ky_k+

-+ 1

σ²_dP_dµ_d . (5.23)

This system of normal equations is independent of the number of observations. The estimated coefficientsdDq, collected in the vectorDd_L up to degreeL, are computed from

Dd_L =* ,

XK

k=1

* ,

1

Dσ²_kA^T_kP_kA_k+

-+ 1

D σ²_dP_d+

-−1 * ,

XK

k=1

* ,

1 D

σ²_kA^T_kP_ky_k+

-+ 1

D

σ²_dP_dµ_d+

- , (5.24)

cf. Eq. (5.18). The componentsdDq have the SI-unit m²/s²of the gravitational potential; the elements of the design matrixAcontain the corresponding conversion factors for the particular functionals.

Variance component estimation

In this work, VCE is used for regulating (1) the relative weighting of heterogeneous observation groups and (2) the amount of regularization by prior information. The variance components are iteratively estimated toDσ²_k andDσ_d², as suggested byKoch and Kusche(2002). The regularization parameter yieldsλσ =Dσ_k²/Dσ²_d. Hereby, the estimates are computed based on stochastic information, i. e. using the residuals of the observations. The iteration starts from an initial value and ends in the point of convergence (the first three digits of two successive

iterations must not change any more).

The reciprocals Dσ_k⁻² and Dσ⁻²_d of the estimated VCs in Eq. (5.24) determine the relative weight of each observation group y_k, and hence, regulate the contribution of each measurement technique to the overall combination result. The influence depends on the accuracy of the observation data, the spectral resolution and the spatial distribution (Lieb et al., 2016). Since low-resolution signal is smoother in the spatial domain, i. e. it shows large scale-variations, and high-resolution signal represents short-scale variations. Further, the spatial distribution and the spectral content of the observations are related to each other, cf. Eq. (4.21). Data which cover wide areas enable to recover medium down to long wavelengths. The spectral content of measurements also depends on the type of the gravitational functional and the aspect of downward continuation. According to the Meissl scheme (Fig. 2.7), different functionals have different spectral sensitivities to high-frequency variations and the downward continuation, e. g. of satellite measurements down to the Earth’s surface, amplifies high frequencies. All those aspects are connected with the need of regularization by prior information µ_d. Details are studied and discussed by means of numerical applications in Chapter 6.

Covariance information of the unknowns

The full covariance information of the vector DdL of estimated coefficients is obtained by applying error propagation on Eq. (5.24).²² The covariance matrix

DDd_L

=Q⁻¹_dd =







v(dD1) c(dD1,dD2) . . . c(dD1,dDQL) c(dD₂,dD₁) v(dD₂) ...

... . . .

c(dDQL,dD1) . . . v(dDQL)







−1

, (5.25)

contains information about the uncertainty of the estimated coefficients (Lieb et al., 2016). The diagonal elements are the variances v(dD_q) of the estimated coefficients. The off-diagonal elements represent the covariances between two coefficients, for instance,c(dD1,dD2)betweendD1anddD2. Q_ddis theQL ×QL normal equation matrix

Q_dd = XK

k=1

* ,

1

σ_k²A^T_kP_kA_k+

-+ 1

σ²_dP_d (5.26)

expressed on the left side of Eq. (5.23), and inverted in Eq. (5.24).

Synthesis: Computation of gravity field functionals

In order to express different functionals of the Earth’s gravity field, the observation equations (5.15) serve as modeling equations in the synthesis. The elements of the vectorDd_Lof estimated coefficients from Eq. (5.24) are again multiplied with adapted scaling functions from Tab. 4.7. However, the functions ˜bl^′(x^c,xq) = Φ˜q,J(x^c,x_q) now are related to a maximum resolution level J, as introduced Sec. 4.3.3, and the function values are computed at locations x = x^c,c = 1, . . . ,C, inside a target area ∂ΩI. Therefore, the SBFs are located, for instance, on a regular grid or along flight tracks in any height on or above the Earth’s surface.

For the example in Fig. 5.1, Blackman SBFs are located on a regular geographical grid in the study area in Northern Germany, multiplied with the estimated coefficients from the analysis, and visualized up toJ =11 in terms of ∆g (lower right). The functional relations in Tab. 4.7 can be further complemented by other functionals, e. g. quasi-geoid heightsζaccording to Eq. (2.63) or deflections of the verticalξ, η according to (2.70).

22In analogy, error propagation can equivalently be applied to the scaling coefficients Eq. (5.18), estimated from one observation group. Since the combination of several observation groups is in the focus of this work, its description was neglected in Sec. 5.2.1.

Choice of basis functions

According to Eq. (4.33), the adapted scaling function yields Φ˜q,J =ΦJ(x^c,x_q) =

lJ

X

l=0

2l+1 4πR²

R r

!l+1

φ_l,J Pl

(r^c)^Tr_q

. (5.27)

Since the components ofDd_Lare estimated in the analysis step by strictly bandlimiting non-smoothing Shannon kernels, it is possible to multiply them in the synthesis with a different kernel function (Schreiner, 1996), cf.

Eqs. (4.11). The use of different functions, defined by φl,j ∈ (

φ^Sha_l,j, φ^Bla_l,j, φ^CuP_l,j )

according to Eqs. (4.29) – (4.31), is analyzed in Sec. 6.1.1.

The bandlimiting maximum degreelJ ≤ Lis chosen smaller than the maximum degree Lof the analysis in order to reduce aliasing errors in the frequency domain, and thus, erroneous side lobes in the spatial domain, cf. Fig. 5.8. Figure. 5.7 visualizes the principle: Modeling a bandlimited signalg(purple curve) according to Eq. (4.15) up to maximum degreeLwithin the analysis (purple box) from a finite set of function values, yields aliasing errors due to frequency folding (yellow) of the originally infinite signal f (Jekeli, 1996). Expanding the series of ˜Φq,J in the synthesis (gray box) up to l_J ≤ L, acting as low-pass filter, here depending on a maximum resolution level J, reduces those aliasing errors (Schmidt et al., 2007).

Modeling equations

In analogy to Eq. (5.20), the system

z_J = B_JDd_L →∆Z_J, (5.28)

of modeling equations can be set up by means ofDd_L, estimated fromK observation groups according to the extended GMM (5.22), andB_J, containing the chosen low-pass filtering SBFs ˜Φq,J, Eq. (5.27). The output functional valueszJ(x^c), collected in theC×1 vectorzJ, contain information up to maximum spectral degree l_J and describe the ”differential signal“∆Z_J w.r.t. the background model. The system of modeling equations is highlighted in red in the synthesis box of Fig. 5.9.

Restoring the previously subtracted background model at the modeling positionsx^cdelivers the ”total signal“

ZJ = Zback+∆ZJ (5.29)

withZ_back=Y[ ˜V_back]. In analogy to computing the differences at locationsx^p, cf. Eq. (5.13), the functional valuesZbackat positionsx^care obtained from series expansion (2.40) in terms of SHs. Note, in this work, ”total signal“ means that the (low up to medium) frequency information of the previously subtracted background model has been restored, in contrast to the ”differential signal“, where the referring spectral content is missing.

Covariance information of the output functionals

Information about the accuracy and the dependencies of the resulting functionals zJ (x^c) from Eq. (5.28) is obtained by computing the variances and covariances at the output positions P(x^c). Applying error propagation on the modeling equation (5.28)²³, theC×Ccovariance matrixD(zJ)reads

D(z_J) =B_JDDd_L B^T_J

=







 v

zJ

x¹ c zJ

x¹ ,zJ

x² . . . c zJ

x¹ ,zJ

x^C c

z_J x²

,z_J

x¹ v

z_J

x² ...

... . . .

c zJ

x^C ,zJ

x¹ . . . v

zJ

x^C









. (5.30)

It contains on the diagonal the variancesv(z_J (x^c)), now indicating the accuracy of the resulting functional zJ(x^c)up to spectral degreelJat pointP^c, cf. Eq. (5.25). The off-diagonal elements represent the covariances

23In analogy, error propagation can be applied to the modeling equation (5.20).

between two values, for instance,c(z_J(x¹),z_J(x²))at the output grid pointsP¹andP². They give a measure for the dependency, i. e. the correlation, between the values of both points.

֒→In the following, all investigations are transferred from single-level to multi-level combination via MRR.

Im Dokument Enhanced regional gravity field modeling from the combination of real data via MRR (Seite 109-113)