Wavelet Based Methods for Improved Wind Profiler Signal Processing

In this paper, we apply wavelet thresholding for removing automatically ground and intermittent clutter (air- plane echoes) from wind profiler radar data. Using the con- cept of discrete multi-resolution analysis and non-parametric estimation theory, we develop wavelet domain thresholding rules, which allow us to identify the coefficients relevant for clutter and to suppress them in order to obtain filtered recon- structions.


Introduction
Radar Wind Profilers (RWP) are versatile tools used to routinely probe the Earth's atmosphere.This technology originally developed for studying the dynamics of the middle atmosphere in the seventies (Hardy and Gage, 1990) is, meanwhile, very prominent in the meteorological research community.Meteorological services started using these systems operationally within the Global Observing System (GOS) (see Monna and Chadwick (1998)).
Most of these RWP employ the Doppler-beam swinging (DBS) method for the determination of the vertical profile of the horizontal wind and, under certain conditions, the vertical wind component.These radars transmit short electromagnetic pulses in a fixed beam direction and sample the small fraction of the electromagnetic field backscattered to the antenna.At least three linear independent beam directions are required to transform the measured 'line-of-sight' radial velocities into the wind vector.Due to the nature of the acting atmospheric scattering processes, the received signal is several orders of magnitude weaker than the transmitted signal.The received signal is Doppler shifted, which is used to determine the velocity component of "the atmosphere" projected onto the beam direction.As the bandwidth B of Correspondence to: V. Lehmann (volker.lehmann@dwd.de)a transmitted electromagnetic pulse of duration τ is much larger (B ∝ 1/τ ≈ 100...1000 kHz) than the Doppler shift (f d ≈ 10...500 Hz), the frequency shift cannot be determined from the processing of a single pulse.Instead, the return of many pulses is evaluated to compute the Doppler frequency from the slowly changing phase of the received signals (Burgess and Ray, 1986).Sampling is done after the receivers quadrature detector (for the in-phase and quadraturephase components of the signal) using sample and hold circuits prior to the A/D conversion.The sampling rate is determined by the pulse repetition period T .The samples at each range gate form a discrete complex time series, which is the raw data of the measurement at this gate.The following digital signal processing has the purpose of extracting the desired atmospheric information from the radar echoes.More details about coherent radar technology and in particular, wind profilers, can be found in standard textbooks (Gossard and Strauch, 1983;Doviak and Zrnić, 1993) and in several review papers, e.g.Röttger and Larsen (1990).
In this paper, we propose a modified signal processing technique for RWPs.It must be noted that signal processing includes all operations that are performed on the radar signal, i.e. analog1 as well as digital processing2 .However, in the following, we will only concentrate on digital signal processing.The incredible development of fast digital processors opens up new opportunities to optimize this latter part of the signal processing chain.The goals of signal processing, as summarized by Keeler and Passarelli (1990), are: -to provide accurate, unbiased estimates of the characteristics of the desired atmospheric echoes; -to estimate the confidence/accuracy of the measurement; -to mitigate effects of interfering signals; -to reduce the data rate.The fundamental base parameters of the atmospheric signal are the reflected power, the radial velocity and the velocity variance (e.g. the first three moments of the Doppler spectrum).Signal processing ends with the estimation of the moments of the Doppler spectrum and further data processing is then performed to finally determine the wind and other meteorological parameters using measurements from all radar beams.This distinction, which goes back originally to Keeler and Passarelli (1990), has become more and more blurred, since some modern algorithms make use of the moments of the Doppler spectrum with the help of continuity and other information (Wilfong et al., 1999b).However, we will refer here to the usually applied and well established "classical" signal processing, as described by Tsuda (1989), Röttger and Larsen (1990), among others.

