Technique for estimating rainfall from a meteorological radar with polarization diversity

[0003] The present invention relates to meteorological techniques for estimating a precipitation rate by means of a radar.

[0004] For cost reasons, the radars of the European network operate in C-band, a frequency which is substantially attenuated on traversing internal rains, this making the retrieval of the precipitation rate much trickier than in S-band where the attenuation effects are more limited.

[0005] In applications to urban hydrology, the X-band is even envisaged, but attenuation by rain then plays a very important role, and it becomes crucial to take it into account.

[0006] In this regard, reference may advantageously be made to their publication: HITSCHFELD-BORDAN, 1954, “Errors inherent in the radar measurement of rainfall at attenuating wavelengths”, J. Meteor, 11, pp. 58-67. Unfortunately, in an actual application where the data are noisy and where one may be in the presence of a calibration error (even a small one), this solution turns out to be numerically unstable.

[0007] To alleviate this drawback, a proposal has already been made to operate under polarization diversity and to utilize the differential reflectivity ZDR, which is the ratio of the reflectivities under H and V polarization.

[0008] As it falls, a raindrop undergoes aerodynamic flattening, the consequence of which is that its scattering cross section is larger for H polarization than for V polarization. In a radar operating at attenuated frequency, the parameter ZDR results from two effects with opposite tendencies: on the one hand, the effect of differential backscattering (which tends to increase ZDR with the mean size of the drops), and on the other hand, the effect of differential attenuation (with the opposite tendency).

[0009] For a recent example of determining a precipitation rate from a dual-polarization radar, reference may advantageously be made to patent application FR-2 742 876, as well as to the publication: “Polarimetric Radar at Attenuated Wavelenghts as a Hydrological sensor”—M. SAUVAGEOT—Journal of etmosyheric and oceanic technology—vol. 13—p. 630-637, 1996.

[0010] However, it is not easy to measure the coefficient ZDR.

[0011] The dynamic range of variation of ZDR is around 1 to 2 and therefore covers a few dB at the very most. Furthermore, this measurement requires cross-calibration of the H and V channels to within a tenth of a dB; integration over numerous independent samples so as to temper the statistical fluctuation of the signal, and a very good performance of the antenna (sidelobes at 30 dB—outward—below the main lobe).

