Analytic characterization of random errors in spectral dual-polarized cloud radar observations

This study presents the first-ever complete characterization of random errors in dual-polarimetric spectral observations of meteorological targets by cloud radars. The characterization is given by means of mathematical equations for joint probability density functions (PDFs) and error covariance matrices. The derived equations are checked for consistency using real radar measurements. One of the main conclusions of the study is that the convenient representation of spectral polarimetric measurements including differential reflectivity Z_DR, correlation coefficient ρ_HV, and differential phase Φ_DP is not suited for the proper characterization of the error covariance matrix. This is because the aforementioned quantities are complex, non-linear functions of the radar raw data, and thus their error covariance matrix is commonly derived using simplified linear relations and by neglecting the correlation of errors. This study formulates the spectral polarimetric measurements in terms of a different set of quantities that allows for a proper analytic treatment of their error covariance matrix. The results given in this study allow for utilization of spectral polarimetric measurements for advanced meteorological applications, among which are variational retrieval techniques, data assimilation, and sensitivity analysis.

1 Introduction

Cloud radars are a major component of state-of-the-art, ground-based observation platforms (Illingworth et al., 2007; Kollias et al., 2020). Their unique capabilities make these instruments extremely valuable for cloud and precipitation research. First, these radars have Doppler capabilities; i.e., they can independently characterize hydrometeors coexisting in the same volume but moving with different speeds relative to the radar (Kollias et al., 2007). Second, the high sensitivity and vast dynamic range make cloud radars capable of measuring return signals from a wide range of particle sizes, which is a challenging task for other instruments like lidars (Bühl et al., 2013). Third, due to relatively low attenuation of microwave signals by liquid water, cloud radars profile clouds up to the top even in the presence of light to moderate rain. These capabilities promote cloud radars for investigation of different formation and development processes throughout the life cycle of clouds. For instance, cloud radars help to characterize initial ice formation and development in mixed-phase clouds (Bühl et al., 2019a, b), improve characterization of pure liquid clouds (Rusli et al., 2017; Acquistapace et al., 2017), estimate rates of aggregation (Kneifel et al., 2015, 2016) and riming (Kalesse et al., 2016; Moisseev et al., 2017; Kneifel and Moisseev, 2020), and quantitatively analyze solid and liquid precipitation (Matrosov, 2005; Matrosov et al., 2006, 2008; Tridon and Battaglia, 2015; Tridon et al., 2017, 2019).

Many cloud radars have dual-polarization capabilities. An interest in polarimetry-based methods in the cloud radar community has been growing, which is indicated by a number of studies during the last decade (Matrosov et al., 2012; Oue et al., 2015; Lu et al., 2015; Myagkov et al., 2016a, b; Matrosov et al., 2017; Oue et al., 2018; Myagkov et al., 2020). Vertically pointed cloud radars often operate in the LDR (linear depolarization ratio) mode; i.e., they transmit a linearly polarized wave (either horizontally or vertically) and receive co- and cross-polarized components of the backscattered signal (e.g., Görsdorf et al., 2015). The LDR mode is efficient for clutter removal and detection of the melting layer and columnar-shaped ice particles. As shown by Matrosov et al. (2001), however, the applicability of the LDR mode at low elevation angles might be limited due to its high sensitivity to the orientation of cloud particles. Therefore, scanning polarimetric cloud radars often have polarimetric modes which are less sensitive to the orientation. One such mode is the hybrid mode (also denoted as the STSR (simultaneous transmission and simultaneous reception) or STAR (simultaneous transmission and reception) mode in the literature). Radars with the hybrid mode emit the horizontal and vertical components of the transmitted wave simultaneously (Myagkov et al., 2015; Bringi and Chandrasekar, 2001, Sect. 4.7). Cloud radars with the hybrid mode allow for adoption of polarimetry-based methods developed during the last several decades for centimeter-wavelength meteorological radars (further denoted as precipitation radars).

Operational precipitation radars are used by weather services to continuously scan the atmosphere, providing polarimetric variables integrated for a scattering volume. In addition to the integrated quantities, cloud radars with the hybrid mode enable spectrally resolved polarimetric observations and, therefore, can provide the same set of polarimetric variables for different types of cloud particles coexisting in the same resolution volume (Oue et al., 2015; Myagkov et al., 2016b, 2020). Spectral observations are in general possible with precipitation radars (Spek et al., 2008; Dufournet and Russchenberg, 2011; Pfitzenmaier et al., 2018). Such measurements, however, are not performed by operational radars due to fast azimuth scanning.

Spectral polarimetry can be used for a development of advanced retrieval methods. For example variational retrievals developed for dual-frequency spectra (Tridon and Battaglia, 2015; Tridon et al., 2017) could be applied also to spectral polarimetry. Moisseev and Chandrasekar (2007) presented first attempts to retrieve profiles of raindrop size distributions using polarimetric spectra from a precipitation radar. This approach, however, has not yet been explored in polarimetric cloud radars.

Recent review studies (Zhang et al., 2019; Morrison et al., 2020; Ryzhkov et al., 2020) demonstrate that polarimetric observations from precipitation radar networks are highly beneficial for the evaluation and development of numerical weather prediction and cloud resolving models. The high value of polarimetric observations is given by their sensitivity to microphysical properties of cloud and precipitation particles such as size, shape, number concentration, state of matter, density, and orientation (Kumjian, 2013). Polarimetric cloud radars are not yet widely used for model improvement. This, however, does not indicate that cloud radar polarimetry is not informative relative to precipitation radars. Conversely, the cloud radar spectral polarimetry can essentially complement available measurements.

The development of both quantitative retrievals and data assimilation algorithms requires the characterization of the systematic and random measurement errors. The former type of errors is solved by a calibration. Calibration aspects of polarimetric quantities have been intensively studied for both precipitation and cloud radars (Chandrasekar et al., 2015) and are out of the scope of this study. In the case of radar observations of meteorological targets, random errors can be characterized from measurements if raw (unaveraged) data are available. Cloud radars, however, rarely store raw data because of the high data rate. Therefore, commonly used approaches to characterize random errors are based on statistical models of the received radar signals. Random errors in radar signals can be represented by a joint probability density function (PDF) of amplitudes and phases in the two orthogonal polarimetric channels. The joint PDF for polarimetric observations obtained for a single pulse can be found in Middleton (1996, chap. 9.2). Single-pulse measurements, however, are rarely used in the radar meteorology because of the low sensitivity and higher requirement for storage space. The observed radar spectra almost always result from the averaging of a number of return pulses. Lee et al. (1994) showed a derivation of a joint probability density function of polarimetric variables for the case of averaging. The authors used a number of assumptions applicable for Earth's surface observations using synthetic-aperture radars. It turns out that the same assumptions are applicable to spectral polarimetric observations of meteorological targets. This allows for using a similar approach in analytic characterization of errors in spectral polarimetric observations.

A number of studies (e.g., Hogan, 2007; Cao et al., 2013; Yoshikawa et al., 2014; Chang et al., 2016; Huang et al., 2020) characterize the joint PDF of polarimetric radar measurements by the error covariance matrix. There are, however, problems with existing approximations of the error covariance matrix for polarimetric observations. First, the elements in the main diagonal of the error covariance matrix – variances of random errors – are found using the first-order Taylor approximation following Bringi and Chandrasekar (2001). Conventional polarimetric variables such as differential reflectivity, correlation coefficient, and differential phase are, however, highly non-linear functions. Therefore, the approximation may lead to biases in the error variance estimates, especially when signal-to-noise ratios (SNRs) and/or the number of averaged samples is low. This problem becomes important for cloud radars collecting polarimetric variables with a high spatial, temporal, and spectral resolution. Second, non-diagonal components of the error covariance matrix are typically set to zero assuming no correlation between errors in measured quantities, but validity and effects of this assumption are not discussed. The information content of measurements is, however, higher when errors are correlated (chap. 3.2.6 in Rodgers, 2000), and therefore non-negligible off-diagonal elements of the covariance matrix should not be ignored.

This study is based on well-known statistical properties of polarimetric radar signals. Using certain simplifications valid for spectral measurements we extend the error model available in the literature and thus derive mathematical expressions characterizing random errors in spectral polarimetric observations of meteorological targets. The study is organized as follows. We review the measurement method of spectral polarimetry with radars operating in the hybrid mode in Sect. 2. In Sect. 3 the likelihood functions of the common polarimetric radar variables are rigorously derived. The error covariance matrix of polarimetric measurements is derived in Sect. 4 by taking into account the correlations among the various measurement random errors. In Sect. 5 the validity of expressions derived for the likelihood functions and error covariance matrix is checked using real raw measurements from a cloud radar.

2 Basics of radar spectral polarimetry in the hybrid mode

This section introduces known relations between a raw cloud radar signal, complex amplitudes, and spectral polarimetric variables for observations of meteorological targets. These relations are based on the same set of assumptions introduced in classical works of Doviak et al. (1979) and Bringi and Chandrasekar (2001) for precipitation radars.

Note that since pulsed radars are currently more common in the meteorological community, we use the term “pulse” to refer to a type of the transmitted radar signal in Sects. 2–4. For radars with frequency-modulated continuous wave (FMCW) signals, however, the term “chirp” should be used. Later, in Sect. 5 we use measurements from a FMCW radar, and therefore the term “chirp” is used there.

2.1 Complex amplitudes of radar measurements

Radar polarimetric measurements are made on an orthogonal measurement basis defined by feeders of the antenna system. In the hybrid mode the measurement basis is typically Cartesian and formed by the horizontal (h) and vertical (v) components. Further this basis is denoted as the h–v basis. Dual-polarimetric cloud radars have two receivers dedicated to the orthogonal polarimetric components of the received signal. For each transmitted pulse the receivers provide range profiles of in-phase I_h,v and quadrature Q_h,v components, where indices h and v denote the polarization state. Note that this study does not cover the radar signal processing to get the I_h,v and Q_h,v profiles. This information can be found in a radar handbook (e.g., Skolnik, 2008, chap. 6). Using N_fft profiles of I_h+iQ_h and I_v+iQ_v, where i is the imaginary unit, the radar calculates complex Doppler spectra in the horizontal and vertical channel, respectively, applying the fast Fourier transformation (FFT) along the time dimension. The complex Doppler spectra are represented by complex amplitudes $\dot{S}$ for each spectral component and each range bin. In the following, ${\dot{S}}_{\mathrm{h}}$ and ${\dot{S}}_{\mathrm{v}}$ denote the measured complex amplitudes of the analyzed spectral component in the horizontal and vertical channels, respectively (the dot hereafter denotes a complex quantity).

2.2 Coherency of complex amplitudes in range and velocity domain

Different range bins as well as different spectral components are often considered to be statistically independent because the corresponding complex amplitudes result from non-coherent scattering of numerous independently moving particles. Some correlation, however, can be expected due to sampling effects and the FFT spectral leakages (e.g., Sect. 5.3 in Marple, 2019). For instance, the power scattered from particles located close to the end of a range bin is distributed between this and the following range bins. These effects depend on filter properties and used FFT windows. It is challenging to give a general analytical solution taking these effects into account. Therefore, these effects are out of the scope of this study. For the sake of simplicity the following analysis is shown only for a single range bin and a single spectral component. Since movements of particles in neighboring range and spectral bins are not related, statistical properties of an individual bin considered in the following are not affected by sampling effects and spectral leakages. The neglection of the dependence of the neighboring bins (due to sampling effects and spectral leakages) leads to an underestimation of the information entropy when a complete spectrum and/or spectral profile is analyzed. This worst case assumption, however, allows for a relatively easy and universal characterization of measurement errors. Future studies may improve the error characterization by considering the sampling and leakage effects.

2.3 Coherency of complex amplitudes in time domain

Unlike precipitation radars which perform rapid azimuth scans, cloud radars are typically pointed to a certain direction or make slow scans to get non-broadened Doppler spectra. Doviak et al. (1979) showed (Eq. 5.2 therein) that the coherency between the adjacent samples depends on the wavelength and the sample repetition period. Cloud radars typically have the pulse repetition frequency on the order of 10 kHz and N_ftt in the range of 128 to 1024. This results in getting a single spectrum every 0.01–0.1 s. For such sampling properties of cloud radars any significant coherency between adjacent samples of a spectral line requires the spectral broadening not exceeding at most a few centimeters per second. The turbulent spectral broadening, however, exceeds a few centimeters per second even in stratiform non-precipitating clouds (Borque et al., 2016). Therefore, consecutive samples of complex amplitudes for a spectral line can be considered to be independent.

2.4 Statistical properties of complex amplitudes

Introduce a measurement column vector

$$\begin{array}{}\text{(1)}& \widehat{\mathit{m}}=[{\widehat{R}}_{\mathrm{h}},{\widehat{J}}_{\mathrm{h}},{\widehat{R}}_{\mathrm{v}},{\widehat{J}}_{\mathrm{v}}{]}^{\mathrm{T}}\end{array}$$

with $\widehat{R}$ and $\widehat{J}$ being real and imaginary parts of a complex amplitude $\dot{S}$, indices h and v denoting the polarization state, and ^T being the transposition sign; the circumflex is used hereafter to emphasize measured quantities. The probability density function (PDF) of $\widehat{\mathit{m}}$, given the true covariance matrix Σ_m of $\widehat{\mathit{m}}$, can be written as follows:

