Effects of mixed scatter on SuperDARN convection maps

The SuperDARN radars map high-latitude ionospheric plasma drift velocities by measuring the Doppler frequency shift of HF signals scattered by decameter electron density irregularities. In many cases the ionospheric returns are contaminated by strong scatter from the ground or sea surface. In this paper we develop and test a two-component fitting algorithm to separate ionospheric and surface scatter components. Application of the technique to sample data reveals that mixed scatter may considerably distort ionospheric convection patterns derived from the radar data and can cause underestimation of the plasma drift velocity.


Introduction
The Super Dual Auroral Radar Network (SuperDARN) consists of pairs of HF radars mapping ionospheric plasma convection over the polar and auroral regions in both hemispheres (Greenwald et al., 1995).To achieve the required ranges (D 3500 km) and Doppler velocities (|V |∼2 km/s), the radars were designed to measure the complex autocorrelation function (ACF) of the returns using a sequence of non-evenly separated pulses (Greenwald et al., 1985;Baker et al., 1995;Ponomarenko and Waters, 2006).In FITACF, the SuperDARN data processing package, velocity V , spectral width W and other characteristics are estimated by fitting analytical functions to the experimental ACF phase (τ ) and power |R(τ )| assuming a single spectral component in the radar echoes (Baker et al., 1995;Ponomarenko and Waters, 2006).The ionospheric plasma convection maps are obtained by fitting a model based on a spherical harmonic expansion Correspondence to: P. V. Ponomarenko (phpp@alinga.newcastle.edu.au) to the velocity data from all available radars (Ruohoniemi and Baker, 1998).
There are two major contributors to radar returns: (i) ionospheric turbulence producing so-called "ionospheric scatter" (IS) and (ii) scatter from the ground or sea surface (surface scatter, SS).The ionospheric component is usually characterised by relatively large values of |V |∼100-1000 m/s and W ≥100 m/s in contrast to |V |, W ≤30-40 m/s typically for SS.While "pure" SS is routinely detected and eliminated using empirical thresholds for |V | and W (G. T. Blanchard and K. Baker, private communication, for details see Ponomarenko et al., 2007), the so-called mixed scatter (MS) returns containing both types of echoes are not easily identified.Application of the conventional FITACF procedures in this case leads to incorrect estimates of V , W and other parameters producing large fitting errors.Chisham and Pinnock (2002) showed that contamination from SS can considerably alter mesoscale features of the fitted convection pattern and developed a technique for eliminating SS and MS from range-velocity distribution of echoes based on empirical identification of HF propagation modes.A similar approach was applied by Ponomarenko et al. (2007), who successfully removed most MS from the data by applying the Blanchard-Baker algorithm with extended empirical thresholds for W and |V |.However, both methods reject potentially useful MS data, which contain information about ionospheric returns.Barthes et al. (1998) demonstrated the possibility of resolving an arbitrary number of modes in multi-component SuperDARN ACFs using the multiple signal classification (MUSIC) method.However, their approach is rather complex and requires extensive additional calculations.In this work we develop and test a simple algorithm for identifying MS returns and deconvolving IS and SS components using a two-component complex fitting function, which can be relatively easily incorporated into the existing SuperDARN software.P. V. Ponomarenko et al.: Mixed scatter in SuperDARN radars Fig. 1.Range-time map of line-of-sight velocity, V , from beam 13 of the TIGER radar (Tasmania) at 06:00-08:00 UT, 1 March 2000 estimated using conventional FITACF package (top) and new twocomponent complex fit (bottom) accompanied by respective analytical ray-tracing for the same distance range (middle).Light green and orange shading in the middle panel show areas of potentially intensive surface and ionospheric scatter, respectively.White rectangles and arrows in the velocity panels correspond to the mixedscatter ACF presented in Fig. 2.

