Abstract

The specific differential phase Kdp is one of the important parameters measured by dual-polarization radar that is being considered for the upgrade of the current Next Generation Weather Radar (NEXRAD) system. Estimation of the specific differential phase requires computing the derivative of range profiles of the differential propagation phase. The existence of possible phase wrapping, noise, and associated fluctuation in the differential propagation phase makes the evaluation of derivatives an unstable numerical process. In this paper, a robust algorithm is presented to estimate the specific differential phase, which is able to work on wrapped phases and keep up with the spatial gradients of rainfall, to provide a high-resolution specific differential phase.

1. Introduction

In dual-linear-polarization bases, the specific differential phase Kdp is defined as the slope of range profiles of the differential propagation phase shift Φdp between horizontal (H) and vertical (V) polarization states (Seliga and Bringi 1978; Jameson 1985; Bringi and Chandrasekar 2001). The specific differential phase is an important parameter for meteorological applications, because it is not affected by propagation attenuation and it is closely related to rain intensity, even in the presence of hail (Aydin et al. 1995; Chandrasekar et al. 2008).

However, the estimation of Kdp involves estimating the slope of Φdp profiles, which is known to be a noisy unstable computation. Evaluation of the derivative is essentially a high-pass filtering, and it expects a smooth and continuous function as the input. The range profile of total differential phase Ψdp contains both Φdp and differential backscatter phase shift δhv, as well as measurement fluctuation. The fluctuation in the estimates of Ψdp will be magnified during the differentiation, resulting in large variance in the estimates of Kdp. Furthermore, phase wrapping may exist in Φdp, and special care is generally required to unfold the wrapped phase profiles. The unambiguous range of Φdp usually is 180° in the alternate H/V transmission mode and 360° in the simultaneous H/V transmission mode. The system phase offset needs to be carefully adjusted to near the lower unambiguous limit for making use of the full unambiguous range. Nevertheless, for a long propagation path in rain or at higher-frequency bands such as X-band, Φdp values can easily exceed the unambiguous range. All these challenges have to be addressed to achieve a satisfactory estimation of Kdp, which can be applied for quantitative applications such as rainfall estimation. The artificial discontinuities introduced by phase wrapping have to be detected first from the noisy phase profiles and to be adjusted. Usually, a piecewise fitting is used to predict the local trend, and any phase sample deviated far from this trend can be attributed to phase wrapping. This solution is neither elegant nor stable because of the statistical fluctuations present in the Φdp field. Indeed, it is a common unfolding technique to estimate the slope first and then recover the intrinsic phase profiles based on the detected discontinuities. In this paper, a completely different approach to estimate Kdp is introduced to avoid the need for unfolding phases. The current techniques to reduce the variance of Kdp estimates include range filtering, linear fitting, or both (Golestani et al. 1989; Hubbert and Bringi 1995). All these techniques reduce the peak Kdp values in the estimation that could introduce biases. In addition, only limited adaptivity can be achieved with filtering or fitting to follow the steep slopes within intense rain cells and reduce the estimation variance in the rest of the segments simultaneously. This paper presents an adaptive scheme to estimate Kdp, which has better range resolution in intense rain cells to capture the small-scale variability. A regularization technique is introduced to control the balance between estimation bias and variance, and it incorporates the adaptivity to keep up with steep change of Φdp. The procedure is developed essentially to improve the robustness in real-time operations. This paper is organized as follows: the conventional estimation of Kdp is summarized first in section 2 with the state-of-the-art Colorado State University–University of Chicago–Illinois State Water Survey (CSU–CHILL) algorithm as an example. Then, the new estimator is described in section 3 to deal with wrapped phases, and an adaptive algorithm for its implementation is developed in section 4. Its application to radar observations from CSU–CHILL S-band radar and the Center for Collaborative Adaptive Sensing of the Atmosphere (CASA) Integrated Project 1 (IP1) X-band radars is presented to demonstrate the performance of the new algorithm.

2. Conventional estimation of specific differential phase

The specific differential phase Kdp is estimated as the range derivative of the differential propagation phase Φdp, following the fundamental definition

 
formula

Radar does not measure the forward propagation parameter Φdp directly. Instead, the total differential phase Ψdp is estimated from copolar covariance that consists of both the phase shifts resulting from forward propagation and backscattering; that is,

 
formula

If the differential backscattering phase δhv is relatively constant or negligible, the profile of Ψdp can be used to estimate Kdp.

Because of different pulse sequencing, separate estimators are used to estimate Ψdp in the alternative transmission mode (VH) and in the simultaneous transmission mode (VHS). In the VH mode, Ψdp can be estimated from the received alternate samples at H and V polarization states as (Mueller 1984; Sachidananda and Zrnić 1989; Bringi and Chandrasekar 2001)

 
formula

where

 
formula
 
formula

where H2i+1 and V2i are the alternately sampled horizontal and vertical polarization returns, respectively. Note that the V polarization is assumed leading the pulse-to-pulse polarization switching on transmission. In the VHS mode, Ψdp is estimated directly from the simultaneous samples Vi and Hi as

 
formula