$$\begin{array}{}\text{(2)}& f_{\mathrm{m}}\left(\widehat{\mathit{m}}|{\mathbf{\Sigma }}_{\mathrm{m}}\right)=\left(\mathrm{2}\mathit{\pi }{)}^{-\mathrm{2}}\det \right({\mathbf{\Sigma }}_{\mathrm{m}}{)}^{-\frac{\mathrm{1}}{\mathrm{2}}}{e}^{-\frac{\mathrm{1}}{\mathrm{2}}{\widehat{\mathit{m}}}^{\mathrm{T}}{\mathbf{\Sigma }}_{\mathbf{m}}\widehat{\mathit{m}}}.\end{array}$$

Note that throughout the study a PDF is a function of measured quantities (e.g., $\widehat{\mathit{m}}$ in Eq. 2) with fixed parameters (e.g., Σ_m in Eq. 2). The same PDF is called a likelihood function if the measured quantities are fixed, and the PDF is viewed as a function of parameters.

Doviak et al. (1979) showed that for meteorological targets I and Q components are jointly normal with zero mean, zero correlation, and equal standard deviation. The authors explain that these properties are due to scattering from a large number of particles moving in an unpredictable way in a scattering volume. Since N_fft is much smaller than the number of particles in a resolution volume, the properties are also valid for relations between ${\widehat{R}}_{\mathrm{h}}$ and ${\widehat{J}}_{\mathrm{h}}$ and between ${\widehat{R}}_{\mathrm{v}}$ and ${\widehat{J}}_{\mathrm{v}}$.

The measured complex amplitudes ${\dot{S}}_{\mathrm{h}}$ and ${\dot{S}}_{\mathrm{v}}$, however, can be correlated. Taking these properties into account, the true covariance matrix Σ_m is defined in the following way (Eq. 5.178 in Bringi and Chandrasekar, 2001):

$$\begin{array}{}\text{(3)}& {\mathrm{\Sigma }}_m=\left(\begin{array}{l}\begin{array}{cccc}{\mathit{\sigma }}_{\mathrm{h}}^{\mathrm{2}}& \mathrm{0}& q{\mathit{\sigma }}_{\mathrm{h}}{\mathit{\sigma }}_{\mathrm{v}}& s{\mathit{\sigma }}_{\mathrm{h}}{\mathit{\sigma }}_{\mathrm{v}}\\ \mathrm{0}& {\mathit{\sigma }}_{\mathrm{h}}^{\mathrm{2}}& -s{\mathit{\sigma }}_{\mathrm{h}}{\mathit{\sigma }}_{\mathrm{v}}& q{\mathit{\sigma }}_{\mathrm{h}}{\mathit{\sigma }}_{\mathrm{v}}\\ q{\mathit{\sigma }}_{\mathrm{h}}{\mathit{\sigma }}_{\mathrm{v}}& -s{\mathit{\sigma }}_{\mathrm{h}}{\mathit{\sigma }}_{\mathrm{v}}& {\mathit{\sigma }}_{\mathrm{v}}^{\mathrm{2}}& \mathrm{0}\\ s{\mathit{\sigma }}_{\mathrm{h}}{\mathit{\sigma }}_{\mathrm{v}}& q{\mathit{\sigma }}_{\mathrm{h}}{\mathit{\sigma }}_{\mathrm{v}}& \mathrm{0}& {\mathit{\sigma }}_{\mathrm{v}}^{\mathrm{2}},\end{array}\end{array}\right),\end{array}$$

where σ_h is the standard deviation of ${\widehat{R}}_{\mathrm{h}}$ and ${\widehat{J}}_{\mathrm{h}}$, σ_v is the standard deviation of ${\widehat{R}}_{\mathrm{v}}$ and ${\widehat{J}}_{\mathrm{v}}$, q is the correlation between ${\widehat{R}}_{\mathrm{h}}$ and ${\widehat{R}}_{\mathrm{v}}$, and s is the correlation between ${\widehat{R}}_{\mathrm{h}}$ and ${\widehat{J}}_{\mathrm{v}}$.

2.5 Polarimetric variables

Since for meteorological targets ${\widehat{R}}_{\mathrm{h}}$ is not correlated with ${\widehat{J}}_{\mathrm{h}}$, and ${\widehat{R}}_{\mathrm{v}}$ is not correlated with ${\widehat{J}}_{\mathrm{v}}$, the absolute phases of ${\dot{S}}_{\mathrm{h}}$ and ${\dot{S}}_{\mathrm{v}}$ are uniformly distributed from 0 to 2π and thus uninformative. Therefore, the polarimetric observations in the hybrid mode can be represented by a 2×2 covariance matrix B (Eq. 4.130 in Bringi and Chandrasekar, 2001) instead of the true covariance matrix Σ_m:

$$\begin{array}{}\text{(4)}& \mathbf{B}=\stackrel{\mathrm{‾}}{\mathit{e}{\mathit{e}}^{\mathrm{T}}}=\begin{array}{c}\left(\begin{array}{lr}{B}_{\mathrm{hh}}& {\dot{B}}_{\mathrm{hv}}\\ {\dot{B}}_{\mathrm{hv}}^{\ast }& B_{\mathrm{vv}}\end{array}\right)\end{array},\end{array}$$

where

$$\begin{array}{}\text{(5)}& \mathit{e}=({\dot{S}}_{\mathrm{h}},{\dot{S}}_{\mathrm{v}}{)}^{\mathrm{T}};\end{array}$$

the overline indicates the expected value; B_hh and B_vv represent total powers of the horizontal and vertical components of the received signal, respectively; ${\dot{B}}_{\mathrm{hv}}$ is the covariance between the horizontal and vertical components of the received signal; and ∗ is the complex conjugation sign. Note that in general B_hh, B_vv, and real and imaginary parts of ${\dot{B}}_{\mathrm{hv}}$ can be calibrated in any quantity that is proportional to the power (watts) received by the radar, e.g., classical radar reflectivity (mm⁶ m⁻³) or even arbitrary units (Myagkov et al., 2016a). Recall that in this study the covariance matrix B corresponds to a single spectral component. Such spectral representation of vector signals was introduced by Wiener (1930).

The elements of B are related to the statistics of the complex amplitudes ${\dot{S}}_{\mathrm{h}}$ and ${\dot{S}}_{\mathrm{v}}$ as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle }& {\displaystyle }{B}_{\mathrm{hh}}=\mathrm{var}\left({\widehat{R}}_{\mathrm{h}}\right)+\mathrm{var}\left({\widehat{J}}_{\mathrm{h}}\right)=\mathrm{2}{\mathit{\sigma }}_{\mathrm{h}}^{\mathrm{2}},\text{(7)}& {\displaystyle }& {\displaystyle }{B}_{\mathrm{vv}}=\mathrm{var}\left({\widehat{R}}_{\mathrm{v}}\right)+\mathrm{var}\left({\widehat{J}}_{\mathrm{v}}\right)=\mathrm{2}{\mathit{\sigma }}_{\mathrm{v}}^{\mathrm{2}},\text{(8)}& {\displaystyle }& {\displaystyle }{\dot{B}}_{\mathrm{hv}}=R_{\mathrm{hv}}+i{J}_{\mathrm{hv}}=(q+js){\mathit{\sigma }}_{\mathrm{h}}{\mathit{\sigma }}_{\mathrm{v}},\end{array}$$

where R_hv and J_hv are real and imaginary parts of ${\dot{B}}_{\mathrm{hv}}$.

In the precipitation radar community, dual-polarized measurements are rarely represented by B. Instead a set of polarimetric variables are used. Therefore, the same polarimetric variables (but spectrally resolved) are introduced in this study. Introduce a vector

$$\begin{array}{}\text{(9)}& \mathit{c}=(B_{\mathrm{hh}},Z_{\mathrm{DR}},{\mathit{\rho }}_{\mathrm{HV}},{\mathrm{\Phi }}_{\mathrm{DP}}{)}^{\mathrm{T}},\end{array}$$

where Z_DR is the differential reflectivity, ρ_HV is the correlation coefficient, and Φ_DP is the differential phase. In this study Z_DR, ρ_HV, and Φ_DP are defined for each spectral line using elements of corresponding B:

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{Z}_{\mathrm{DR}}={\displaystyle \frac{B_{\mathrm{hh}}}{B_{\mathrm{vv}}}},\text{(11)}& {\displaystyle }& {\displaystyle }{\mathit{\rho }}_{\mathrm{HV}}=\sqrt{{\displaystyle \frac{R_{\mathrm{hv}}^{\mathrm{2}}+J_{\mathrm{hv}}^{\mathrm{2}}}{B_{\mathrm{hh}}{B}_{\mathrm{vv}}}}},\text{(12)}& {\displaystyle }& {\displaystyle }{\mathrm{\Phi }}_{\mathrm{DP}}=\mathrm{atan}\left(-{\displaystyle \frac{J_{\mathrm{hv}}}{R_{\mathrm{hv}}}}\right).\end{array}$$

Note that elements of the matrix B are in general affected by noise. The noise in both polarimetric channels is not known exactly. Typically, it is estimated from spectra using, for example, the algorithm from Hildebrand and Sekhon (1974). A subtraction of noise levels from corresponding diagonal terms of the covariance matrix B to get an estimate of signal-only powers leads to occasions when the covariance matrix is no longer positively semi-definite. In this case, ρ_HV calculated from the noise-corrected covariance matrix can exceed 1, which is beyond the range of valid values. In order to avoid this problem, we characterize radar measurements without noise subtraction. A further advantage of this approach is that spectral lines containing noise only can also be correctly characterized.

3 Likelihood of elements of the covariance matrix B

Any measurement is affected by inherent uncertainty. As many other measurement devices, radars also attempt to reduce uncertainty in measurements by means of an average over multiple independent samples. The result of the average maximizes the likelihood of the measurements, while the characterization of the distribution of the observations yields an estimate of the uncertainty in the measurements.

Assume the following problem. The state of the atmosphere is represented by the state vector x. A forward model F maps x into a vector

$$\begin{array}{}\text{(13)}& F\left(\mathit{x}\right)=\mathit{b}=(B_{\mathrm{hh}},R_{\mathrm{hv}},J_{\mathrm{hv}},B_{\mathrm{vv}}{)}^{\mathrm{T}}\end{array}$$

in the space of observations. The actual measurement vector is

$$\begin{array}{}\text{(14)}& \widehat{\mathit{b}}=({\widehat{B}}_{\mathrm{hh}},{\widehat{R}}_{\mathrm{hv}},{\widehat{J}}_{\mathrm{hv}},{\widehat{B}}_{\mathrm{vv}}{)}^{\mathrm{T}}=\mathit{b}+\mathit{ϵ},\end{array}$$

where

$$\begin{array}{}\text{(15)}& {\displaystyle }& {\displaystyle }{\widehat{B}}_{\mathrm{hh}}=\langle {\dot{S}}_{\mathrm{h}}{\dot{S}}_{\mathrm{h}}^{*}\rangle ,\text{(16)}& {\displaystyle }& {\displaystyle }{\widehat{R}}_{\mathrm{hv}}=\mathrm{Re}\left(\langle {\dot{S}}_{\mathrm{h}}{\dot{S}}_{\mathrm{v}}^{\ast }\rangle \right),\text{(17)}& {\displaystyle }& {\displaystyle }{\widehat{J}}_{\mathrm{hv}}=\mathrm{Im}\left(\langle {\dot{S}}_{\mathrm{h}}{\dot{S}}_{\mathrm{v}}^{*}\rangle \right),\text{and}\text{(18)}& {\displaystyle }& {\displaystyle }{\widehat{B}}_{\mathrm{vv}}=\langle {\dot{S}}_{\mathrm{v}}{\dot{S}}_{\mathrm{v}}^{*}\rangle \end{array}$$

are constituents of the measured covariance matrix $\widehat{\mathbf{B}}$, and ϵ represents the vector of measurement random errors in each component of $\widehat{\mathit{b}}$. In Eqs. (15)–(18) Re and Im are the real and imaginary parts of a complex number; $<>$ denotes averaging over N_s independent complex spectra calculated from non-overlapping time sequences. The estimators Eqs. (15)–(18) are the same as given in Bringi and Chandrasekar (2001, chap. 6.4.5). The only difference is that within this work the variables are calculated using complex amplitudes for a spectral line instead of using in-phase and quadrature components (I/Q hereafter) as is done by precipitation radars. What is the likelihood of $\widehat{\mathit{b}}$ given the state vector x? In the case that the forward model provides a unique and accurate relation between x and b, the problem is equivalent to finding $f_b\left(\widehat{\mathit{b}}\right|\mathit{b},N_{\mathrm{s}})$ – the likelihood of $\widehat{\mathit{b}}$ – given the true vector of measurements b and the number of averaged spectra N_s.

In a general case, elements of the vector $\widehat{\mathit{b}}$ can be correlated. In this case the derivation of the likelihood function $f_b\left(\widehat{\mathit{b}}\right|\mathit{b},N_{\mathrm{s}})$ is challenging. In order to simplify the derivation, we follow an approach identical to the one demonstrated in Rodgers (2000, Sect. 2.3.1 therein). The author considers a multivariate PDF with correlated errors. He transforms the coordinate system in such a way that orthogonal components of the error vector are independent (uncorrelated). On this basis, the joint PDF can be represented by the product of independent, univariate PDFs for each individual component. Thus, following a similar approach, the derivation of $f_b\left(\widehat{\mathit{b}}\right|\mathit{b},N_{\mathrm{s}})$ provided in this section is developed in four steps.