Example of mixed scatter data
To illustrate the effects of mixed scatter, in Fig. 1 we show a range-time map of V for beam 13 of the TIGER (Tasmania) SuperDARN radar during 06:00-08:00 UT on 1 March 2000.The radar was operating in common mode consecutively scanning all 16 beams with the integration time 7 s/beam (number of averaged pulse sequences, N a 70) at the working frequency f 0 12.2 MHz.We present data only with signal-to-noise ratio (SNR) ≥6 dB.The measured ve- locity (top panel) exhibits two distinct components with relatively low (olive green, range D 1000-1700 km) and high (orange-yellow, D 1300-2000 km) values of V , which apparently overlap at 07:00-07:30 UT.
The middle panel shows results of analytical ray tracing covering the same distance range and based on expressions for a parabolic ionospheric layer and spherical geometry (Kelso, 1964, p. 254-257).The layer parameters were provided by the IRI model 1 for the interval shown in the top panel, and we plotted only those rays reflected by the ionosphere.The multi-hop HF propagation results in periodic spatial variations of the wave field strength due to focusing effects.The "bright spots" are located at D D skip ×n on the ground (green shading) and D D skip ×(n−0.5) in the ionosphere (pink shading), where D skip is the skip zone distance, and n is hop number.Ionospheric and ground regions partially overlap in the range (time delay) domain so that some SS and IS echoes may arrive at the receiver at the same times producing mixed scatter returns.By comparing the ray-tracing results with the data from the top panel it is reasonable to assume that most of the low-velocity echoes from D 1000-1700 km belong to the single-hop SS mode, while the majority of the high-speed returns from D 1500-1800 km may be interpreted as 1.5-hop IS returns.Some signatures of 0.5, 2, and 2.5-hop modes can also be seen after 07:15 UT.
In Fig. 2 we show an example of mixed scatter ACF power and phase for the time and range identified by the white rectangle and arrow in Fig. 1.Both |R(τ )| and (τ ) exhibit pronounced oscillations so that the FITACF single-component fit (dashed lines) becomes meaningless producing large errors in estimating V , W and other echo parameters and illustrating the need for a more sophisticated fitting function.Importantly, this ACF was identified by FITACF as "ground scatter" and would normally be excluded from further analysis.

