Geophysicae Amplitude domain analysis of strong range and Doppler spread radar echos

We present a novel method for analyzing range and Doppler spread targets in the amplitude domain using linear statistical inversion. The result of the analysis is an estimate of the range dependent amplitude behaviour of the target backscatter during the time that the transmission passes the target. A meteor head echo and strong backscatter from artificially heated regions of the ionosphere are used to demonstrate this novel analysis method. Plans to apply amplitude-domain radar target estimation methods to more complicated noisy underdetermined targets are also briefly discussed.


Introduction
Incoherent scatter radar targets are usually analyzed in the power domain using lag-profile or correlation based methods.For example, Virtanen et al. (2008a) discusses autocorrelation function estimation of range and Doppler spread ionospheric targets through statistical inversion.On the other hand, for targets that are spread only in range, matched filters or range sidelobe-free inverse filters have been used to analyze targets in the amplitude domain (e.g.Sulzer, 1989;Ruprecht, 1989).
Lag-profile analysis usually implies pre-defined integration times, range gates and lags to be estimated.These settings do not necessarily preserve all the information of the target.Also, lag-profile analysis inherently implies that the target backscatter is modelled as a stationary stochastic process -an assumption which is not always true.
While filter based amplitude domain decoding methods are fast and well proven, they are not suitable for all situ-Correspondence to: J. Vierinen (juha.vierinen@iki.fi)ations.For example, the matched filter suffers from range ambiguities and has an underlying assumption of a point-like target.The sidelobe-free inverse filter on the other hand does not have range ambiguity problems (Lehtinen et al., 2004;Vierinen et al., 2006), but just like the matched filter, there is an assumption that the target scattering coefficient (being defined as the ratio of target backscatter to the complex amplitude of the transmission) stays constant while the transmission pulse travels through the target.In reality this assumption is often violated.A good example is the F-region heating that is discussed later.
In this study we present a novel method for estimation of the target backscattering in the amplitude domain.To do this, we model the time evolution of the reflection amplitude for each range gate using a parametric model.For a wide target that is also Doppler spread, this results in a difficult underdetermined problem with many more parameters than measurements.But when the target is sufficiently narrow in both range and Doppler spread, the problem becomes an overdetermined linear statistical inverse problem which can be solved.We describe this analysis procedure and as an example we show how to get high spatial and temporal resolution amplitude estimates of narrow and strong radar targets, even with transmissions that are coded with bauds longer than the range resolution.The fundamental limit is set by the sample rate used to measure the echo.
The strong artificial ionospheric heating effects shown in this study were seen at the EISCAT Tromsø site on 18 October 2007 with O-mode heating during an experiment that was mainly intended for D-region studies.The heating was pointed in vertical direction with a 10 s on 10 s off modulation.The heater was operating at 5.4 MHz with an effective radiated power of 600 MW.Strong backscatter was often seen during the heater on period.The radar experiment was designed to also probe ranges up to 1100 km unambiguously by use of uneven inter pulse periods, which enabled us to also see strong heating effects in the F-region with three out of four echos.In addition to the strong F-region heating effects, we also saw a strong sporadic E-layer heating, although it was less frequent and often of much shorter duration.The heating effects were seen on both UHF and VHF radars.The short transmission pulse length of 150 µs, while necessary for D-region studies, prevented us from forming a high resolution spectrum of the target, but this could be easily remedied by using a longer transmission pulse in future experiments.