The argument operation in (3) and (5) leads to finite unambiguous ranges of 180° and 360° in the VH and VHS modes, respectively. As a cumulative value, Φdp is a continuous range function. However, discontinuities usually appear in the estimated Ψdp because of phase wrapping, statistical fluctuations in estimation, and gate-to-gate variation of δhv. The artifact introduced by phase wrapping needs to be identified and corrected for, whereas the statistical errors and gate-to-gate variation should be then suppressed using various smoothing techniques. The error model is presented to describe the estimated Ψdp in term of Kdp as

 
formula

where Φdp(0) represents the system differential phase and ɛ represents fluctuations resulting from statistical errors and the variation of δhv. In the presence of significant δhv over a short range such as a “bump,” Hubbert and Bringi (1995) introduced an iterative technique to detect and remove the extra change resulting from the δhv bumps.

a. Phase unfolding

Phase wrapping can occur, depending on the system differential phase Φdp(0), even if the total Φdp accumulation is within the unambiguous range. For this reason, the system differential phase usually is intentionally set near the lower bound of the unambiguous range. Nevertheless, this system phase can drift, and Φdp can exceed the unambiguous range. Therefore, phase unfolding needs to be done first at the very beginning during the evaluation. In CSU–CHILL, a set of straightforward but complex logic is adopted to detect the wrapped phase and accordingly unfold the range profile of Φdp. The procedure is described in the following.

The first step is to determine the location and phase offset where the meaningful Φdp profile begins. Within the path free of weather echoes, low signal-to-noise ratio (SNR) and low copolar correlation result in large phase fluctuation, which should not be misjudged as phase wrapping. This decision is made based on the standard deviation of Φdp (σϕ) and the copolar correlation coefficient ρhv. A consensus of small σϕ and high ρhv over consecutive gates is used to ensure robustness in the detection of weather echoes.

Subsequently, σϕ and the slope of Φdp are continuously computed. A segment of profile with large σϕ is labeled as nonweather echoes. Otherwise, the estimated slope is accumulated into a smoothed Φdp profile for reference. Based on this knowledge, the current Φdp value is compared against the reference value of preceding Φdp to determine whether the current one needs to be unfolded.

The flowchart describing the current CSU–CHILL algorithm is shown in Fig. 1. The computation for σϕ and the slope is done to prevent false phase unfolding resulting from the contamination by noise. It should be noted that hard thresholds are used throughout the procedure that are obtained from the system metrics of the CSU–CHILL radar. It is also assumed that Φdp would not be wrapped at the beginning range.

Fig. 1.

Flowchart of the phase unfolding algorithm used in the CSU–CHILL radar. The judgments and adjustments are valid in VH mode.

Fig. 1.

Flowchart of the phase unfolding algorithm used in the CSU–CHILL radar. The judgments and adjustments are valid in VH mode.

b. Range filtering

The differentiation process is equivalent to applying high-pass filtering on the range profiles of Φdp, as seen in the frequency domain

 
formula

where Kdp(ω) and Φdp(ω) are Fourier transforms of range profiles of Kdp(r) and Φdp(r), respectively. Therefore, statistical errors in the estimated Ψdp will be magnified in the estimated Kdp profile. If Kdp is estimated using finite difference as

 
formula

then the standard deviation of the estimate of Kdp is related to the phase variation through (Bringi and Chandrasekar 2001)

 
formula

assuming uniform differential phase variance across range. However, the differential phase variance depends on Doppler spectrum, ρhv, and δhv. Low-pass range filtering can be further employed to reduce the estimation variance. To gain insight into this filtering process, the error term in (6) is assumed as a white, ergodic, wide-sense stationary process. Using a filter as specified by F(ω), the variance of the estimates of Ψdp after filtering can be written as

 
formula

The variance reduction is approximately proportional to the normalized cutoff frequency of the filter. The cutoff frequency can be translated to a range scale. At gate spacing Δr, the Nyquist frequency can be written as 1/(2Δr). Assuming the filter is expected to allow the trend of Φdp down to a range scale of rc, the normalized cutoff frequency for the passband of the filter must be (Oppenheim and Schafer 1989)

 
formula

Then, the fluctuations at scales smaller than rc will be suppressed. A good design of the range filter should suppress gate-to-gate fluctuations but preserve the underlying trend of Φdp (Hubbert et al. 1993). In the CSU–CHILL algorithm, the cutoff scale is chosen at 3 km, and iterative filtering is used to adaptively compensate for the persistence of large δhv. Note that at different gate spacing the filter needs to be redesigned as shown in appendix A.

c. Least square fitting

Polynomial regression can be applied as an alternate estimator to finite difference (Bringi and Chandrasekar 2001), where the range profile of Ψdp is approximated by a polynomial function as

 
formula

As the order of this polynomial function increases, the regression can better capture the intrinsic change of Φdp. However, high-order regression will also increase the computational complexity and the estimation sensitivity. In practice, low-order piecewise polynomial fitting is rather performed.