Separation algorithm
Following Barthes et al. (1998), we separate different propagation modes in MS by analytically representing a mixedscatter ACF as a sum of exponentially decaying harmonic oscillations but using only two terms describing IS and SS: where τ is time lag, R i,s (0), τ i,s and f i,s =1/T i,s are ACF magnitude, decay time and frequency shift, respectively, indices i and s depict IS and SS, respectively, and j = √ −1.The above parameters relate to SuperDARN velocity and spectral width as V i,s =λf i,s /2 and W i,s =λ/(2πτ i,s ), where λ is radar wavelength, and W is measured at the half-power level and expressed in ms −1 .
To avoid dealing with 2π phase skips, the new algorithm uses a single procedure to directly fit the model in Eq. (1) to a complex ACF.
For actual fitting we used a free IDL package MPFIT2 , which applies the Levenberg-Marquardt technique to solve the least-squares problem.Initial values for varying the fitting parameters were: R i (0)=R s (0)=0.5R(0),τ i =20 ms, T i =30 ms, T s =∞ (V s =0).To avoid "swapping" IS and SS or obtaining two components with essentially the same τ and T , the SS decay time was fixed at a large value τ s =100 s.No other limitations have been imposed on the input or output parameters.
Prior to fitting experimental data, the new technique was thoroughly tested against simulated ACFs emulating a sum of IS, SS and background noise.All three components were generated using the SuperDARN mlutipulse algorithm (Baker et al., 1995;Ponomarenko and Waters, 2006) applied to the following simulated timeseries: -IS: sum of signals scattered by a large number of irregularities randomly distributed within a range gate and characterized by exponential decay time τ i (Lorentzian ACF) and background velocity V i ; Applying different weighting to the above components simulates the relative SS/IS contributions and SNR.We used the same number of averaged pulse sequences as in Fig. 2, N=78, so that the resulting statistical variations resembled those of the experimental data analyzed in this paper.To avoid distortion of the shape of the ACF power around τ =0, prior to fitting the background noise power was extracted from R(0), so that this component only contributed to the statistical fluctuation level (for more details see Ponomarenko and Waters, 2006).
Intuitively, one would expect the best results from the new algorithm when it is applied to a mixed scatter ACF with a pronounced interference pattern similar to that presented in Fig. 2. The IS and SS powers are comparable, R s (0) R i (0), and there is at least one full oscillation in the interference  pattern over the duration of the ACF, T i ≤ min[τ max , τ i ] (here τ max 44 ms is the maximum ACF time lag).Deviation from these conditions should lead to a less reliable fit.These effects are illustrated in Fig. 3 and 4, which show simulated examples of "good" and "bad" ACFs obtained for the SS/IS power ratio ρ=R s (0)/R i (0)=2 and SNR=0.The "good" ACF covers 1 1 2 oscillation periods, and in this case the new algorithm comprehensively deconvolves IS and SS components and reliably reproduces the input parameters for both components, while the conventional technique fails.For the "bad" ACF covering only 1 4 of the oscillation period, while not performing as good as for V i =500 m/s, the new algorithm still provides a better estimate of IS parameters compared with FITACF.
To establish statistically reliable limitations of the new technique, we applied it to a number of simulated datasets with different values of V i and ρ.Each dataset contained n=1000 ACFs with constant values of τ i =30 ms, V s =10.33 m/s and SNR=6 dB, which is the lowest SNR value used in our experimental data analysis (here the "signal" presents a sum of IS and SS).Results of the statistical analysis are presented in Fig. 5, where we show dependence of V i (a) and R i (0) (b) on ρ.The dashed lines in both pan- els represent simulated (input) parameters while the symbols depict median values for the fitting results (output), and the error bars correspond to one standard deviation.Median output velocity and power follow closely the input values for ρ≤10−20.The velocity error bars increase with decreasing V i , but the algorithm is reliable even for the lowest velocity magnitude, V i =30 m/s.For ρ≥10−20 the fitting results become unreliable exhibiting strong bias from the input values and unacceptably large fluctuations.
In order to correctly apply Eq. ( 1) it is necessary to design a procedure for detecting the presence or absence of MS in a single ACF.Previously, André et al. (2002) proposed identification of multi-component ionospheric echoes recognising that a single-component model (FITACF) applied to multi-component ACFs generally produces larger fitting errors, δ fit , than those obtained for single-mode echoes.They applied empirical thresholds for the standard deviation of the experimental data from the fitted functions for both power, δ fit R /R(0)≥0.15, and phase, δ fit ≥0.3 rad.The disadvantage of using phase error δ fit is that it relies on correctly accounting for 2π skips.Using δ fit R seems to be more reliable, but requires a more rigorous approach when choosing the cut-off value.
The natural candidate for a physically justified cutoff is a statistical fluctuation level for ACF power, σ R .If both R i (0) and R s (0) exceed σ R , then MS should be considered present in this particular ACF.The theoretical estimate for the fluctuation level, σ R =R(0)/ √ N, is used in FITACF to reject ACF lags with low power (Ponomarenko and Waters, 2006).However, MS presents a combination of a quickly varying IS component (τ i ≤τ max ) and an essentially regular SS (τ s τ max ).In this case the statistical fluctuations are determined predominantly by IS, and their relative magnitude decreases with increasing contribution from SS ∝1/(1+ρ).Also, for large values of ρ the contribution from the background noise becomes important.To account for all possible sources of fluctuations it is more practical to use fitting errors defined as , where R(τ ) is the experimental complex ACF, R fit (τ ) is the fitted model and . . .means averaging over τ .We assume that the IS or SS component is present in a particular experimental ACF if the respective lag 0 power exceeds the fitting error by at least 100%, i.e.R i (0), R s (0)≥2δ fit R .To illustrate the validity of the selection criterion, Fig. 5b shows the dependence of 2δ fit R on ρ (red solid lines).The dash-dotted red line corresponds to the "conventional" statistical fluctuation level, σ R .Notably, the fitting error exhibits little variation with changing V i .For small values of ρ the MS threshold is close in magnitude to σ R , but with increasing contributions from SS it decreases until ρ 8−10 and then starts to saturate because R i (0) becomes comparable to the external noise (−6 dB).Increasing SNR leads to respective increase in the effective threshold in ρ, which is illustrated by the dotted red line showing average 2δ fit R for SNR=10.Therefore, according to our criterion, in Fig. 5b the IS component is considered to be present in the data wherever the power curves are located above the error curves.Comparison with the velocity data in Fig. 5a shows that our algorithm effectively rejects unreliable data with large velocity and power fluctuations.