In Sect. 3.1 we change the basis from h–v to the one on which elements of the vector $\widehat{\mathit{b}}$ become independent. On the new basis, the joint multivariate likelihood function can be represented by the product of the likelihood functions of each independent element. The likelihood of a single independent element is relatively simple to describe analytically. In Sect. 3.2 a formal derivation of the likelihood function on this new basis is provided. The solution for $f_b\left(\widehat{\mathit{b}}\right|\mathit{b},N_{\mathrm{s}})$ is given in Sect. 3.3 converting back to the original space and applying the rule of change in variables. As mentioned above, the radar observations are often represented by the vector c. Therefore, Sect. 3.3 also provides the likelihood $f_c\left(\widehat{\mathit{c}}\right|\mathit{b},N_{\mathrm{s}})$ for the conventional representation of polarimetric measurements.

Note that in this section we keep only equations required to understand the principle of the derivation. The extensive calculus required to prove the formulas used in the section is provided in the Appendix.

3.1 Step 1: change in basis and diagonalization of the covariance matrix B

As previously mentioned, ${\dot{S}}_{\mathrm{h}}$ and ${\dot{S}}_{\mathrm{v}}$ are, in general, correlated. There is, however, always a basis on which the projections of ${\dot{S}}_{\mathrm{h}}$ and ${\dot{S}}_{\mathrm{v}}$ become completely uncorrelated. This basis is further denoted as the c–x (co-polar and cross-polar) basis. The conversion of the vector e on the h–v basis to the vector e_D on c–x basis is made using the unitary operator Q:

$$\begin{array}{}\text{(19)}& {\displaystyle }{\mathit{e}}_D=\begin{array}{c}\left(\begin{array}{l}{\dot{S}}_c\\ {\dot{S}}_x\end{array}\right)\end{array}=\mathbf{Q}\mathit{e}.\end{array}$$

The calculation of the matrix Q is given in Appendix A. Real and imaginary parts of ${\dot{S}}_c$ are jointly distributed normally with the zero mean, zero correlation, and standard deviation σ_c. Real and imaginary parts of ${\dot{S}}_x$ are also jointly distributed normally with zero mean and zero correlation but have, in general, a different standard deviation σ_x.

A transformation from the basis h–v to the basis c–x also changes the covariance matrix of the measurements. The covariance matrix D of measurements on the c–x basis is diagonal and can be found as follows:

$$\begin{array}{}\text{(20)}& \mathbf{D}=\begin{array}{c}\left(\begin{array}{lr}{D}_{cc}& \mathrm{0}\\ \mathrm{0}& D_{xx}\end{array}\right)\end{array}={\mathbf{Q}}^{\mathit{†}}\mathbf{BQ}.\end{array}$$

In Eq. (20) † is the Hermitian conjugate. Zero off-diagonal terms in D indicate that there is no correlation between the orthogonal components, i.e., ${\dot{S}}_c$ and ${\dot{S}}_x$. Expanding Eq. (20), the elements of the matrix D can be found as follows:

$$\begin{array}{}\text{(21)}& {\displaystyle }& {\displaystyle }{D}_{cc}=q_{\mathrm{11}}^{\mathrm{2}}{B}_{\mathrm{hh}}+{\left|{\dot{q}}_{\mathrm{12}}\right|}^{\mathrm{2}}{B}_{\mathrm{vv}}-\mathrm{2}{q}_{\mathrm{11}}\left(R_{\mathrm{12}}{R}_{\mathrm{hv}}+J_{\mathrm{12}}{J}_{\mathrm{hv}}\right)\text{(22)}& {\displaystyle }& {\displaystyle }{D}_{xx}={\left|{\dot{q}}_{\mathrm{12}}\right|}^{\mathrm{2}}{B}_{\mathrm{hh}}+q_{\mathrm{11}}^{\mathrm{2}}{B}_{\mathrm{vv}}+\mathrm{2}{q}_{\mathrm{11}}\left(R_{\mathrm{12}}{R}_{\mathrm{hv}}+J_{\mathrm{12}}{J}_{\mathrm{hv}}\right),\end{array}$$

where ${\dot{q}}_{nm}$ represents elements of Q, with n and m being indices of row and column, respectively;

$$\begin{array}{}\text{(23)}& {\dot{q}}_{\mathrm{12}}=R_{\mathrm{12}}+i{J}_{\mathrm{12}}.\end{array}$$

Similar to relations between the powers and the standard deviations given in Eqs. (6) and (7), σ₁ and σ₂ are related to D_cc and D_xx, respectively:

$$\begin{array}{}\text{(24)}& {\displaystyle }& {\displaystyle }{D}_{cc}=\mathrm{var}\left(R_c\right)+\mathrm{var}\left(J_c\right)=\mathrm{2}{\mathit{\sigma }}_c^{\mathrm{2}}\text{(25)}& {\displaystyle }& {\displaystyle }{D}_{xx}=\mathrm{var}\left(R_x\right)+\mathrm{var}\left(J_x\right)=\mathrm{2}{\mathit{\sigma }}_x^{\mathrm{2}}.\end{array}$$

The measured values ${\widehat{D}}_{cc}$,

$$\begin{array}{}\text{(26)}& {\widehat{D}}_{cx}={\widehat{R}}_{cx}+i{\widehat{J}}_{cx},\end{array}$$

and ${\widehat{D}}_{xx}$ represent elements of the matrix $\widehat{\mathbf{D}}$:

$$\begin{array}{}\text{(27)}& \widehat{\mathbf{D}}={\mathbf{Q}}^{\mathit{†}}\widehat{\mathbf{B}}\mathbf{Q}.\end{array}$$

Note that the operator Q here is the same as in Eq. (20) and not recalculated using $\widehat{\mathbf{B}}$.

3.2 Step 2: likelihood function of the measurements on the c–x basis

In the previous step, measurements were represented on a new – c–x – basis on which the orthogonal components of the measurement vector are independent. On this basis the joint multivariate likelihood function can be represented as a product of likelihood functions with a single element as an argument. This allows for a relatively easy mathematical description of the likelihood function on the c–x basis.

By definition, the off-diagonal elements of the covariance matrix D are zeros (see Eq. 20). This implies no correlation between ${\dot{S}}_c$ and ${\dot{S}}_x$. In this case, the likelihood function $f_d\left(\widehat{\mathit{d}}\right|\mathit{b},N_{\mathrm{s}})$, where

$$\begin{array}{}\text{(28)}& \widehat{\mathit{d}}={\left({\widehat{D}}_{cc},{\widehat{R}}_{cx},{\widehat{J}}_{cx},{\widehat{D}}_{xx}\right)}^{\mathrm{T}},\end{array}$$

can be written as a multiplication of likelihood functions of individual components:

$$\begin{array}{}\text{(29)}& \begin{array}{rl}{f}_d\left(\widehat{\mathit{d}}\right|\mathit{b},N_{\mathrm{s}})=& f\left({\widehat{D}}_{cc}|\mathit{b},N_{\mathrm{s}}\right)f\left({\widehat{R}}_{cx}|\mathit{b},N_{\mathrm{s}}\right)\\ \times & f\left({\widehat{J}}_{cx}|\mathit{b},N_{\mathrm{s}}\right)f\left({\widehat{D}}_{xx}|\mathit{b},N_{\mathrm{s}}\right).\end{array}\end{array}$$

The derivation of the formulas for the calculation of the likelihood functions is tedious and provided in full in Appendix B for the interested reader. The likelihoods of the individual components can be computed as follows:

$$\begin{array}{}\text{(30)}& {\displaystyle }& {\displaystyle }f\left({\widehat{D}}_{cc}|\mathit{b},N_{\mathrm{s}}\right)={\displaystyle \frac{N_{\mathrm{s}}}{{\mathit{\sigma }}_c^{\mathrm{2}}}}{\mathit{\chi }}_{\mathrm{2}{N}_{\mathrm{s}}}^{\mathrm{2}}\left({\displaystyle \frac{N_{\mathrm{s}}}{{\mathit{\sigma }}_c^{\mathrm{2}}}}{\widehat{D}}_{cc}\right),\text{(31)}& {\displaystyle }& {\displaystyle }f\left({\widehat{D}}_{xx}|\mathit{b},N_{\mathrm{s}}\right)={\displaystyle \frac{N_{\mathrm{s}}}{{\mathit{\sigma }}_x^{\mathrm{2}}}}{\mathit{\chi }}_{\mathrm{2}{N}_{\mathrm{s}}}^{\mathrm{2}}\left({\displaystyle \frac{N_{\mathrm{s}}}{{\mathit{\sigma }}_x^{\mathrm{2}}}}{\widehat{D}}_{xx}\right),\end{array}$$

$$\begin{array}{}\text{(32)}& {\displaystyle }\begin{array}{rl}f\left({\widehat{R}}_{cx}|\mathit{b},N_{\mathrm{s}}\right)=& {\displaystyle \frac{\left(\mathrm{2}{N}_{\mathrm{s}}{)}^a\right|{\widehat{R}}_{cx}{|}^{-b}}{\sqrt{\mathit{\pi }}{\mathrm{2}}^{\mathrm{2}{N}_{\mathrm{s}}}\left({\mathit{\sigma }}_c{\mathit{\sigma }}_x{)}^a\mathrm{\Gamma }\right(N_{\mathrm{s}})}}{K}_b\\ \times & \left({\displaystyle \frac{\left|\mathrm{2}{N}_{\mathrm{s}}{\widehat{R}}_{cx}\right|}{{\mathit{\sigma }}_c{\mathit{\sigma }}_x}}\right),\end{array}\text{(33)}& {\displaystyle }\begin{array}{rl}f\left({\widehat{J}}_{cx}|\mathit{b},N_{\mathrm{s}}\right)=& {\displaystyle \frac{\left(\mathrm{2}{N}_{\mathrm{s}}{)}^a\right|{\widehat{J}}_{cx}{|}^{-b}}{\sqrt{\mathit{\pi }}{\mathrm{2}}^{\mathrm{2}{N}_{\mathrm{s}}}\left({\mathit{\sigma }}_c{\mathit{\sigma }}_x{)}^a\mathrm{\Gamma }\right(N_{\mathrm{s}})}}{K}_b\\ \times & \left({\displaystyle \frac{\left|\mathrm{2}{N}_{\mathrm{s}}{\widehat{J}}_{cx}\right|}{{\mathit{\sigma }}_c{\mathit{\sigma }}_x}}\right),\end{array}\end{array}$$

where ${\mathit{\chi }}_k^{\mathrm{2}}$ is the chi-squared distribution with k degrees of freedom,

$$\begin{array}{}\text{(34)}& {\displaystyle }a& {\displaystyle }=(\mathrm{2}{N}_{\mathrm{s}}+\mathrm{1})/\mathrm{2},\text{(35)}& {\displaystyle }b& {\displaystyle }=(\mathrm{1}-\mathrm{2}{N}_{\mathrm{s}})/\mathrm{2},\end{array}$$

Γ is the gamma function, and K_μ is the Bessel function of the second kind of order μ. Recall that σ_c and σ_x in Eqs. (30)–(33) are derived from the elements of b using Eqs. (21)–(22) and Eqs. (24) and (25). Derivation and Monte Carlo evaluation of Eqs. (30)–(33) are given in Appendix B. Appendix B3 shows how to handle Eqs. (32) and (33) when ${\widehat{R}}_{cx}$ and ${\widehat{J}}_{cx}$ are close to 0.

3.3 Step 3: likelihood function on the h–v basis

In the previous step, the likelihood function of measurements represented on the c–x basis was derived. In this subsection we perform a transformation back from the c–x basis to the original h–v basis that allows for the comparison of radar measurements in a common orthogonal reference frame.

Applying the rule of changing variables in a multivariate PDF (e.g.,  Walpole et al., 2012, Theorem 7.4), $f_b\left(\widehat{\mathit{b}}\right|\mathit{b},N_{\mathrm{s}})$ can be found from Eq. (29) as follows:

$$\begin{array}{}\text{(36)}& f_b\left(\widehat{\mathit{b}}\right|\mathit{b},N_{\mathrm{s}})=|{\mathbf{J}}_{bd}\left|f_d\right(\widehat{\mathit{d}}|\mathit{b},N_{\mathrm{s}}).\end{array}$$

As shown in Appendix E, the determinant of the Jacobian J_bd of the transformation from $\widehat{\mathit{b}}$ to $\widehat{\mathit{d}}$ is equal to 1.

3.4 Step 4: likelihood for the conventional representation of polarimetric measurements

As already mentioned, the polarimetric measurements are commonly described by means of a set of quantities (Z_DR, ρ_HV, Φ_DP) that are non-linear functions of the elements of B. It is therefore interesting to derive the likelihood function of those quantities.

Likelihood $f_c\left(\widehat{\mathit{c}}\right|\mathit{b},N_{\mathrm{s}})$ of a vector

$$\begin{array}{}\text{(37)}& \widehat{\mathit{c}}=({\widehat{B}}_{\mathrm{hh}},{\widehat{Z}}_{\mathrm{DR}},{\widehat{\mathit{\rho }}}_{HV},{\widehat{\mathrm{\Phi }}}_{\mathrm{DP}})\end{array}$$