[0012] Alternatively, there has also been a proposal, in particular in the publication: “Differential propagation phase shift and rainfall rate estimation”—M. SACHINANDA, D.S. ZRNIC—Radio science—vol. 21, no. 2, p. 235-247, March-April 1986. to utilize another parameter, namely the differential phase (denoted Φdp). This is because the flattening of the raindrop does not affect only the backscattering cross section inducing the asymmetry σ_H>σ_V(where σ_Hand σ_Vrepresent the cross sections for H and V polar, respectively), but also the propagation of the radar wave whose wave vector is affected as much in respect of its imaginary part (the specific attenuation is greater for H polarization than for V polarization, differential attenuation effect mentioned earlier), as in respect of its real part (where the asymmetry between H and V translates into a differential variation of the phase along the path. In actual fact, the derivative of Φdp along the radius (rate of variation of Φdp denoted Kdp) is almost proportional to the precipitation rate (and hence Φdp, to the integrated precipitation rate), which explains the benefit of measuring it.

[0013] The estimator of Kdp exhibits numerous advantages: it is insensitive to attenuation along the path; it is insensitive to a radar calibration error; it is much less affected than ZDR by the sidelobes of the antenna; it is entirely unaffected by partial masking of the antenna beam (which occurs routinely when operating at low elevation).

[0014] However, it exhibits the major drawback of being greatly affected by noise since it results from differentiating Φdp along the beam. The measurement of Kdp therefore requires a long integration time, which is incompatible with the hydrological application which demands fast scanning of the radar beam.

[0015] For its part, the invention proposes a technique which makes it possible to alleviate the drawbacks of the prior techniques and which implements simple processing which is reliable and robust, in particular with regard to the statistical measurement noise (thereby allowing fast scanning).

[0016] More particularly, the invention proposes a process for estimating a precipitation rate by means of a bipolar radar, characterized by the following various steps:

[0017] the differential phase Φdp and the attenuated reflectivity Z according to at least one of the polarizations H or V are measured by means of said bipolar radar, over a given interval [r₁, r₀] of path radius r with respect to said radar,

[0018] an estimate of the value K(r₀) of the attenuation at r₀is determined from the attenuated reflectivity profile thus measured, as well as from the difference in the differential phase between r₀and r₁;

[0019] an estimate K(r) of the specific attenuation at r as a function of the attenuation K(r₀) thus determined and of the attenuated reflectivity profile Z(r) is determined;

[0020] the rate of precipitation R(r) is determined knowing K(r).

[0021] Advantageously, an estimate of the value K(r₀) of the attenuation r₀is determined from the equation $K\left(r_0\right){\int }_{r_1}^{r_0}\frac{Z^b\left(s\right)}{Z^b\left(r_0\right)+K\left(r_0\right)·I\left(s,r_0\right)}ds=\gamma ·\left({\phi }_{\mathrm{dp}}\left(r_0\right)-{\phi }_{\mathrm{dp}}\left(r_1\right)\right)$

[0022] where:

[0023] $I\left(s,r_0\right)=0.46b{\int }_{r_1}^{r_0}{Z}^b\left(u\right)du$

[0024] and b is the exponent of the power relation

[0025] K(r)=aZ(r)^b

[0026] and where γ is the parameter of proportionality between the specific attenuation and the differential rate of variation of the phase.

[0027] Preferably, an estimate K(r) of the specific attenuation is determined as a function of r from the equation $\frac{Z^b\left(r\right)}{K\left(r\right)}-\frac{Z^b\left(r_0\right)}{K\left(r_0\right)}=I\left(r,r_0\right).$

[0028] Advantageously, the attenuated reflectivity Z is measured for both of the two polarizations H or V and the precipitation rate R(r) for a given path radius r is determined for both of these two polarizations.

[0029] The invention also relates to a device for estimating a precipitation rate comprising a bipolar radar, as well as processing means, characterized in that said radar comprises means for measuring the differential phase Φdp and the attenuated reflectivity Z according to at least one of the polarizations H or V and in that the processing means implement the various processing steps of the process according to claim 1.

[0030] Other characteristics and advantages of the invention will emerge further from the description which follows, which is purely illustrative and non-limiting and should be read in conjunction with the appended drawings in which:

[0031] FIGS. 1aand 1bare graphs on which are plotted scattergrams illustrating the relevance of a parametric representation used in one mode of implementation of the invention;

[0032] FIG. 2 illustrates a device for implementing the process in accordance with the invention;

[0033] FIGS. 3ato 3c, 4ato 4c, 5ato 5care graphs, in which it has been assumed that the radar beam intercepts two rain cells of gaussian form and on which have been plotted the parameters Z, K and R for various values of the parameter No (defined later in the description) of dimensional distribution of the raindrops.

[0034] FIG. 6 is a graph on which has been plotted a simulation of variations for the coefficient N₀* as a function of the path radius r

[0035] FIGS. 7. to 7care graphs on which have been plotted, for the case where N₀* is variable, the values obtained for the parameters Z, K and R;

[0036] FIGS. 8ato 8cand FIGS. 9ato 9care graphs on which have been plotted the parameters Z, K and R obtained for absolute biases of 2dB and 5 dB respectively and for a value of N₀and a capture equal to 8000;

[0037] FIGS. 10ato 10cand 11ato 11care figures similar to those of FIGS. 8ato 8cand 9ato 9cin the case of relative biases taken equal to 0.2 dB and 0.5 dB respectively;

[0038] FIGS. 12ato 12cand 13ato 13care graphs on which have been plotted the values obtained for the parameters Z, K and R in the case of retrieval from a number of samples equal to 60 and from a number of samples equal to 10;

[0039] FIGS. 14ato 14care graphs on which have been plotted the values obtained for the parameters Z, K and R in the case of a maximum rainfall rate of 200 mm/hour and of a noisy signal;

[0040] FIGS. 15ato 15dare histograms illustrating a possible calibration technique;

[0041] FIGS. 16ato 16care scattergrams illustrating another possible calibration technique.

[0042] The process which will now be described implements a processing on the differential phase Φdp and the attenuated reflectivity Z_H(or Z_V) according to the H (or V) polarization (reflectivity Z in the remainder of the text).

[0043] In a first step, an interval [r₁, r₀] in which the data Z and Φdp are continuously available, is determined for the radial distance r from the radar antenna.

[0044] It is known that the Hitschfeld-Bordan equation may be written: $\begin{array}{cc}{u}^{\phi -1}\frac{du}{dr}=-0.46a{Z}^b& \left(1\right)\end{array}$

[0045] with u=Z/Z₀where Z₀is the unattenuated reflectivity and Z the apparent reflectivity, and where a and b are the two parameters of the power law K=a Z₀^b, in which K is the specific attenuation. It is recalled that a varies strongly, whereas b is stable.

[0046] The solution of this equation (1) is, when integrating between two limits r₁and r₂:

[0047] [u^b]_r^r=−0.46ab∫_r^rZ^bds  (2)

[0048] By expressing (2) between r and r₀, it follows that: $\begin{array}{cc}\frac{Z^b\left(r\right)}{K\left(r\right)}-\frac{Z^b\left(r_0\right)}{K\left(r_0\right)}=I\left(r,r_0\right)& \left(3\right)\end{array}$

[0049] where

[0050] I(r,r₀)=0.46b∫_r^rZ^bds  (4)

[0051] and where K=a Z₀^b

[0052] Equation (3) provides an estimator of K(r) which we will seek to constrain through the Φdp information while adjusting K(r⁰). To do this, (3) is transformed into $\begin{array}{cc}K\left(r\right)=K\left(r_0\right)\frac{Z^b\left(r\right)}{Z\left(r_0\right)+K\left(r_0\right)·I\left(r,r_0\right)}& \left(5\right)\end{array}$

[0053] Now, the derivative of the phase Φdp with respect to the path r, or rate of variation of the differential phase Kdp, is a quasi-linear function of the specific attenuation K (that is to say K_Hor K_V, depending on the polarization with which one is working).

[0054] Consequently,

[0055] ∫_r^rK(s)ds=γ(Φ_dp(r₀)−Φ_dp(r₁))  (6)

[0056] where γ is a known proportionality parameter i.e.: $\begin{array}{cc}K\left(r_0\right){\int }_{r_1}^{r_0}\frac{Z^b\left(s\right)}{Z\left(r_0\right)+K\left(r_0\right)·I\left(s,r_0\right)}ds=\gamma ·\left({\phi }_{\mathrm{dp}}\left(r_0\right)-{\phi }_{\mathrm{dp}}\left(r_1\right)\right)& \left(7\right)\end{array}$

[0057] i.e., after integration: $\begin{array}{cc}K\left(r_0\right)=\frac{Z^b\left(r_0\right)}{I\left(r_1,r_0\right)}\left\{{10}^{0.2b\gamma ·\mathrm{\Delta \phi }}-1\right\}& \left(8\right)\end{array}$

[0058] where ΔΦ=Φ_dp(r₀)−dp(r1).

[0059] From this estimate of K(r₀), we determine K(r) at every point of the interval [r₁, r₀] from equation (3).

[0060] We thus find: $\begin{array}{cc}K\left(r\right)=\frac{Z^b\left(r\right)·\left\{{10}^{0.2b\gamma ·\mathrm{\Delta \phi }}-1\right\}}{I\left(r,r_0\right)+\left\{{10}^{0.2b\gamma ·\mathrm{\Delta \phi }}-1\right\}·I\left(r,r_0\right)}.& \left(9\right)\end{array}$

[0061] We note that the benefit of this formulation is on the one hand that it provides a solution K(r) which is independent of the parameter a which is precisely the one which fluctuates the most, and, on the other hand, that it harnesses a differential measurement of the differential phase, which is much easier to achieve than an absolute measurement.

[0062] This makes it possible to estimate the value of the precipitation rate R(r), by using for example for this purpose charts giving R(r) as a function of K(r) for various values of the parameter a, which can itself be determined from the relation:

[0063] K(r₀)=a Z₀(R₀)^b

[0064] By implementing the processing just described for both of the two polarizations H and V, we end up with two estimates of R which it is possible to amalgamate in order to validate a final value.

[0065] In an implementational variant allowing still greater fineness of analysis, the “universal” parametric relations between the true reflectivity (unattenuated) Z₀(mm⁶.m⁻³), the specific attenuation K (db.km⁻¹), the specific differential phase K_DP(deg.km⁻¹), and the precipitation rate R (mm.hr⁻¹) are used; These are written:

[0066] K=a[N₀*]^(1−d).Z₀^b

[0067] R=c[N₀*]^(1−d).K^d

[0068] K=γ/2K_DP

[0069] where the parameters a, b, c, d, γ are known (for a given temperature), and where N₀* is a parameter characteristic of the dimensional distribution N(D) [number density of raindrops per interval of diameter in m⁻⁴] of the raindrops which one is proposing to determine by the algorithm N₀* (m⁻⁴) is expressed as a function of the liquid water content W (g/m⁻³) and of the median volumic diameter D₀as:

[0070] N₀*=[181/(πρ_w]W/D₀

[0071] (where ρ_wis the density of water in g/m³).

[0072] The relevance of this parametric representation is illustrated by FIGS. 1aand 1bon which have respectively been represented on the one hand the scattergram of the R-Z relation obtained from 11853 spectra N(D) measured by in-situ probe during the TOGA-COARE airborne experiment and on the other hand the same scattergram after normalizing each of the parameters by N₀* (the TOGA-COARE experiment took place in the west Pacific between November 1992 and February 1993, the PMS sensor being mounted on the NCAR's Electra aircraft. The primary processing of the data was done by R. BLACK from the NOAA (HRD Miami). The data are available on the HRD Miami server).

[0073] An estimator of N₀* is given by: $\begin{array}{cc}{N}_0^{*}={\left[a^{-1}\frac{K\left(r_0\right)}{Z^b\left(r_0\right)+K\left(r_0\right)·I\left(r_1,r_0\right)}\right]}^{\frac{1}{1-b}}& \left(10\right)\end{array}$

[0074] We then calculate R(r) between r₁and r₀by substituting K(r) and the value of N₀* into the relation:

[0075] R=c[N₀*]^(1−d).K^d  (11)

[0076] The processing just described has been tested with a rig of the type illustrated in FIG. 2 which comprises:

[0077] a transmitter 1;

[0078] a 3dB power divider 2 receiving the output signal from the transmitter 1;

[0079] a mode (or orthomode) extractor 3 linked by two channels H and V to the outputs of a power divider 2;

[0080] a horn antenna 4 transmitting or receiving the two H and V polarizations;

[0081] circulators 5 (filtering tee) interposed, on the H and V channels, between the mode extractor 3 and the power divider 2;

[0082] a phase shifter 6 interposed, on the V channel, between the power divider 2 and the mode extractor 3;

[0083] reception means 7 linked, for the H and V channels, to the circulators 5;

[0084] an oscillator 9 which synchronizes the transmitter 1 and the reception means 9;

[0085] means 8 for sampling, and then for processing the I and Q outputs from these reception means 7, these processing means 8 allowing the determination from the parameters ZH, ZV and ΦDP, which in complex notation satisfy:

[0086] Z_H,V=[I_II,V²+Q_H,V²]^½

[0087] and

[0088] ΦDP=1 arg [(I_H+i_QH)·(I_V−iQ_v)]

[0089] where I_H,vand Q_H,vare the values obtained by sampling the I and Q outputs of the H and V channels.

[0090] It is these processing means 8 which implement the determination processing described hereinabove.

[0091] The operating parameters of this device were as follows:

[0092] detection threshold of the radar at 0 dB SNR and at 10 km: 0 dBZ;

[0093] number of independent samples N_iin the resolution cell: 10 or 60

[0094] radial resolution: 75 metres

[0095] noise in the measurement of Φ_DP: ±1 degree (for N_i=60), degrees ±3 degrees (for N_i=10)

[0096] dimensional distribution of the drops exponential: N(D)=N₀exp (−ΛD) [where D is the mean diameter of the drops and where diameter Λ is a parameter], with N₀constant along the radial (with the following values: 800, 8000 [Marshall and Palmer value], 80,000, or continuously variable proportionally to R.

[0097] Under these conditions, the following various results were obtained:

[0098] 1) N₀* Constant

[0099] FIGS. 3ato 3c, 4ato 4c, 5ato 5cdepict the retrieved rainfall rates for constant N₀* over the radial and equal to 800, 8000, 80,000 respectively. In these three cases, it may be seen that the retrieved rainfall rates are very similar to the simulated rainfall rates, and that the retrieved value of the parameter N₀* is on each occasion relatively similar to that used for the simulations.

[0100] 2) N₀* Varying with Rainfall Rate

[0101] The inventors have also tested the processing on simulations which are still realistic but carried out with an N₀* which depends on the rainfall rate, this complying fairly well with reality but rendering the retrieval of the precipitation rate “fairly perilous”. The relation used between N₀* and the precipitation rate is the following: $N_0^{*}=\frac{{10}^{\log \left(R\right)+66}}{1000}$

[0102] For rainfall rates varying from 0 to 80 mm/h, it may be seen that the parameter N₀* sweeps a very wide range, covering more than one decade. By way of indication, presented in FIG. 6 are the variations of N₀* over the simulated radial. The retrievals obtained with the processing described earlier are depicted in FIGS. 6ato 6c. The maximum deviation observed is of the order of 20%, which seems to be reasonable in the case of “such unfavorable” simulations.

[0103] 3) Effects of Relative and Absolute Bias

[0104] FIGS. 8ato 8cand 9 to 9cdepict the results obtained with an absolute bias of 2 dB and 5 dB respectively. A value of absolute bias of 2 dB induces a non-negligible error in the estimation of the parameter N₀*, but does not penalize the retrieval since the retrieved rainfall rates differ by only a few per cent from those which were simulated. On the other hand, an absolute bias of 5 dB seems to be crippling. The results obtained in the case of relative biases of 0.2 dB and 0.5 dB, and which are depicted in FIGS. 10ato 10cand 11 to 11c, show the weak sensitivity of this algorithm to a relative calibration error between the two polarizations.

[0105] 4) Effects of Noise.

