AD-AI0A 577 PACIFIC-SIERRA RESEARCH CORP SANTA MONICA CA F/6 4/1

CALCULATION OF IONOSPHERIC CONDUCTIVITY PROFILES BY INVERTING V--ETCIU)

OCT Al R E WARREN, E C FIELD, C R WARBER F1962R 7R C-008A

NL '

,IEIIIIIIIIII.
UNCLASSIFIED PSR-114 RADC-TR- 81286

IIEEEEEEIIIIIE
III IIIffllflllllllllf
 4 11 J 111...

1.1
1111IL25 ±6

.25
MICROCOPY .4 I?CHART

" *il MI('ROL01'Y fRI SOLUflON If S1 ( I-tARI

I
 CALQDLT ION Of IONOmERic,
COUCTIVITY
:. PROFIES BY
to INVERTING VLF/LF REFLCTION
LODATA RSOROCWPAu
o PtdA.lem nw& cbCp

. Wounni
EL
L C. Field
Q. L Webm 1

APWovm XW uuwc MMEASt. wsYVMWDN UNUMIWD

ROME1- AIR MINOMT~ COO=

Systems Commeuid
AlrfFgr5
Al 41

L11 11
 This report has -been rev'Ied by the RADC Public Affairs Office ( md "
is releasable to the National Tecbical Infomation Service (NTIS). At,8 ,
it will be releasable to the general public, Including foreign nations.

RADC-TR-81-286 has been reviewed and is approved for publication.,

APPROVED:
TERRENCE J. ELKINS
Chief, Propagation Branch
Electromagnetic Sciences Division

APPROVED:

ALLAN C. SC4LL
Chief, Electromagnetic Sciences Division

FOR THE MADR

JOHN P. HUSS
Acting Chief, Plans Office

If your address has changed ot if you wish to be rmwed from the VADC
sailing list, or if the addressee Is no -loagr aloyed by ys erosiatim,
please notify RADC (EEPL) Hmscom APB VA 01731. ..This wiLl asit us la
maintaining a current miling list.

Do not return copies of this report unless contractual obligltim or notices

on a specific document requires that it be returned.
 UNCLASSIFIED