Amplitude model of an incoherent scatter target
Using discrete time and range, the direct theory for a signal measured from a radar receiver can be expressed as a sum of range lagged transmission envelopes multiplied by the target backscatter amplitude (1) Here m(t)∈C is the measured baseband signal (the notation ∈C means that the signal is complex-valued), (t)∈C is the transmission modulation envelope, ζ (r, t)∈C is the range and time dependent target scattering coefficient and ξ(t)∈C is measurement noise consisting of thermal noise and skynoise from cosmic radio sources.The measurement noise is assumed to be a zero mean complex Gausian white noise with variance E ξ(t i ) ξ(t j )=δ i,j σ 2 .Ranges r j are defined in round-trip time at one sample intervals and t i denotes time as samples.
There are many possible ways to model ζ (r, t).One possibility is to use a Fourier series in time, so our model parameters will consist of k terms of a Fourier series representation of the target scattering coefficient for each range of interest.This has the advantage that we can define the frequency characteristics that we expect to see in a target, as it is often the spectral properties that are of interest.Thus, we can express ζ (r, t) using coefficients c j,k ∈C of the series with frequency parameters ω k selected so that the frequency domain characteristics can be determined from the data.The backscatter amplitude of the target can thus be modelled using the parameter set θ={c j,k }, which has N r ×N f parameters, where N r is the number of ranges and N f is the number of elements in the Fourier series representation of the target amplitude.Thus, θ contains the parameters that we will attempt to infer based on the measurements.
We are left with a simple statistical parameter estimation problem, with parameters in set θ, which can be solved using statistical inversion.Using Eqs. ( 1) and (2), we can then write our direct theory z(t i , θ ) using the model as: We can write a likelihood function as a product of independent complex Gaussian densities, as our measurements are assumed to be distributed this way.Here D represents the set of measurements D={m(t 1 ), ..., m(t N )}: Normally, if the target range extent is wide, we would need many more parameters in θ than there are measurements.In this case it would be necessary either to use prior information or instead of backscatter coefficients, estimate the second order statistical properties of the target backscatter coefficients: r, t+τ ).This is what is done in traditional analysis using lagged product data m(t) m(t+τ ) to determine σ (r, τ ) without estimating ζ (r, t).
If we are interested in a narrow region only, as depicted in Fig. 1, we can leave out all parameters that are not from ranges that are interesting to us, assuming that the backscatter from these ranges merely adds to the measurement noise.If the range we are interested in has a very strong signal compared to the surrounding ranges, this is a good assumption to make.In this case, the problem becomes easy to solve as we have more measurements than model parameters.This study focuses on narrow strong targets that fulfill this criteria.

Numerical details
Assuming that we know the white noise variance σ 2 , our problem is a linear statistical inverse problem (Kaipio and Somersalo, 2004).We can find the maximum a posteriori parameters θ MAP1 using linear algebra if we write Eqs. ( 1) and (3) in the form where the measurements and parameters are vectors and the theory is expressed as a matrix.The measurement vector is m=[m(t 1 ), ..., m(t N )] T and the number of measurements N=N r +l−1 is a sum of target ranges and transmission envelope length l.The parameter vector is θ=[c 1,1 , c 1,2 , ..., c N r ,N f ] T , which has N r ×N f elements.Errors are uncorrelated so ξ ∼N(0, ), with =diag(σ2 , ..., σ 2 ).The theory matrix A can be expressed using Eq.(3).
To solve this problem efficiently, we used a software package called FLIPS 2 (Orispää and Lehtinen, 2008 3 ).The library uses QR-factorization via Givens rotations to solve the system of overdetermined linear equations.FLIPS can also be used to evaluate the posterior distribution of the parameters, which can be used to express errors associated with the parameters.

