## Abstract

The specific differential phase *K*_{dp} 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 *K*_{dp} 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 *K*_{dp} 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 *K*_{dp}. 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 *K*_{dp}, 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 *K*_{dp} is introduced to avoid the need for unfolding phases. The current techniques to reduce the variance of *K*_{dp} estimates include range filtering, linear fitting, or both (Golestani et al. 1989; Hubbert and Bringi 1995). All these techniques reduce the peak *K*_{dp} 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 *K*_{dp}, 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 *K*_{dp} 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 *K*_{dp} is estimated as the range derivative of the differential propagation phase Φ_{dp}, following the fundamental definition

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,

If the differential backscattering phase *δ*_{hv} is relatively constant or negligible, the profile of Ψ_{dp} can be used to estimate *K*_{dp}.

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)

where

where *H*_{2i+1} and *V*_{2i} 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 *V _{i}* and

*H*as

_{i}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 *K*_{dp} as

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.

### 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

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

then the standard deviation of the estimate of *K*_{dp} is related to the phase variation through (Bringi and Chandrasekar 2001)

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

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 *r _{c}*, the normalized cutoff frequency for the passband of the filter must be (Oppenheim and Schafer 1989)

Then, the fluctuations at scales smaller than *r _{c}* 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

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 *K*_{dp} is given by *β*_{1}. The standard deviation of estimated *K*_{dp} in this case can be derived as (Gorgucci et al. 1999)

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 *K*_{dp} using adaptive length *N*, which is adjusted according to three reflectivity levels, such as

This adaptive adjustment is meant to follow the steep phase changes in intensive precipitation regions while keeping the variation of estimated *K*_{dp} 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 *K*_{dp}. 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 *K*_{dp} 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

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

The only difference between (15) and (16) is the factor 2, which needs to be applied back to the final estimates of *K*_{dp}. 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

Therefore, the estimate of *K*_{dp} can be written as

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 *K*_{dp} 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 *K*_{dp}.

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

respectively, which satisfy the following properties under a constant shift by *ϕ*_{0}:

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

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 *K*_{dp}. Figure 2 illustrates the simple logic for data masking with use of the dispersion parameter.

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 *K*_{dp} 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 *K*_{dp} estimates will be filtering angular Φ_{dp} profiles based on angular Ψ_{dp} statistics and estimating *K*_{dp} using cubic spline. Thus, phase unfolding is avoided, and the accuracy is also ensured. An adaptive estimation of *K*_{dp} is developed based on the principles described and is presented in the next section.

## 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):

where *f* (*r _{k}*) 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 *k*th interval between nodes *r _{k}* and

*r*

_{k+1}can be written as

Normally, *K*_{dp} 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 *d _{k}* corresponds to the smoothed angular Φ

_{dp}and the coefficient

*c*corresponds to its range derivatives. Note that only the derivatives at the nodes are of our interest to estimate

_{k}*K*

_{dp}. Therefore, the estimates can be simply obtained as

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 *K*_{dp} region to track the peak *K*_{dp}. The corresponding adaptivity is enabled by introducing varying weights in the minimization (24) as

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., *K*_{dp}); 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 *d _{k}* should be equal to the mean of

*f*(

*r*). Therefore,

_{k}*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

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

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

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

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 K_{dp}

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

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

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 *K*_{dp} 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 *K*_{dp} can be assumed in each range interval. At different *K*_{dp} levels, the smoothness is quantitatively assessed and shown in Table 2 by allowing a 10% variation in *K*_{dp}. The last column can be used as a baseline for choosing *λ*. The existence of extra *K*_{dp} in (33) requires a different smoothing factor for different peak *K*_{dp} values. For example, targeted at *K*_{dp} of 30° km^{−1}, the smoothing factor should be around 1.1Δ*r*.

Figure 3 shows a simulated *K*_{dp} profile, composed of a Gaussian function at 15 km, to demonstrate an infinitely differentiable profile, and a triangular function at 35 km, to demonstrate *K*_{dp} discontinuities. The profile exhibits sharp gradient on *K*_{dp}, 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 *K*_{dp} 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 *K*_{dp} 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 *K*_{dp} 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.

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 *K*_{dp} estimation. The scaling as shown in (33) cannot be applied directly in practice. However, the *K*_{dp} in (33) serves only as a relative weighting factor and does not need to be exact. An approximation of *K*_{dp} can be first translated from the measured reflectivity using the following power-law relation:

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 *K*_{dp} 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 *K*_{dp} 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 *K*_{dp} profile and the resulting scaling functions very well. This two-step iteration technique is adopted in this paper to estimate *K*_{dp} 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 *K*_{dp} of up to 7° km^{−1}. Figures 4a–c present the radar reflectivity, the estimated *K*_{dp} using the conventional filtering and fitting approach, and the estimated *K*_{dp} using the new adaptive approach, respectively. Both the conventional approach and the adaptive approach result in a similar *K*_{dp} field, with the adaptive approach showing less “range streaks.” The adaptive approach presents better resolution, which captures the *K*_{dp} peaks as well as the storm texture that is exhibited in the radar reflectivity. In addition, the *K*_{dp} 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 *K*_{dp} and the intrinsic *K*_{dp} is shown in Fig. 4d. If the intrinsic *K*_{dp} is small, the new estimates are more stable; if the intrinsic *K*_{dp} 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.

### 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 *K*_{dp} estimates in the overall pattern by comparing Figs. 5b,c. The *K*_{dp} peaks are well estimated by the adaptive estimator in the storm core, whereas the *K*_{dp} 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 *K*_{dp} 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.

### 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 *K*_{dp} from the principal Φ_{dp} is shown in Fig. 6b, and the new *K*_{dp} estimates are shown in Fig. 6c. Visually, these two estimators give a comparable *K*_{dp} 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 *K*_{dp} observations than the case of lower range resolutions. However, the *K*_{dp} field in Fig. 6c matches the storm structure much better, and the negative *K*_{dp} in Fig. 6b is largely eliminated. The improved quality of the new adaptive *K*_{dp} is also shown in the *K*_{dp}–*Z _{h}* scatterplot in Fig. 6d.

## 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 *K*_{dp} 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 *K*_{dp}.

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).

## 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

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

### 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,

where

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

respectively.

The following define vectors and matrices as

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

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

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,

Therefore, the coefficients **b** and **d** can be solved as

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,

## Footnotes

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