SiCUNITY CLASSIPICATION or TmtS PAGE (owen Dea. Exeorod)
REPORT DOCUMENTATION PAGE READ INSTRUCTIONS
BEFORE COMPLETING FORM
I. RE[PORtT NUMSq n2. QOVT ACCtSSION NO. S. RECIPIENT'S CATALOG NUMBER

RADC-TR-81-286
4. TITLE (dad Subtile) S. TYPE OF REPORT & PERIGO COVEREO

CALCULATION OF IONOSPHERIC CONDUCTIVITY Final Technical Report

VLF/LF REFLECTION
INVERTING BYRINRERTING
PROFILES BY PROFILES VLFRLPORTFLECTEO Oct 79 - May 81
DATA: I. ISOTROPIC PROPAGATION .PR
PSR Repor
1.
111
REPORT trut11

7. AuTNORW(a S. CONTRACT OR GRANT MU9El 'S)

R. E. Warren
E. C. Field F19628-78-C-0088
C. R. Warber
9. PERFORMING ORGANIZATION NAME AND AOORESS 10. PROGRAM ELEMNT. PROJECT. TASK

AREA II WORK UNIT NUMBERS

Pacific-Sierra Research Corporation
1456 Cloverfield Blvd 61102F
Santa Monica CA 90404 2304J322
It. CONTROLLING OPICE NAMC ANO ADORESS 12. REPORT OAT%

Deputy for Electronic Technology (RADC/EEP) October 1981

13. NUMBER Of PAGES
Hanscom AFB MA 01731
14. MONITORING AGENCY NAME A AOORIESS(It diffeen from Co.ngr.Iln Office) IS. SECURITY CLASS. (o ti report)

Same UNCLASSIFIED
ISa. OECLASS PI CATIOW/ O'WNGRAOING
N/AscEOU
1. OISTRIBUTION STATEMENT (of thle Repot)

Approved for public release; distribution unlimited.

17. OSTRIBUTION STATEMENT (ofthe obettrt a.ered in Block 20, It different

frm Report)

Same

r ~It. SUPPLEMENtTARY NOTES

RADC Project Engineer: Paul A. Kossey (EEP)

It. KEY WORDS (ConghWO On,fwaO eide Itne0Ce.WP B

ed ift.(ilf
by block ne"Obe)
Lower Ionosphere
VLF/LF reflection
VLF/LF inverse problem
lonosounding
20. AOSTRItACT (Con~tinue en Poereree &#defIt
necoeev and IdenIF~ br block nutoher)
This report presents an iterative method for inverting VLF/LF
ionosounding data to obtain ionospheric conductivity profiles. The
relationship between the detail that can be realized in the calculated
profile and the quality and quantity of the reflection data is given
explicitly. Only altitudes below about 70 km where the propagation can
be assumed isotropic are considered, although the method can be extended
to anisotropic propagation. The method is demonstrated by applying it

DD | 1473 EDITION OF 1 0ov 65 IS OBSOLETE UNCLASSIFIED

SECURITY CLASSIFICATION OF TMllS PAGE (When De nteredI4)
 UNCLASSIFIED
SeCUMIVI CLA5MICATIO14, F THIS PAGS(Ubat f e".
ol d)
..

to two qualitatively different ionospheric models--a strong SPE

disturbance and a C-layer exhibiting a well-defined conductivity peak.
Ground-level reflection data contain information only about those heights
from which significant energy is returned. The calculated profiles
agree well with the true profiles over the height range where the
strongest reflections occur, and break down outside of that range. The
solutions converge toward the true profiles despite the use of initial
estimates that contain no information about the true ionospheric
structure. The height resolution actually achieved for the two
examples is much better than the theoretically predicted limit.

UNCLASSIFIED
ISCUAITY CLASSIPI1CAYO OP , PAGS&Vflf n Does gnhtedl

- - ---.
 SUMMARY

This report presents an iterative method for inverting VLF/LF

ionosound data to obtain ionospheric conductivity profiles. The re-
lationship between the detail that can be realized in the calculated
profile and the quality and quantity of the reflection data is given
explicitly by our approach. We consider only altitudes below about
70 km where the propagation can be assumed isotropic, although the
method can be extended to anisotropic propagation.
We demonstrate the method by applying it to two qualitatively
different ionospheric models--a strong SPE disturbance and a C-layer
exhibiting a well-defined conductivity peak. Ground-level reflection
data contain information about only those heights from which signifi-
cant energy is returned. Our calculated profiles agree well with
the true profiles over the height range where the strongest ref lec-
tions occur, and break down outside of this range. The solutions
converge toward the true profiles despite our use of initial esti-
mates that contain no information about the true ionospheric struc-
ture. The height resolution actually achieved for the two examples
is much better than the theoretically predicted limit.
 PREFACE

This report describes a method for inverting VLF/LF reflection

data to obtain ionospheric conductivity profiles undet conditions
where the propagation can be assumed isotropic. We demonstrate the
method by inverting computer-generated reflection coefficients for
two generic model ionospheres. In addition, we present a means of
choosing the best compromise between the deleterious effects of (1)
noisy data and (2) loss of uniqueness due to incomplete data.
 -vii-

CONTENTS

SUMMARY .. . . . . . . . . . . . . . . . . . . . . . . . . iii

PREFACE...................................................... V

FIGURES......................................................... ix

TABLES.......................................................... xi

Section
1. INTRODUCTION ............................................ 1

II. INVERSE THEORY APPLIED TO TONOSOUNDING DATA ..... 3

Fr~chet kernels for isotropic ionospheric
conductivity.......................................3
Construction of the inverse solution...................7
Uniqueness of the inverse solution.....................9
III. NUMERICAL EXAMPLES......................................14
Strong SPE ............................................ 14
Inverse solution....................................16
Sensitivity to noise and limited data................16
Physical interpretation..............................23
Spacing of frequencies...............................27
Height-gain function................................28
Model C-layer.........................................29

IV. CONCLUSIONS.............................................35

REFERENCES....................................................37

APPENDIX: THE SPECTRAL EXPANSION METHOD ...................... 39

Fr~~L ~ ~~~m
jjua. j1JAjED
F
 -ix-

FIGURES

1. Conductivity height profile for a strong SPE .......... 15

2. Calculated conductivity profiles (solid curves) aiud

true profile (dashed curve) for the strong SPE
example ............................................. 17

3. Calculated and true profiles after 15 iterations for

the strong SPE example .............................. 18

Magnitude of the averaging function A for the strong

4.
SPE example ......................................... 21 )
5. Tradeoff between height resolution and noise-induced
variance for the strong SPE example ................. 22

6. Reflection per unit height for the strong SPE example 25

7. Normalized height-gain function for (a) true and

(b) calculated profiles ............................. 30

8. Calculated and true profiles after 20 iterations for

the C-layer example ................................. 31

9. Tradeoff between height resolution and noise-induced

variance for the C-layer example .................... 33

... ,......L n
 -xi-

TABLES

I. Assumed uncertainties in reflection coefficients for

the strong SPE example .............................. 19

2. Effects of noise for strong SPE example ............... 20

3. Reflection regions for strong SPE example ............. 26

4. Approximate peak reflection heights for strong SPE

example ............................................. 28

5. Assumed uncertainties in reflection coefficients for

the C-layer example ................................. 32

- s-

1iLZ
A blAi Z 1
 -1-

I. INTRODUCTION

Radio-sounding at VLF (3 to 30 kHz) and LF (30 to 300 kHz) is

often used to sample the ionosphere at heights below a1-out 100 km
where high-frequency (HF) ionosounds provide little lnformation.
Further, at altitudes below about 70 km the electron density is often
too small to be sensed by rocket-borne probes, and VLF/LF sounding
is the only sampling means available. Once ground-level VLF/LF re-
flection coefficients have been measured, they must be inverted to
infer ionospheric structure. This "inverse" problem is much more
difficult than the "forward" problem of calculating reflection coef-
ficients from a specified model ionosphere.
Inversion of ionosound data is often performed by trial and
error, with reflection coefficients being calculated from a number of
model ionospheres, and the model giving the coefficients that--in the
judgment of the analyst--most closely agree with experiment being deemed
correct. This report describes a method ior inverting VLF/LF reflec-
tion data to obtain ionospheric conductivity profiles without using
trial and error. We consider only altitudes below about 70 km where
the propagation can be assumed isotropic, although the method can be
extended to anisotropic propagation. This restriction on altitude
limits the applicability to disturbed conditions--such as solar-
proton events (SPE)--and some ambient daytime situations.
Our procedure extends the spectral method described by Parker
[19771 for implementing the Backus and Gilbert theory [1967, 1968,
1970, 1971], which was developed for inverting seismic or electro-
magnetic soundings of the earth's interior. We find that method to
be well suited to the ionosphere as well. Chuang and Yeh [1977]
applied the Backus and Gilbert theory to HF ionosound data. However,
they used geometric optics, whereas VLF/LF propagation must usually
be described by full-wave theory. Shellman [1973] and Morfitt and
Shellman [1977] inverted VLF/LF ionospheric reflection data, including
geomagnetic anisotropy, by performing a least-squares optimization of
the electron density and collision frequency. That approach, although
 -2-

direct, requires that the number and location of the sampling heights
be selected beforehand. The Backus-Gilbert method requires no such
preselection and imposes fewer constraints on the functional form of
the inverse solution. Moreover, it gives the relationship between the
detail that can be realized in the calculated conductivity height
profile and the quality and quantity of the reflection data.
Section II describes the mathematical basis of our inversion
method. Section III applies the method to actual and computer-
synthesized reflection-coefficient data and evaluates uncertainties
in the inferred profiles caused by noise and incomplete data. Sec-
tion III also presents a means of estimating the height range over
which waves at the sounding frequencies interact strongly with the
ionosphere; and, therefore, the height range over which the sounding
data can provide information on ionospheric structure. Section IV
gives the conclusions, and the Appendix presents certain mathemati-
cal details.
 -3-

II. INVERSE THEORY APPLIED TO IONOSOUNDING DATA

Inferring ionospheric height profiles from ground-based reflec-

tion data is one of a class of problems treated by inverse theory, the
object of which is to use limited and usually noisy data to determine
the maximum information about geophysical structure. Inverse problem-
are typically much more difficult than forward problems, which require
calculation of the signal emerging from specified propagation media.
This difficulty is due primarily to the nonuniqueness of the inferred
profile. Backus and Gilbert [1967] showed that--under fairly general
conditions--a finite data set is consistent with an infinite number
of possible profiles, if it is consistent with any one. In addition,
in many cases (including the one addressed here), the measured signal
has a nonlinear dependence on the propagation medium. This nonlinearity
could introduce nonuniqueness, because the unknown profile must be
constructed iteratively from a starting estimate and different start-
ing estimates can, in principle, produce different profiles.
Despite these difficulties, the importance of obtaining informa-
tion about the earth or its environment from remote probes has led to
at least a partial resolution of questions pertaining to uniqueness
[Backus and Gilbert, 1967, 1968, 1970; Oldenburg, 1978]. The non-
uniqueness due to incomplete data is fairly well understood and char-
acterized. That due to nonlinear dependence on the medium is still
unresolved, but seems to cause no difficulty for the numerical examples
that we give in Sec. II.

FRECHET KERNELS FOR ISOTROPIC IONOSPHERIC CONDUCTIVITY

The first step in the inverse solution is to characterize the
dependence of the reflection coefficients on the structure of the
ionosphere. For TM (vertically polarized) waves the reflection coef-
ficient R below the ionosphere can be defined as

R = H(d)/H(u) (1)
y ,y
 -4-

where H (u)u and H (d)

d
are incident and reflected magnetic intensities
y Y
measured at the same point, and y is the cartesian coordinate normal
to the plane of incidence. The equation governing R for plane waves
reflected from an isotropic ionosphere is [Budden, 1961]

2
2i dR2dR=Cn2(1
2 - R)2
2 _I_
2 (I + R) 2 (2)
k dz 2C

2 2 2
where n(z) is the complex refractive index, q = n - S , S and C
are the sine and cosine of the incident wave in free space, and k is )
the free-space wave number. An analogous equation exists for the TE
(horizontally polarized) reflection coefficient [Budden, 1961], but
we will consider only vertically polarized waves.
The general refractive index expression, including geomagnetic
effects, is

2e2 NY(z)(3
n= We_0 m (z)[W - iv Wz ± Wc() 'j

where the sum is over all charged species (electrons, positive and
negative ions), w is the angular frequency, e is the electronic charge,
FO is the permittivity of free space, and Ny, my, vy, and wCy are the
height-dependent number density, mass, collision frequency, and gyro-
frequency of the yth species. The isotropic equation (2) is valid
provided all important reflections occur below altitudes of 70 to
<
75 km, where v > wC" Moreover,
2C since << w C V at VLF/LF in this
low-altitude regime, n assumes the simplified form

n2(z) i O(z) (4)

WE
0

where the total conductivity is given by

a(z) - -
N -
e2 Y (5)

YY

- -- r --. -- 1~-~-
 -5-

Heavy ions make important contributions to the conductivity only at

altitudes below 50 to 55 km above which the conductivity may be con-
sidered entirely due to electrons. At altitudes where Eqs. (2) and
(4) are valid, reflection data can provide information only about the
total conductivity. Auxiliary information is needed to separate the
components appearing in Eq. (5).
Equation (2)--with n 2 given by Eq. (4)--can be integrated numer-
ically for a given w to obtain R(O) at the ground. We start the inte-
gration from the initial estimate

2C
R(z ) nC
m n 2C + q

where zm is chosen large enough for R to be slowly varying, and Eq.

(6) results from setting dR/dz = 0 in Eq. (2). This solution for the
reflection coefficient R(O), corresponding to a given conductivity
profile a(z), comprises the forward problem.
The inverse problem is to obtain the best possible estimate of
a(z) from R(O), specified at a finite set of frequencies. We con-
struct a by computing successive corrections to an initial estimate
a0(z), each correction being determined by relating a small change
6a(z) in conductivity at height z to a change 6Ri(0) in the ground-
level reflection coefficient at frequency f.
1
= w./27T.
1
We find the
necessary relationship by expanding an equation analogous to Eq. (2)
for R'(z) = R(z) + 6R(z) in terms of a'(z) = a(z) + 6a(z), retaining
only terms linear in 6R and 6a, and using Eq. (2). The resulting
linear first-order differential equation for 6R(z) is

d 6R + (z)6R(z) -Q(z) (7)

dz

where

P(z) [Cn2(1
C - R) + (I + R) (8)
n 2C
 -6-

Sro . ,2 2
E (I+R)] 1 ()2 (9)
Q(z) = -- 3 R)" 4.+ R

to give
Equation (7) can be integrated

dz Q(z) ep dz' P(z'

6R(O) -

0 f

+ 6R(zm) exp dz P(z)

term is negligible
component; so the second
P,conLtains a negative real
at the ground, and

6R(O) f dz G(z)6o(z) , (11)

in con-
which relates a small change
where C is the Fr~chet kernel,
coefficient at the
at altitude z to a change in reflection
ductivity
ground:

(;(Z,) j.i- (I + R'

2 n

R) + 2 I R)jI (12)j
exISj dz{ Cfl(I-
I If n2
 -7-

CONSTRUCTION OF THE INVERSE SOLUTION

As emphasized by Backus and Gilbert 11967, 1968, 1970], the esti-
mation of an unknown propagation medium from a finite set of field
measurements requires two steps. First, we must construct a conduc-
tivity profile that reproduces the data. Second, and iqually important,
we must characterize the class of solutions to which our profile and--
presumably--the true profile belong. Here we address the first step.
The uniqueness of our inverse solution is considered in the following
subsection.
To construct the inverse solution we must find a conductivity
profile that matches a set of reflection coefficient data at a number
(n) of frequencies. Following Oldenburg [19781, we pick an initial
0 0
estimate a , compute R.i from Eq. (2) at each frequency f.' and set

6R =I R - R (0) i = 1, ... ' n , (13)

i 1 i

where superscripts denote the order of the correction (i.e., the

number of iterations), and subscripts index the frequencies. We then
O 1
express the first-order correction to 0 , 6a , as

n
6a (z) = cxG (z)
( , (14)
j=l

where the a. are found by solving the n × n set of equations

n
r 6 (15)
= 1

and

z
m
= f dz Gi(z)G (z) . (16)
0
 -8-

Next, R1 is computed from Eq. (2) using the new conductivity estimate

a (z) = Y (z) + 6a (z), and the process is successively repeated. At

each iteration N we use the rms relative error

n 21 /2
N = SR.N/Rata) (17)

i= 1

as a measure of convergence.

In practice, we use tne following modifications to the above

procedure:

1. Because the data and Fr~chet kernels are expressed in com-

plex form, the differential conduccivities 6o I produced by

Eq. (14) can have imaginary parts. We have found it suffi-

cient to use the fact that the conductivity in Eq. (14) must

be real and simply suppress the imaginary part of the cor-

rection at each iteration. A more rigorous--and complicated--

approach would be to express all quantities in real form at

the outset.
2. As discussed by Oldenburg, direct use of the output correc-

tions to the conductivity at each iteration tends to pro-

duce oscillatory profiles, and a nonphysical build-up of

conductivity at low heights. These pathologies require that

each correction be substantially smoothed and filtered before

being used in the next iteration. In addition, we do not

allow the total conductivity at any iteration to become

negative but redefine it as zero, if required, before pro-

ceeding. These steps slow the convergence, but appear

necessary to produce realistic profiles.

3. Following Parker [19771 we diagonalize the matrix F before

inversion, so the expansion JEq. (14)] is actually performed

in the eigenfunctions of the r matrix. This step allows each

term in the diagonalized basis to be associated with a

LAL_I'.- -
 -9-

relative uncertainty in the data, and components that intro-

duce large statistical uncertainties can be omitted. The
Appendix gives the details of this spectral expansion method.

UNIQUENESS OF THE INVERSE SOLUTION

The above algorithm is iterated until the relatve rms error
[Eq. (17)] asymptotically approaches a minimum. We have then done as
well as possible in constructing a conductivity profile with the data
and starting estimate used. The remaining task is to assess the
validity of the profile at various heights. As discussed above there
are three sources of nonuniqueness:

1. Noise in the reflection coeificient data causes uncertainties.

2. The data set is finite, and therefore incomplete.
3. The conductivity profile and reflection coefficient data are
nonlinearly related.

The first error source is easily assessed with multivariate

error analysis. The spectral expansion method (see Appendix) provides
a simple means of estimating the effect of measurement uncertainties
on the accuracy of the inferred conductivity profile. More important,
variances in the reflection data are included in the solution in a
manner that deemphasizes frequencies at which measurements are ex-
ceptionally noisy.
The second type of nonuniqueness has received the most attention
in the literature, since it is present regardless of the care taken
to reduce experimental uncertainty. It is therefore the ultimate
limitation on using a finite data set to characterize an unknown func-
tion. Before specifically addressing this type of nonuniqueness, we
consider a somewhat similar situation, which is helpful in understand-
ing inaccuracies caused by incomplete data.
The Fourier series representation of an arbitrary periodic func-
tion O(x), 0 < x < 27r, is given by
 -10-

-
W= C e (18a)

where

21r
C= dx $(x) e- . (18b)

It is well known that--subject to some weak requirements--the i'zit

set of coefficients {C n} completely defines $. If, however, only a
jin-c set -N : n ! N of Fourier coefficients is given, the high-
frequency content of (pwill be lost; and any ' that differs from

only in these neglected components will have the same coefficients

Cn , for Inl r N.
To quantify the uncertainty caused by omitting the coefficients
of index Inl > N, we let ($)denote the estimate for 1 based on the
incomplete set of Cn, and write

N
Iinx
OW(x) = - C e (19)
n= - N

By using Eq. (18b), we find

2Tr

OW(x)) = J dx' O(x')AN(x, x') , (20a)

0

where

N
1')
N(X, x') =2Tr ein(x - x') (20b)

n--N
AN is the Dirichlet kernel, and is a measure of how closely 0 can be
resolved by the finite harmonic set. The width of AN is proportional
to 1/N, so the resolution improves as the number of "data" points
becomes larger. In the limit N - -, (u) - , and the resolution be-
comes perfect. This behavior occurs because AN approaclies a Dirac
delta function for large N.
The above approach can be used in the ionospheric inverse problem
to quantify the loss of resolution due to having incomplete reflection
data. For a given iteration, Eq. (14) can be written

j=l

where we place ) around 6a to emphasize that only a finite set of

frequencies is available, and our estimate can therefore differ from
the "true" correction 6S. By using Eqs. (15) and (11), we find

n
aj= r-l' 6Rk(0)

k=l

n zm
klr-1
- dz0(22
dz' Gk(z')Sa(z') (22)

k-l 0

and substitution into Eq. (21) gives

(6cl (z)) =f dz' 6o(z')A(z, z') (23)

with
 -12-

A(z, z') P 1 Gj (Z)Gk(Z') (24)

j,k-l

The averaging function (Eq. (24)] is the analog of the Dirichlet kernel
for Fourier series. Variations of the conductivity on a scale shorter
than the width of A(z, z') about a given height will not be resolvable
by the data set used to generate the profile. Even though we con-
structed nte dveraging function for a small variation in conductivity,
Backus and Gilbert [1968] prove that--to order j6o(z) 12--Eq. (24) pro-
videz i measi,1re of the height resolution of the total conductivity.
show in the Appendix that Eq. (24) can be rewritten in essen-
t.'iii, ,"icform

kzcn
A(z, z') H(z)H(z') , (25)

i= 1

where the functions Hi(z) are produced by orthonormalizing the set of

Fr~chet kernels Gi(z). As for the orthonormal Fourier series basis
function discussed above, the H. functions can be ranked in order of
1
increasingly oscillatory behavior. The diagonal form of the averag-
ing function given by Eq. (25) is convenient for evaluating the im-
portance of retaining the more oscillatory components. In general,
a better fit to the data is achieved by increasing the number of com-
ponents used in the construction of 6((z). However, the resulting
conductivity profile can exhibit nonphysical oscillations if too many
terms are retained. The optimum number of components H. depends
partly on the estimated noise in the data, and partly on the analyst's
intuition regarding profile smoothness. The relation of resolution

tTbis situation is analogous to that which can be encountered in

fitting a high-order polynomial to a large number of data points. The
polynomial can match the data exactly at each point, but exhibit
spurious oscillations between the points.

- -I
 -13-

to noise is treated in the next section and the Appendix. Its rela-
tion to the analyst's intuition is difficult to quantify, however.
The remaining source of nonuniqueness is the nonlinear relation
between conductivity and the ground-level reflection coefficient.
This source is by far the most difficult to characterize in general,
because there is no guarantee that other "globally" iifferent conduc-
tivity profiles would not be obtained with different initial esti-
0
mates a . The most practical means of exploting this possibility is
to choose widely different starting estimates and see if the calcu-
lated profiles agree, aside from minor variations. We have found that
our algorithm produces iterative solutions even with a null initial
estimate (a 0 = 0). Moreover, calculations made from the same data
sets--but with various nonzero estimates--gave virtually identical
profiles. These results strongly suggest (but do not prove) that--
in practice--the nonlinear dependence of the reflection coefficients
on ionospheric structure causes no loss of uniqueness.
 -14-

III. NUMERICAL EXAMPLES

This section applies the inversion method developed in Sec. II

and the Appendix. Recall that our assumption of isotropic propagation
restricts this application to conditions where all important ref lec-
tions occur below 70 to 75 km. We use the examples given below to
test the ability of our inversion method to obtain correct conduc-
tivity height profiles from a set of ground-level reflection coef-
ficients; and to assess the effects of measurement noise and incom-
plete data on the accuracy and height resolution of our results. The
examples were chosen to illustrate the salient aspects of our approach
and are intentionally simple. our method should work equally well for
more complicated situations.

I STRONG SPE

Our first example is

based on high-latitude satellite measure-
ments of incident particle fluxes at 1508 UT time on 4 August 1972.
These measurements were made near the peak of one of the strongest
SPEs on record. Height profiles of electron and ion densities were
calculated from the measured fluxes with Lockheed's air-chemistry code,
and provided to the authors by Reagan [1978]. These densities, along
with nominal height profiles of electron and ion collision frequencies
and masses were used with Eq. (5) to obtain the conductivity height
profile shown (solid lime) in Fig. 1. This profile is an example of
an extremely severe natural disturbance.
Here we illustrate application of our method rather than perform
detailed calculations for specific situations. Therefore, for conve-
nience, we approximate the strong SPE profile with the often-used form

a)=0 0 e ,(26)

-15 -l
where .94 x 10
1o= mhos/m and 0.39 km- . Figure 1 shows the
agreement between the approximation [Eq. (26)] and the actual conduc-
tivity profile.
 -profile from meuured flux
-- excponential fit (Eq.(26) I

lore

30 0/ s
/egt(m
/rn
Fi.ICndciiyhegtpofl o P
 -ib-

Inverse Solution
To provide a starting point for our inverse solution, we have
calculated TM reflection coefficients for the profile given by Eq.
(26) by integrating Eq. (2) for an incidence angle of 65 deg at fre-
quencies of 4, 6, 9, 14, 21, 30, and 48 kHz. Our reason for choosing
this set of frequencies is discussed below. These coefficients serve
as artificial "data," which are analogous to measured data that would
be available from ionosounders. To suppress the unphysical build-up
of conductivity below 35 km, we use the weighting technique of Olden-

burg [1978] (see Appendix), and choose

0 z-b
w(z) = (27)
1 - exp I-y(z - b)] , z ; b

=
for the weighting function w(z) with b = 35 kin, y 0.5 km-1 . We
compute successive corrections to the conductivity using the spectral
0
expansion method and a null (a = 0) starting estimate. Therefore,
our initial profile contains no information about the true profile.

Oscillations in the corrections 6a for each iteration are controlled

by our use of a weighted average of the correction at height z (weight
0.5) and at 0.7 km above and below z (weight 0.25).
Figure 2 shows the computed conductivity profile following the
I.;', 5th, 100/,, and 15th iterations, the true profile (dashed lines),
and the relative rms error c. The computed profile was still slowly
increasing near 50 km after the 1502 iteration, and more iterations
presumably would have given closer agreement with the true profile
between 49 and 50 km. Figure 3 gives an expanded plot of the calcu-
lated and true profiles after the 15th iteration. The agreement is
good at altitudes between about 4() and 49 km. Outside of this alti-
tude range the calculated conductivity is much smaller than the true
conductivity. We give the physical reason for this disagreement below.

Sensttivity to Noise and Limited hata

To evaluate the effect of measurement noise on the inferred con-
ductivity, we include hypothetical uncertainties S1 in the reflection
 -17--

00

co oo CM4
0c co (

9- . 4-
E,~ 0"

- IL,

CL

4J

oo - 0cc

AI!Ao 0WS4
npuo
 10-5
/
I

true profile/I
I
/
I
I

/
-6
0
E

0 calculated profile

10-8-/
/
I
I

/
I
I
/

109-
30
I I
40 50 60
d
Height (kin)

Fig. 3--Calculated and true profiles after 15 iterations for the strong SPE
example. Vertical bars show uncertainty caused by hypothetical noise
in the reflection coefficients. Horizontal bars show uncertainty
caused by incompleteness of the reflection data.

.&L . ..
 -19-

coefficients used in our calculations. The values for S i are derived

from typical measurement uncertainties (Rasmussen [19811) and are
listed in Table 1.

Table 1

ASSUMED UNCERTAINTIES IN REFLECTION COEFFICIENTS

FOR THE STRONG SPE EXAMPLE

Frequency Is/
(kHz) i
i/R(O)

4 0.3
6 0.3
9 0.3
14 0.25
21 0.1
30 0.1
48 0.15

By including Si. we can examine sensitivity of the profile at various

heights to me;-urement noise, even though the computer-generated "data"
are noiseless.
The variance in conductivity is defined as

2 2
Var [W(z)] (2(z))
W - (a(z)) (28)

where ( ) denotes averages. The Appendix shows that

kcn
Var [o(z)] = [w(z)]2 . 1 2 (29)
i=I '

where the sum is over the positive-definite contributions from each

diagonalized component H.(z). The eigenvalues X for the 15th iter-
ation are listed in Table 2, as is the estimated variance at height
44 km formed by truncating Eq. (29) at different values of k. Table
2 shows that retention of fewer terms in Eq. (29) causes a decreased
 -20-

Table 2

EFFECTS OF NOISE FOR STRONG SPE EXAMPLE

Cumulative
[Var (0)1/2
Eigenvalue (mhos/m)

I = 1 1.81 . 10-6 k = 1 1.70 x 10 - 9

7 = - 9
1= 2 4.32 x 10- k 2 3.89 x 10

I = 3 9.46 x 10 k 3
= 1.69 x 10 8
-9 - 8
i = 4 1.72 x 10 k = 4 2.38 x 10
- I0 - 8
i = 5 4.51 x 10 k = 5 9.19 x 10
- 7
1 = 6
-
1.95 x 10 11 k = 6 2.35 x 10
- 13 - 6
i = 7 2.65 x 10 k = 7 5.40 x 10

variance. If variance were the only consideration, retention of only

the first spectral component would provide the "best" solution. How-
ever, as discussed below, there is a tradeoff between noise variance
and height resolution, and the solution retaining only i = 1 would
provide poor resolution.
We evaluate resolution of the inferred profile by the general-
ized Dirichlet kernel method described in Sec. II and the Appendix.
As for the variance, the achievable height resolution depends on the
number of spectral components retained. Figure 4 shows the modulus
of the averaging function A at height 44 km formed by truncating the
Dirichlet representation JEq. (25)] at various values of k. The best
resolution--i.e., the narrowest averaging function--is achieved by
retaining all seven terms. Unfortunately, as discussed above, the
variance increases with the number of terms retained, and a compromise
is needed.
Figure 5 illustrates the tradeoff between resolution and noise-
induced variance in the conductivity height profile computed for the
strong S11E. The resolution--defined as the full width at half-peak

This behavior occurs because higher-order terms tend to be

oscillatory and, therefore, exhibit a larger variance.
 -21-

12, K=5
20

16
811
12
I I
4 I8

I 4

0 30 01
40 so 60 70 30 401 50 60 70
1

20 I K 2 24 K=6

16 20

12 16
IAI 12
8
8
4 4

30 40 50 60 70 030 40 50 60 70

IAI
[K =3 35 -K =7

16 25
12
, 20
15
010 I

30 50 60 70 30 40 50 60 70
Height (km)

20
16.

12

4
0I
30 40 50 s0 70

Height (kin)

Fig. 4--Magnitude of the averaging function A for the strong

SPE example (height 44 km)
 -22-

10-67

I
6
0
E 7

4
3

K 1
10- 1 1 I
0 2 4 6 8 10 12 14
Resolution (kin)

Fig. 5--Tradeoff between height resolution and noise-induced variance

for the strong SPE example
 -23-

of A--decreases as the number k of retained terms is increased, but

at the expense of a greater variance. The vertical and horizontal
1/2
bars in Fig. 3 show the resolution and [Var (a)] , respectively
for k = 5, which we chose as a good compromise. If we had retained
the maximum number of terms (k = n = 7), the resolutior would improve

from 5.1 to 3.5 km, but [Var ()]1/2 would increase by nearly two
orders of magnitude. The optimum number of terms depends in the num-
ber of frequencies at which reflection data are available, and tie
amount of measurement noise. All terms (k = n) should be retained
if noise is negligible, but only a few retained if the data are noisy.
The optimum value for k must be selected on a case-by-case basis,
using tradeoff graphs such as Fig. 5.

Physical Interpretation
We next show physically why the calculated profile shown in
Fig. 3 agrees with the true profile only at heights between about 40
and 50 km. In addition, we give diagnostic methods for estimating

the height range over which a set of reflection data can yield accu-
rate conductivity profiles.
Ground-level reflection coefficients contain information about
only those heights from which significant energy is returned. To
estimate the heights at which important interactions between the in-
cident field and the ionosphere occur, we use the following approxi-
mate expression (Field and Lewinstein 119781) for the TM reflection
coefficient:

R(0) = - dz [Q exp [-i2w (z)/c] , (30)

0

where

J(z)
f
0
z
dz' q(z') , (31)
 -24-

and Q = q/n 2 . Equation (30) is the second-order WKB approximation,

and is strictly valid only for weak reflections. Nonetheless, we
will show that it provides useful insight into the reflection process.
The integrand 5 of Eq. (30) is the contribution of a thin stratum
at height z to the ground-level reflection coefficient. The quantity
in brackets represents conversion of upgoing to downgoing waves in
the stratum; the exponential term represents the round-trip transmis-
sion loss suffered in propagating to and from height z.
If we set

n 2 (z) = I - ia(z) , (32)

the reflection-per-unit height Ji becomes

- dag I 2C 2 + *
id 1(I22 exp [-i2wO(z)/c]I , (33)
S dz - iC)(C - ia)

which has magnitude

4 +i
44221/2 1/+_
=() I- ( 1 -a4 2S2 l12 exp [2wIm *(z)cI

(34)

where

Im *(Z) - ILf (_C4 + a2)1 2 dz' .(35)

Figure 6 plots for the exponential conductivity profile

given by Eq. (26) at the frequencies used in the inversion, and Table
3 summarizes information from the figure. The quantity zR denotes
the height at which 121 attains its maximum value, and zB and zT de-
note heights below and above zR, respectively, where IiJis one-half
its peak value.

L .- - • - .
 -25-

00
1S%

C,

4-

L-

-WVed uilooiet
004 ~un
 - ---
----
--
--

-26-

Table 3

REFLECTION REGIONS FOR STRONG SPE EXAMPLE

Frequency ZR ZB ZT Wdh(T B
(kHz) (kin) (kmn) (kmn) (kmn)

4 44.1 38.8 50.0 11.2

6 44.1 39.6 49.6 10.0
9 44.4 40.3 49.3 9.0
14 44.8 41.1 49.1 8.0
21 45.5 41.7 49.0 7.2
30 45.8 42.2 49.0 6.7
48 46.2 42.7 49.0 6.3

The maximum contribution to the ground-level reflection co-

eff icient comes from heights near z R9 whereas relatively little con-
tribution comes from heights below zB or above zT' Therefore, we
expect the ground-level reflection coefficient to provide the best
information on ionospheric conductivity in the height range bounded
by zB and zTV and only poor information outside of this range.
Table 3 shows that zRis a slowly increasing function of fre-
quency, but remains near the center of the 40 to 50 kmn height range
over the entire band considered. The lower height zB increases from
39 to 43 km as the frequency is raised from 4 to 48 kHz, but the upper
height zT remains between 49 and 50 km. From Table 3 and Fig. 6 we
conclude that, for the frequencies used, the reflection coefficients
give little information about heights below 40 or above 50 kmn; and,
even within the 40 to 50 km height range, the information is concen-
trated between about 44 and 50 km. This conclusion is consistent
with Fig. 3, which shows excellent agreement between the calculated
and true profiles between 44 and 49 kin; fair agreement between 40 and
44 kmn; and a total breakdown below 40 or above 50 km.
The insensitivity of z7 to frequency indicates that the height
range sampled by VLF/LF sounding could not be expanded substantirlly
by using higher frequencies; and the use of frequencies much lower
than 4 kHz is probably not practical. Thus, for the strong SPE ex-
ample, VLF/LF sounding would appear incapable of determining the
 -27-

state of the ionosphere below 40 or above 50 km. Finally, the fact

that zT is much lower than 70 km provides a posteriori justification
of our neglect of geomagnetic effects.

Spacing of Frequencies
We have found that the best results are obtained by choosing an
unequal spacing for the frequencies used. On physical grounds we
expect that the optimum choice of frequencies would correspond to a
uniform spacing of the heights of peak reflection z We therefore
need a mapping between reflection height and frequency. Unfortunately,
Eq. (34) is inconvenient for this purpose. However, Field and Engel
[1965] showed that--subject to certain restrictions--the major reflec-
tions take place within a few kilometers of the altitude where

f C2 (36)

Specifically, by examining the second-order correction to the WKB

approximation, they showed that substantial reflection can occur only
at altitudes where

--
' 2k-- (37)
n
32 2

where we have used Eq. (26). Provided that q2 << n 2 , the peak value
of the left side of Eq. (37) is easily shown to occur at the height
where a is given by Eq. (36). By using Eqs. (4), (25), (32), and
(37), we find
1r - C2 27f 01
zR 2 C J (38)

Table 4 shows the values of zR given by Eq. (38) for the strong SPE
example. Note that the logarithmically-spaced frequencies correspond
to evenly-spaced reflection heights. The differences between the
approximate zR in Table 4 and the more accurate ones in Table 3 do
not appear to strongly affect the optimum choice of frequencies.
 -28-

Table 4

APPROXIMATE PEAK REFLECTION HEIGHTS

FOR STRONG SPE EXAMPLE

Frequency zR [Eq. (38)]

(kHz) (kin)

4 40.7
6 41.8
9 42.8
14 43.9
21 45.0
30 45.9
48 47.1

In practice, we selected the frequency spacing for the strong

SPE model by (1) using arbitrarily spaced frequencies to construct
a preliminary profile, (2) making an exponential fit to this prelim-

inary profile, and (3) using Eq. (38) to select logarithmically spaced
frequencies corresponding to evenly spaced zR.

Height-Gain Function

It is desirable to have a method of assessing profiles as they

are being constructed. We have found the height-gain function H(z)
to be useful in this regard. Here we normalize H so that--below the
ionosphere--it equals the upgoing wave, which is normalized to unity
at the ground. The resulting expression is (Field, Lewinstein, and
Dore [19761)

= -
H I + R(O) exp -i dz',n2W (39)
1 + R(z) c A(z')

where R(z) is the reflection coefficient and A(z) is the admittance

A(z) = I + R(z) 1
1 - R(z) C
l -29-

Figure 7a shows H(z) for frequencies of 4 and 48 kHz and the

exponential fit to the strong SPE IE'q. (26)1. Comparison with Fig. 6
shows that most reflections occur in the height range where 0.5 < H(z)

1. At lower altitudes the ionosphere is too rarefied to reflect the

wave strongly; at higher altitudes the signal is so wck that tile re-
turn is small even if the reflectivity is high. Thu shaded region
thus indicates where both the signal -,nml the ref ectivity are strong.
Figure 7b shows H for the calculated profile, which we have assumLd
retains its 49-km value at greater altitudes. Again, we see that
0.5 < H < 1 provides a measure of the height range over which the
calculated profile is valid; viz., between about 40 and 50 km. Com-
puting H during profile construction gives a criterion for terminating
the iterations. When The major contribution for an iteration falls
outside the height range defined by 0.5 < H < 1, we have done as well
as possible and accept that profile as the final one.

MODEL C-LAYER
The above example pertains to conditions where the conductivity
increases monotonically with height. We now consider a qualitatively
different situation, where the lower ionosphere exhibits a well-defined
layer. Specifically, we used the model given in Fig. 8 (dashed lines),
which was shown hv Rasmussen, Kossey, and Lewis 1.19801 to give close
agreement between measured and calculated ionosound returns. This
-7
profile consists of a sharply bounded layer of conductivity 1.8 x lO
mhos/m extending from 63 to 69 km.
The magnitudes of the reflection coefficients were measured at
an incidence angle of about 64 deg b, Rasmussen, Kossey, and Lewis, and
were input to our algorithm at 10 evenly spaced frequencies between
5 and 50 kHz. Since no phase measurements are available, we had to
manufacture phase "data" by integrating Eq. (2) to find Arg R(O) for
the "true" profile in Fig. 8. We used the weighting function given
by Eq. (26) with y - 0.5 km' i and b = 50 kin; and, as before, we used
a null (a0 = 0) initial profile and a three-point smoothing filter.

The logarithmic spacing discussed above is not suitable for the

type of profile being considered here.

-- L---,--- I-
 - 3(-

weak reflection - main weak signal -t

2interaction

- 1

0.5-

4 kHz
-- 48 kHz

0.11
0 10 20 30 40 50 60 70
Height (kin)

(a) Exponential fit [Eq.(26)]

weak reflection main weak signal

interaction
2
N

0.5

0.1 0 A
10 20 3,0 40 50 60 70

Height (km)

(b) Calculated profile (Fig. 3)

Fig. 7--Normalized height-gain function for (a) true

and (b) calculated profiles

• _ li _ . ...
-. -
.- " .... .. . V. . . . . .,
 2.2 -

2.0-

1.8- r - - - true profile

calculated
1.6 - profile-

E 1.4-

a I
._•
0 1.0 -I

o.8-
E

0.6-

0.4-
0
0.2 - I

58 62 66 70 74

Height (kin)

Fig. 8--Calculated and true profiles after 20 iterations for the C-layer
example. Vertical bars show uncertainty caused by hypothetical noise
in the reflection coefficients. Horizontal bars show uncertainty
caused by incompleteness of the reflection data.
 -32-

Figure 8 shows the profile calculated from the above reflection

data after 20 iterations. The calculated and true profiles agree well,
particularly with regard to the peak conductivity and location of the
layer. The failure of the inverse solution to reproduce the sharp
boundaries is caused by the limited height resolution achievable with
incomplete data.
To examine the effects of measurement noise on the variance of
the computed profile, we use the uncertainties in reflection coef-
ficients given in Table 5. As before, these values are representa-
tive of actual measurement uncertainties. The tradeoff between vari-
ance and height resolution at a height of 66 km is shown in Fig. 9 for
the calculated profile after 20 iterations. The averaging functions
for k < 4 did not display well-defined peaks, so their widths could
not be determined.

Table 5

ASSUMED UNCERTAINTIES IN REFLECTION

COEFFICIENTS FOR THE
C-LAYER EXAMPLE

Frequency i/Ri(0)l
(klz)

5 0.3
10 0.25
15 0.18
20 0.15
25 0.1
30 0.1
35 0.1
40 0.1
45 0.15
50 0.15

The vertical and horizontal bars on Fig. 8 correspond to the

1/2
uncertainties caused by noise ([Var ()] ) and incomplete data
(height rogolution) with k = 8. Note that the agreement between the
calculateu and true profiles is much closer than indicated by the
 10-6

10-7

8
~h 7

10-8
6

K-4

l 9 I I '
6 7 8 9 10 11 12 13 14

Spread (kn)

Fig. 9--Tradeoff between height resolution and noise-induced variance

for the C-layer example

I
 -34-

horizontal bars. This agreement suggests that the height resolution

is, in fact, far better than the width of the averaging function A.
Shellman [1973] also concluded that this theoretical estimate of
height resolution is, in many cases, much too pessimistic. In prac-
tice, the achievable resolution is probably best determined by exer-
cising our inversion algorithm for a number of generic profiles, such
as the two examples given above.

- ,
 -35-

IV. CONCLUSIONS

Our extension of the Backus-Gilbert theory appears well suited

for inferring ionospheric properties from ground-level VLF/LF iono-
sound data. It permits height profiles of ionospheric parameters to
be determined without resorting to trial and error and explicitly
gives the relationship between uncertainties caused by noisy and in-
complete data. The analyst can therefore rhoose the best compromise
between variance and height resolution of calculated profiles.
We consider only altitudes below about 70 kmn, where the propa-
gation can be assumed isotropic. This altitude restriction limits
application to certain ambient daytime or disturbed conditions. More-
over, only the profile for total conductivity can be found from the
reflection data. Determination of profiles for the density and col-
lision frequency of each charged constituent requires auxiliary in-
formation. The method can be extended to include geomagnetic effects.
The versatility of our approach is evidenced by its satisfactory
performance for two qualitatively different examples--a strong SPE
where the conductivity increases quasi-exponentially with altitude
and a model C-layer exhibiting a well-defined conductivity peak. In
each example, the iterative solution converges toward the true pro-
file despite our use of initial estimates that contain no information
about ionospheric structure. The nonlinear relationship between re-
flection coefficients and ionospheric conductivity apparently causes
no difficulties with uniqueness.
We found the height resolution of the calculated profiles for
the two examples to be much better than predicted by theory. Other
investigators also have reported that the theoretical limit on height
resolution is often too pessimistic. The resolution achievable in
practice is best determined empirically by inverting computer-generated
reflection coefficients for model conductivity profiles that corre-
spond to generic ionospheric conditions.
Ground-level reflection data contain information about only those
heights from which significant energy is returned. Our calculated

ALL-~ -v--------- ~iT T '

 -36-

profiles therefore agree with the true profiles only over a limited
height range, which corresponds to the region where the most signifi-
cant reflections occur. This reflection-height region provides a
physical basis for judging where (a) the calculated profiles are
reliable, or (b) the interaction between the ionosound signal and
the ionosphere is too weak to permit a valid inverse solution.
 -37-

REFERENCES

Backus, G., and F. Gilbert, "Numerical Applications of a Formalism

for Geophysical Inverse Problems," Geophys. J. R. Astr. Soc.,
Vol. 13, 1967, pp. 247-276.

"The Resolving Power of Gross Earth Data," Geophys. J. R. Astr.

Soc., Vol. 16, 1968, pp. 169-205.

-----, "Uniqueness in the Inversion of Inaccurate Gross Earth Data,"

Phil. Trans. R. Soc. A, Vol. 266, 1970, pp. 123-192.

Budden, K. G., Radio Waves in the Ionosphere, Cambridge University

Press, 1961.

Chuang, S. L., and K. C. Yeh, "A Method for Inverting Oblique Sound-
ing Data in the Ionosphere," Radio Sci., Vol. 12, 1977, pp. 135-140.

Field, E. C., and R. D. Engel, "The Detection of Daytime Nuclear

Bursts Below 150 km by Prompt VLF Phase Anomalies," PROC. IEEE,
Vol. 53, 1965, pp. 2009-2017.

Field, E. C., M. Lewinstein, and M. A. Dore, Effects of Antenna Ele-

vation md Inclination of VLF/LF Signal Structure, Pacific-Sierra
Research Corporation, Report RADC-TR-76-375, December 1976.

Field, E. C., and X. Lewinstein, Analysis of Electron Density-Profiles

in the Lower Ionosphere, Pacific-Sierra Research Corporation,
Report RADC-TR-78-79, March 1978.

Gilbert, F., "Ranking and Winnowing Gross Earth Data for Inversion
and Resolution," Geophys. J. R. Astr. Soc., Vol. 23, 1971
pp. 125-128.

Jackson, D. D., "Interpretation of Inaccurate, Insufficient and

Inconsistent Data," Geophys. J. R. Astr. Soc., Vol. 28, 1972,
pp. 97-109.

Jupp, D. L. B., and K. Vozoff, "Stable Iterative Methods for the

Inversion of Geophysical Data," Geophys. J. R. Astr. Soc., Vol. 42,
1975, pp. 957-976.

Morfitt, D. G., and C. H. Shellman, A Technique for Obtaining D-Region

Electron Density Profiles from VLF Reflection Coefficients, Naval
Ocean Systems Center, Report 781, November 1977.

Oldenburg, D. W., "The Interpretation of Direct Current Resistivity

Measurements," Goophys., Vol. 43, 1978, pp. 610-625.
 -38-

Parker, R. L., "Understanding Inverse Theory," Ann. Rell. Earth Planet.

:'(.., Vol. 5, 1977, pp. 35-64.

Rasmussen, .1. E., private communication, 1981.

Ramussen, J. E., P. A. Kossey, and E. A. Lewis, "Evidence of an Ion-

ospheric Reflecting Layer Below the Classical D Region," J. Geophys.
?. Vol. 85, 1980, pp. 3037-3044.

Reagan, J., Lockheed Palo Alto Research Laboratory, private communica-

tion, 17 March 1978.

Shellman, C. It., Datemination of D-Region Electron-Density Distribu-

"onsq from Rzdio I'ropagation Data, Naval Electronics Laboratory
Center, NELC/TR 1856, January 1973.

Wiggins, R. A., "The General Linear Inverse Problem: Implication of

Surface Waves and Free Oscillations for Earth Structure," Ret.
;,-opy.s. .,'r1zc Phyr., Vol. 10, 1972, pp. 251-285.

.
 -39-

Appendix

THE SPECTRAL EXPANSION METHOD

Although the inverse solution and its interpretation can be

carried out entirely in terms of Fr~chet kernels G.(z) as discussed
1
in Sec. 11, significant advantages result from using the set of basis
functions formed by orthogonalizing the set {G.(z)}. These advan-
1
tages include (I) improved numerical stability in the iterative inverse
solution, (2) the ability to construct resolution functions that are
direct generalizations of Dirichlet kernels for Fourier series, and
(3) explicit presentation of the tradeoff between achieving high
resolution and minimum noise-induced errors. Our analysis has been
guided by Parker's [1977] review, which was based in turn on work of
Gilbert [1971], Jackson 119721, Wiggins [1972], and Jupp and Vozoff
[19751.
The first-order variation in the reflection coefficient 6R.(O)
at frequency fi, due to a small change in 60, is given by Eq. (11),
from which it follows that

6R.(o) m i(z)

S f dz . 6o(z) , (A.1)
0

where S i represents the measurement uncertainty in R1 (O). Following

Oldenburg [19781 and Parker [1977] we slightly generalize the matrix

r given by Eq. (16) by including the weighting function w(z) and the

error estimates Si:

Z

fm Gi(z) C1 (z) 2
1' dz i~-* Iw(z)I (A.2)
ij Si S

The use of w(z) allows us to emphasize height regions where significant

conductivity is anticipated, and suppress the nonphysical build-up of

ALL
 -40-

conductivity at low heights. By construction r is a Hermitian,

positive-definite matrix and can therefore be diagonalized by a
unitary matrix U as

St U = A, (A.3)

where (A)i. XjS.." Next, we define the function Hi(z) as

-j j. i

H.(z X-112 : G W(z) (A.4)

The set {1i(z)} is orthonormal since

dz Ht(z)H (z) = t1 / 2 X-1/2 (Utik

0k,P
z
m
f Gk.(Z) Gg(z) ]2

X dz [w(z) 1JU
Sk S*
f

Ai
112 12 : U t ik
A kkY
k,9.

Xi-1/2 X-/2 X16 j 6ij (A.5)

The orthonormality allows the conductivity correction 6o(z) to be

expanded as

n
&o(z) aiHi(z)w(z) + 6&(z) , (A.6)
i-1

- - ... , _ 7 7- ... . . . ..-

 -41-

where

f zin
a = dz H.(z)
.i w(z)
0

= 1/2iij (Ut) fs. dz 60(z)

j 0

= x-I/2 (U5 ij S.-- (A.7)

and 6a(z) is any function orthogonal to the finite set {H.}. Equa-
tions (A.6) and (A.7) provide the means of computing the correction
to the conductivity at each iteration by taking

N data N-i
6R = i Ri (A.8)

One advantage to finding 6 in this diagonalized basis is that pos-

sible numerical instabilities in computing r - , which is needed in
the nondiagonal solution, are avoided.
The spectral expansion leads to a simple method of estimating
the effect of uncertainties in the measured reflection coefficients
on the inferred conductivity profile. Equations (A.6) and (A.7) re-
late the effect of small variations in the ground-level reflection
coefficients on the conductivity at height z. Equivalently these
equations may be interpreted as the uncertainty in conductivity due
to uncertainties in the data. The variance in conductivity then
becomes

Var fa(z)I tw(z) 1 (<a i a*)Hi(z)H j ( z ) , (A.9)

1,3

- - -----------
. : , ------- ~ -k---- I
 -42-

with

(a. a.) = 1/ -1/2 (U )iU k (A.1O)

j i
I ~
k,2
''** k SI
S£

where the ( ) denote ensemble averages. Because S.

1
are the standard
deviations in the data and the noise is assumed to be uncorrelated,
the covariance matrix for the data is simply

_ 6(A.11)
Rk Rj\

and the unitarity of U gives

j
(ai a*) = i (A.12)
j ii

Substitution into Eq. (A.9) produces

Var [(z)] = [w(z)12 l Hi(z) (A.13)

1
i

Equation (A.12) states that the expansion coefficients ai are

s;tatistically independent and have variance 1/Xi, where Xi are the
positive-definite eigenvalues of the r matrix. The spectral expan-
sion [Eq. (A.6)] can therefore be interpreted as a superposition of
terms having progressively greater statistical uncertainty. The
digonalization process recombines the original Frichet kernels Gi
to produce components Hi, which allow the conductivity features that
can most reliably be inferred from the data to be identified. In
practice, a large reduction in Var (o) can be achieved by suppressing
Few of the smallest elgenvalue terms from Eq. (A.6). Because eigen-
Value components with large X1 tend to be smoother than those with
small Xi. this truncation removes the highly oscillatory terms, which

--------- y.------ _-
 -43-

are typically the noisiest in the data set. This filtering capability
is not possible in the undiagonalized basis set.
The resolution analysis treated in Sec. I can be assily carried
out in the diagonalized basis. Denoting by (65(z)) the estimate of
60(z) produced by retaining k :<n terms in the spectral expansion
[Eq. (A.6)], we find

k
(60(z)> = w(z) aiH (z)
i=l1

z
m k
dz' 6a(z') w(z) H(z)H(z)

f dz' 6c(z')A(z, z') (A.14)

The averaging function A has the same form as the Fourier series
Dirichlet kernel discussed in Sec. It. Its width is a measure of
the maximum height range over which the conductivity profile could--
in principle--be deformed and still be consistent with the reflection-
coefficient data set. in practice, we have found that the height
resolution is better than indicated by this conservative criterion.
In general, retaining fewer terms in Eq. (A.14) has the effect
of increasing the width of the averaging function. Because this
process also decreases the variance in the conductivity [see Eq. (A.13)]
there is a tradeoff between high resolution and statistical reliability.
Our experience indicates that dropping the smallest elgenvalue terms
produces a significant improvement in estimated variance, while caus-
ing only a negligible degradation in height resolution. This behavior
is illustrated by the examples in Sec. III.
 MISSION
Of
Rom Air Development Center
RAVC Ptam1 a~nd exWLuteA tAewch, devetopnext, te~t &#d
setea~.ed aquiait pxoghuf in 6UWo'%t 0j Comuand, Contot
conui&-w~tona and inetgm ic 3ri aw u% TthacaZ t5
and engine~ 4uppo~t uVkin ate,". oj teaija copeteme
46. P'Lovided to EsP PuLgu o66ice, (poa) and otkeot ESP
eteinenta. The p'clnepat tedinZnizL ia.6ion #Aeah qAq
Woreflw9W oaaf, etee wmagnetu gcu~aue and eontwto, 6U&-
vegftance oj gJ4oud and da40 pace objecitA, ite ene dafta
cottee.tion and handting, in6ohmotion 60ateJR tecJuw gy.
iOno,6pheA&i p/Lopagation, 6otid 6&UL .enoceA, R&~AMVe
phy6Zeh and ete~tni.c xetiabitity, maiz nab.iUt ad

Poetud by s.Fw
Hm~COMM
 DAT

FILM

DTI