Example: F-region heating effect
During our 18 October 2007 daytime D-region heating experiment there was a sporadic E region during most of the experiment.In addition to this, we saw many strong O-mode heating related backscatter enhancements from the F-region The transmitted signal is in the beginning of the signal, followed by ground clutter and the ionospheric echo.At around 1500 µs one can see the F-region heating related echo.The time on the x-axis is in samples which are 0.5 µs long.The red and blue represent real and complex parts of the baseband signal.
and the sporadic E region on both VHF and UHF radars.By looking at the raw echos, shown in Fig. 2, it is evident that the heating effect was very strong and concentrated in a narrow region.By looking at the individual echos it was clear that the target was not completely coherent because the strong echo was not even close to an exact copy of the transmission pulse.An example of a transmission and the corresponding echo from the heated F-region is shown in Fig. 3.
To examine the amplitude of the F-region heating, we modelled 12 ranges 1 µs apart.Our coding, described in Virtanen et al. (2008b), used four 150 µ pulses with 10 µs bauds.The transmission envelope (t) was sampled directly from the waveguide.We modelled the range dependent amplitude using seven Fourier series parameters ω k 6.667 kHz apart within a ±20 kHz spectral area.The number of parameters was chosen so that the fit was good, while still giving residuals of correct magnitude.The signal was strong enough for us to be able to construct a decent estimate for each separate echo.Figure 4 shows the modulus of the parameters c j,k for each of the modelled ranges as a function of time during the 10 s heating period.This parameter plot can also be interpreted as a dynamic spectrum of the range dependent backscatter amplitude.The modelled backscatter amplitude at ranges 199.65-201.3km during the first 100 ms of heating is shown in Fig. 5 199.65 199.8 199.95 200.1 200.25 200.4 200.55 200.7 200.85 201 201.15 201.3 Fig. 4. Heating in the F-region at 150 m resolution from one heating period starting from 10:45:20.The temporal resolution is approximately 2.5 ms (uneven IPPs).The figure contains the modulus of one set of spectral parameters for each transmit pulse, the values are in linear scale.The values are the statistically most probable values given the measurements.Each range gate is represented with a ±20 kHz spectrum at a 6.67 kHz frequency resolution.The spectrum is dominated by one central peak.The heated layer is completely contained within a 1.8 km range interval and most of it is within a 600 m region.After recovering from the strong overshoot in the beginning, the heated region moves down at about 45 m/s during a single 10 s heating period.The spectrum seems to broaden and strengthen slightly towards the end of the heating period.
The results are similar to the ones obtained by Djuth et al. (2004), except that we have slightly worse frequency resolution due to the shorter transmission pulse.But we are able to obtain much better temporal resolution.During this experiment, we did not record plasma lines, but this same method is applicable for analyzing them, provided the plasma line bands are sampled.

Example: sporadic E-layer heating effect
The sporadic E region heating effect was analyzed using the same Fourier series parameters ω k as the F-region heating in the previous section.The combined results of one 10 s heater on period are shown in Fig. 6.Compared to the F-region heating effect, both the spectrum and the layer itself is very narrow.The spectrum is so narrow that only the zero frequency component of the Fourier series has significant power -the change in amplitude is only apparent when inspecting  The estimate is based on the statistically most probable value of parameter vector θ given the measurements.The red and blue lines represent the real and complex parts of the baseband signal in linear scale.This F-region heating event is the same as the one in Fig. 4. The amplitude is modelled for 150 µs, which is the time that the transmission pulse travels through the range gate.Discontinuities in the figure are greater than they appear in, they are determined by the inter-pulse period, which is approximately 2.5 ms in this case.
the signal on a pulse to pulse basis, where the slow changes begin to appear.
The results from the sporadic E layer heating show variation in backscatter power during the on-period.The heating effect is mostly contained in one 150 m range gate, with a weak signal in the neighboring gates in the beginning of heating (the measurements could not be explained with only one range without causing worse residuals, which is an indication that these additional ranges are needed in the model).There is certain similarity to heating effects reported by Rietveld et al. (2002), with the exception that the ion-line spectrum obtained here is very narrow, less than 10 Hz.
In this case, the amplitude mostly contained very slow changes and one can easily see the main Doppler shift of −8 Hz by inspecting the estimated amplitude data.During the first echo received after heating on, there is a strong overshoot, which is not there any more during the next echo.In addition to this there were at least three detectable harmonics of 50 Hz, with 50 Hz the strongest of them, only approximately 10 dB lower than the main peak centered at −8 Hz.It is unclear what causes these harmonics, but we have ruled out the EISCAT VHF transmitter by inspecting the transmitter envelope sampled from the wave guide.The receiver chain also seems to be free of any of these components, as e.g. the ground clutter does not contain any of these compo-nents.Two feasible alternatives could be the heater RF or direct power transmission line modulation of the sporadic E region in the ionosphere.