Statement of the problem
Before we discuss the problems that are associated with the "classical" processing, let us briefly repeat the steps as visualized in Fig. 1.In particular, we refer to the signal processing as it is implemented in the RWP, whose data are used in this study.
Digital signal processing in a system using an analog receiver3 starts with the sampling of the in-and quadraturephase components of the received signal at a rate that is determined by the pulse repetition period T .To reduce the data rate for further processing, hardware adder circuits perform a so-called coherent integration (Barth et al., 1994;Carter et al., 1995;Wilfong et al., 1999a), adding some N (typically ten to hundred) complex samples together.Mathematically, this operation can be seen as a combination of a digital boxcar filtering, followed by an undersampling at a rate of N T (Schmidt et al., 1979;Farley, 1985).If the radar system uses pulse compression techniques (e.g.phase coding using complementary sequences), then the next step is decoding (Schmidt et al., 1979;Sulzer and Woodman, 1984;Farley, 1985;Ghebrebrhan and Crochet, 1992;Spano and Ghebrebrhan, 1996).The coherently averaged and decoded samples are then used to compute the Doppler spectrum using the Windowed Fourier Transform (FFT) and the Periodogram method (see Keeler and Passarelli, 1990).In our system, a Fourier transformed Hanning-window is convolved with the result of the FFT.A number (typically some ten) of individual Doppler spectra is then incoherently averaged to improve the detectability of the signal (Tsuda, 1989, see).Finally, the noise level is estimated with the method proposed by Hildebrand and Sekhon (1974), and the moments of the maximum signal in the spectrum are computed over the range where the signal is above the noise level (May and Strauch, 1989).
The problem with this type of signal processing is the underlying assumption that the signal consists of only two parts: the signal, that is produced by one atmospheric scattering process, and noise (different sources, mainly thermal electronic noise and cosmic noise).This is certainly not true, especially at UHF, where the desired atmospheric signal itself is often the result of two distinct scattering processes, namely scattering at inhomogenities of the refractive index (Bragg scattering) and scattering at particles, such as droplets or ice crystals (Rayleigh scattering) (see, for instance, Gossard, 1979;Gossard andStrauch, 1981, 1983;Ralph et al., 1995Ralph et al., , 1996;;Gage et al., 1999).Therefore, even the desired atmospheric signal may have different characteristics.But, as experience shows us, the most serious problems are caused by the following contributions to the signal: Ground Clutter.Echo returns from the ground surrounding the site, which emerge from antenna's sidelobes; Intermittent Clutter.Returns from unwanted targets, such as airplanes or birds, from both the antenna'ss main lobe and the sidelobes; Radio Frequency Interference (RFI).RFI can emerge from external radio-frequency transmissions within the passband of the receiver (matched filter), or it can be generated internally due to imperfections of the radar hardware.
Recently, much work has and continues to be done to develop frequency domain processing algorithms, i.e. to improve the process of moment estimation.The purpose of these methods is to select the "true" atmospheric signal in the Doppler spectrum even in the presence of severe contamination.Only this signal will then be used for the determination of the wind vector.Several criteria are used to make an "intelligent" selection of the signal (Clothiaux et al., 1994;Gossard, 1997;Griesser, 1998;Cornman et al., 1998;Schumann et al., 1999;Wilfong et al., 1999b;Morse et al., 2000).The emphasis on frequency domain processing was probably caused by the fact that it is much easier to handle spectral data, as the data volume is significantly reduced due to the data compression effect of the periodogram computation and the spectral integration.Some of these "multiple moment estimation" algorithms additionally assign a quality indicator to the computed wind values, which does not only depends on the quality of the moment estimation, but also uses continuity criteria (Wilfong et al., 1999b) and the testing of assumptions that are inherent to the DBS method (Goodrich et al., 2000).First evaluations have indeed shown a very promising improvement of those new algorithms (Cohn et al., 2000), but no long-term evaluation against independent measurement systems, such as the Rawinsonde, has been performed so far.
Modified time domain processing has been proposed to reduce the problems caused by contaminating signals.One problem emerges from the fact that the receiver filter of the radar is matched to the transmitted pulse in order to optimize the single pulse signal-to-noise ratio for improved signal detection in the presence of noise (Tsuda, 1989;Papoulis, 1991;Doviak and Zrnić, 1993).This implies a receiver bandwidth of B ∝ 1/τ .Yet, the sampling is done at a rate of T (without coherent averaging) or even N T (with coherent averaging).Thus, the Nyquist frequency, after coherent averaging, is severely smaller than the frequency that a received signal might have.Of course, it is true that the desired atmospheric signal is band-limited by a sufficiently long coherence time of the scattering process, so that this undersampling has no consequences (aside from some modification due to the filtering characteristics of the coherent integration).If there is, however, some artificial signal, such as RFI present, whose spectrum falls into the receiver passband, the complex I/Qtimeseries then represents a process that is only band-limited by the receiver hardware.The consequence of undersampling is frequency aliasing of higher frequency components into the atmospheric band of interest.This problem is especially critical in the U.S., where profilers at 449 MHz operate simultaneously with amateur radios.Although the problem of principal undersampling cannot be solved due to the fact that 1/τ 1/T , Wilfong et al. (1999a) achieved an improved time domain filtering using a four-term Blackman-Harris filter (Harris, 1978), instead of the usually applied boxcar filter of coherent averaging.While this kind of digital filtering helps to reject RFI, it is not helpful in the presence of ground and intermittent clutter.Those clutter signals fall well into the region of the desired atmospheric signal.May and Strauch (1998) proposed the use of linear convolution filters (digital FIR 4 filters) with a band rejection characteristic around zero Doppler shift (DC).This requires, however, a long filter sequence and also does not protect against intermittent clutter signals, which can occur at any frequency.Additionally, the transfer characteristic is fixed for a set of given filter coefficients.For that reason, wavelet domain filtering of ground and intermittent clutter (Jordan et al., 1997;Boisse et al., 1999) has been proposed.The main purpose of all these time domain operations is the filtering aspect, i.e. the intention is to "clean" the raw data from contaminating signals while leaving the desired atmospheric contribution ideally intact.In the following, we will concentrate on the clutter problem and investigate the properties of these 4 Finite Impulse Response signals and the possibility of applying nonlinear wavelet filters.There are not many investigations about the properties of RWP raw data.Normally, using statistical arguments, one assumes simply a Gaussian signal characteristic for atmospheric and clutter signals, as well as for noise (Doviak and Zrnić, 1993;Petitdidier et al., 1997).Recently, Muschinski et al. (1999) used data from a large-eddy simulation to derive I/Q signals for clear air scattering, and Capsoni and D'Amico (1998) presented a software-based radar simulator for generating time series from a synthetic distribution of hydrometeors.For our purpose, we assume that the Gaussian model describes sufficiently well both the atmospheric scattering component and the ground clutter signal.Intermittent clutter returns can be described by the simple model given by Boisse et al. (1999), with their main property being the transient character.
The 482 MHz wind profiler, whose data are used in this study, was installed at the Meteorological Observatory Lindenberg during the summer of 1996.The system is the prototype for three additional profilers to be installed in Germany in the future to supplement the operational aerological network of the Deutscher Wetterdienst (DWD).A summary of the main characteristics of the system is given in Table 1.For a more detailed description, the reader is referred to Steinhagen et al. (1998).The system is operated quasi continuously using a five beam configuration.All the main system parameters can be freely programmed which eases special investigations, such as the investigation of the detrimental ground clutter signal that was present in the system's The gap in the data was caused by this detailed investigation, as the radar was programmed to store time series data for about 30 minutes, only in the East beam, thus, no wind computations were possible for that period of time.
East beam from the 30 November to the 1 December 1999.During this event, the profiler was operated for a short period using only this beam (and low mode), while the huge amount of time series data was stored for further investigation (namely, the wavelet filtering).We now substantiate the radar parameter settings that were used in collecting the radar raw data: from the table of the radar's parameter settings, we find that the spacing of the time series data is t = NT = 8.784 ms.This corresponds to a Nyquist frequency of f N = 1/2 t = 56.92Hz, which gives, in turn, the maximum resolvable radial velocity v R = λf d /2 = 17.6 m/s.Clearly visible in Fig. 2 is the detrimental impact of ground clutter at the heights around 1400 m and 3000 m.The computed winds are obviously wrong and we will, therefore, look in detail into the problem.A more detailed, exemplary look into the raw data (coherently integrated I/Q-Time series) of gates 11 and 17, and the resulting power spectra (Fig. 3), immediately reveals that advanced signal processing for RWP is necessary to increase the accuracy of wind vector reconstruction.The time series at gate 11 shows the typical signature of a slowly fading, large amplitude ground clut-ter signal component, which corresponds to the narrow spike centered around point 1024 (zero Doppler shift) in the resulting power spectrum (compare also May and Strauch, 1998).Additionally, the time series at gate 17 shows a strong transient component in the last quarter.Such a signature is quite typical for a flier echo, as was shown by Boisse et al. (1999).This transient almost completely covers up any atmospheric signal in the power spectrum.

