Mitigation of Wind Turbine Clutter With Digital Phased Array Radar

Wind turbine clutter (WTC) contamination affects polarimetric meteorological observations and automatic high-impact weather detection algorithms in several ways, such as, misidentification of precipitation echoes, false mesocyclone detections, and incorrect storm cell identification. With the dramatic growth of the wind power industry, this would only aggravate in the future. Since fully digital Phased Array Radar (PAR) is a promising candidate technology for the next generation of weather radars, evolutionary signal processing algorithms that make use of their capabilities should be investigated to mitigate WTC contamination. This article investigates the use of Space-Time Adaptive Processing (STAP) to mitigate WTC contamination with ground-based polarimetric PAR. First, a flexible wind turbine time-series signal simulator is developed to characterize the contamination signatures in the space-time domain. Then, a STAP algorithm that removes the contamination is presented and demonstrated on simulated data. Two digital radar back-end architectures are considered to evaluate the performance of the proposed algorithm. One with digital sub-array outputs the other one with fully digital outputs (i.e., element-level digital). Results indicate the WTC spectrum has a characteristic structure in the space-time domain, and that biases induced in polarimetric weather variables can be more effectively mitigated using an fully digital PAR.


I. INTRODUCTION
The use of alternative energy sources to produce electricity is growing dramatically in the United States. Driven by the U.S. Department of Energy initiative to increase wind energy production to 20% nationwide by 2030, and to 35% by 2050 [1], [2], wind-turbine farm installations are increasing rapidly. In addition to the quick growth in the total installed capacity, the size of individual wind turbines has exponentially increased over the last 30 years [3] to obtain a reduced price per generated kilowatt hour. For instance, modern turbines have rotor diameters and tower heights on the order of 145 m and 180 m [3]. As the turbines get larger, The associate editor coordinating the review of this manuscript and approving it for publication was Hasan S. Mir. they are seen by weather radars at farther ranges in their line of sight.
Radar signals received from wind turbines can contaminate meteorological measurements impacting data quality. These are commonly referred to as wind turbine clutter (WTC), since these come from non-meteorological scatterers. Clutter filtering algorithms conventionally used in the signal processor are effective in removing signals from non-meteorological scatterers with narrow Doppler spectra and mean velocities near zero ms −1 . However, it is well-known that wind turbines produce clutter signals that are challenging to identify and remove from the Doppler spectrum using conventional spectral-filtering techniques, and can result in significant biases on radar-variable estimates [4], [5]. The WTC contamination is especially difficult to remove because despite the fact that wind turbines are non-moving scatterers, FIGURE 1. Simulated Doppler spectra from stationary ground clutter (red), weather (blue), and wind turbines (yellow). It can be seen that due to the high spectral width of the WTC spectra, and its overlap with that of the weather signal, it can be challenging to remove. their blade rotation produces non-zero Doppler velocities that result in relatively wide spectra and could overlap with Doppler velocities from hydrometeors being measured. Consequently, conventional Doppler-spectrum-based clutter filters are usually ineffective at removing it. Simulated Doppler spectra from stationary ground clutter (red), weather (blue), and wind turbines (yellow), are shown in in Fig.1. It can be seen that due to the high spectral width of the WTC spectra, and its overlap with that of the weather signal, it can be challenging to remove. This can cause false mesocyclone detections in automatic weather-data algorithms used by forecasters [6]. With the target date for increased wind energy production rapidly approaching, WTC contamination could potentially impact the U.S. National Weather Service (NWS) operations in the future [6], [7].
Many research efforts have been carried out over the years to characterize, detect, and remove WTC from weather radar observations. Reference [8] studied the radar cross section characteristics of wind turbines as a function of the azimuthal orientation of the turbine with respect to the radar using a three-dimensional model. Reference [9] characterized WTC contamination as seen by radars of the Weather Surveillance Radar -1988 Doppler (WSR-88D) network and proposed a multiquadric interpolation clutter-removal technique. The technique relies on a neighborhood of uncontaminated weather-variable estimates and information regarding the exact locations of the wind turbines. Reference [10] measured the radar cross-section of a scaled wind turbine model and suggested the use of telemetry information to retrieve blade velocities to model and attempt to mitigate the contamination. Reference [11] developed an automatic fuzzy-logic-based WTC contamination detection algorithm capable of flagging contaminated resolution volumes with a probability of detection of approximately 90% and a probability of false alarm lower than 1%. Nevertheless, their algorithm does not remove WTC contamination and they suggest using a sophisticated signal processing technique on range gates with positive detections. Finally, [12] and kong2012wind explored the use of adaptive spectral processing to remove WTC contamination by means of a band-pass filter applied to the Doppler spectrum. None of these considered the use of modern radar technology with advanced beamforming capabilities.
Phased Array Radar (PAR) technology offers advanced capabilities that could support new signal processing methods to improve the quality of weather radar observations, and is being considered for the next generation of weather radars [14], [15]. Although PAR technology can provide advanced capabilities, important questions regarding the PAR architecture and level of signal digitization need to be investigated to inform a future NWS radar-procurement process. Recent research efforts indicate that the single-face rotating PAR architecture is capable of meeting temporal resolution needs for future polarimetric observations [16]. Nevertheless, these studies have not considered the number of digital channels needed. More system flexibility is achieved with higher number of digital channels, at the price of increased cost and data throughput.
Digital PAR signals received at multiple phase centers could provide additional information used by a signal processing technique to mitigate contamination from moving clutter scatterers [17]. Traditional signal processing techniques like Fourier beamforming and Doppler processing are normally applied independently. This limits the signal processing operation to one domain at a time: spatial for beamforming, temporal for Doppler processing. Researchers have reported that this separation of functions may be sacrificing capabilities because the digital signal processor is not making use of all available degrees of freedom [18]. To make full use of the radar's resources, joint space-time processing should be considered. However, a high number of digital channels is needed to minimize the effect of antenna grating lobes, which could create spectral replicas of the clutter spectrum in the spatial domain, making the filtering operation more challenging.
Of interest here is to explore the joint use of digital beamforming (DBF) and Doppler processing to characterize the WTC space-time spectra, and evaluate potential ways to remove it. With DBF, digital signals received at different array locations are combined to spatially filter (i.e., through custom antenna beam pattern shapes) undesired returns. Similarly, Doppler processing is traditionally used to filter undesired signals in the temporal frequency domain. The combined used of these methods is known as space-time adaptive processing (STAP). The STAP concept has not been applied to data from ground-based weather surveillance radars because traditional parabolic-reflector radars (like the WSR-88D) have a single phase center. Additionally, since ground-based weather radars are stationary, ground clutter returns have no Doppler shift and can thus be removed by a narrow one-dimensional Doppler notch filter (considering internal clutter motion) [19], [20]. Nevertheless, contamination from moving clutter scatterers, such as blades from a wind turbine, are difficult to distinguish and remove from weather signals.
Ongoing research efforts at the University of Oklahoma's (OU) Advanced Radar Research Center (ARRC) are aimed at demonstrating unique capabilities of digital PARs for polarimetric weather observations [21]. In particular, the fully digital Horus radar is a dual polarization, mobile, S-band system recently completed at the ARRC with support from the National Oceanic and Atmospheric Administration (NOAA). It is a unique platform to demonstrate the benefits of digital PAR technology for meteorological observations [22]. Exploiting the advanced capabilities of a highly-digital PAR, the quality of mission-critical data used by NWS forecasters and automatic algorithms could be improved.
The development of a STAP-based algorithm to mitigate WTC contamination is the focus of this article. Being PAR the candidate system to replace the polarimetric, reflectorbased WSR-88D, leveraging their unique capabilities to mitigate WTC contamination is fundamental. With access to a fully digital dual-polarization PAR system like Horus, these evolutionary signal processing techniques could be demonstrated on a real system. Providing data with higher quality to forecasters and automatic algorithms could lead to improved severe weather warnings and weather forecasts. Therefore, it is important to investigate the feasibility of using STAP techniques on PARs to reduce WTC contamination on radar-variable estimates. The article is structured as follows. In Section II, a short review of the STAP theory is provided. A WTC time-series simulation for digital PARs is developed in Section III to evaluate WTC spectra in the spacetime domain. In Section IV, a WTC-mitigation algorithm capable of recovering meteorological signals of interest is described. Experimental results are presented in Section V. These include an evaluation of the proposed algorithm on a sub-digital radar architecture (NOAA's Advanced Technology Demonstrator), a fully digital architecture (Horus), and an emulation using actual digital array weather data and simulated WTC signals. Conclusions are provided in Section VI.

II. SPACE-TIME ADAPTIVE PROCESSING
The development of the active PAR technology motivated [23] to introduce STAP techniques to the airborne radar community in the early 1970's. It consists of a two-dimensional space-time filter designed adaptively to maximize signals of interest and filter unwanted signals, herein referred to as clutter signals. This promising signal processing technique inspired the airborne radar researchers who then vigorously established and demonstrated the theory of STAP to improve the quality of airborne radar detections [24], [25], [26], [27]. Nevertheless, most of the STAP theory has been developed for airborne or spaceborne radars to detect and track point targets, with little-to-none for ground-based weather PARs.
Airborne surveillance radars are typically equipped with PAR systems to detect and track targets of interest. The need for combined space-time processing in airborne radar comes from the inherent problem of platform motion-induced Doppler shift. Because of the platform motion, stationary ground clutter has none-zero velocity relative to the radar. Clutter targets ahead of the aircraft have positive Doppler shifts, and clutter targets behind the aircraft have negative shifts. Thus, the signature of ground clutter echoes in the space-time domain appears as a diagonal line where the slope is determined by the radar frequency, the platform velocity, and the inter-element spacing [28].
Extensive research on STAP has been published for airborne systems. After [23] introduced the theory of adaptive array processing, researchers started investigating the characteristics of ground clutter targets in the space-time domain. Reference [29] developed theoretical clutter models for STAP, characterized its power spectra for different scenarios, and studied the eigenvalues of the clutter covariance matrix. He concluded that to suppress ground clutter effectively, the radar should have as many degrees of freedom as dominant clutter eigenvalues in the clutter covariance matrix. Adaptive estimation of the clutter and interference covariance matrices is a challenging problem in STAP because it usually involves access to large sets of training data (i.e., real clutter or interference observations) not readily available in practical scenarios [30]. An important milestone for the establishment of STAP was the Mountaintop Program, which was initiated in 1990 to support the evaluation of STAP techniques to meet mission requirements for next generation airborne early warning systems [31]. The radar consisted of a 14-element, horizontally-polarized, ground-based PAR system with independent receive channels for each element. The inverse displaced phase center array technique was implemented in a separate co-located transmit array to emulate platform motion and generate aircraft-like clutter returns. Sets of Mountaintop data were made available to the public, which further motivated the research community. There are however drawbacks associated with STAP, the most serious ones being the intensive computational load imposed on the signal processor, and the limited amount of training clutter data available in practical scenarios.

A. STAP SIGNAL MODEL
Consider a unit-amplitude electromagnetic plane wave impinging on an N -element uniform linear array (ULA) with inter-element spacing d. A one-dimensional ULA is adopted here for simplicity, but the derivations can be extended to planar arrays without loss of generality. Define the plane wave angle of arrival θ s as relative to boresight (array normal) and designate the first antenna element as the phase reference. Then, the received spatial signal vector is x s = a s s s (f sp ), being a s a random complex voltage and is the spatial steering vector, with f sp being the spatial frequency given by Assume now that we transmit a periodic train of M pulses at a Pulse Repetition Time (PRT) of T s . The time delay of the n th pulse due to the pulse-to-pulse motion of r of a point scatterer initially located at a range of r 0 is where c is the speed of light. The pulse-to-pulse phase is then with f being the radar frequency, ω is the angular frequency, and λ is the wavelength. From (3), the Doppler frequency is Then, the Doppler shift observed is proportional to the radial velocity v r of the point scatterer. The received sample-time signal vector of a point scatterer is then x t = a t s t (f d ), being a t a random complex voltage and is defined as the temporal steering vector. The space-time signal vector corresponding to the echoes returned from a point scatterer with spatial frequency f sp and Doppler frequency f d is [28] where a s−t is a random complex voltage, s s−t (f sp , f d ) is the space-time steering vector, and ⊗ represents the Kronecker product. Notice that x s−t is a NM -long space-time signal vector or snapshot, that represents the signals received at each of the N antenna elements for each of the M transmitted pulses. This is illustrated in Fig. 2 using the radar data cube concept, where the dimensions of the cube are: spatial for phase centers, temporal for sample-time pulses, and range along the radial. The space-time processor linearly combines the elements of the data snapshot, yielding the scalar output In (8), 'H' denotes Hermitian transposition and w is a NM -length weight vector. The goal of the digital signal processor is to maximize the signal-to-interference-plusnoise ratio (SINR) by carefully selecting the space-time weight vector w. From [26], the SINR power at the filter output is being R the covariance matrix of the undesired signals (e.g., clutter, interference, white noise), and σ 2 s is the average signal power per channel per pulse. The optimum w that maximizes the SINR as posed in (3) That is, to design an effective STAP filter, the signal processor needs a good estimate of the clutter covariance matrix R.
Signals received at the radar generally consist of the sum of signals from clutter, interference, noise, and scatterers of interest. Therefore, it is challenging to know the clutter covariance matrix R a priori. Instead, one can determine which spectral components come from clutter, and estimate the covariance matrix from these. For the question of interest in this article, it is important to characterize the space-time spectra of wind-turbine signals, therefore, a time-series simulation is developed.

III. WTC TIME-SERIES SIGNAL SIMULATION
The first step to investigate the potential of applying STAP to mitigate WTC contamination using a digital PAR is to characterize WTC signals in the space-time domain. Since no data of this nature are available, a time-series simulation framework is developed.

A. SIMULATION OF WTC SIGNALS
Simulations provide a mechanism to evaluate the performance of clutter mitigation algorithms, since the simulation parameters are known. To this end, and without loss of generality, assume a ULA with N radiating elements. The far-field antenna pattern produced by the i th element at a point (x, y, z) in space can be expressed as the product of the antenna element electric-field pattern function f (θ, φ) and a propagating plane wave [32] where λ is the wavelength, k = 2π/λ is the free space wave number, (θ, φ) are the azimuth and elevation angles of the point were the field is being measured, and being (x i , y i , z i ) the coordinates of the i th antenna element. It can be shown that the complex voltage received at the i th element from the field backscattered by a point scatterer is proportional to E i (r i , θ, φ) 2 , and the backscattering cross-section of the scatterer σ b (θ, φ). Then, applying the superposition principle we can add the far-field contributions received from L scattering centers at the i th antenna element, Using (12), signals received at each antenna element of the array can be computed for a large collection of scattering centers. These complex values determine the magnitudes and phases of signals received at each antenna element of the array. Then, to generate time-series data for a train of M periodic pulses, we apply (12) to M sets of scattering centers sampled every T s seconds.
To simulate time-series data from a wind turbine, a structural model of a commercial-grade wind turbine (Vestas 4MW) was imported into the Altair Feko computational electromagnetics software [33]. An illustration of the wind turbine model is shown in Fig. 3. Back-scattering crosssections for horizontally and vertically polarized incident waves as a function of angle (θ l , φ l ) were derived using Feko, with a sampling of 2 • in both dimensions. The model was discretized using a triangular mesh, which represents the wind turbine as a collection of scattering centers at time t 0 , and imported into Matlab. In this article, only the horizontal polarization is considered for simplicity.
The complex voltage for the first sample is calculated by adding the back-scattered signals from every point in the discretized wind-turbine model, using (12), and for every element in the array. To simulate the remaining M − 1 samples, a rotation matrix is applied to the scattering centers representing the blades. The rotation angle is θ WT = ω WT T s , where ω WT is the angular velocity of the blades in rad. s −1 . An illustration of the blade rotation is shown in Fig. 4, with ω WT = 90 • , and T s = 0.1 s (these unusual parameters are used here only for ease of visualisation). Finally, M sets of scattering centers are obtained as the turbine's blades rotate, and using (12) for each time-step, the signals for each antenna element of the array can be simulated. Additive white Gaussian noise is added to set the desired SNR level and account for effects in the transmit/receive chains (e.g., receiver noise power, waveforms). Table 1 shows the main simulation parameters and benchmark values used.
This time-series simulation framework serves to investigate characteristics of WTC in the space-time domain, as well as the evaluation of candidate clutter removal algorithms. Future improvements proposed for the simulator include: considering multi-path reflections, incorporating mutual coupling and diffraction effects of finite arrays, using actual element patterns (co-and cross-polarization).

B. CHARACTERIZATION OF WTC IN THE SPACE-TIME DOMAIN
With a flexible simulation framework able to produce timeseries in-phase and quadrature (IQ) data as a function of the wind turbine's physical characteristics and the radar acquisition parameters (i.e., M , T s , and polarization), WTC signatures are explored in the space-time domain. This would reveal distinctive signatures produced by wind turbines that an automatic algorithm could recognize and attempt to remove.
Wind turbines consist of three main structural components: tower, hub, and blades. Because the tower is a large stationary structure typically made of steel, it is expected to produce strong clutter returns near the center of the Doppler spectrum (a velocity of 0 ms −1 ), where a narrow spectrum width caused by slow structural motion could be observed. The hub is at   Table 1, (b) with parameters shown in Table 1, except a turbine orientation of 20 • . the top of the tower and holds the blades together. It rotates at the same angular velocity as the blades but its linear velocity is low because it is at the rotation center. Researchers have recognized this slowly moving structures near the hub that cause the ''hub contamination'' [11]. Finally, blade velocities are largest at the tips and decrease linearly as their structure gets closer to the hub (i.e., v = r ω WT ).
Assuming a radar beam is pointed directly at the wind turbine (i.e., turbine azimuth = 0 • as in Table 1), and considering the blade motion described, it is expected that the WTC contamination will produce signatures in the shape of straight lines in the space-time domain. In other words, the velocity of points along the blades' structure decreases linearly as they get closer to the hub and so does the angular position of these points as seen from the radar's phase centers. The signature depends on several factors including the angular velocity of the blades, the location (i.e., azimuth and range) and orientation of the turbine with respect to the radar. It is expected that as the range of the turbine increases its angular width as seen by the phase centers will decrease and so will the slope of the WTC contamination lines in the spectral domain.
The simulator is used to evaluate the previous intuitive analysis. A ULA with 64 digital elements separated by d = λ/2 = 5 cm observing a wind turbine is simulated. The turbine is located at a range of 500 m with its are blades  rotating at an angular velocity of 60 • s −1 . The number of pulses is M = 64, the PRT is 1µs, transmitted waves are polarized horizontally, and the SNR is set to 20 dB. Two turbine orientations are simulated, one with the blades directly facing radar, the other with a 20 • rotation about its z axis. The simulated WTC spectra in the space-time domain are illustrated in Fig. 5. This simulation framework can characterize WTC signatures with different turbine configurations and radar acquisition parameters. Note that the WTC spectral components will generally appear as strongly aligned ridges in the spectrum, each corresponding to a turbine blade. However, their locations and widths may be different depending on several parameters related to the radar and the wind-turbines themselves. On the radar side, the width of WTC ridges will be dictated by the antenna half power beamwidth in the spatial axis, and the number of samples (M ) in the Doppler axis. Further, if the PRT is relatively long, WTC spectral components may be aliased and wrap around the Doppler axis. On the turbine side, the location of the spectral components will change depending on the orientation of the wind turbine with respect to the radar, their spatial extent will change as a function of blade size and range (i.e., wider at closer ranges and narrower at farther ranges), and their Doppler extent will change as a function of blade rotation speed. Nevertheless, simulations show straight lines with strong WTC components will always results, albeit in different locations and/or different widths.. This is now used to develop an adaptive algorithm capable of mitigating WTC on meteorological signals.

IV. MITIGATION OF WTC CONTAMINATION USING STAP
A fundamental contribution in this article is the evaluation of removing WTC contamination from weather signals using digital PARs. WTC contamination detection algorithms could be used to detect and censor this clutter when no weather signals are present. The challenge is to remove WTC contamination from weather signals to accurately estimate the spectral moments and polarimetric weather variables [34]. At the current pace wind farm deployment is growing, and with their physical size also increasing, this could impact automatic Quantitative Precipitation VOLUME 11, 2023

Estimation (QPE) algorithms and consequently the NWS operations in the future.
To simulate polarimetric weather signals, the time-series simulation technique outlined in [35] is employed. The simulation produces M dual-polarization weatherlike time-series samples for a PRT of T s , and with desired radar-variable characteristics: signal power (S), radial velocity (v r ), spectrum width (σ v ), differential reflectivity (Z DR ), differential phase ( DP ), and copolar correlation coefficient (ρ hv ). Even though this simulation technique does not produce time-series data for multiple phase centers, its extension only requires adding the appropriate phase shift for each phase center corresponding to the element locations in the digital array. Nevertheless, because weather echoes are generally uniform over large volumes (i.e., with respect to the radar resolution volume), simulations in this article assume that weather signals with the same characteristics are received on all phase centers. In other words, independent realizations are carried out for each phase center using the same weathersignal characteristics.
Simulations of weather signals contaminated by WTC are produced using the techniques described. First, timeseries data for a wind turbine located at 500 m from the radar with a blade angular velocity of ω WT = 60 deg s −1 are simulated, as shown in Fig. 6(a). Then, weather signals covering the scanned sector are simulated for the horizontal polarization with SNR = 20 dB, v r = 5 m/s, σ v = 2 m/s, for M = 64 pulses and T s = 1 ms, as shown in Fig. 6(b). The WTC and weather signals are added to generate the contaminated weather spectra. The space-time domain representation of the WTC-contaminated weather signals is shown in Fig. 6(c). Assuming the WTC covariance matrix R WT is known, the SINR-optimum STAP filter w = R WT −1 s s−t (f sp , f d ) can be applied to accurately recover the weather signals, as shown in Fig. 6(d). The recovered weather space-time spectra are nearly identical to the uncontaminated one in Fig. 6(b), with the exception that some residual WTC power is observed. However, the clutter covariance matrix is unknown in practical scenarios and estimating it is not trivial. In the next section, we present a possible way of estimating R WT from measurements.

A. AUTOMATIC WTC CONTAMINATION REMOVAL USING STAP
Perfect knowledge of the clutter covariance matrix is not feasible in practical scenarios. Typical ways of estimating the clutter covariance include using clutter models, from the actual data (i.e., gates with clutter returns only), or estimating it from an eigen-spectra analysis [36]. However, a mismatch between the real and estimated covariance matrices can severely degrade the performance of the STAP filter to recover signals of interest.  An automatic WTC contamination filtering algorithm is developed next. The WTC covariance matrix will be estimated using the simulated received contaminated weather signals and knowledge gained from modeling WTC. Through space-time simulations in Section III-B, it is apparent that space-time WTC signatures are composed of one or more linearly oriented spectral components. These depend on the turbine's location and blade velocity with respect to the radar, and the radar acquisition parameters.
Gabor filters are used to determine the orientation of WTC ridges in the spectral space-time domain. These filters have been used extensively for feature discrimination in image processing applications [37], [38] and provide a way to identify WTC contamination in the space-time domain. A two-dimensional Gabor filter is a sinusoidal plane wave of a certain spatial frequency and orientation within a  Gaussian envelope. By applying a bank of Gabor filters with various orientations to the covariance matrix of the received signal (weather and WTC), the orientation of the WTC contamination ridges can be retrieved. They are determined by applying each filter (e.g., going from 0 • to 180 • in steps of 10 • ), and extracting the normalized outputs that exceed 0.7 (outputs are normalized to 1 using the maximum value). Therefore, spectral components in the Gabor filtered space-time spectrum exceeding 0.7 (filter output normalized to 1) are considered WTC returns. This is done independently for each Gabor filter (e.g., different orientations) and the final output is composed of all the components classified as WTC. Filter outputs with a spatial width narrower than 0.1 in the normalized angle domain (d/λ sin(d)) or a Doppler width narrower than 1 ms −1 are discarded to increase robustness. Once the orientations are known, the widths of the spectral WTC ridges can be estimated from the data. VOLUME 11, 2023 FIGURE 12. Simulation of WTC using the fully-digital Horus S-band radar architecture. The panels are similar to those in Fig. 6, for a wind turbine with blades rotating at 40 • s −1 , a turbine orientation of 60 • with respect to the radar, a range of 0.5 km (due to coarse beamwidth), a CSR of 20 dB, and with weather spectra at -10 m s −1 (v /2v a ∼ -0.2). Other parameters are the same as in Table 1 (e.g., M = 60 samples, f = 3 GHz).
The WTC covariance matrix R WT is then estimated assuming WTC spectra are composed of lines with widths (in space/time) as determined by the Gabor filtering process. Then, the total clutter covariance matrix is computed as R C = R WT + R N , where R N is the noise covariance matrix. A simplified block diagram of the algorithm is shown in Fig. 7. The algorithm described can be applied to each polarization (horizontal and vertical) independently to improve the quality of polarimetric-variable estimates. The filter designed using this algorithm is now applied to remove the WTC from the weather signals. The same characteristics as the previous example are used for the radar, the weather echoes, and the wind turbine (default parameters in Table 1). Notice that the range, orientation, or angular blade velocity of the turbine are unknown to the algorithm. The recovered weather space-time spectrum is shown in Fig. 8.
The algorithm is effective in removing most of the WTC contamination preserving the weather signals. The preliminary results show that it is feasible to remove WTC contamination using STAP. While there are several STAP-based algorithms for other applications (as those described in Section II), these were designed for different applications, and they leverage some knowledge of the specific clutter target of interest to estimate the covariance matrix. Considering that there are no previous studies characterizing the space-time spectral signatures of wind turbines using digital phased arrays, or the use of adaptive array processing methods to mitigate them, it is challenging to compare this method to others.
There are additional challenges that should be considered in future research. For example, a challenging situation occurs when the blades' radial velocities are larger than the maximum unambiguous velocity or Nyquist interval (v a ), then WTC contamination aliases and becomes more difficult to mitigate. It is also important to study the clutter signatures produced by wind farms, where several turbines are located within close proximity. A follow-up study systematically quantifying the performance of the algorithm over a large FIGURE 13. Similar to Fig. 12, for a wind turbine at 10 km from the radar. It can be seen that as the turbine is narrower with respect to the angular resolution (i.e., beamwidth) the space-time signature follows the normalized angle of 0 • , which can be easily mitigated.
parameter space (e.g., M , PRT, f , number of antenna elements, turbine orientation and range, etc.) is needed.

V. EXPERIMENTAL RESULTS
The algorithm presented in the previous section is capable of automatically removing WTC contamination using STAP on active PARs. However, several assumptions made to produce these simulations may not hold for actual radar systems. Characterization of WTC contamination for a sub-array digital PAR and a fully digital PAR are presented in this section.

A. SUB-ARRAY DIGITAL PAR
The Advanced Technology Demonstrator (ATD) radar system is a full-size, active, S-band, planar, dual-polarization PAR installed in Norman, OK. It was funded jointly by the NOAA and the Federal Aviation Administration (FAA) [39], [40]. The ATD antenna is composed of 76 panels, where each panel consists of an 8 × 8 set of radiating patch-antenna elements with dual linear polarization (H and V), for a total of 4,864 elements.On receive, the antenna is partitioned into 24 overlapped sub-arrays (consisting of 8 panels each, 2 in y by 4 in z) per polarization to produce lower sidelobes and partially mitigate grating lobes of the sub-array pattern outside of the main beam [41]. The radar's aperture is illustrated in Fig. 9.
The WTC simulation framework presented was extended to simulate the ATD radar. To this end, signals for each of the 4,864 antenna elements are simulated and added coherently according to the sub-array distribution to produce time-series data for 24 digital phase centers. With a broadside half-power beamwidth of 1.58 • , and its phase centers separated by approximately 4λ, some WTC can be mitigated, although spectral replicas of the WTC are present due to the grating lobes in the sub-array pattern (i.e., d = 4λ > λ/2). A simulation of WTC as would be seen by the ATD is presented in Fig. 10. In this case, the simulated wind turbine is located 20 km from the radar, facing it directly (i.e., turbine orientation = 0 • ), and the blade angular speed is 40 • /s. VOLUME 11, 2023 The algorithm proposed in this article is intended for fully digital PARs, and was not tested for suppressing spectral replicas. Further, it is likely that for other configurations of the wind turbine (e.g., oriented at an angle with respect to the radar), these spectral replicas will overlap significantly, making the WTC mitigation a very challenging problem. Therefore, mitigation of WTC using sub-array digital PARs is not addressed in this article.

B. FULLY DIGITAL PAR
The WTC simulation is now used to simulate the fully digital Horus radar. Horus is a polarimetric, mobile, S-band radar developed by the ARRC with support from NOAA [21]. The aperture has 32 × 32 active digital elements, as shown in Fig. 11, which results in a broadside half-power beamwdith of about 4 • . For simplicity, signals for 32 antenna elements are simulated (one dimensional beamforming) and used to produce the space-time spectrum.
The first case is for a wind turbine with blades rotating at 40 • s −1 , a turbine orientation of 60 • with respect to the radar, at 0.5 km from the radar (due to coarse beamwidth), with a CSR of 20 dB, and with weather spectra at −10 m s −1 (v/2v a ∼ −0.2). Other parameters are as described in Table 1 (e.g., M = 60 samples, f = 3 GHz). The simulation output is presented in Fig. 12. Panel (a) shows that the WTC has three linearly oriented spectral components, one for each blade. Further, a ground clutter component is observed near the 0 m s −1 Doppler, corresponding to the tower and the hub (stationary parts). Panel (b) shows the space-time spectra of the simulated weather signals. Panel (c) shows the contaminated received echo, obtained by adding the simulated complex echo voltages from the WTC and weather signals. Lastly, the algorithm output is presented in panel (d), where three deep notches (below −40 dB) are filtering the WTC, and an additional narrow filter is placed around the 0 m s −1 Doppler velocities to remove ground clutter returns. This indicates the algorithm is capable of identifying and removing space-time WTC components with minimal distortion on weather spectra.
The previous case was for wind turbines at close range. Next, we investigate the effect of wind turbines at 10 km from the radar. In this case, due to the angular resolution of the Horus system (∼4 • ) and the size of the simulated wind turbine, WTC contamination becomes narrow in the spatial domain. This is shown in Fig. 13, where most of the contamination from the wind turbine is near the normalized angle of 0. Panels are organized as in Fig. 12. It can be seen that although this proof-of-concept radar does not have very fine angular resolution, WTC can still be effectively mitigated using its fully digital capability. That is, even though the simulated wind turbine is relatively far with respect to the beamwidth, their contamination can still be removed from the received signal, as long as there are no spatial replicas. This key benefit of fully digital PAR technology will be critical in the future, as clean and renewable wind energy continues to grow. A performance analysis of this algorithm is carried out by systematically increasing the power of the WTC signals while keeping the power of the weather signals constant. The average power of the recovered weather signals is calculated and compared to the simulated (true) value. The estimated signal power as a function of CSR is shown in Fig. 14. It is noted that there is a small power bias (∼0.01 dB) for the CSR range from −10 dB to approximately 22 dB due to the application of the filter on the weather signal space-time spectrum. That is, the filter not only removes the wind-turbine contamination, but also part of the weather signal spectrum. Although the bias in weather-signal power estimates is nearly negligible up to about 22 dB of CSR, estimates produced by the proposed algorithm would still be acceptable at 40 dB CSR since the requirement for reflectivity bias (derived from signal power) set forth by NOAA is 1 dB (see page 35, row 23 in [42]). A study on the performance of the proposed algorithm on polarimetric estimates is left for future work.

VI. CONCLUSION
This article presents the first evaluation of applying Space-Time Adaptive Processing (STAP) on digital Phased Array Radars (PAR) to reduce the unwanted contamination introduced by wind turbines. This contamination affects weather radar signal processing algorithms (e.g., ground clutter filter) and QPE algorithms (e.g., rain rate estimation) in several ways such as misidentification of precipitation echoes, false mesocyclone detections, and incorrect storm cell identification. With the dramatic growth of the wind power industry in the United States, this could impact NWS forecasters' ability to issue accurate severe weather warnings and weather forecasts in the future. Furthermore, it would increase the likelihood that mobile radar field deployments contain WTC contamination, hindering research efforts that may rely on careful data analysis and interpretation.
Being PAR the candidate technology for the future weather radar network, evolutionary signal processing algorithms that make full use of their capabilities should be investigated to mitigate WTC contamination. In this article, a flexible wind turbine time-series signal simulator is developed to characterize the contamination in a controlled environment. An algorithm to remove the contamination is proposed and illustrated using simulated data. It was designed using knowledge gained from simulations on the structure of WTC space-time signatures. The mitigation algorithm exploits Gabor filters to detect the orientations from WTC-ridges in the space-time domain, and estimates the clutter covariance matrix. Different WTC configurations were simulated and used to evaluate the algorithm's performance. It was shown that the algorithm can effectively remove most WTC contamination from the time-series IQ data, with minimal impact on weather signals.
Two radar configurations were considered, one based on sub-array digital outputs and the other on fully digital elements. Results show that it is challenging to use space-time techniques on PARs based on sub-array digital architectures, because the physical spacing between digital phase centers exceeds half wavelength, creating grating lobes in the subarray pattern. These grating lobes create replicas of the WTC spectrum that may be harder to mitigate from weather signals. Since there are no grating lobes in fully digital arrays, this configuration is more beneficial for suppressing WTC returns and retrieving weather signals of interest. The algorithm performance was quantified as a function of CSR, showing a small bias (≤ 0.1 dB) for CSRs from −5 up to approximately 25 dB, and a bias of 1 dB at 40 dB CSR.
The simulation developed is for a single wind turbine in the resolution volume and does not consider challenges associated with real-time signal processing. Future efforts could extend this simulator to include multiple wind turbines in the resolution volume, a more practical situation for existing wind farms. The extension would involve replicating the wind turbine structure and placing copies of it (with slightly different parameters) in the near vicinity. Then, signals from the collection of scattering centers from all wind turbines can be added to produce signals received at the radars' digital elements. Real-time processing of STAP-based method is challenging to implement, especially as the size of the aperture increases. If wind turbines are sufficiently far from the radar, the mitigation algorithm could be simplified to deterministically remove clutter along the normalized zeroth angle and Doppler lines. This would take care of stationary ground clutter scatterers (as conventionally done in weather radars) and also spatially narrow moving scatterers, such as the blades of a wind turbine.
This article provides a new perspective on the mitigation of WTC capitalizing on the high flexibility of fully digital PARs. Although no actual digital PAR data were available at the time of this writing, future efforts include evaluating the proposed algorithm with actual observations. Time-series data will be collected using the mobile Horus radar once the system is ready for recording element-level IQ data. This will support refining the technique presented in this article, and establishing the first STAP-based mitigation technique to provide more accurate polarimetric weather estimates to forecasters and QPE algorithms. The long-term outcomes of this research shall provide information regarding the opportunities and challenges of fully digital PARs to improve the quality of polarimetric weather observations.