can be found by multiplying $f_b\left(\widehat{\mathit{b}}\right|\mathit{b},N_{\mathrm{s}})$ by $\left|{\mathbf{J}}_{cb}\right|$, with

$$\begin{array}{}\text{(38)}& {\mathbf{J}}_{cb}=-B_{\mathrm{hh}}^{\mathrm{3}}{Z}_{\mathrm{DR}}^{-\mathrm{3}}{\mathit{\rho }}_{\mathrm{HV}}\end{array}$$

being the Jacobian of the transformation from $\widehat{\mathit{c}}$ to $\widehat{\mathit{b}}$ (see Appendix F):

$$\begin{array}{}\text{(39)}& f_c\left(\widehat{\mathit{c}}\right|\mathit{b},N_{\mathrm{s}})=B_{\mathrm{hh}}^{\mathrm{3}}{Z}_{\mathrm{DR}}^{-\mathrm{3}}{\mathit{\rho }}_{\mathrm{HV}}{f}_b(\widehat{\mathit{b}}|\mathit{b},N_{\mathrm{s}}).\end{array}$$

The final results of this section – Eqs. (36) and (39) – can be used for the maximum likelihood optimization and Bayesian inference methods. Ready-to-use MATLAB implementations of these equations are provided in the Supplement.

4 Error covariance matrices

In the previous section we derived mathematical expressions for the likelihood for polarimetric radar observations. A number of scientific studies, however, require the numerical computation of the covariance matrix of the measurement errors. For instance, optimal estimation, data assimilation, and sensitivity analysis are often performed using error covariance matrices. Unfortunately, an analytical integration of Eqs. (29), (36), and (39) required for the statistical moment calculation is challenging. In this section we therefore follow a different and more viable way to calculate elements of the error covariance matrix. This is done by going back to the representation of the measurements on the convenient c–x basis and applying well-known rules for the calculation of the covariance matrix after a linear transformation.

4.1 Error covariance matrix of b

In this section we start from the representation of measurements on the c–x basis because the measurement errors are independent in this case. Recall that Eq. (27) relates the covariance matrix $\widehat{\mathbf{D}}$ and $\widehat{\mathbf{B}}$. This equation thus can be used to find relations between elements of the vector $\widehat{\mathit{b}}$ on the original h–v basis and elements of the vector $\widehat{\mathit{d}}$ on the c–x basis.

The covariance matrix $\widehat{\mathbf{B}}$ estimated from measurements is related to the matrix $\widehat{\mathbf{D}}$ as follows:

$$\begin{array}{}\text{(40)}& \widehat{\mathbf{B}}=\mathbf{Q}\widehat{\mathbf{D}}{\mathbf{Q}}^{\mathit{†}}.\end{array}$$

After expanding Eq. (40) it can be seen that the elements of the vector $\widehat{\mathit{b}}$ can be found as linear combinations of the elements of the vector $\widehat{\mathit{d}}$:

$$\begin{array}{}\text{(41)}& {\displaystyle }\begin{array}{rl}{\widehat{B}}_{\mathrm{hh}}& =q_{\mathrm{11}}^{\mathrm{2}}{\widehat{D}}_{cc}+{\left|{\dot{q}}_{\mathrm{12}}\right|}^{\mathrm{2}}{\widehat{D}}_{xx}+\mathrm{2}{q}_{\mathrm{11}}\left(R_{\mathrm{12}}{\widehat{R}}_{cx}\right.\\ & \left.+J_{\mathrm{12}}{\widehat{J}}_{cx}\right),\end{array}\text{(42)}& {\displaystyle }\begin{array}{rl}{\widehat{R}}_{\mathrm{hv}}& =q_{\mathrm{11}}{R}_{\mathrm{12}}\left({\widehat{D}}_{xx}-{\widehat{D}}_{cc}\right)+\left(q_{\mathrm{11}}^{\mathrm{2}}-R_{\mathrm{12}}^{\mathrm{2}}+J_{\mathrm{12}}^{\mathrm{2}}\right){\widehat{R}}_{cx}\\ & -\mathrm{2}{R}_{\mathrm{12}}{J}_{\mathrm{12}}{\widehat{J}}_{cx}\end{array}\end{array}$$

$$\begin{array}{}\text{(43)}& {\displaystyle }\begin{array}{rl}{\widehat{J}}_{\mathrm{hv}}& =q_{\mathrm{11}}{J}_{\mathrm{12}}\left({\widehat{D}}_{xx}-{\widehat{D}}_{cc}\right)+\left(q_{\mathrm{11}}^{\mathrm{2}}+R_{\mathrm{12}}^{\mathrm{2}}-J_{\mathrm{12}}^{\mathrm{2}}\right){\widehat{J}}_{cx}\\ & -\mathrm{2}{R}_{\mathrm{12}}{J}_{\mathrm{12}}{\widehat{R}}_{cx}\end{array}\text{(44)}& {\displaystyle }\begin{array}{rl}{\widehat{D}}_{\mathrm{vv}}& ={\left|{\dot{q}}_{\mathrm{12}}\right|}^{\mathrm{2}}{\widehat{D}}_{cc}+q_{\mathrm{11}}^{\mathrm{2}}{\widehat{D}}_{xx}-\mathrm{2}{q}_{\mathrm{11}}\left(R_{\mathrm{12}}{\widehat{R}}_{cx}\right.\\ & \left.+J_{\mathrm{12}}{\widehat{J}}_{cx}\right).\end{array}\end{array}$$

Or they can be found in matrix form:

$$\begin{array}{}\text{(45)}& \begin{array}{rl}\widehat{\mathit{b}}& =\left(\begin{array}{cccc}{q}_{\mathrm{11}}^{\mathrm{2}}& \mathrm{2}{q}_{\mathrm{11}}{R}_{\mathrm{12}}& \mathrm{2}{q}_{\mathrm{11}}{J}_{\mathrm{12}}& {\left|{\dot{q}}_{\mathrm{12}}\right|}^{\mathrm{2}}\\ -q_{\mathrm{11}}{R}_{\mathrm{12}}& q_{\mathrm{11}}^{\mathrm{2}}-R_{\mathrm{12}}^{\mathrm{2}}+J_{\mathrm{12}}^{\mathrm{2}}& -\mathrm{2}{R}_{\mathrm{12}}{J}_{\mathrm{12}}& q_{\mathrm{11}}{R}_{\mathrm{12}}\\ -q_{\mathrm{11}}{J}_{\mathrm{12}}& -\mathrm{2}{R}_{\mathrm{12}}{J}_{\mathrm{12}}& q_{\mathrm{11}}^{\mathrm{2}}+R_{\mathrm{12}}^{\mathrm{2}}-J_{\mathrm{12}}^{\mathrm{2}}& q_{\mathrm{11}}{J}_{\mathrm{12}}\\ {\left|{\dot{q}}_{\mathrm{12}}\right|}^{\mathrm{2}}& -\mathrm{2}{q}_{\mathrm{11}}{R}_{\mathrm{12}}& -\mathrm{2}{q}_{\mathrm{11}}{J}_{\mathrm{12}}& q_{\mathrm{11}}^{\mathrm{2}}\end{array}\right)\widehat{\mathit{d}}\\ & =\mathbf{M}\widehat{\mathit{d}}.\end{array}\end{array}$$

In this case, as shown in Wilks (chap. 10.4.3), the error covariance matrix Σ_b of $\widehat{\mathit{b}}$ can be calculated from the error covariance matrix Σ_d of $\widehat{\mathit{d}}$:

$$\begin{array}{}\text{(46)}& {\mathrm{\Sigma }}_b=\mathbf{M}{\mathrm{\Sigma }}_d{\mathbf{M}}^{\mathrm{T}},\end{array}$$

where