For short paths, linear approximation can be used, and the estimate of Kdp is given by β1. The standard deviation of estimated Kdp in this case can be derived as (Gorgucci et al. 1999)

 
formula

where N is the number of consecutive gates for the piecewise linear regression. The estimation variance can by reduced by increasing pathlength (i.e., NΔr) or decreasing the gate spacing. In the current CSU–CHILL algorithm, piecewise linear regression is done to estimate Kdp using adaptive length N, which is adjusted according to three reflectivity levels, such as

 
formula

This adaptive adjustment is meant to follow the steep phase changes in intensive precipitation regions while keeping the variation of estimated Kdp low in light precipitation regions. The regression length is tailored for the CSU–CHILL radar, which has a range resolution of 150 m. Therefore, the estimation can keep up with a small-scale variability up to 1.5 km in heavy rain and a large-scale variability up to 4.5 km in light rain.

3. New estimator for specific differential phase

A practical phase unfolding approach has been described in the previous section. One feature of the procedure is somewhat complicated in the sense that a slope is first approximated for unfolding the phase profiles and then a more accurate slope is estimated for Kdp. Its performance also relies on the choice of several hard thresholds. Nevertheless, wrapped phase profiles have to be unfolded first before the subsequent filtering and fitting, because the discontinuities are not “physical” features. In this section, a new estimator will be proposed for estimating Kdp based on the complex-valued range profiles of the differential phase shift. To clarify the difference from the conventional description, hereafter the differential phase shift in the complex domain will be referred as the angular Φdp, whereas its counterpart in the real domain (e.g., from 0 to 2π) will be referred as the principal Φdp.

For example, Φdp repeats in periods of 2π in the VHS mode. When the profile crosses a full period, even though the principal phase wraps back with an abrupt jump, locally the profile is still continuous in the complex domain. The angular Φdp profile can be obtained by defining a profile of unit vector as

 
formula

In the VH mode, the period is only π and the angular Φdp in this case is known as an axial variable, which can be transformed to circular variable by magnifying the angle by a factor of 2. Again, the phase profile can be vectorized as

 
formula

The only difference between (15) and (16) is the factor 2, which needs to be applied back to the final estimates of Kdp. Without loss of generality, the following discussion will be based on the profile as defined in (15). The same algorithm works for the VH mode as well.

Taking the derivative of (15) with respect to r, we have

 
formula

Therefore, the estimate of Kdp can be written as

 
formula

This shows that the quotient ϑ′(r)/ϑ(r) is an imaginary number, provided the continuity in the complex-valued phase profile and its imaginary part corresponds to Kdp with an extra term associated with estimation error. On the other hand, the existence of phase fluctuations can disturb the continuity. In addition, the numerical errors in derivatives make the quotient likely not a pure imaginary number. If the phase fluctuation is sufficiently small, as implied earlier in (8), the imaginary part of the quotient can be taken as the estimates of Kdp.

The statistical fluctuations of the differential phase shift can be also studied through the angular Φdp. Its mean and its dispersion can be defined as

 
formula
 
formula

respectively, which satisfy the following properties under a constant shift by ϕ0:

 
formula
 
formula

It can be shown that the mean defined in (19) is the same as the mean of unfolded principal Φdp if it follows a symmetric distribution. The dispersion defined in (20) decreases as the standard deviation of Φdp increases, and it approaches unity if Φdp is highly concentrated, whereas it approaches 0 if Φdp is uniformly distributed between 0 and 2π. Specially, if Φdp is a narrowly distributed Gaussian variable, it can be shown that

 
formula

Under a Gaussian distribution assumption, the dispersion can be used to examine the degree of fluctuation in the estimated Ψdp and predict the data quality for the estimates of Kdp. Figure 2 illustrates the simple logic for data masking with use of the dispersion parameter.

Fig. 2.

Flowchart of decision for rain cells and good data mask, where MG and MB are constants set for consecutive “good” and “bad” samples, respectively. Usually, MB is smaller than MG.

Fig. 2.

Flowchart of decision for rain cells and good data mask, where MG and MB are constants set for consecutive “good” and “bad” samples, respectively. Usually, MB is smaller than MG.

Numerical derivative errors can be reduced by using high-order estimation techniques. Linear finite difference is a genuine slope estimator for a linear function. However, the definition in (15) and (16) does require higher-order numerical differentiation because of the exponential term, unless the sampling interval is sufficiently small. As a simple example, Table 1 shows the numerical derivative errors in the angular Φdp profiles using different solutions. The principal Φdp profile is assumed linear from which all these estimators can accurately retrieve the true slope. However, as the slope increases, finite difference gives larger errors in the estimates of Kdp from the angular Φdp profiles. Numerical differentiation based on interpolation polynomials of a higher degree can substantially reduce the errors, simply because the high-degree polynomials can better approximate the underlying function. Note that cubic spline is also a piecewise interpolation polynomial with only extra constraint of continuities at the knots. It can be seen that cubic spline performs uniformly better than a four-degree Lagrange interpolation polynomial. Therefore, a feasible solution for Kdp estimates will be filtering angular Φdp profiles based on angular Ψdp statistics and estimating Kdp using cubic spline. Thus, phase unfolding is avoided, and the accuracy is also ensured. An adaptive estimation of Kdp is developed based on the principles described and is presented in the next section.