[0106] The inventors have tested the effect of noise on the retrievals in the case of a number of samples of 60 (FIGS. 12ato 12c) and in the less favorable case of 10 samples (FIGS. 13ato 13c). In these two cases, the measurement the measurement of the integrated phase exhibits noise of 3 degrees. Retrieval is effected so long as the signal-to-noise ratio is sufficient. For a maximum rainfall rate of 80 mm/h, such is the case up to a radial distance of around 19 km. The results obtained confirm the very stable behavior of the KDP algorithm. Under these conditions, correlation coefficients of greater than 0.95 and biases of less than 4 mm/h are entirely satisfactory.

[0107] 5) Case of Very High Rainfall Rates

[0108] The inventors have verified the robustness of the algorithm in the case of high rainfall rates, reaching 200 mm/h and a noisy signal. In this case, the signal is fairly rapidly swamped by the noise, and it obviously becomes impossible to retrieve the second precipitation cell on account of the very strong attenuation of the radar signal by the first. However, the rainfall rates obtained up to around 10 km are consistent and the whole of the first precipitation cell is very adequately retrieved (FIGS. 14ato 14c).

[0109] The processing has also been tested in bands other than the X-band and in particular in the C-band and S-band. The results obtained are just as satisfactory as in the X-band.

[0110] Thus, with the processing just described, one has a precipitation rate estimator which is particularly stable with regard to noise and/or with regard to the bias in the reflectivities and in the phase. This processing succeeds, even in difficult cases (high rainfall rate, No varying rapidly with rainfall rate) in retrieving the rainfall rate in a very encouraging manner.