$$\begin{array}{}\text{(47)}& {\mathrm{\Sigma }}_d=\left(\begin{array}{cccc}\mathrm{4}{\mathit{\sigma }}_c^{\mathrm{4}}/N_{\mathrm{s}}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\mathit{\sigma }}_c^{\mathrm{2}}{\mathit{\sigma }}_x^{\mathrm{2}}/N_{\mathrm{s}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& {\mathit{\sigma }}_c^{\mathrm{2}}{\mathit{\sigma }}_x^{\mathrm{2}}/N_{\mathrm{s}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{4}{\mathit{\sigma }}_x^{\mathrm{4}}/N_{\mathrm{s}}\end{array}\right).\end{array}$$

Recall that the off-diagonal terms of Σ_d are set to 0 taking into account that the elements of $\widehat{\mathit{d}}$ are not correlated. The derivation of diagonal terms – variances of elements of $\widehat{\mathit{d}}$ – is given in Appendix C.

The main result of this subsection – Eq. (46) – was implemented as a ready-to-use MATLAB function that is available in the Supplement.

In the next subsection we also consider the error covariance matrix of the vector $\widehat{\mathit{c}}$ – the conventional representation of polarimetric measurements. Note however that, as is shown in Sect. 5, the approximation of the error covariance matrix of $\widehat{\mathit{c}}$ has issues which may limit its applicability.

Table 1 Operational setting of the used W-band radar.

4.2 Error covariance matrix of the conventional measurement vector c

As is shown in Sect. 4.1, the error covariance matrix Σ_b can be used to characterize uncertainties in spectral radar observations. By analogy to what is done in Sect. 3.4, one might think about applying again the rules of linear transformation to obtain the error covariance of the vector $\widehat{\mathit{c}}$. In this section we do that by means of a linearization of the formulas that define the components of c (Eqs. 10, 11, and 12). It is further demonstrated in Sect. 5 that such representation of measurement uncertainties for $\widehat{\mathit{c}}$ is deficient.

Recall that the calculation of $\widehat{\mathit{c}}$ includes highly non-linear functions (see Sect. 2.5). Therefore, the error covariance matrix Σ_c of the vector $\widehat{\mathit{c}}$ is estimated using the first-order Taylor approximation:

$$\begin{array}{}\text{(48)}& {\mathrm{\Sigma }}_c=\mathbf{S}{\mathrm{\Sigma }}_b{\mathbf{S}}^{\mathrm{T}},\end{array}$$

where S is the sensitivity matrix:

$$\begin{array}{}\text{(49)}& \mathbf{S}=\left(\begin{array}{cccc}\frac{\partial B_{\mathrm{hh}}}{\partial B_{\mathrm{hh}}}& \frac{\partial B_{\mathrm{hh}}}{\partial R_{\mathrm{hv}}}& \frac{\partial B_{\mathrm{hh}}}{\partial J_{\mathrm{hv}}}& \frac{\partial B_{\mathrm{hh}}}{\partial B_{\mathrm{vv}}}\\ \frac{\partial Z_{\mathrm{DR}}}{\partial B_{\mathrm{hh}}}& \frac{\partial Z_{\mathrm{DR}}}{\partial R_{\mathrm{hv}}}& \frac{\partial Z_{\mathrm{DR}}}{\partial J_{\mathrm{hv}}}& \frac{\partial Z_{\mathrm{DR}}}{\partial B_{\mathrm{vv}}}\\ \frac{\partial {\mathit{\rho }}_{\mathrm{HV}}}{\partial B_{\mathrm{hh}}}& \frac{\partial {\mathit{\rho }}_{\mathrm{HV}}}{\partial R_{\mathrm{hv}}}& \frac{\partial {\mathit{\rho }}_{\mathrm{HV}}}{\partial J_{\mathrm{hv}}}& \frac{\partial {\mathit{\rho }}_{\mathrm{HV}}}{\partial B_{\mathrm{vv}}}\\ \frac{\partial {\mathrm{\Phi }}_{\mathrm{DP}}}{\partial B_{\mathrm{hh}}}& \frac{\partial {\mathrm{\Phi }}_{\mathrm{DP}}}{\partial R_{\mathrm{hv}}}& \frac{\partial {\mathrm{\Phi }}_{\mathrm{DP}}}{\partial J_{\mathrm{hv}}}& \frac{\partial {\mathrm{\Phi }}_{\mathrm{DP}}}{\partial B_{\mathrm{vv}}}\end{array}\right).\end{array}$$

Note that the utilization of the first-order Taylor approximation for variances of polarimetric variables was proposed in the classical book of Bringi and Chandrasekar (2001). Equation (D1) in Appendix D shows the complete matrix S in terms of B_hh, B_vv, ${\dot{B}}_{\mathrm{hv}}$, R_hv, and J_hv.

A ready-to-use MATLAB implementation of Eq. (48) is provided in the Supplement. As is shown in the next section, the error covariance matrix Σ_c does not always reflects the true statistical properties of polarimetric observations. Therefore, this approximation is provided only for demonstration purposes, and it is not recommended.

5 Consistency checks on radar observations

In order to check consistency of Eqs. (36), (39), (46), and (48) with radar measurements, I/Q data collected with a W-band cloud radar with the hybrid polarimetric mode were used (Myagkov and Unal, 2021). The radar is a part of a dual-frequency system owned and operated by the Technical University of Delft in Cabauw, the Netherlands. Technical specifications of the radar can be found in Myagkov et al. (2020). The radar uses frequency-modulated continuous signals. Küchler et al. (2017) explain the operation principle and show that the radar profiles the atmosphere using several chirp types. Each chirp type is dedicated to a certain distance range. During measurements chirp types are switched consequently. For each chirp type a number of chirps (chirp sequence hereafter) are processed continuously. Operational settings used during I/Q measurements are listed in Table 1.

Measurements were made during a rain event on 21 June 2021 at 7:44 UTC. I/Q measurements provide a high data rate of about 900 MB min⁻¹. Therefore, about 3 min of I/Q measurements were collected for the analysis. The radar was pointed to 45^∘ elevation. Since different chirp types have different properties, in the following only I/Q data collected with the first chirp type are used. Since the first chirp sequence covers the lowest part of the atmosphere, the analyzed data correspond to rain. As explained in Sect. 2, no noise subtraction is required to describe the statistics of the measurements. We therefore use all available spectral lines, including those containing noise only. A total of 90 % of spectral noise power was from 0.2–1.3$\times {\mathrm{10}}^{-\mathrm{3}}$ a.u (arbitrary units). Signal-to-noise ratio (defined here as a ratio of signal power in a spectral line divided by the mean spectral noise power in the same range bin) specified in linear units was from 0 (no signal) to 10⁶. We would like to emphasize that no filtering based on signal-to-noise ratio was applied. Taking into account that the first chirp type has 37 range bins, in total 2.2×10³ chirp sequences (15.9×10⁶ chirps) are available in each polarimetric channel.

5.1 Processing

All I/Q measurements within a chirp sequence in every polarimetric channel are split into 224 continuous blocks. Each block contains 32 I/Q pairs. The FFT with the Blackman weighting window is applied to each block to get complex Doppler spectra. Then the 224 blocks are split into 28 sub-blocks with 8 spectra in each sub-block. Within each sub-block elements of the vector $\widehat{\mathit{b}}$ are calculated according to Eqs. (15)–(18) with N_s=8 for every spectral line. For each $\widehat{\mathit{b}}$ the vector $\widehat{\mathit{c}}$ is obtained. Note that for this, Eqs. (10)–(12) were applied to elements of $\widehat{\mathit{b}}$ instead of b. Using vectors $\widehat{\mathit{b}}$ and $\widehat{\mathit{c}}$ within a sequence, the error covariance matrices ${\widehat{\mathrm{\Sigma }}}_b$ and ${\widehat{\mathrm{\Sigma }}}_c$ are calculated numerically. The circumflex here indicates that the error covariance matrices are estimated from measurements.

Figure 1 Schematic illustration of the error covariance matrix calculation.

The calculation of the likelihood functions using Eqs. (36) and (39) requires b. The approximation of covariance matrices using Eqs. (46) and (48) requires the matrix B. In order to estimate b and B, elements of the vector $\widehat{\mathit{b}}$ are averaged over 28 sub-blocks available within a single chirp sequence. These averaged values are assumed to be elements of the vector b from which the matrix B is obtained. Using B and N_s=8, Σ_b and Σ_c are calculated for each chirp sequence as shown in Fig. 1.

5.2 Filtering

Figure 2 Distributions of the ratio of mean power over the power standard deviation for the horizontal (blue line) and vertical (yellow line) channels. The expected distribution is shown with the red line. The vertical black line indicates the threshold corresponding to the 5th percentile of the distribution for the randomly generated complex numbers.

The random error analysis provided in this study is only applicable to volume-distributed scattering and noise. As discussed in Sect. 2, in this case ${\widehat{R}}_{\mathrm{h}}$ is not correlated with ${\widehat{J}}_{\mathrm{h}}$, and ${\widehat{R}}_{\mathrm{v}}$ is not correlated with ${\widehat{J}}_{\mathrm{v}}$. However, radar observations in general contain scattering from atmospheric plankton, ground clutter, and coherent receiver noise, which do not fulfill the assumption. In order to filter out spectral lines with correlated real and imaginary parts, a simple filtering rule was applied. It is known that for a signal with uncorrelated in-phase and quadrature components, its mean power and power standard deviation are related to each other (Eq. 5.193 in Bringi and Chandrasekar, 2001). Figure 2 shows distributions of the mean power over the power standard deviation calculated in the horizontal and vertical polarization channels shown by blue and yellow lines, respectively. It can be seen that the mode of the distributions is close to the theoretical value of $\sqrt{N_{\mathrm{s}}}=\mathrm{2.8}$. The distributions, however, have a considerable tail on the left side. These small values of the ratio are expected for correlated in-phase and quadrature components. Thus, a threshold in the ratio of the mean power over the standard deviation of power can be used to filter out unwanted spectral lines. In order to specify the threshold, the Monte Carlo approach was used. A total of 15.9×10⁶ random complex values with normal distribution, zero mean, and a standard deviation of 1 were generated. The same processing as for measured I/Q data was applied to the generated complex values. The distribution of the ratio of the mean power over the power standard deviation for the generated data (denoted as expected distribution) is shown in Fig. 2 by the red line. The expected distribution has a much smaller tail on the left side relative to the ones of the measured distributions. The threshold of 2.3 used for filtering is chosen as the 5th percentile of the expected distribution. Vectors $\widehat{\mathit{b}}$ and $\widehat{\mathit{c}}$ are excluded from the analysis if for the corresponding spectral component within a chirp sequence the ratio of the mean power over the power standard deviation is below the threshold in at least one of the polarimetric channels. Around 18 % of the data are excluded.

5.3 Evaluation of $f_b\left(\widehat{\mathit{b}}\right|\mathit{b},N_{\mathrm{s}})$ and $f_b\left(\widehat{\mathit{c}}\right|\mathit{b},N_{\mathrm{s}})$

Recall that b is estimated from measurements by averaging all available sub-blocks within a chirp sequence; b, however, can also be estimated by maximization of the likelihood functions given in Eqs. (36) and (39). In this case, an optimization algorithm needs to be employed to find a set of elements of b corresponding to the global maximum in either Eq. (36) or Eq. (39). This study uses a derivative-free optimization method available by default in MATLAB (Lagarias et al., 1998). Since the optimization method minimizes a function, the likelihood functions were not used directly. Instead, the following cost functions were used for the minimization:

$$\begin{array}{}\text{(50)}& {\displaystyle }& {\displaystyle }{C}_b=-\sum _{l=\mathrm{1}}^{\mathrm{28}}{\log }_{\mathrm{10}}\left(f_b\left(\widehat{\mathit{b}}|\mathit{b},N_{\mathrm{s}}\right)\right),\text{(51)}& {\displaystyle }& {\displaystyle }{C}_c=-\sum _{l=\mathrm{1}}^{\mathrm{28}}{\log }_{\mathrm{10}}\left(f_c\left(\widehat{\mathit{c}}|\mathit{b},N_{\mathrm{s}}\right)\right).\end{array}$$

Here the index l runs over 28 sub-blocks within a chirp sequence. Equations (50) and (51) take into account that the consecutive $\widehat{\mathit{b}}$ vectors are not correlated. In this case the total likelihood of 28 vector $\widehat{\mathit{b}}$'s is a product of likelihood of each individual $\widehat{\mathit{b}}$. In order to avoid an overflow of double numbers, the logarithm was used. In this case the logarithm of the product is replaced by the sum of logarithms. The logarithm is a monotonically increasing function, and therefore it does not change the position of the maximum of the likelihood function. Finally, the minus sign was introduced to have a smaller value of a cost function corresponding to a higher value of the likelihood. For the evaluation, 1000 chirp sequences were chosen randomly for the maximum likelihood estimation using $f_b\left(\widehat{\mathit{b}}\right|\mathit{b},N_{\mathrm{s}})$. In each chirp sequence a single spectral line was randomly chosen for the analysis. Thus, there are 28 vector $\widehat{\mathit{b}}$'s available in each of the 1000 chirp sequences. For each sequence, the optimization algorithm requires an initial guess of b. In order to avoid local minima, five different initial guesses were used, which are a coefficient P multiplied by the first $\widehat{\mathit{b}}$ in the analyzed chirp sequence. The values of P were 0.5, 0.75, 1, 1.25, and 1.5. The solution giving the lowest cost function out of the five outcomes was chosen as the result. Similarly the maximum likelihood estimation using $f_c\left(\widehat{\mathit{c}}\right|\mathit{b},N_{\mathrm{s}})$ was done using independently chosen 1000 chirp sequences. Figure 3 shows a comparison of elements of b estimated by averaging over 28 sub-blocks and those estimated by the maximum likelihood approach. All panels show a good agreement indicated by the close-to-unity slope of the linear regression. Both $f_b\left(\widehat{\mathit{b}}\right|\mathit{b},N_{\mathrm{s}})$ (results in the first row of Fig. 3) and $f_c\left(\widehat{\mathit{c}}\right|\mathit{b},N_{\mathrm{s}})$ (results in the second row of Fig. 3) show the same level of agreement and, therefore, can be used with no difference.

Figure 3 Comparison of elements of b estimated by averaging over 28 sub-blocks (x axis) with those estimated by the maximum likelihood approach (y axis); $f_b\left(\widehat{\mathit{b}}\right|\mathit{b},N_{\mathrm{s}})$ was used for (a)–(d). $f_c\left(\widehat{\mathit{c}}\right|\mathit{b},N_{\mathrm{s}})$ was used for (e)–(h). Each panel contains 1000 points described in the text. Linear regressions are shown by solid red lines. Each panel has a text box with the slope of the corresponding linear regression. Uncertainties in the slopes were estimated using bootstrapping. Note that units are not critical for the evaluation of the correctness of the derived likelihood functions. Therefore, arbitrary units (a.u.) are used.

5.4 Evaluation of Σ_c

Figure 4 Comparison of variances of (a) ${\widehat{B}}_{\mathrm{hh}}$, (b) ${\widehat{Z}}_{\mathrm{DR}}$, (c) ${\widehat{\mathit{\rho }}}_{HV}$, and (d) $\widehat{{\mathrm{\Phi }}_{\mathrm{DP}}}$. Approximations developed in this study are on the x axis. Approximations from Bringi and Chandrasekar (2001) are on the y axis; ρ_HV is color-coded in (c) and (d) to illustrate the values of ρ_HV at which approximations lead to erroneous values (see details in text). Note that units are not critical for the evaluation of the derived equations. Therefore, arbitrary units (a.u.) are used in (a).

Diagonal elements of Σ_c – variances of ${\widehat{B}}_{\mathrm{hh}}$, ${\widehat{Z}}_{\mathrm{DR}}$, ${\widehat{\mathit{\rho }}}_{HV}$, and ${\widehat{\mathrm{\Phi }}}_{\mathrm{DP}}$ – were checked against those calculated using Eqs. (6.139a), (6.141), (6.144), and (6.143) in Bringi and Chandrasekar (2001), respectively. Taking into account that samples for a spectral line are not correlated, approximations for variances of ${\widehat{B}}_{\mathrm{hh}}$, ${\widehat{Z}}_{\mathrm{DR}}$, ${\widehat{\mathit{\rho }}}_{HV}$, and ${\widehat{\mathrm{\Phi }}}_{\mathrm{DP}}$ based on the equations in Bringi and Chandrasekar (2001) are

$$\begin{array}{}\text{(52)}& {\displaystyle }& {\displaystyle }{\mathrm{VAR}}_{b\mathrm{hh}}={\displaystyle \frac{B_{\mathrm{hh}}^{\mathrm{2}}}{N_{\mathrm{s}}}},\text{(53)}& {\displaystyle }& {\displaystyle }{\mathrm{VAR}}_{z\mathrm{dr}}={\displaystyle \frac{\mathrm{2}{Z}_{\mathrm{DR}}^{\mathrm{2}}(\mathrm{1}-{\mathit{\rho }}_{\mathrm{HV}}^{\mathrm{2}})}{N_{\mathrm{s}}}},\text{(54)}& {\displaystyle }& {\displaystyle }{\mathrm{VAR}}_{\mathit{\rho }}={\displaystyle \frac{(\mathrm{1}-{\mathit{\rho }}_{\mathrm{HV}}^{\mathrm{2}}{)}^{\mathrm{2}}}{\mathrm{2}{N}_{\mathrm{s}}{\mathit{\rho }}_{\mathrm{HV}}^{\mathrm{2}}}},\text{and}\text{(55)}& {\displaystyle }& {\displaystyle }{\mathrm{VAR}}_{\mathrm{\Phi }}={\displaystyle \frac{(\mathrm{1}-{\mathit{\rho }}_{\mathrm{HV}}^{\mathrm{2}})}{\mathrm{2}{N}_{\mathrm{s}}{\mathit{\rho }}_{\mathrm{HV}}^{\mathrm{2}}}},\end{array}$$

respectively.

Figure 4 shows that VAR_bhh, VAR_zdr, and VAR_Φ match exactly Σ_c(1,1), Σ_c(2,2), and Σ_c(4,4), respectively. VAR_ρ, however, agrees with Σ_c(3,3) only at values of ρ_HV>0.95. Below this value VAR_ρ overestimates the variance of ${\widehat{\mathit{\rho }}}_{HV}$. At values of ρ_HV close to 0, VAR_ρ has unrealistically high values, which result from ρ_HV in the denominator of Eq. (54).

Figure 4d also shows unrealistic values with both approximations of the ${\widehat{\mathrm{\Phi }}}_{\mathrm{DP}}$ variance. Taking into account that ${\widehat{\mathrm{\Phi }}}_{\mathrm{DP}}$ can take values within the range of 0 to 2π rad, the variance of ${\widehat{\mathrm{\Phi }}}_{\mathrm{DP}}$ exceeding 10³ rad² is definitely erroneous. The high variance of ${\widehat{\mathrm{\Phi }}}_{\mathrm{DP}}$ corresponds to values of ρ_HV<0.3. This effect results from the first-order Taylor approximation of Eq. (12), which is a highly non-linear function.

Figure 5 Comparison of ${\widehat{\mathrm{\Sigma }}}_c$ estimated from the radar measurements with Σ_c obtained from Eq. (48). Elements of Σ_c are given on the x axes. Elements of ${\widehat{\mathrm{\Sigma }}}_c$ are given on the y axes. The first and the second numbers in brackets indicate the row and the column of the corresponding matrix, respectively. Linear regressions are shown by red lines. Slopes of the linear regressions and Pearson correlations are given in boxes in each panel. Uncertainties in the slope and the correlation are represented by ±1 standard deviation of the corresponding parameter. The standard deviations are obtained using bootstrapping. Panels without linear regressions show elements for which Eq. (48) gives only near-zero values. Note that units are not critical for the evaluation of the derived equations. Therefore, arbitrary units (a.u.) are used. Also note that only values on the x and y axes in an individual panel should be compared. Values in different panels should not be compared.

A comparison of the error covariance matrices ${\widehat{\mathrm{\Sigma }}}_c$ with the calculated one Σ_c is shown in Fig. 5. Figure 5f, k, and p indicate considerable differences caused by the first-order Taylor approximation in variances of ${\widehat{Z}}_{\mathrm{DR}}$, ${\widehat{\mathit{\rho }}}_{HV}$, and ${\widehat{\mathrm{\Phi }}}_{\mathrm{DP}}$, respectively. The results also reveal that the first-order Taylor approximation cannot adequately represent most of the non-diagonal components of the error covariance matrix.

5.5 Evaluation of Σ_b

Figure 6 Comparison of ${\widehat{\mathrm{\Sigma }}}_b$ estimated from the radar measurements with Σ_b obtained from Eq. (46). Elements of Σ_b are given on the x axes. Elements of ${\widehat{\mathrm{\Sigma }}}_b$ are given on the y axes. The first and the second numbers in brackets indicate the row and the column of the corresponding matrix, respectively. Linear regressions are shown by red lines. Slopes of the linear regressions and Pearson correlations are given in boxes in each panel. Uncertainties in the slope and the correlation are represented by ±1 standard deviation of the corresponding parameter. The standard deviations are obtained using bootstrapping. Note that units are not critical for the evaluation of the derived equations. Therefore, arbitrary units (a.u.) are used.

Figure 6 shows a comparison of elements of error covariance matrices ${\widehat{\mathrm{\Sigma }}}_b$ estimated from the radar measurements with those calculated using Eq. (46). Estimated and calculated elements are in a good agreement. Linear regressions shown in the panels by red lines have slopes close to 1. Pearson correlations between estimated and calculated elements exceed 0.96. These results indicate an agreement of the theoretical calculation with measurements and thus confirm the correctness of Eq. (46). As expected, Figs. 6a and 5a show equivalent results. This is because the co-polar signal B_hh is effectively the same in both measurement representations and highlights the relevance of the present study only for dual-polarimetric quantities.

It is thus concluded that any application of spectral polarimetric measurements which requires the estimate of the error covariance matrix (e.g., variational retrievals, data assimilation, and sensitivity analysis) should be performed in the space of observations $\widehat{\mathit{b}}$ rather than $\widehat{\mathit{c}}$.

6 Summary and outlook

Spectral and polarimetric cloud radar observations have a great potential in the cloud science (Kollias et al., 2020). Decades of such measurements have been already collected by, for example, the ARM (Atmospheric Radiation Measurement) and CLOUDNET communities. An advanced application of these vast datasets requires an accurate characterization of measurement uncertainties. Systematic errors in moment radar data and polarimetric variables have been discussed in many studies. Random measurement errors, in contrast, are rarely considered in the literature. There are three main problems in existing random-error-characterization methods in meteorological studies, namely (1) a lack of joint PDFs for averaged spectral polarimetric measurements, (2) neglection of non-diagonal components of the error covariance matrix, and (3) inaccuracy of the first-order approximation in variances of polarimetric variables. This study thus aims to provide solutions for these three problems.

Equations provided in Sect. 3 give an exact mathematical solution for the joint PDFs of spectral polarimetric observations. The PDFs are given for two equivalent representations of the measurements: (1) $\mathit{b}={\left(B_{\mathrm{hh}},R_{\mathrm{hv}},J_{\mathrm{hv}},B_{\mathrm{vv}}\right)}^{\mathrm{T}}$ and (2) $\mathit{c}={\left(B_{\mathrm{hh}},Z_{\mathrm{DR}},{\mathit{\rho }}_{\mathrm{HV}},{\mathrm{\Phi }}_{\mathrm{DP}}\right)}^{\mathrm{T}}$. The obtained equations take into account non-coherent averaging of spectra, which is applied by a majority of cloud radars to improve the sensitivity. Maximum likelihood estimators of b based on Eqs. (36) and (39) were compared with the estimator based on longer averaging. The comparison was based on dual-polarimetric cloud radar observations. The comparison showed a good agreement. Both PDFs can be equivalently used for methods based on the maximum likelihood and Bayesian inference.

Section 4 is focused on the error covariance matrix required for a number of applications such as data assimilation, sensitivity analysis, and variational retrievals. The error covariance matrices Σ_b and Σ_c for b and c, respectively, are obtained using the characteristic functions of the PDFs described in Sect. 3. Since the calculation of the c includes highly non-linear functions, Σ_c was derived using the first-order Taylor approximation. The same approach was used by Bringi and Chandrasekar (2001) to get equations for variances of polarimetric observations.

The error covariance matrices were evaluated using I/Q observations from a polarimetric W-band radar. It is illustrated that elements of Σ_c have considerable differences from those estimated from the measurements. First, we found differences in variances of Z_DR, ρ_HV, and Φ_DP of up to a factor of 10, 5, and 100, respectively. Second, the calculated variance of Φ_DP shows unrealistically high values by far exceeding the range of possible values. Third, most of the off-diagonal terms of Σ_c are not correlated with corresponding values estimated from observations. We relate the differences to the first-order Taylor approximation. The Taylor approximation assumes linear relations between elements of the vector b and the elements of the vector c, while the relations include highly non-linear functions. In contrast, Σ_b agrees well with the observations. The correlation between calculated elements of Σ_b with those estimated from the observations exceeds 0.965.

Thus, based on the results found within this study, it is recommended to use the vector b to represent polarimetric cloud radar observations for applications requiring the error covariance matrix. This representation has a better characterization of random errors in comparison with widely used representation c. When the signal-to-noise ratio is high (>35 dB), however, the variances are quite low, and the Taylor approximation may give reasonable results. We would like to emphasize that there is no additional processing required to get the vector b. Elements of the vector b are an intermediate processing step on the way from I/Q data to conventional spectral polarimetric variables and thus have been already calculated by Doppler cloud radars with the hybrid mode.

In order to demonstrate a practical application of the developed characterization of the measurements errors, a few retrieval techniques are currently being developed. The first one is an improvement of the ice-shape retrieval described in Myagkov et al. (2016a). Another one is an adoption of the drop size distribution retrieval from Tridon and Battaglia (2015) for dual-polarimetric cloud radar observations.

Appendix A: Diagonalization matrix Q

The operator Q, which is used to diagonalize the covariance matrix B in Eq. (20), is calculated as follows (Kanareykin et al., 1968, chap.  2.5):

$$\begin{array}{}\text{(A1)}& \mathbf{Q}=\begin{array}{c}\left(\begin{array}{lr}{q}_{\mathrm{11}}& {\dot{q}}_{\mathrm{12}}\\ -{\dot{q}}_{\mathrm{12}}^{*}& q_{\mathrm{11}}\end{array}\right)\end{array},\end{array}$$

where

$$\begin{array}{}\text{(A2)}& {\displaystyle }& {\displaystyle }{q}_{\mathrm{11}}={\left(\mathrm{1}+{\left|\dot{d}\right|}^{\mathrm{2}}\right)}^{-\mathrm{0.5}},\text{(A3)}& {\displaystyle }& {\displaystyle }{\dot{q}}_{\mathrm{12}}=-{\dot{d}}^{*}{q}_{\mathrm{11}},\text{ and}\text{(A4)}& {\displaystyle }& {\displaystyle }d={\displaystyle \frac{{\dot{B}}_{\mathrm{hv}}^{*}}{\mathrm{0.5}\left[\mathrm{Tr}\mathbf{B}+\sqrt{\mathrm{Tr}{\mathbf{B}}^{\mathrm{2}}-\mathrm{4}\det \left(\mathbf{B}\right)}\right]-B_{\mathrm{vv}}}}.\end{array}$$

In Eq. (A4) Tr is the matrix trace.

Appendix B: Derivation of likelihood functions

B1 Change in variables in a PDF

Consider a vector a with n random variables a_1…n. Assume the joint PDF f_a(a) of the variables is known. The joint PDF f_y(y) of a vector

$$\begin{array}{}\text{(B1)}& \mathit{y}=G\left(\mathit{a}\right)\end{array}$$

can be found by changing the variables in f_a(a):

$$\begin{array}{}\text{(B2)}& f_y\left(\mathit{y}\right)=\left|\mathbf{J}\right|f_a\left[G^{-\mathrm{1}}\left(\mathit{y}\right)\right],\end{array}$$

where G⁻¹ is the reverse transformation from y to a, and J is the determinant of the Jacobian of the transformation $\mathit{a}=G^{-\mathrm{1}}\left(\mathit{y}\right)$.

B2 Likelihood functions for D_cc and D_xx

It is known that the PDF of z_s being a sum of squares of independent standard normal samples (i.e., distributed normally with a mean of 0 and standard deviation of 1) is the chi-squared distribution ${\mathit{\chi }}_k^{\mathrm{2}}\left(z_{\mathrm{s}}\right)$, where the degree of freedom k shows how many samples have been summed. Taking into account that

$$\begin{array}{}\text{(B3)}& \begin{array}{rl}{\widehat{D}}_{cc}& =N_{\mathrm{s}}^{-\mathrm{1}}{\mathit{\sigma }}_c^{\mathrm{2}}\left\{{\mathit{\sigma }}_c^{-\mathrm{2}}\sum _{l=\mathrm{1}}^{N_{\mathrm{s}}}\mathrm{Re}{\left({\dot{S}}_c\right)}_l^{\mathrm{2}}\right.\\ & +\left.{\mathit{\sigma }}_c^{-\mathrm{2}}\sum _{l=\mathrm{1}}^{N_{\mathrm{s}}}\mathrm{Im}{\left({\dot{S}}_c\right)}_l^{\mathrm{2}}\right\},\end{array}\end{array}$$

where the first and the second summed terms in the curly brackets are sums of squares of independent standard normal samples, the likelihood function $f\left({\widehat{D}}_{cc}|{\mathit{\sigma }}_c,N_{\mathrm{s}}\right)$ can be found by changing the variable z_s to $N_{\mathrm{s}}{\mathit{\sigma }}_c^{-\mathrm{2}}{\widehat{D}}_c$:

$$\begin{array}{}\text{(B4)}& f\left({\widehat{D}}_{cc}|{\mathit{\sigma }}_c,N_{\mathrm{s}}\right)={\displaystyle \frac{N_{\mathrm{s}}}{{\mathit{\sigma }}_c^{\mathrm{2}}}}{\mathit{\chi }}_{\mathrm{2}{N}_{\mathrm{s}}}^{\mathrm{2}}\left({\displaystyle \frac{N_{\mathrm{s}}}{{\mathit{\sigma }}_c^{\mathrm{2}}}}{\widehat{D}}_{cc}\right).\end{array}$$

The factor of 2 in the degree of freedom is because there are 2N_s summed components in the curly brackets in Eq. (B3). The equation for ${\widehat{D}}_{xx}$ is derived in a similar manner as for ${\widehat{D}}_{cc}$, resulting in

$$\begin{array}{}\text{(B5)}& f\left({\widehat{D}}_{xx}|{\mathit{\sigma }}_x,N_{\mathrm{s}}\right)={\displaystyle \frac{N_{\mathrm{s}}}{{\mathit{\sigma }}_x^{\mathrm{2}}}}{\mathit{\chi }}_{\mathrm{2}{N}_{\mathrm{s}}}^{\mathrm{2}}\left({\displaystyle \frac{N_{\mathrm{s}}}{{\mathit{\sigma }}_x^{\mathrm{2}}}}{\widehat{D}}_{xx}\right).\end{array}$$

B3 Likelihood functions for R_cx and J_cx

Nadarajah and Pogány (2016) provide a solution for the PDF of an averaged multiplication z_m of two standard normal variables. For two uncorrelated variables the PDF is defined as follows:

$$\begin{array}{}\text{(B6)}& \begin{array}{rl}{f}_z\left(z_{\mathrm{m}}\right)=& {\displaystyle \frac{n^{(n+\mathrm{1})/\mathrm{2}}{\mathrm{2}}^{(\mathrm{1}-n)/\mathrm{2}}|z_{\mathrm{m}}{|}^{(n-\mathrm{1})/\mathrm{2}}}{\sqrt{\mathit{\pi }}\mathrm{\Gamma }(n/\mathrm{2})}}\\ \times & K_{(\mathrm{1}-n)/\mathrm{2}}\left(n\left|z_{\mathrm{m}}\right|\right),\end{array}\end{array}$$

where n is the number of averaged multiplications, Γ is the gamma function, and K_μ is the Bessel function of the second kind of order μ.

${\widehat{R}}_{cx}$ is calculated as follows:

$$\begin{array}{}\text{(B7)}& \begin{array}{rl}{\widehat{R}}_{cx}& =\mathrm{2}{\mathit{\sigma }}_c{\mathit{\sigma }}_x\left\{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{N}_{\mathrm{s}}{\mathit{\sigma }}_c{\mathit{\sigma }}_x}}\left[\sum _{l=\mathrm{1}}^{N_{\mathrm{s}}}\mathrm{Re}\left({\dot{S}}_c\right)\mathrm{Re}\left({\dot{S}}_x\right)\right.\right.\\ & \left.\left.+\sum _{l=\mathrm{1}}^{N_{\mathrm{s}}}\mathrm{Im}\left({\dot{S}}_c\right)\mathrm{Im}\left({\dot{S}}_x\right)\right]\right\},\end{array}\end{array}$$