Table 1.

Numerical derivative errors in Kdp estimates from angular Φdp profile. The range spacing is assumed at 0.5 km, and the intrinsic Kdp is assumed constant. The derivative using higher-order polynomials is computed from piecewise Lagrangian interpolation polynomials.

Numerical derivative errors in Kdp estimates from angular Φdp profile. The range spacing is assumed at 0.5 km, and the intrinsic Kdp is assumed constant. The derivative using higher-order polynomials is computed from piecewise Lagrangian interpolation polynomials.
Numerical derivative errors in Kdp estimates from angular Φdp profile. The range spacing is assumed at 0.5 km, and the intrinsic Kdp is assumed constant. The derivative using higher-order polynomials is computed from piecewise Lagrangian interpolation polynomials.

4. Adaptive estimation of specific differential phase

The instability in numerical differentiation of noisy data can be well controlled in a regularization framework. To reach stable solutions, the differentiation of an arbitrary function f (·) (which is the range profile of angular Φdp in this paper) is to find a function s(·) to minimize both the smoothness of the underlying regression function and the regression errors (Reinsch 1967):

 
formula

where f (rk) represents the observations and λ is the Lagrangian parameter. The first term in (24) evaluates the degree of smoothness. An extreme smoothness occurs when λ = 0, in which case a linear continuous function is the solution to the minimization. The second term in (24) assesses the degree of data fidelity. As λ → ∞, the regression function is dictated to pass all the measurement. It can be shown that the cubic spline function has the greatest smoothness over all second-order continuous functions that pass the nodes (Pollock 1999). Other intermediate values of λ control the trade-off between the expected model smoothness and the allowed regression errors. The minimization will impose more emphasis on the degree of smoothness if small λ is used. Through the choice of λ, the solution balances between the bias and variance of the estimates.

Assuming cubic splines as the regression function, the polynomial function for the kth interval between nodes rk and rk+1 can be written as

 
formula

Normally, Kdp can be assumed as constant at both ends of a radial profile over a rain cell, which means that the second-order derivatives at end nodes should be zero. This end condition imposes the complete determination of all the coefficients and leads to “natural cubic splines.” On the cubic splines solved, the coefficient dk corresponds to the smoothed angular Φdp and the coefficient ck corresponds to its range derivatives. Note that only the derivatives at the nodes are of our interest to estimate Kdp. Therefore, the estimates can be simply obtained as

 
formula

Thus, the conventional filtering process is avoided in this regularization framework, and stable numerical derivative can be obtained with an appropriate λ. However, cubic smoothing spline still is largely a global smoothing procedure. Smoothing the segments of high variation at a satisfactory level leads to an oversmoothed profile in the segments of low variation. In addition, oversmoothness needs to be avoided in high Kdp region to track the peak Kdp. The corresponding adaptivity is enabled by introducing varying weights in the minimization (24) as

 
formula

Substituting (25) into (27), the objective function for minimization can be rewritten as

 
formula

Its solution is provided in appendix B. The new introduced parameters q and w are both real-valued scaling factors. With them, the balance between bias and variance can be selectively controlled for different segments. The former one will be used to compensate for the variation in slope (i.e., Kdp); the latter one will be used to compensate for the statistical fluctuation of the input profiles.

a. Adaptivity for statistical fluctuation

To minimize the second term in (28), the smoothed profile dk should be equal to the mean of f (rk). Therefore, w can be chosen as the inverse of the standard deviation of the angular Φdp such that the variation is normalized. The variance of angular exp{ jΦdp} can be evaluated as

 
formula

The maximum variation occurs when noise dominates, in which case Φdp distributes uniformly between −π and π. It then follows that

 
formula

When the variation in the angular Φdp is small, Taylor expansion can be used to approximate (29) as

 
formula

This shows that the variance of the angular Φdp is approximately equal to the variance of the principal Φdp. The variance σϕ2 can be evaluated from the Doppler spectrum parameters (Bringi and Chandrasekar 2001). It depends both on the Doppler spread and the copolar correlation coefficient. For simplicity, it can be approximated for large ρhv as

 
formula

Note that the approximations in (31) and (32) are only valid for small variation or large ρhv. However, the weighting function (32) captures well the overall relation between the variation of angular Φdp and the copolar correlation coefficient, noting that it also satisfies the maximum variation condition (30) as the correlation becomes poor. Nonnegligible δhv can result from melting hydrometeors (Hubbert and Bringi 1995). In this case, the radar signatures will have lower ρhv values and the adaptive method will smooth the differential phase profiles more over the bumped segments. It should also be noted that the iterative filtering technique developed in Hubbert and Bringi (1995) can be additionally applied, if the user desires.

b. Adaptivity for Kdp