Applying multiresolution analysis and statistical estimations
For the problem at hand, the goal of the signal processing should be signal component separation, i.e. an automatic, reliable and stable extraction of the different contributions to the signal (noise, clutter, interference).Motivated by Daubechies (1992); Vetterli and Kovacević (1995); Louis et al. (1998); Meyer (1993) and Holschneider (1995), our purpose was to embed the filtering procedure into the known mathematical theory of wavelets.In general, mathematical experience concerning problems related to contamination removal or denoising shows that usually more than time domain filtering and Fourier domain filtering techniques are required to obtain optimum results.Often, most of the existing and implemented methods are insufficient.The main reasons for the particular effectiveness of wavelet analysis can be summarized as follows: -The fact that contamination appears often instationary or transient, and with a priori unknown scale structure, favors the superior localization properties of the wavelets.
A wavelet expansion may allow the separation of signal components that overlap both in time and frequency (Burrus et al., 1998).
-In order to effectively localize clutter components, one can use a great variety of wavelet filters (Daubechies, 1992;Dahlke et al., 2000;Teschke, 1998;Stark, 1992;Dahlke and Maaß, 1995).To choose a certain wavelet that especially suits the desired signal component, one can determine the properties of the clutter signal; otherwise, one can select a wavelet empirically.
-The wavelet expansion coefficients, β j k , drop off rapidly for a large class of signals, which makes the expansion very efficient (Burrus et al., 1998).
-The fast wavelet transform has a computationally complexity that is lesser than or equal to the fast Fourier transform; the algorithm is recursive (Kaiser, 1994;Louis et al., 1998;Burrus et al., 1998).This allows for an efficient implementation on digital computers.
Thus, the application of wavelet techniques to our particular problem seems to be promising.Before we start, let us briefly repeat the basics of multi-resolution analysis.
Let L 2 (R) be the space of functions of finite energy.Let φ be some function in L 2 (R), such that the family of translates of φ form an orthonormal system.We define Further, we define linear spaces by Assuming that φ is chosen in such a way that the spaces are nested: then the sequence {V j , j ∈ Z} is called a multi-resolution analysis.This concept was introduced by Mallat (1989); Meyer (1993).φ is called the father wavelet.Furthermore, one may define subspaces W j by V j +1 := V j ⊕ W j and iterating this we have Assuming that our data may be described by some f ∈ L 2 (R), we can represent the signal as a series where {ψ j k }, k ∈ Z is an orthonormal basis in W j .The function ψ is called mother wavelet.This expansion is a special kind of orthogonal series.Hence, it would be useful to search in the framework of nonparametric statistical estimation theory for an applicable method to solve our problem (Donoho and Johnstone, 1992).In case of orthogonal series estimation, the idea of reconstructing the desired atmospheric signal is simple.Basically, we replace the unknown wavelet coefficients in the wavelet expansion by estimates which are based on observed data.For that, we need a selection procedure to choose relevant coefficients since the main emphasis of performing wavelet domain filtering is to create a suitable, i.e. problem matched, coefficient selecting procedure.To separate the atmospheric signal component, we apply statistical estimation theory.A side effect of using statistics is to obtain a measure of reconstruction quality.A typical quality measure is a loss function/estimation error.Minimizing the error function reveals an objective evaluation and a self-acting filter algorithm.
The following sub-section describes the construction of our atmospheric-signal-estimator.In advance, we briefly remark that in the following section, we assume that our signal belongs to some Besov space, i.e. a generalized mathematical function space.One special example is the previously introduced function space L 2 (R).But sometimes it makes more sense to suppose that the derivatives of our signal are of finite energy as well.In this and other situations, the framework of Besov spaces is an adequate mathematical tool for our application.A Besov space, denoted by B s pq , depends on three parameters: s smoothness, the number of bounded derivatives and p, q which describe the underlying function space L q (l p ).In the following, we make use of some well-known facts of estimation theory, which are valid for almost all Besov spaces (Donoho and Johnstone, 1992;Donoho et al., 1993;Johnstone and Silverman, 1995;v. Sachs and MacGibbon, 1998;Dahlhaus et al., 1998;Härdle et al., 1998).If our signal is an element of one of these spaces (which is true for all practical signals), we can adapt wavelet threshold estimators.The main advantage of this framework is that we can use existing rules for evaluating bounds and rates of convergence for our loss function, which describes the quality of our reconstructed atmospheric signal component.By optimizing bounds and rates of convergence, we obtain self acting algorithms.
For our purpose, we only need the following characterization of Besov spaces: A function f belongs to B s pq if (2 j (s+1/2−1/p) β j. l p ) q 1/q < ∞.
We are looking for optimal reconstructions of functions belonging to some subset F s pq (M) = {f ∈ B s pq : J s pq < M}.For our calculations, we assume that the function is in L 2 (R) and s is small.
From given measurements (Y 1 , . . ., Y n ), we want to estimate the function f in the simple model We assume that we have the X i on a regular grid and ε is a random variable (a stochastic process which describes all non-atmospheric components).The basic idea is to replace the wavelet coefficients in the series expansion by empirical estimates where the X i are time stamps and the Y i are observations.A straightforward linear estimation is given by the projection onto a subspace V j 1 To appraise this estimator, it is known that one may measure the expected loss or the risk (in L 2 sense) E fj 1 − f 2 2 .This measure is the so-called MISE (mean integrated squared error).To determine the MISE, one may decompose it into E fj 1 − E fj 1 2 2 (stochastic contribution) and E E fj 1 − f 2 2 (deterministic contribution).Under certain conditions, one may find bounds for MISE: sup A minimum of the sum is given by sup furthermore, one can generalize this result for p > 2 sup 2s+1) .
This gives us an upper bound for the maximum risk.This bound becomes small if the number of observation increases and if j 1 (n) is determined in such a way that the bias and stochastic bound are balanced (for detailed computations of bounds, see v. Sachs and MacGibbon, 1998;Donoho et al., 1993;Donoho and Johnstone, 1992;Dahlhaus et al., 1998).
Obviously, this kind of linear estimation includes oscillating components, in particular, the clutter components.This phenomenon occurs because we have taken the whole set of wavelet coefficients up to scale j 1 , i.e. we have not performed any filtering step thus far.In the following, we need a suitable selection procedure for the coefficients in order to perform the necessary filtering step.We apply a so-called hard thresholding and soft thresholding, respectively.This methodology was introduced and adapted to several problems by Donoho and Johnstone (1992); Donoho et al. (1993).It is based on taking the discrete wavelet transform (using a multiresolution analysis), passing the transform through a threshold (actually, the expansion coefficients are thresholded) and then taking the inverse DWT to obtain a filtered reconstruction.Note that this type of thresholding is usually applied in a different way, by removing coefficients below a certain threshold in order to "de-noise" the data (Burrus et al., 1998, see Fig. 4).The functions for hard and soft thresholding are defined by and the modified functions used here for hard and soft thresholding are given by the rule η * (u) = u − θ * (u): Here, λ is an adequate threshold.Applying this rule to our linear wavelet estimator, we obtain a nonlinear estimator where η * is η s or η h , respectively.
If the threshold λ is specified according to the asymptotic distribution of the empirical coefficients, then only those coefficients remain which are supposed to carry significant signal information.These are finally used for the reconstruction by the inverse wavelet transform.For the correct level of significance, an appropriate choice of the threshold λ is needed.In general, this does not only depend on the sample size n, but also on the resolution scale j , and location k of the coefficients.In the case of regression with nonstationary errors, we have to use both a level and location dependent threshold rule (v. Sachs and MacGibbon, 1998).The resulting non-linear estimator does not only provide local smoothers, but, in many situations, achieves the nearminimax L 2 -rate for the risk of estimation, i.e. v. Sachs and MacGibbon (1998) for (random) thresholds λ j k satisfying n for any positive constant C: where σ j k is the variance and M j denotes the number of the coefficients used in the nonlinear estimator.The optimal threshold rate (1/n) 2s/(2s+1) is attained only for the ideal threshold.However, in practice, this is unknown.Therefore, we have to replace σ j k by some estimation σjk , which results in random thresholds λjk = σjk 2 log M j .Hence, the logterm has to be understood as the price for some data-driven threshold rule, and it originates due to the estimation of the unknown variance σ 2 j k = Var( βjk ).We conclude that we may adapt an estimation rule for our desired atmospheric signal component where the quality is measurable in the sense of L 2 -risk.This means the procedure used displays bounds for our reconstruction, and we may easily determine the rate of convergence.The calculation of the wavelet coefficients can be done by using the fast wavelet algorithm which is easily implemented.