Fitting experimental data
First, we applied the new technique to the example shown in Fig. 2. Resembling the simulation results in Fig. 3, the new two-component fitting function (dash-dot) closely reproduces the experimental power and phase (diamonds), in contrast to the single-component FITACF (dashed lines).The two-component analysis reveals that in this case the ionospheric echoes characterized by high V i 300 m/s and W 200 m/s values are effectively masked by SS, and the ACF is mis-identified by FITACF as surface scatter with low V s 3.5 m/s and W 50 m/s.
The bottom panel in Fig. 1 shows ionospheric velocities V i calculated using the new fitting technique.While the new processing algorithm closely reproduces major features of the ionospheric scatter obtained by the FITACF (top panel), it also reveals additional details hidden by SS, such as a highspeed patch centered at D 1500 km during 07:00-07:30 UT.
Next, we applied the new fitting procedure to a full 16beam scan 07:12-07:14 UT corresponding to the white rectangle in Fig. 1. Figure 6 shows spatial maps of line-of-sight velocities calculated by FITACF (left) and using the ionospheric component of the new fit (right).Both plots exhibit a clear change from V 1000 m/s in the eastern part of the radar's field of view to V −600 m/s in the western part, which is a characteristic response of the radar's lineof-sight velocity to the regularly observed westward flow of the ionospheric plasma.However, in the left panel at D≤1600 km this pattern is almost completely "blanketed" by the 1-hop SS component (olive green), while the secondhop SS echoes are affecting scatter from D≥2300 km.In contrast, the two-component fit effectively separates radar echoes into IS (right panel) and SS (not shown) components and expands the effective spatial coverage revealing further details of the V i distribution over the "blanketed" distance range, D 1300−1600 km.
Figure 7 shows respective convection patterns and electric potential contours obtained by fitting the conventional APL model to all available radar data from the Southern Hemisphere, for the interval considered in Fig. 2.Only three radars, TIGER, Syowa East and Sanae, had sufficient scatter during that period.In contrast to TIGER, the data from Sanae and Syowa East had virtually no SS component, and their reprocessing with the new technique did not show much difference in spatial velocity distribution (not shown) compared with FITACF.When calculating the convection maps, we applied the spherical harmonic expansion of 8th order (Ruohoniemi and Baker, 1998), while the IMF components obtained from the ACE satellite at the second Lagrangian point were B x =6.2 nT, B y =3.9 nT, B z =−5.1 nT accounting for the propagation time 48 min.
As expected, the new algorithm produces an extended coverage clearly seen over the TIGER field-of-view on the lefthand side of the map.It clarifies the latitudinal extent of the high-speed flow around 70 • MLAT.While the total cross-cap potential remains virtually unchanged, the new procedure results in substantially higher local fitted velocities and a noticeable change in the configuration of equipotential lines at the mesoscale level, in agreement with Chisham and Pinnock (2002).

Summary and conclusions
Based on numerical model studies, we have developed physically justified selection criteria for detecting SuperDARN ACFs containing both ionospheric and ground/sea scatter modes (mixed scatter).These criteria were implemented in a new two-component fitting algorithm, which allows us to successfully deconvolve ionospheric and surface scatter modes and to obtain reliable ionospheric velocity estimates from the mixed returns.
Applying the new technique to sample data from TIGER SuperDARN radar, we have illustrated the effects of mixed scatter on radar echoes.When sea/ground scatter is mixed with the ionospheric component, it may "blanket" ionospheric returns.Processing the data using the conventional single-component algorithm may considerably decrease the effective radar coverage because the mixed-scatter echoes are either misidentified as SS and excluded from analysis, or V and W attributed to the ionospheric component are generally underestimated.Based on real data we show that the latter may result in distortion of the convection pattern at mesoscale level and systematically lower convection plasma velocities leading to underestimation of electric potential gradient magnitudes.
The new technique may also be applied to deconvolving multiple IS modes, provided that a suitable set of modified selection/separation criteria is designed.

Fig. 3 .-
Fig. 3. Simulated example of a "good" ACF (open diamonds) covering 1 12 periods of the interference pattern between IS and SS.The dashed line corresponds to the FITACF fit, while the dash-dot line shows the two-component fit.

Fig. 4 .
Fig. 4. Simulated example of a "bad" model ACF covering 1 4 of the interference period presented in the same format as Fig. 3.

Fig. 5 .
Fig. 5. (a) Dependence of median IS velocity estimate by twocomponent model on IS/SS power ratio.(b) Same dependence for IS power and doubled fitting error illustrating the selection criterion developed in this paper.The dashed lines in both panels represent simulated (input) parameters while the symbols depict median fitting results, and the error bars correspond to one standard deviation.

Fig. 6 .
Fig. 6.Spatial distributions of line-of-sight velocity for 07:12-07:14 UT 1 March 2000.The left panel shows results obtained for the conventional single-component fitting procedure.The right panel represent the same parameter but estimated for the ionospheric component of the two-component fitting.

Fig. 7 .
Fig. 7. Fitted model plasma convection velocities (coloured vectors) and electric field potential contours over the Southern Hemisphere at 07:12-07:14 UT 1 March 2000.As in Fig. 6, the left panel shows results obtained using the conventional FITACF procedure, while the right panel represent the results based on the new technique.