To follow the profiles in large Kdp, larger variation shall be allowed such that bias is minimized to avoid excessive smoothness. Therefore, let q−1 = 2Kdp, where the factor 2 is in place to compensate the two-way Φdp profile. With this adjustment, the smoothness requirement at large Kdp will be less emphasized. On the contrary, the smoothness at small Kdp will be heavily weighted. The scaling factor q−1 cuts off at a Kdp value of 0.1° km−1 for numerical stability.

Based on the angular Φdp profile, the smoothness term in (27) can be evaluated as

 
formula

One should note that this term (33) is not unitless after the scaling factor q is applied. It is still related to the actual value of slope Kdp owing to the translation from the principal Φdp to the angular Φdp.

c. Choice of the Lagrangian parameter

The adaptive approach is only left with the freedom in choosing the smoothing factor λ. Cross-validation can be used to find the optimal λ from a dataset (Craven and Wahba 1979). Such determination does not suit this application because of excessive ray-by-ray variation in radar observations. However, with the scaling scheme, a simple procedure is used to make the evaluation of smoothness and that of data fidelity equally important. The variation of angular Φdp has been normalized using w; therefore, a reasonable choice of λ should be close to the expected smoothness by evaluating (33). As an approximation, constant Kdp can be assumed in each range interval. At different Kdp levels, the smoothness is quantitatively assessed and shown in Table 2 by allowing a 10% variation in Kdp. The last column can be used as a baseline for choosing λ. The existence of extra Kdp in (33) requires a different smoothing factor for different peak Kdp values. For example, targeted at Kdp of 30° km−1, the smoothing factor should be around 1.1Δr.

Table 2.

Approximate quantitative evaluation of the smoothness of the cubic regression splines in single interval.

Approximate quantitative evaluation of the smoothness of the cubic regression splines in single interval.
Approximate quantitative evaluation of the smoothness of the cubic regression splines in single interval.

Figure 3 shows a simulated Kdp profile, composed of a Gaussian function at 15 km, to demonstrate an infinitely differentiable profile, and a triangular function at 35 km, to demonstrate Kdp discontinuities. The profile exhibits sharp gradient on Kdp, which changes from 0° to 30° km−1 within a 3-km range. Statistical fluctuation is introduced to Φdp by setting the intrinsic ρhv at 0.99, the normalized spectral width at 0.12, and the integration length at 64 samples. The Lagrangian factor is chosen between 0.1Δr and 1.1Δr for comparison. The conventional Kdp estimation is also plotted from the filtered Φdp profiles with a six-order finite impulse response (FIR) filter (designed to keep the variation at 3-km scale). Near the Kdp peaks, it can be seen that at λ of 1.1Δr the estimation biases are almost completely removed and the standard deviation of the estimation increases but less than 1° km−1. In the nonstorm regions, both the biases and the standard deviation of the estimation are reduced. With smaller λ, oversmoothed Kdp estimates are obtained that exhibit large biases. It should be noted that the λ value is not guaranteed to be the optimal, but Table 2 can be used as a baseline to empirically pick an appropriate λ. The value at 1.1Δr is used for all the case studies that are presented in the next section.

Fig. 3.

Kdp estimation from principal Φdp (black dashed line) using simple finite difference and from angular Φdp using adaptive cubic smoothing spline: (a) the mean and (b) standard deviation of estimated Kdp based on 100 realizations. The intrinsic Kdp profile (thick gray line) is simulated as (left) a narrow Gaussian function and (right) a triangular function, at gate spacing of 0.5 km. Fluctuation resulting from noise, Doppler spread, and copolar correlation is simulated.

Fig. 3.

Kdp estimation from principal Φdp (black dashed line) using simple finite difference and from angular Φdp using adaptive cubic smoothing spline: (a) the mean and (b) standard deviation of estimated Kdp based on 100 realizations. The intrinsic Kdp profile (thick gray line) is simulated as (left) a narrow Gaussian function and (right) a triangular function, at gate spacing of 0.5 km. Fluctuation resulting from noise, Doppler spread, and copolar correlation is simulated.

The algorithm can be easily implemented in a real-time environment. Even though matrix inversion and matrix multiplication are involved in the solution, the computation complexity is well within the capability of modern-day computation power because all the matrices are band limited. Most importantly, the calculation for phase unfolding is completely avoided.

5. Evaluation

In this section, a simulated storm and actual radar observations are used to evaluate and demonstrate the new Kdp estimation. The scaling as shown in (33) cannot be applied directly in practice. However, the Kdp in (33) serves only as a relative weighting factor and does not need to be exact. An approximation of Kdp can be first translated from the measured reflectivity using the following power-law relation:

 
formula

which can be obtained assuming a given raindrop size and shape distribution. This conversion generally works very well at S band. At shorter wavelengths, rain attenuation on the measured reflectivity would limit its applicability. Hail contamination will be another factor violating the conversion.

To solve these problems, a two-step iteration is adopted, in which a rough estimation of Kdp is made first using the nonadaptive cubic splines derived from (24) and then it is used in the adaptive cubic splines to perform the final Kdp estimation. This introduces a slight increase of computation time without any perceivable impact on real-time implementation. During the first step, a smaller Lagrangian parameter fixed at 0.1 is used to obtain the overall trend of the Kdp profile and the resulting scaling functions very well. This two-step iteration technique is adopted in this paper to estimate Kdp from radar observations.