Example: meteor echo
Meteor head echos are also one example of strong point-like radar targets.Two meteor head echos are shown in Fig. 2 below the F-region heating effect.Meteor head echos are routinely measured with high power large aperature radars such as EISCAT or Arecibo radars (e.g.Mathews et al., 1997;Pellinen-Wannberg, 2005).These measurements are usually modelled with a delayed transmission envelope multiplied by a complex sinusoid The meteor velocity and range are then determined by finding the best fitting parameters ζ , r and ω.This is actually a good model, but it cannot describe arbitrary amplitude behaviour, and moreover it cannot be used to model range dependence very well.Typically, there is an underlying assumption of a point-like target, which results in range ambiguities for a spread target.Fig. 6.Backscatter amplitude estimate of the sporadic E-region heating effect.The time is relative to the start of the 10 s O-mode heater on period.The figure above shows the modulus of the seven ±20 kHz Fourier series coefficients c j,k used to estimate the amplitude of each range gate during the 150 µs that the transmission pulse passes each range.The heating effect is mostly concentrated in only one range gate, with slight hints of power on the neighboring range gates, which cannot be explained by a model with only one range gate.The figure on the lower left depicts the amplitude behaviour of the first 200 ms of heating at 107.55 km.Blue is real and red is the imaginary part of the signal, the black line is the modulus of amplitude.The first echo is stronger and at a different phase than the rest of the backscattered waveform.On the lower right is the low frequency spectrum of the reflection amplitude from 107.55 km estimated over the whole 10 s heater on period.The main Doppler shift is centered at around −8 Hz.The blue vertical lines depict 50 Hz harmonics shifted by −8 Hz.It is unclear why the 50 Hz harmonics are in the received backscatter signal -but it does not seem to be caused by the EISCAT VHF TX or RX receive path.Possibilities include heater modulation or a direct modulation by ground based power transmission lines.
To demonstrate amplitude domain analysis of meteor head echos, we modelled ζ (r j , t i ) at 9 ranges using 9 Fourier series coefficients centered around 40 kHz, which was approximately the Doppler shift of the meteor head echo.The raw voltage data was sampled at 8 MHz bandwidth.The modulus of the coefficients c j,k for one meteor head echo is shown in Fig. 7.The code length was 104 µs with 2 µs bauds.The backscatter amplitude is concentrated in a 100 m region with a backscatter magnitude decreasing with range.This could be a signature of the quickly vanishing trailing edge of the meteor head echo, but a more rigorous analysis would be required to verify this.