where the term in the curly brackets is an average over 2N_s multiplications of independent standard normal samples. In this case, the likelihood function $f\left({\widehat{R}}_{cx}|{\mathit{\sigma }}_c,{\mathit{\sigma }}_x,N_{\mathrm{s}}\right)$ can be found by changing z_m by $(\mathrm{2}{\mathit{\sigma }}_c{\mathit{\sigma }}_x{)}^{-\mathrm{1}}{\widehat{R}}_{cx}$:

$$\begin{array}{}\text{(B8)}& \begin{array}{rl}f\left({\widehat{R}}_{cx}|{\mathit{\sigma }}_c,{\mathit{\sigma }}_x,N_{\mathrm{s}}\right)=& {\displaystyle \frac{\left(\mathrm{2}{N}_{\mathrm{s}}{)}^a\right|{\widehat{R}}_{cx}{|}^{-b}}{\sqrt{\mathit{\pi }}{\mathrm{2}}^{\mathrm{2}{N}_{\mathrm{s}}}\left({\mathit{\sigma }}_c{\mathit{\sigma }}_x{)}^a\mathrm{\Gamma }\right(N_{\mathrm{s}})}}{K}_b\\ \times & \left({\displaystyle \frac{N_{\mathrm{s}}\left|{\widehat{R}}_{cx}\right|}{{\mathit{\sigma }}_c{\mathit{\sigma }}_x}}\right),\end{array}\end{array}$$