a. Simulated storm

Chandrasekar et al. (2006) have developed a procedure to simulated X-band radar observation from S-band radar data with the natural variation of microphysical parameters preserved. X-band radar observations are simulated using the CSU–CHILL radar observations of a thunderstorm during the Severe Thunderstorm Electrification and Precipitation Study (STEPS). The simulated X-band observations are shown in Fig. 4, with Kdp of up to 7° km−1. Figures 4a–c present the radar reflectivity, the estimated Kdp using the conventional filtering and fitting approach, and the estimated Kdp using the new adaptive approach, respectively. Both the conventional approach and the adaptive approach result in a similar Kdp field, with the adaptive approach showing less “range streaks.” The adaptive approach presents better resolution, which captures the Kdp peaks as well as the storm texture that is exhibited in the radar reflectivity. In addition, the Kdp field is much cleaner with smaller variation in the regions of low radar reflectivities. Because the “true” radar parameters are available in the simulation, the scatterplot between the estimated Kdp and the intrinsic Kdp is shown in Fig. 4d. If the intrinsic Kdp is small, the new estimates are more stable; if the intrinsic Kdp is large, the new estimates have a larger standard deviation, but the mean bias is smaller. The results of Fig. 4 clearly illustrate the advantage of the new adaptive approach.

Fig. 4.

Simulated X-band radar observations at elevation angle of 0.5° and gate spacing of 0.1 km. (a) PPI plot of the radar reflectivity; (b) PPI plot of the estimated Kdp using conventional approach on the principal Φdp; (c) PPI plot of the estimated Kdp using the adaptive estimation on the angular Φdp; and (d) the scattering plot between the estimated Kdp and the intrinsic Kdp, where the black points illustrate the conventional estimates and the gray points illustrate the adaptive estimates.

Fig. 4.

Simulated X-band radar observations at elevation angle of 0.5° and gate spacing of 0.1 km. (a) PPI plot of the radar reflectivity; (b) PPI plot of the estimated Kdp using conventional approach on the principal Φdp; (c) PPI plot of the estimated Kdp using the adaptive estimation on the angular Φdp; and (d) the scattering plot between the estimated Kdp and the intrinsic Kdp, where the black points illustrate the conventional estimates and the gray points illustrate the adaptive estimates.

b. CSU–CHILL S-band radar observation

Figure 5 shows the S-band observations of a heavy rain event on 31 May 2003 UTC from the CSU–CHILL radar. A large area of thunderstorms developed southeast of the CSU–CHILL radar and some of the beams were continuously filled with heavy rain over range segments of nearly 50 km. The value of Φdp is accumulated over 200° around an azimuth of 128°, where phase is wrapped. Again, it can be found that both approaches give similar Kdp estimates in the overall pattern by comparing Figs. 5b,c. The Kdp peaks are well estimated by the adaptive estimator in the storm core, whereas the Kdp estimates outside of storm core are fairly smooth compared to these by the conventional estimator. Note that the storm is located at range of 50 km and beyond. The spatial resolution there is much lower than the previous simulated cases. The scatterplot between the estimated Kdp and radar reflectivity is presented in Fig. 5d. It can be seen once again that the large number of extraneous points from the conventional algorithm are eliminated.

Fig. 5.

CSU–CHILL S-band radar observations at elevation angle of 0.5° and gate spacing of 0.15 km, at 0059 UTC 31 May 2003, Greeley, Colorado. (a) PPI plot of the attenuated reflectivity; (b) PPI plot of the estimated Kdp using conventional approach on the principal Φdp; (c) PPI plot of the estimated Kdp using the adaptive estimation on the angular Φdp; and (d) the scattering plot between the estimated Kdp and the attenuation-corrected Zh, where the black points illustrate the conventional estimates and the gray points illustrate the adaptive estimates.

Fig. 5.

CSU–CHILL S-band radar observations at elevation angle of 0.5° and gate spacing of 0.15 km, at 0059 UTC 31 May 2003, Greeley, Colorado. (a) PPI plot of the attenuated reflectivity; (b) PPI plot of the estimated Kdp using conventional approach on the principal Φdp; (c) PPI plot of the estimated Kdp using the adaptive estimation on the angular Φdp; and (d) the scattering plot between the estimated Kdp and the attenuation-corrected Zh, where the black points illustrate the conventional estimates and the gray points illustrate the adaptive estimates.

c. CASA X-band radar observation

On 14 June 2007, an intense rainfall event developed over the CASA IP1 test bed in Oklahoma (Chandrasekar et al. 2007) for which a flash flood warning was issued by the National Weather Service (NWS). Reflectivity observations by the X-band IP1 radar in Rush Springs are shown in Fig. 6a. The estimated Kdp from the principal Φdp is shown in Fig. 6b, and the new Kdp estimates are shown in Fig. 6c. Visually, these two estimators give a comparable Kdp field in terms of the overall pattern and intensity. For this case, the radar has a range resolution of 48 m, and the conventional approach has better estimates for large Kdp observations than the case of lower range resolutions. However, the Kdp field in Fig. 6c matches the storm structure much better, and the negative Kdp in Fig. 6b is largely eliminated. The improved quality of the new adaptive Kdp is also shown in the KdpZh scatterplot in Fig. 6d.