Removing clutter
In this section, we will demonstrate the performance of nonlinear wavelet filtering.This is done both with simulated and with real data.For a better understanding, we particularize Fig. 1 to see where we have inserted the wavelet tool.To apply our procedure, a more substantiated algorithm flow diagram is shown in Fig. 5.
Following the first box in the algorithm flow diagram, one has first to determine the analyzing wavelet (high and low pass filter coefficients).Usually, the decomposition of a signal in a basis (i.e. a wavelet series) has the goal of highlighting particular properties of the signal (Mallat, 1999).In the problem of wind profiler signal filtering, the desired atmospheric signal component can be contaminated with spurious signal components.The ultimate goal is obviously to find a wavelet basis, which would allow a separation of the desired and the unwanted parts of the signal, i.e. which would have the ability to approximate the unwanted signal components (ground clutter, intermittent clutter) with only a few non-zero wavelet coefficients.In other words, the wavelet ψ has to be chosen in such a way that a maximum number of wavelet coefficients, β j k , are close to zero.This depends primarily on the regularity of the (contaminating) signal f , the number of vanishing moments of the wavelet ψ, and the size of the wavelets support.If f is regular and ψ has enough vanishing moments, then the coefficients β j k are guaranteed to be small for small scales.If, however, the signal f contains isolated singularities, the strategy to have a maximum number of small wavelet coefficients would be to reduce the support size of the wavelet.Unfortunately, there is a tradeoff between both properties for orthogonal wavelets: if ψ has p vanishing moments, then its support size is at least 2p − 1.The best compromise between those two requirements are Daubechies wavelets, which are optimal in the sense that they have minimum support for a given number of vanishing moments.
There have been no detailed investigations thus far about the regularity properties of contaminating wind profiler signals, but there is evidence that these can be both "quite regular" (ground clutter) or "not so regular" (intermittent clutter).Thus, the Daubechies family was selected.The order of the Daubechies wavelet was chosen according to the regularity condition, which we have conservatively chosen to be rather small (s ≤ 1).To approximate correctly a function of B s pq , we need to select an analyzing wavelet of regularity [s] + 1.A wavelet with regularity of the order of s = 2 and minimal compact support is the Daubechies-2-wavelet; hence, we have chosen this one for our calculations.Mathematically, it is no problem to increase the wavelet order (regularity), but the wavelet support size and the number of filter coefficients also increases, and this will decelerate the algorithm.Finally, we note, in passing, that we have concentrated on the fast wavelet transform (multiresolution analysis), which is a special case of the discrete wavelet transform.Obviously, for an online algorithm, the number of operations per data point is limited.The fast wavelet transform is, therefore, the best choice, since it has the highest numerical efficiency (i.e. it is faster than the fast Fourier transform).This, of course, restricts the possible choices of the underlying basis wavelet.The number of decomposition scales is determined by balancing the stochastic and the deterministic part of the MISE.Thus, the optimal scale may be evaluated automatically by the rule 2 j 1 (n) n 1/(2s+1) .After fixing the main parameters, one may start the wavelet decomposition of the in-phase and the quadrature-phase time series.To separate the atmospheric component, the algorithm calculates for each decomposition level the local thresholds λjk .
Additionally, one may use histogram information, which displays the empirical distribution of the coefficients αk and βjk .In particular, if the signal was contaminated by an airplane echo, the main part of the observations is concentrated in a small neighborhood around zero.If there is no airplane echo, the coefficients are exponentially distributed (see Fig. 6).
The histogram methods acts as follows: we denote by h j (k), the histogram function of the coefficient sequence of scale j , and by H j (z), the connected empirical distribution function.We know that H j is monotonic increasing, continuous from the right and a step function.If H j (z) is given then all values z i may be recognized completely; this means that one can detect the smallest value q 0 with H j (q 0 ) ≤ c j,α (socalled empirical α-quantile).If we now determine a lower bound for the number of coefficients we want to have available for reconstructing, we may easily evaluate q j,α by solving q j,α −q j,α h j (k)dk = c j,α .
We define the histogram-based threshold by λ hist j k := q j,α .
Since there is only empirically information and no model about the characteristics of intermittent clutter echoes, we assume that the histogram method should remove a maximum of 15 percent of the observations and hence, we have chosen c j,α = 0.85.This is, of course, just a heuristic value.For our dataset, this value has given the best results for the loss function.We are quite confident that the rule is robust if a larger percentage of the dwell time is contaminated with flier echos.However, more research about the properties of the contaminating signals and the distribution of the wavelet coefficients is certainly needed.
In the case of having only ground clutter, λjk and λ hist j k are almost equal.In the case of intermittent clutter, we choose our data driven threshold, λjk , by taking the minimum of λjk and λ hist j k .Hence, in case of having bird or airplane echoes, it may occur that the resulting threshold underestimates the threshold evaluated by minimizing the MISE.Yet from extensive test calculations, we know that no problems accrue.In addition, by using this simple histogram rule, the accuracy of the thresholding step increases.But the price of applying the rule is an enlargement of the number of calculations per data point.
To observe how this algorithm works, we start by simulating one easy test sample.Using the statistical-stochastic approach of Zrnić (1975) to generate I/Q-timeseries, we first generate an atmospheric signal with Gaussian characteristics in the frequency domain.We choose the Doppler frequency of the atmospheric signal close to zero to force the separation problem.Now we add a noise variable and a ground clutter peak, which is generated by a narrow Gaussian.The order of the ground clutter amplitude is much higher than the atmospheric signal amplitude.Since the algorithm removes the ground clutter completely, the reconstructed signal consists only of the atmospheric part (and some noise).This demonstrates impressively the difference between the nonlinear wavelet filtering method and the Fourier methods and digital filtering: the spectra of clutter and atmospheric signal can overlap as much as they want; nonetheless, we can still separate the two components.The different amplitude of both signals allows for the discrimination.
For intermittent clutter, one of the advantages of waveletbased techniques is certainly the ability to describe a transient signal with only a few wavelet coefficients.This is caused by the finite support of the wavelet basis, in contrast to the basis functions (the so-called windowed Fourier atoms) e iωt g(t − u) of the windowed Fourier transform.It is the localizing properties of wavelets (Burrus et al., 1998) that makes the wavelet transform especially suited for filtering of transient signals (e.g.intermittent clutter).
To expose how the routine is acting on measured RWP time series, we eventually go back to the presented "real life" problem (example Fig. 3) in order to demonstrate the robustness of the method.The problem was that the signal at gate 17 was contaminated by intermittent clutter (aircraft echo) and the signal at gate 11 was contaminated by persistent ground clutter.The spectra obtained with the standard signal processing were severely affected by clutter contributions to the received signal and thus, the moment estimation and finally the wind vector determination were significantly biased.Figure 7 shows exemplarily how wavelet thresholding was realized in decomposition sequences α k and β 1k of gates 11 and 17.The dotted lines may be identified with the threshold.It can be observed that in both cases, the clutter components could be completely removed.