[0111] The processing of the data is possible as long as the signal-to-noise ratio is greater than 1. Under these conditions, the bias in the estimator of R does not exceed 5% and its standard deviation 10%. The performance of the estimator is independent of N₀* (when it is assumed to be constant along the radial), and is only slightly impaired under the conditions of variable N₀* . Absolute calibration of the radar to ±2 or 3 dBZ is sufficient to ensure such performance. This calibration is not necessary which to refine the R-K relation from which one estimates R. No specification is necessary as regards the relative calibration of the H and V channels.

[0112] Several calibration techniques may be envisaged.

[0113] According to one of them, the histogram of the N₀* observed by spectropluviometer on the ground is compared with the histogram of the N₀* retrieved by the radar over a large number of radii. In the presence of a reflectivity calibration error δC, we then measure δC.Z believing that we are measuring Z), the estimator of N₀* given by (10) is then multiplied by the factor (δC)^−b/(1−b).

[0114] By correlating the histogram measured on the ground and that measured by the radar, we determine δC. If no ground spectropluviometer is available, fairly precise calibration can be obtained by demanding that the histogram of the N₀*, as measured by the radar, be centered on the Marshall and Palmer value (N₀*=0.8 10⁷m⁻⁴) after re-calibration. An illustration of this calibration technique is given in FIGS. 15ato 15cwhich represent the histogram of N₀* which have been retrieved from a volume scanning of the Darwin radar (polarimetric C-band radar), (a) leaving the measured reflectivities as is, (b), raising them by 1 dB, (c), lowering them by 1 dB.