Fig. 6.

CASA X-band radar observations at elevation angle of 2.0° and gate spacing of 0.048 km, at 0708 UTC 14 Jun 2007, in Rush Springs, Oklahoma. (a) PPI plot of the attenuated reflectivity; (b) PPI plot of the estimated Kdp using conventional approach on the principal Φdp; (c) PPI plot of the estimated Kdp using the adaptive estimation on the angular Φdp; and (d) the scattering plot between the estimated Kdp and the attenuation-corrected Zh, where the black points illustrate the conventional estimates and the gray points illustrate the adaptive estimates.

Fig. 6.

CASA X-band radar observations at elevation angle of 2.0° and gate spacing of 0.048 km, at 0708 UTC 14 Jun 2007, in Rush Springs, Oklahoma. (a) PPI plot of the attenuated reflectivity; (b) PPI plot of the estimated Kdp using conventional approach on the principal Φdp; (c) PPI plot of the estimated Kdp using the adaptive estimation on the angular Φdp; and (d) the scattering plot between the estimated Kdp and the attenuation-corrected Zh, where the black points illustrate the conventional estimates and the gray points illustrate the adaptive estimates.

6. Discussion and summary

Challenges in the estimation of the specific differential propagation phase lie in the existence of phase wrapping and statistical fluctuations in the range profiles of the differential propagation phase, because differentiation is involved. Even though practical unfolding approach exists, the current implementation involves several ad hoc adjustments. Nevertheless, this procedure has been very useful so far and has been extensively used at CHILL. Phase wrapping is in fact an artifact coming from mapping an angular variable along unit circle into the real axis. In this paper, the estimation is shifted to the complex-valued range profile of the differential propagation phase exponentials. This yields a simpler estimation that is inherently immune to phase wrapping problem. Its numerical evaluation requires better accuracy as a result of the exponential form function. To this end, a regularization framework was developed, and a cubic smoothing spline was adopted to fit a better regression function.

The specific differential phase is an important parameter derived from dual-polarization radars. Using the specific differential phase can improve the accuracy of quantitative rainfall estimation. However, its estimates are highly susceptible to large variance in light and moderate rain, where Kdp is intrinsically small; on the other hand, its estimates are likely underestimated in heavy rain, where it is intrinsically large. These two factors could deter its wider use in quantitative rainfall mapping. A new algorithm is introduced to incorporate adaptivity through scaling the regression errors for estimation of Kdp.

Bias and variance are both typically a dilemma in inferring estimation from noisy data. During filtering, the bias and variance can be controlled through the cutoff frequency of the filter. As the sampling interval changes, new filter design needs to be in place. In the new adaptive algorithm, the bias and variance are controlled through a single Lagrangian parameter, or smoothing factor, which is also the only control parameter bounded. A simple method is proposed based on equilibrium principle to choose the smoothing factor. The proposed algorithm is well suitable for real-time implementation, primarily because of the banded matrix structure in the governing equations. It takes the same order of computation time as the conventional filtering or fitting procedure. The adaptive technique has been implemented in the CSU–CHILL radar, and it is expected to find importance in operational applications (Chandrasekar et al. 2008).

Fig. A1. Frequency responses of FIR range filters at various gate spacings.

Fig. A1. Frequency responses of FIR range filters at various gate spacings.

Table A1. Specification of FIR range filters to pass Φdp trend at 3-km scale.

Table A1. Specification of FIR range filters to pass Φdp trend at 3-km scale.
Table A1. Specification of FIR range filters to pass Φdp trend at 3-km scale.

Acknowledgments

This research is supported by NASA via the Precipitation Measurement Mission (PMM) program and by Colorado State University. The authors acknowledge the opportunity to test the algorithm for CASA radar data. The CSU-CHILL radar facility is supported by the National Science Foundation through ATM-0118021.

REFERENCES