Discussion
We have demonstrated a method that gives very good temporal and spatial resolution for decoding strong sufficiently narrow targets.The method works with many types of radar transmissions, and can thus be run as a secondary analysis for situations where strong echos are observed.The method, although very promising, is still new and thus there remains work to be done with testing, parametrization, estimation errors, transmission code optimality, and numerical solution methods.
In this study we used a Fourier series to model the target, as it was the most straightforward one and it resulted in a linear model.Because of the small number of parameters in the series, there will certainly be some artifacts caused by this parametrization.The most notable one is that the amplitude behaviour tends to be periodic at the ends of the estimation interval, which is visible in Fig. 5.In cases where the target backscatter amplitude is sufficiently narrow band, a better parametrization for target backscatter amplitude would more likely be a complex sinusoid multiplied by a cubic spline.This would also be more suitable for meteor head echos, as this would allow more precise determination of the Doppler shift.This approach results in a non-linear statistical inverse problem which can be solved, e.g. by using MCMC (Hastings, 1970).
In addition to the examples presented in this study, there are also many other possible applications for this method.In the case of strong targets, our method will be directly Ann.Geophys., 26, 2419Geophys., 26, -2426Geophys., 26, , 2008 www.ann-geophys.net/26/2419/2008/q q q q q q q q q kHz 0 50 100.22km q q q q q q q q q kHz Mod(res$spec[((r − 1) * nfreq + 1):(nfreq * r)]) 0 50 100.24km q q q q q q q q q kHz Mod(res$spec[((r − 1) * nfreq + 1):(nfreq * r)]) 0 50 100.26km q q q q q q q q q kHz Mod(res$spec[((r − 1) * nfreq + 1):(nfreq * r)]) 0 50 100.28km q q q q q q q q q kHz Mod(res$spec[((r − 1) * nfreq + 1):(nfreq * r)]) 0 50 100.29 km q q q q q q q q q kHz Mod(res$spec[((r − 1) * nfreq + 1):(nfreq * r)]) 0 50 100.31km q q q q q q q q q kHz Mod(res$spec[((r − 1) * nfreq + 1):(nfreq * r)]) 0 50 100.33km q q q q q q q q q kHz Mod(res$spec[((r − 1) * nfreq + 1):(nfreq * r)]) 0 50 100.35km q q q q q q q q q kHz Mod(res$spec[((r − 1) * nfreq + 1):(nfreq * r)]) 0 50 100.37km Fig. 7.The modulus of the Fourier series coefficients c j,k .The meteor head echo is about 100 m wide and the main 40 kHz Doppler shift dominates the model.applicable.For weak targets it is yet unknown how our method will perform, and it is a topic of future work.For example, near earth asteroids are an example of narrow radar targets that are fairly weak.Range-Doppler measurements of near earth asteroids (Hudson, 1993) are routinely used to determine the shape of the target.The range-amplitude measurement presented in this paper would offer more information than a traditional power domain range-Doppler estimate.The reason for this is that it is possible the reduce the rangeamplitude estimate into a power domain range-Doppler estimate by simply taking the modulus of the Fourier domain representation of the amplitude, but this would mean discarding the phase information.
In this paper, we have not covered transmission code optimality in terms of amplitude domain inversion.In order to optimize a transmission code for a certain kind of target, one needs to calculate the covariance matrix for the parameter vector θ.In matrix form this is where A is the theory matrix from Eq. ( 5) and σ 2 is the measurement noise variance.This matrix contains the transmission envelope (t) and the Fourier series terms e iω k t i .Here A H is the complex conjugated transpose of the theory matrix A. There are several aspects of the covariance matrix that one can optimize, but in general the errors of the parameters should be small and as independent as possible.This leads to several different code optimization criteria, such as minimization of the determinant of the error covariance matrix.Code optimization is a topic of future work.Also, we have not yet visualized the estimation errors properly with the results, although this is pretty straightforward to do, as the problem gives a well defined Gaussian posterior covariance.This will be important when infering physical parameters from amplitude domain estimates.
It would be interesting to repeat the experiments shown in the examples with a longer transmit pulse and a higher sample rate to achieve better frequency and height resolution.The plasma lines should also be measured and analyzed using the method described here.Preferably the data should be sampled at a large enough rate to fit the whole signal.
In this work we have used a discrete time and range model.Some improvements in estimation accuracy can be expected if the model would include proper range ambiguities that also take into account the impulse response of the receiver chain.
While we have only applied this method to strong overdetermined targets, there might also be a possibility to extend this method to analyze underdetermined and weak incoherent scatter targets.It is not yet completely clear how this would be carried out.However, we plan to develop methods for a calculus of singular distributions for the target scattering coefficients, which could then be used in a further step of analysis modeling the scattering autocorrelation function as an unknown instead.

Fig. 1 .
Fig. 1.Simplified range-time diagram of backscatter from a strong narrow region.In this example there are two transmit samples and three ranges that cause backscatter.The red lines visualize the changing amplitude of backscatter at each range.The gray area represents the area where the backscatter of one sample originates from (assuming boxcar impulse response).

Fig. 2 .
Fig.2.Modulus of raw VHF measurements from a strong narrow layer in the F-region.Two meteor head echos can also be seen below the F-region backscatter.The strong echos below are ground clutter echos.Origin of time is the end of the TX pulse.

Fig. 3 .
Fig.3.Example transmission and echo from a point-like heating effect in the F-region.In this case from the EISCAT UHF signal.The transmitted signal is in the beginning of the signal, followed by ground clutter and the ionospheric echo.At around 1500 µs one can see the F-region heating related echo.The time on the x-axis is in samples which are 0.5 µs long.The red and blue represent real and complex parts of the baseband signal.

Fig. 5 .
Fig.5.Estimated backscatter amplitude from three ranges during the first 100 ms of heating.The estimate is based on the statistically most probable value of parameter vector θ given the measurements.The red and blue lines represent the real and complex parts of the baseband signal in linear scale.This F-region heating event is the same as the one in Fig.4.The amplitude is modelled for 150 µs, which is the time that the transmission pulse travels through the range gate.Discontinuities in the figure are greater than they appear in, they are determined by the inter-pulse period, which is approximately 2.5 ms in this case.