[0115] By comparing these histograms with the “actual” histogram of N₀*, i.e. as determined in the same region from TOGA-COARE airborne data (see FIG. 15d), it may be seen that the 1 dB lowering enables the histogram obtained with the processing proposed by the invention to be matched to the actual histogram.

[0116] According to another calibration technique, which advantageously supplements the techniques just described, the following other “universal” relations are used:

[0117] (i) between R, K and ZDR (where ZDR=ZH/ZV)

[0118] R=K.(ZDR)^f  (12)

[0119] (ii) between the differential attenuation K_H−K_Vand the attenuation under H polarization, K_H:

[0120] K_H−K_V=α[N₀*]^1−βK_H^β  (13)

[0121] Where α, β and f are known parameters.

[0122] We calculate ZDR from $\begin{array}{cc}\mathrm{ZDR}={\mathrm{ZDR}}_a·{10}^{{\int }_{r_1}^{r_0}\alpha V_{*}^{-1\beta -Z^b}ds}& \left(14\right)\end{array}$

[0123] then we calculate R from (12).

[0124] Two independent estimators of the precipitation rate (the estimator (12) and that given by (11)) are then compared statistically (in a scattergram).

[0125] We then adjust δC to obtain optimal consistency. For the same set of data the scattergrams given in FIGS. 16ato 16cillustrate this technique.

[0126] Plotted in FIG. 16ais a scattergram for the scatter between the estimator (11) [deduced from K alone], and the estimator (12) [deduced from K and from ZDR]. Here the measured reflectivities are used as is, and it is found that (12) overestimates by 12% relative to (11).

[0127] Plotted in FIG. 16bis a scattergram of the scatter between the estimator (11) [deduced from K alone] and the estimator (12) [deduced from K and ZDR]. Here, the measured reflectivities are raised by 1 dB, and it is found that (12) overestimates by 21% relative to (11).

[0128] Plotted in FIG. 16cis a scattergram of the scatter between the estimator (11) [deduced from K alone] and the estimator (12) [deduced from K and ZDR]. Here, the measured reflectivities are lowered by 1 dB, and it is found that there is perfect agreement (to within 2%) between (11) and (12).

[0129] Thus, the best consistency is again obtained by considering a 1 dB lowering of the reflectivities.

[0130] The rainfall rate estimator just described is particularly advantageous in the case of estimating rainfall over drainage basins for the monitoring of highwater floods and the management of water resources. This monitoring and this management nowadays rely on runoff/infiltration models whose main limitation as regards predictive performance stems from poor knowledge of the water sheet entering the system.