REFERENCES
Aydin
,
K.
,
V. N.
Bringi
, and
L.
Liu
,
1995
:
Rain-rate estimation in the presence of hail using S-band specific differential phase and other radar parameters.
J. Appl. Meteor.
,
34
,
404
410
.
Bringi
,
V. N.
, and
V.
Chandrasekar
,
2001
:
Polarimetric Doppler Weather Radar: Principles and Applications.
Cambridge University Press, 662 pp
.
Chandrasekar
,
V.
,
S.
Lim
, and
E.
Gorgucci
,
2006
:
Simulation of X-band rainfall observations from S-band radar data.
J. Atmos. Oceanic Technol.
,
23
,
1195
1205
.
Chandrasekar
,
V.
,
D. J.
McLaughlin
,
J.
Brotzge
,
M.
Zink
,
B.
Philips
, and
Y.
Wang
,
2007
:
Distributed collaborative adaptive radar network: The CASA IP-1 network and tornado observations.
Preprints, 33rd Conf. Radar Meteor., Cairns, Australia, Amer. Meteor. Soc., 13A.3A
.
Chandrasekar
,
V.
,
A.
Hou
,
E.
Smith
,
V. N.
Bringi
,
S. A.
Rutledge
,
E.
Gorgucci
,
W. A.
Petersen
, and
G. S.
Jackson
,
2008
:
Potential role of dual- polarization radar in the validation of satellite precipitation measurements: Rationale and opportunities.
Bull. Amer. Meteor. Soc.
,
89
,
1127
1145
.
Craven
,
P.
, and
G.
Wahba
,
1979
:
Smoothing noisy data with spline functions.
Numer. Math.
,
31
,
377
403
.
Golestani
,
Y.
,
V.
Chandrasekar
, and
V. N.
Bringi
,
1989
:
Intercomparison of multiparameter radar measurements.
Preprints, 24th Conf. Radar Meteorology, Tallahassee, FL, Amer. Meteor. Soc., 309–314
.
Gorgucci
,
E.
,
G.
Scarchilli
, and
V.
Chandrasekar
,
1999
:
Specific differential phase estimation in the presence of nonuniform rainfall medium along the path.
J. Atmos. Oceanic Technol.
,
16
,
1690
1697
.
Hubbert
,
J.
, and
V. N.
Bringi
,
1995
:
An iterative filtering technique for the analysis of copolar differential phase and dual-frequency radar measurements.
J. Atmos. Oceanic Technol.
,
12
,
643
648
.
Hubbert
,
J.
,
V. N.
Bringi
, and
V.
Chandrasekar
,
1993
:
Processing and interpretation of coherent dual-polarized radar measurements.
J. Atmos. Oceanic Technol.
,
10
,
155
164
.
Jameson
,
A.
,
1985
:
Microphysical interpretation of multiparameter radar measurements in rain. Part III: Interpretation and measurement of propagation differential phase shift between orthogonal linear polarizations.
J. Atmos. Sci.
,
42
,
607
614
.
Mueller
,
E. A.
,
1984
:
Calculation procedure for differential propagation phase shift.
Preprints, 22nd Conference on Radar Meteorology, Zurich, Switzerland, Amer. Meteor. Soc., 10–13
.
Oppenheim
,
A. V.
, and
R. W.
Schafer
,
1989
:
Discrete-Time Signal Processing.
Prentice Hall, 870 pp
.
Pollock
,
D. S. G.
,
1999
:
A Handbook of Time-Series Analysis, Signal Processing and Dynamics.
Academic Press, 782 pp
.
Reinsch
,
C. H.
,
1967
:
Smoothing by spline functions.
Numer. Math.
,
10
,
177
183
.
Sachidananda
,
M.
, and
D. S.
Zrnić
,
1989
:
Efficient processing of alternately polarized radar signals.
J. Atmos. Oceanic Technol.
,
6
,
173
181
.
Seliga
,
T. A.
, and
V. N.
Bringi
,
1978
:
Differential reflectivity and differential phase shift: Applications in radar meteorology.
Radio Sci.
,
13
,
271
275
.

APPENDIX A

Range Filter Design for Differential Propagation Phase

The Nyquist sampling rate in range is directly related to the along-range gate spacing. As shown in (11), the cutoff frequency depends on the gate spacing and the required cutoff scale of range. Because of these changes, the range filter needs to be redesigned for different gate spacings. As Δr decreases, the normalized passband width becomes narrower and more samples are available for the same pathlength. Using an FIR filter, its order can be correspondingly increased. The specifications of the FIR filters used in this paper are listed in Table A1, and their frequency responses are plotted in Fig. A1.

APPENDIX B

Solution of the Cubic Smoothing Spline

The cubic splines are required to be continuous in the intermediate nodes; therefore,

 
formula

where

 
formula

Additionally, their first-order and second-order derivatives are also required to be continuous, which leads to

 
formula
 
formula

respectively.

The following define vectors and matrices as

 
formula
 
formula
 
formula
 
formula
 
formula
 
formula

The continuity constraints in (B1), (B3), and (B4) can be combined in matrix form as (Pollock 1999)

 
formula

where 𝗠 is an invertible square matrix. Similarly, the scaled objective function (28) used for adaptive smoothing can be written as

 
formula

Substitute (B11) into (B12) to express the objection function in term of d, and the minimization can be solved by setting the derivative with respect to d to zero; that is,

 
formula

Using (B11) in (B13), we can further derive

 
formula
 
formula

Therefore, the coefficients b and d can be solved as

 
formula
 
formula

The solution to the cubic smoothing spline as defined in (24) can be easily found by setting 𝗪 as an identity matrix and 𝗠q = 𝗠; that is,

 
formula
 
formula

With b and d available, a can be solved from (B4). Replacing a in (B1) using (B4) arrives the solution to c as

 
formula

Footnotes

Corresponding author address: Yanting Wang, 1373 Campus Delivery, Colorado State University, Fort Collins, CO 80523. Email: ytwang@engr.colostate.edu