where $a=\left(\mathrm{2}{N}_{\mathrm{s}}+\mathrm{1}\right)/\mathrm{2}$, $b=\left(\mathrm{1}-\mathrm{2}{N}_{\mathrm{s}}\right)/\mathrm{2}$, Γ is the gamma function, and K_μ is the Bessel function of the second kind of order μ. When ${\widehat{R}}_{cx}\to \mathrm{0}$, the modified Bessel function $K_b\left(N_{\mathrm{s}}\left|{\widehat{R}}_{cx}\right|({\mathit{\sigma }}_c{\mathit{\sigma }}_x{)}^{-\mathrm{1}}\right)\to \mathrm{\infty }$. Therefore, for ${\widehat{R}}_{cx}$ close to 0, the following approximation based on Eqs. (9.6.6) and (9.6.8) from Abramowitz and Stegun (1972) should be used:

$$\begin{array}{}\text{(B9)}& f\left({\widehat{R}}_{cx}|{\mathit{\sigma }}_c,{\mathit{\sigma }}_x,N_{\mathrm{s}}\right)\approx {\displaystyle \frac{N_{\mathrm{s}}\mathrm{\Gamma }(-b)}{\mathrm{2}\sqrt{\mathit{\pi }}{\mathit{\sigma }}_c{\mathit{\sigma }}_x\mathrm{\Gamma }\left(N_{\mathrm{s}}\right)}}.\end{array}$$

Formulas for ${\widehat{J}}_{cx}$ are defined in a similar manner:

$$\begin{array}{}\text{(B10)}& \begin{array}{rl}f\left({\widehat{J}}_{cx}|{\mathit{\sigma }}_c,{\mathit{\sigma }}_x,N_{\mathrm{s}}\right)=& {\displaystyle \frac{\left(\mathrm{2}{N}_{\mathrm{s}}{)}^a\right|{\widehat{J}}_{cx}{|}^{-b}}{\sqrt{\mathit{\pi }}{\mathrm{2}}^{\mathrm{2}{N}_{\mathrm{s}}}\left({\mathit{\sigma }}_c{\mathit{\sigma }}_x{)}^a\mathrm{\Gamma }\right(N_{\mathrm{s}})}}{K}_b\\ \times & \left({\displaystyle \frac{N_{\mathrm{s}}\left|{\widehat{J}}_{cx}\right|}{{\mathit{\sigma }}_c{\mathit{\sigma }}_x}}\right).\end{array}\end{array}$$

The approximation for ${\widehat{J}}_{cx}$ is close to 0:

$$\begin{array}{}\text{(B11)}& f\left({\widehat{J}}_{cx}|{\mathit{\sigma }}_c,{\mathit{\sigma }}_x,N_{\mathrm{s}}\right)\approx {\displaystyle \frac{N_{\mathrm{s}}\mathrm{\Gamma }(-b)}{\mathrm{2}\sqrt{\mathit{\pi }}{\mathit{\sigma }}_c{\mathit{\sigma }}_x\mathrm{\Gamma }\left(N_{\mathrm{s}}\right)}}.\end{array}$$

B4 Monte Carlo evaluation of Eqs. (B4), (B5), (B8), and (B10)

Table B1 Percentage of test-statistic values exceeding critical values for different significance levels. Percentages are given in percent. The names of the four columns on the right side of the table indicate the distribution for which a percentage is given.

For the equation evaluation a simulated dataset was generated. In total 1000 sets of distributions were simulated using the Monte Carlo approach. A single set included distributions of ${\widehat{B}}_{cc}$, ${\widehat{B}}_{xx}$, ${\widehat{R}}_{cx}$, and ${\widehat{J}}_{cx}$. For a single set 10⁵ vector $\widehat{\mathit{b}}$'s were generated. A single vector $\widehat{\mathit{b}}$ resulted from N_s randomly generated vector m's. For a single set of distributions a single covariance matrix B was taken. The elements of the covariance matrix B and N_s were randomly generated according to the following rules (values have linear arbitrary units):

1. B_hh is a sum of mean powers of signal P_sh and noise P_nh.

2. B_vv is a sum of mean powers of signal P_sv and noise P_nv.

3. $P_{\mathrm{nh}}=P_{\mathrm{nv}}=\mathrm{1}$.

4. P_sh and P_sv were randomly and independently generated using the uniform distribution from 1 to 5.

5. ${\dot{B}}_{\mathrm{hv}}$ was calculated as ${\mathit{\rho }}_{\mathrm{HV}}{e}^{i{\mathrm{\Phi }}_{\mathrm{DP}}}\sqrt{P_{\mathrm{sh}}{P}_{\mathrm{sv}}}$.

6. ρ_HV was chosen randomly using the uniform distribution from 0 to 1.

7. Φ_DP was chosen randomly using the uniform distribution from 0 to 2π.

8. N_s was chosen as a random integer number in the range of 2 to 80.

From the covariance matrix B the true covariance matrix Σ_m was obtained. A total of 10⁵×N_s vector m's were generated according to the PDF given in Eq. (2). Then, 10⁵ elements of the $\widehat{\mathit{b}}$ were calculated according to Eqs. (15)–(18). Elements of the vector $\widehat{\mathit{d}}$ were derived from the vector $\widehat{\mathit{b}}$'s using Eq. (27).

Using the 10⁵ vector $\widehat{\mathit{d}}$'s individual histograms for each of the variables ${\widehat{B}}_{cc}$, ${\widehat{B}}_{xx}$, ${\widehat{R}}_{cx}$, and ${\widehat{J}}_{cx}$ are derived. A histogram has 10 bins covering the range from the minimum to maximum values of the corresponding variable. Widths of bins were adjusted to have 10 000 samples in each bin. For the same bins the expected number of samples is calculated using the corresponding PDF. Since integration of Eqs. (B4), (B5), (B8), and (B10) is challenging, the integration is done numerically. Then the Pearson's chi-squared test is applied. The same procedure is repeated for all 1000 sets of distributions. Thus, for each PDF (Eqs. B4, B5, B8, and B10) 1000 test-statistic values were obtained.

The Pearson's chi-squared test implies a comparison of the test-statistic values with critical values for a given level of significance. A test-statistic value exceeding the critical value would indicate that there is a chance (equal to the significance level) that the data significantly differ from the PDF. There is, however, a small chance that the conclusion that the data differ from the PDF is erroneous. Table B1 shows the percentage of the test-statistic values exceeding critical values. It can be seen that the number of test-statistic values exceeding corresponding critical values is very close to the theoretical values, i.e., 5, 2.5, and 1 % at 0.95, 0.975, and 0.99 significance levels, respectively. This confirms the validity of the obtained PDFs.

Appendix C: Variances of elements of the vector $\widehat{\mathit{d}}$

To derive solutions for the mean and variances of elements of $\widehat{\mathit{d}}$, the distribution of the elements is represented by characteristic functions. A γth raw statistical moment M_γ of a random variable with a characteristic function ϕ(t) can be found as follows:

$$\begin{array}{}\text{(C1)}& M_{\mathit{\gamma }}=i^{-\mathit{\gamma }}{\displaystyle \frac{d^{\mathit{\gamma }}\mathit{\varphi }\left(t\right)}{\mathrm{d}{t}^{\mathit{\gamma }}}}{|}_{t=\mathrm{0}}.\end{array}$$

The calculation of derivatives of the characteristic functions is in general easier to obtain than integration of the corresponding PDFs.

The characteristic function of the chi-squared distribution ${\mathit{\chi }}_k^{\mathrm{2}}\left(z_{\mathrm{s}}\right)$ is

$$\begin{array}{}\text{(C2)}& {\mathit{\varphi }}_s\left(t\right)={\left(\mathrm{1}-\mathrm{2}it\right)}^{-k/\mathrm{2}}.\end{array}$$

Therefore, the characteristic function for ${\widehat{D}}_{cc}$ for a given σ_c and N_s can be written in the following way:

$$\begin{array}{}\text{(C3)}& {\mathit{\varphi }}_{cc}\left(t\right)={\left(\mathrm{1}-{\displaystyle \frac{\mathrm{2}i{\mathit{\sigma }}_c^{\mathrm{2}}t}{N_{\mathrm{s}}}}\right)}^{-{\mathrm{N}}_{\mathrm{s}}}.\end{array}$$

The mean value and variance of ${\widehat{D}}_{cc}$ are calculated as follows:

$$\begin{array}{}\text{(C4)}& {\displaystyle }\stackrel{\mathrm{‾}}{{\widehat{D}}_{cc}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{i}}{\displaystyle \frac{d{\mathit{\varphi }}_{cc}\left(t\right)}{\mathrm{d}t}}{|}_{t=\mathrm{0}}=\mathrm{2}{\mathit{\sigma }}_c^{\mathrm{2}},\text{(C5)}& {\displaystyle }\mathrm{var}\left({\widehat{D}}_{cc}\right)& {\displaystyle }=-{\displaystyle \frac{d^{\mathrm{2}}{\mathit{\varphi }}_{cc}\left(t\right)}{\mathrm{d}{t}^{\mathrm{2}}}}{|}_{t=\mathrm{0}}-{\stackrel{\mathrm{‾}}{{\widehat{D}}_{cc}}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{4}{\mathit{\sigma }}_c^{\mathrm{4}}}{N_{\mathrm{s}}}}.\end{array}$$

Similarly,

$$\begin{array}{}\text{(C6)}& {\displaystyle }\stackrel{\mathrm{‾}}{{\widehat{D}}_{xx}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{i}}{\displaystyle \frac{d{\mathit{\varphi }}_{xx}\left(t\right)}{\mathrm{d}t}}{|}_{t=\mathrm{0}}=\mathrm{2}{\mathit{\sigma }}_x^{\mathrm{2}},\text{ and}\text{(C7)}& {\displaystyle }\mathrm{var}\left({\widehat{D}}_{xx}\right)& {\displaystyle }={\displaystyle \frac{\mathrm{4}{\mathit{\sigma }}_x^{\mathrm{4}}}{N_{\mathrm{s}}}}.\end{array}$$

Based on Nadarajah and Pogány (2016) the characteristic function corresponding to f_z(z_m) is

$$\begin{array}{}\text{(C8)}& {\mathit{\varphi }}_z\left(t\right)={\left(\mathrm{1}+{\displaystyle \frac{t^{\mathrm{2}}}{n^{\mathrm{2}}}}\right)}^{-n/\mathrm{2}}.\end{array}$$

Therefore, the characteristic function for ${\widehat{R}}_{cx}$ and ${\widehat{J}}_{cx}$ for given σ_c, σ_x, and N_s is as follows:

$$\begin{array}{}\text{(C9)}& {\mathit{\varphi }}_{cx}\left(t\right)={\left(\mathrm{1}+{\displaystyle \frac{{\mathit{\sigma }}_c^{\mathrm{2}}{\mathit{\sigma }}_x^{\mathrm{2}}{t}^{\mathrm{2}}}{N_{\mathrm{s}}^{\mathrm{2}}}}\right)}^{-{\mathrm{N}}_{\mathrm{s}}/\mathrm{2}}.\end{array}$$

As expected for a multiplication of two uncorrelated variables, the mean values of ${\widehat{R}}_{cx}$ and ${\widehat{J}}_{cx}$ are as follows:

$$\begin{array}{}\text{(C10)}& \stackrel{\mathrm{‾}}{{\widehat{R}}_{cx}}=\stackrel{\mathrm{‾}}{{\widehat{J}}_{cx}}={\displaystyle \frac{\mathrm{1}}{i}}{\displaystyle \frac{d{\mathit{\varphi }}_{cx}\left(t\right)}{\mathrm{d}t}}{|}_{t=\mathrm{0}}=\mathrm{0}.\end{array}$$

The variance of ${\widehat{R}}_{cx}$ and ${\widehat{J}}_{cx}$ can be found as follows:

$$\begin{array}{}\text{(C11)}& \mathrm{var}\left({\widehat{R}}_{cx}\right)=\mathrm{var}\left({\widehat{J}}_{cx}\right)=-{\displaystyle \frac{d^{\mathrm{2}}{\mathit{\varphi }}_{cx}\left(t\right)}{\mathrm{d}{t}^{\mathrm{2}}}}{|}_{t=\mathrm{0}}={\displaystyle \frac{{\mathit{\sigma }}_c^{\mathrm{2}}{\mathit{\sigma }}_x^{\mathrm{2}}}{\mathrm{2}{N}_{\mathrm{s}}}}.\end{array}$$

Appendix D: Sensitivity S

$$\begin{array}{rl}\mathbf{S}& =\left(\begin{array}{cccc}\frac{\partial B_{\mathrm{hh}}}{\partial B_{\mathrm{hh}}}& \frac{\partial B_{\mathrm{hh}}}{\partial R_{\mathrm{hv}}}& \frac{\partial B_{\mathrm{hh}}}{\partial J_{\mathrm{hv}}}& \frac{\partial B_{\mathrm{hh}}}{\partial B_{\mathrm{vv}}}\\ \frac{\partial Z_{\mathrm{DR}}}{\partial B_{\mathrm{hh}}}& \frac{\partial Z_{\mathrm{DR}}}{\partial R_{\mathrm{hv}}}& \frac{\partial Z_{\mathrm{DR}}}{\partial J_{\mathrm{hv}}}& \frac{\partial Z_{\mathrm{DR}}}{\partial B_{\mathrm{vv}}}\\ \frac{\partial {\mathit{\rho }}_{\mathrm{HV}}}{\partial B_{\mathrm{hh}}}& \frac{\partial {\mathit{\rho }}_{\mathrm{HV}}}{\partial R_{\mathrm{hv}}}& \frac{\partial {\mathit{\rho }}_{\mathrm{HV}}}{\partial J_{\mathrm{hv}}}& \frac{\partial {\mathit{\rho }}_{\mathrm{HV}}}{\partial B_{\mathrm{vv}}}\\ \frac{\partial {\mathrm{\Phi }}_{\mathrm{DP}}}{\partial B_{\mathrm{hh}}}& \frac{\partial {\mathrm{\Phi }}_{\mathrm{DP}}}{\partial R_{\mathrm{hv}}}& \frac{\partial {\mathrm{\Phi }}_{\mathrm{DP}}}{\partial J_{\mathrm{hv}}}& \frac{\partial {\mathrm{\Phi }}_{\mathrm{DP}}}{\partial B_{\mathrm{vv}}}\end{array}\right)\\ & =\left(\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ B_{vv}^{-\mathrm{1}}& \mathrm{0}& \mathrm{0}& -B_{hh}{B}_{vv}^{-\mathrm{2}}\\ -\mathrm{0.5}\left|{\dot{B}}_{\mathrm{hv}}\right|B_{vv}^{-\mathrm{0.5}}{B}_{hh}^{-\mathrm{1.5}}& R_{\mathrm{hv}}\left|{\dot{B}}_{\mathrm{hv}}{|}^{-\mathrm{1}}\right(B_{hh}{B}_{vv}{)}^{-\mathrm{0.5}}& J_{\mathrm{hv}}\left|{\dot{B}}_{\mathrm{hv}}{|}^{-\mathrm{1}}\right(B_{hh}{B}_{vv}{)}^{-\mathrm{0.5}}& -\mathrm{0.5}\left|{\dot{B}}_{\mathrm{hv}}\right|B_{hh}^{-\mathrm{0.5}}{B}_{vv}^{-\mathrm{1.5}}\\ \mathrm{0}& -J_{\mathrm{hv}}|{\dot{B}}_{\mathrm{hv}}{|}^{-\mathrm{2}}& R_{\mathrm{hv}}|{\dot{B}}_{\mathrm{hv}}{|}^{-\mathrm{2}}& \mathrm{0}\end{array}\right)\text{(D1)}\end{array}$$

Appendix E: Jacobian J_bd of the transformation from $\widehat{\mathit{b}}$ to $\widehat{\mathit{d}}$

Using Eqs. (21)–(22) J_bd can be written as follows:

$$\begin{array}{rlr}{\mathbf{J}}_{bd}& =\left|\begin{array}{cccc}\frac{\partial D_{cc}}{\partial B_{\mathrm{hh}}}& \frac{\partial D_{cc}}{\partial R_{\mathrm{hv}}}& \frac{\partial D_{cc}}{\partial J_{\mathrm{hv}}}& \frac{\partial D_{cc}}{\partial B_{\mathrm{vv}}}\\ \frac{\partial R_{cx}}{\partial B_{\mathrm{hh}}}& \frac{\partial R_{cx}}{\partial R_{\mathrm{hv}}}& \frac{\partial R_{cx}}{\partial J_{\mathrm{hv}}}& \frac{\partial R_{cx}}{\partial B_{\mathrm{vv}}}\\ \frac{\partial J_{cx}}{\partial B_{\mathrm{hh}}}& \frac{\partial J_{cx}}{\partial R_{\mathrm{hv}}}& \frac{\partial J_{cx}}{\partial J_{\mathrm{hv}}}& \frac{\partial J_{cx}}{\partial B_{\mathrm{vv}}}\\ \frac{\partial D_{xx}}{\partial B_{\mathrm{hh}}}& \frac{\partial D_{xx}}{\partial R_{\mathrm{hv}}}& \frac{\partial D_{xx}}{\partial J_{\mathrm{hv}}}& \frac{\partial D_{xx}}{\partial B_{\mathrm{vv}}}\end{array}\right|& =\left|\begin{array}{cccc}{q}_{\mathrm{11}}^{\mathrm{2}}& |{\dot{q}}_{\mathrm{12}}{|}^{\mathrm{2}}& -\mathrm{2}{q}_{\mathrm{11}}{R}_{\mathrm{12}}& -\mathrm{2}{q}_{\mathrm{11}}{J}_{\mathrm{12}}\\ |{\dot{q}}_{\mathrm{12}}{|}^{\mathrm{2}}& q_{\mathrm{11}}^{\mathrm{2}}& \mathrm{2}{q}_{\mathrm{11}}{R}_{\mathrm{12}}& \mathrm{2}{q}_{\mathrm{11}}{J}_{\mathrm{12}}\\ q_{\mathrm{11}}{R}_{\mathrm{12}}& -q_{\mathrm{11}}{R}_{\mathrm{12}}& q_{\mathrm{11}}^{\mathrm{2}}-R_{\mathrm{12}}^{\mathrm{2}}+J_{\mathrm{12}}^{\mathrm{2}}& -\mathrm{2}{R}_{\mathrm{12}}{J}_{\mathrm{12}}\\ q_{\mathrm{11}}{J}_{\mathrm{12}}& -q_{\mathrm{11}}{J}_{\mathrm{12}}& -\mathrm{2}{R}_{\mathrm{12}}{J}_{\mathrm{12}}& q_{\mathrm{11}}^{\mathrm{2}}+R_{\mathrm{12}}^{\mathrm{2}}-J_{\mathrm{12}}^{\mathrm{2}}\end{array}\right|\\ & =(q_{\mathrm{11}}^{\mathrm{2}}+|q_{\mathrm{12}}{|}^{\mathrm{2}}{)}^{\mathrm{4}}\end{array}\text{(E1)}$$

Taking into account Eqs. (A2) and (A3), J_bd=1.

Appendix F: Jacobian J_cb of the transformation from $\widehat{\mathit{c}}$ to $\widehat{\mathit{b}}$

Using Eqs. (10)–(12) J_cb can be written as follows:

$$\begin{array}{rl}{\mathbf{J}}_{cb}& =\left|\begin{array}{cccc}\frac{\partial B_{\mathrm{hh}}}{\partial B_{\mathrm{hh}}}& \frac{\partial B_{\mathrm{hh}}}{\partial Z_{\mathrm{DR}}}& \frac{\partial B_{\mathrm{hh}}}{\partial {\mathit{\rho }}_{\mathrm{HV}}}& \frac{\partial B_{\mathrm{hh}}}{\partial {\mathrm{\Phi }}_{\mathrm{DP}}}\\ \frac{\partial R_{\mathrm{hv}}}{\partial B_{\mathrm{hh}}}& \frac{\partial R_{\mathrm{hv}}}{\partial Z_{\mathrm{DR}}}& \frac{\partial R_{\mathrm{hv}}}{\partial {\mathit{\rho }}_{\mathrm{HV}}}& \frac{\partial R_{\mathrm{hv}}}{\partial {\mathrm{\Phi }}_{\mathrm{DP}}}\\ \frac{\partial J_{\mathrm{hv}}}{\partial B_{\mathrm{hh}}}& \frac{\partial J_{\mathrm{hv}}}{\partial Z_{\mathrm{DR}}}& \frac{\partial J_{\mathrm{hv}}}{\partial {\mathit{\rho }}_{\mathrm{HV}}}& \frac{\partial J_{\mathrm{hv}}}{\partial {\mathrm{\Phi }}_{\mathrm{DP}}}\\ \frac{\partial B_{\mathrm{vv}}}{\partial B_{\mathrm{hh}}}& \frac{\partial B_{\mathrm{vv}}}{\partial R_{\mathrm{hv}}}& \frac{\partial B_{\mathrm{vv}}}{\partial J_{\mathrm{hv}}}& \frac{\partial B_{\mathrm{vv}}}{\partial B_{\mathrm{vv}}}\end{array}\right|\\ & =\left|\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ {\mathit{\rho }}_{\mathrm{HV}}\cos \left({\mathrm{\Phi }}_{\mathrm{DP}}\right)Z_{\mathrm{DR}}^{-\mathrm{0.5}}& -\mathrm{0.5}{B}_{\mathrm{hh}}{\mathit{\rho }}_{\mathrm{HV}}\cos \left({\mathrm{\Phi }}_{\mathrm{DP}}\right)Z_{\mathrm{DR}}^{-\mathrm{1.5}}& B_{\mathrm{hh}}\cos \left({\mathrm{\Phi }}_{\mathrm{DP}}\right)Z_{\mathrm{DR}}^{-\mathrm{0.5}}& -B_{\mathrm{hh}}{\mathit{\rho }}_{\mathrm{HV}}\sin \left({\mathrm{\Phi }}_{\mathrm{DP}}\right)Z_{\mathrm{DR}}^{-\mathrm{0.5}}\\ {\mathit{\rho }}_{\mathrm{HV}}\sin \left({\mathrm{\Phi }}_{\mathrm{DP}}\right)Z_{\mathrm{DR}}^{-\mathrm{0.5}}& -\mathrm{0.5}{B}_{\mathrm{hh}}{\mathit{\rho }}_{\mathrm{HV}}\sin \left({\mathrm{\Phi }}_{\mathrm{DP}}\right)Z_{\mathrm{DR}}^{-\mathrm{1.5}}& B_{\mathrm{hh}}\sin \left({\mathrm{\Phi }}_{\mathrm{DP}}\right)Z_{\mathrm{DR}}^{-\mathrm{0.5}}& -B_{\mathrm{hh}}{\mathit{\rho }}_{\mathrm{HV}}\cos \left({\mathrm{\Phi }}_{\mathrm{DP}}\right)Z_{\mathrm{DR}}^{-\mathrm{0.5}}\\ Z_{\mathrm{DR}}^{-\mathrm{1}}& -B_{\mathrm{hh}}{Z}_{\mathrm{DR}}^{-\mathrm{2}}& \mathrm{0}& \mathrm{0}\end{array}\right|\\ & =-B_{\mathrm{hh}}^{\mathrm{3}}{Z}_{\mathrm{DR}}^{-\mathrm{3}}{\mathit{\rho }}_{\mathrm{HV}}\end{array}\text{(F1)}$$

Appendix G: Table of symbols

Table G1 Main symbols used throughout the study. The overdot indicates a complex number. Indices h and v indicate the polarization of the receiver channel. The circumflex indicates a measured quantity.

Bringi and Chandrasekar (2001)Bringi and Chandrasekar (2001)Bringi and Chandrasekar (2001)Bringi and Chandrasekar (2001)