Conclusions
This paper discusses an algorithm that employs discrete multiresolution analysis and nonlinear estimation theory to separate the atmospheric Doppler signal in RWP measurements in the presence of contaminating signals.Using simulated and real wind profiler data, we have demonstrated that wavelet thresholding is effective in removing ground and intermittent clutter (airplane echoes) from the RWP raw data (I/Q timeseries).The presented wavelet based filtering technique is self-acting and, therewith, a step toward an automatic algorithm for clutter removal in Doppler spectra.Real time implementation in profiler systems is required to test the new method with a substantially longer dataset, preferably in parallel with the standard processing (comparison), and to demonstrate its use for operational applications.

Fig. 1 .
Fig. 1.The figure shows the flow diagram of 'classical' digital signal processing.

Fig. 2 .
Fig. 2.The figure shows the final result of the measurement with the 482 MHz RWP at Lindenberg (Germany) on the 30 November and 1 December 1999.The wind barbs are color-coded according to the wind speed.Note the effects of the persistent ground clutter around the 1500 m and 3000 m heights.The gap in the data was caused by this detailed investigation, as the radar was programmed to store time series data for about 30 minutes, only in the East beam, thus, no wind computations were possible for that period of time.

Fig. 3 .
Fig.3.The left part shows the "stacked spectrum" plot, i.e. the Doppler spectra for each range gate, for the radar dwell at 08:53:38 UTC on 1 December 1999.The right figures give a detailed look into the raw data (I/Q-timeseries) and the Doppler spectrum for the gates 11 and 17 (whose data will be wavelet processed).The black arrows indicate the estimated first moment (i.e. the radial velocity).

Fig. 5 .
Fig. 5. Left: The flow diagram extended by the wavelet tool.Right: The wavelet algorithm flow diagram.

Fig. 6 .
Fig. 6.This figure shows typical histograms of the wavelet coefficients (see text).The upper histogram represents an in-phase series without an airplane echo and the lower histogram represents an inphase series with an airplane reflection.

Fig. 7 .
Fig.7.Decomposition (α k und β 1k ), reconstruction and Fourier power spectrum of gate 17 (top) and gate 11 (below).The black curves in the power spectra representations display the decontaminated spectra.One clearly recognizes are the differences of moment estimations; see the computed first moment before (gray arrow) and after (black arrow) the filtering step.

Fig. 8 .
Fig.8.Top left: Simulated Fourierpower-spectrum with strong ground clutter influence and with an atmospheric signal overlapping the ground clutter peak.Lower left: I/Q series derived from the simulated Fourierpower-spectrum using theZrnić (1975) method.Lower right: I/Q time series after applying the nonlinear wavelet filter.Top right: Resulting Fourier-powerspectrum based on the reconstructed (filtered) signal.

Table 1 .
Technical specification of RWP and radar operating parameters