weatherTimeSeries

Simulate I/Q signals for weather returns using a Monte Carlo approach

Since R2024a

Syntax

[Vh,Vv] = weatherTimeSeries(freq,Pr,Vr,SW,ZDR,Rho,Phi)

[Vh,Vv] = weatherTimeSeries(freq,Pr,Vr,SW,ZDR,Rho,Phi,Name=Value)

Description

This function generates I/Q signals for monostatic polarimetric weather radar systems. Each weatherTimeSeries simulation produces an independent radar return from a particular resolution volume, or cell, that is modeled as a complex stationary Gaussian random process using Monte Carlo method (see Algorithms).

[Vh,Vv] = weatherTimeSeries(freq,Pr,Vr,SW,ZDR,Rho,Phi) returns I/Q signals as complex voltages for horizontal and vertical polarizations for a monostatic polarimetric weather radar system. Resolution cells are characterized by polarimetric variables specified by Pr, Vr, SW, ZDR, Rho, and Phi.

[Vh,Vv] = weatherTimeSeries(freq,Pr,Vr,SW,ZDR,Rho,Phi,Name=Value) returns I/Q signals for a monostatic polarimetric weather radar with additional options specified using one or more name-value arguments. For example, NumPulses sets the number of pulses and NumTrials controls the number of trials.

Examples

weatherTimeSeries IQ Simulation Example

This example generates weatherTimeSeries simulations and displays the radar estimates alongside the truth values.

Simulate Weather Time Series I/Q Data

Simulate 100 trials of time series I/Q signals for a resolution volume collected using a polarimetric weather radar operating at a frequency of 3 GHz with a PRF of 1000 Hz and 64 pulses in the dwell time. Assume the resolution volume is comprised of pure rain which has a signal power of -80 dBW, a radial velocity of 10 m/s, a spectrum width of 4 m/s, a differential reflectivity of 2 dB, a correlation coefficient of 0.99, and a differential phase of 60 degrees.

% Define radar and weather parameters
freq = 3e9;                     % Operating frequency (Hz)
Pr = -80;                       % Signal power (dBW)
Vr = 10;                        % Radial velocity (m/s)
SW = 4;                         % Spectrum width (m/s)
ZDR = 2;                        % Differential reflectivity (dB)
Rho = 0.99;                     % Correlation coefficient
Phi = 60;                       % Differential phase (degree)
prf = 1000;                     % PRF (Hz)
numPulses = 64;                 % Number of pulses per trial
numTrials = 100;                % Number of trials
    
% Generate time series
[Vh,Vv] = weatherTimeSeries(freq,Pr,Vr,SW,ZDR,Rho,Phi,...
    PRF=prf,NumPulses=numPulses,NumTrials=numTrials); 

Verify Simulation Results

Verify that the radar estimates derived from the simulated I/Q are consistent with the truth.

% Calculate covariance matrix
R0_hh = mean(abs(Vh.^2)); 
R0_vv = mean(abs(Vv.^2));
C0_hv = mean(conj(Vh).*Vv); 
R1_hh = mean(conj(Vh(1:numPulses-1,:)).*(Vh(2:numPulses,:))); 
 
% Pulse Pair Processing
lambda = freq2wavelen(freq);
Pr1 = pow2db(R0_hh);
Vr1 = -lambda*prf/(4*pi)*angle(R1_hh); 
SW1 = abs(lambda*prf/(4*pi).*(2*log(abs(R0_hh)./abs(R1_hh))).^0.5); 
ZDR1 = pow2db(R0_hh./R0_vv);
Rho1 = abs(C0_hv)./(R0_hh).^0.5./(R0_vv).^0.5;
Phi1 = mod(angle(C0_hv)*180/pi,360);
 
% Display estimation result and compare with truth
VariableName = {'Pr';'Vr';'SW';'ZDR';'Rho';'Phi'};
Truth = [Pr;Vr;SW;ZDR;Rho;Phi];
Estimate = [round(mean(Pr1),2);round(mean(Vr1),2);round(mean(SW1),2);...
    round(mean(ZDR1),2);round(mean(Rho1),3);round(mean(Phi1),2)];
Unit = {'dBW';'m/s';'m/s';'dB';'dimensionless';'degree'};
T = table(VariableName,Truth,Estimate,Unit)

T=6×4 table
    VariableName    Truth    Estimate          Unit       
    ____________    _____    ________    _________________

      {'Pr' }        -80      -80.1      {'dBW'          }
      {'Vr' }         10       9.85      {'m/s'          }
      {'SW' }          4       3.96      {'m/s'          }
      {'ZDR'}          2       1.99      {'dB'           }
      {'Rho'}       0.99      0.991      {'dimensionless'}
      {'Phi'}         60      59.97      {'degree'       }

Input Arguments

freq — Radar operating frequency

positive scalar

Radar operating frequency, specified as a positive scalar. Units are in hertz (Hz).

Data Types: double

Pr — Signal power

scalar | length-N vector

Signal power, specified as a scalar or length-N vector where N is the number of resolution cells. Units are in decibel-watt (dBW). Pr is related to the radar constant C, target range R, and reflectivity factor Z by Pr = Z -20log₁₀(R) + C, where C is given in dB, R is in km, and Z is in dBZ [1].

Data Types: double

Vr — Radial velocity

scalar | length-N vector

Radial velocity of the weather returns relative to the radar system, specified as a scalar or length-N vector where N is the number of resolution cells. Units are in meters per second (m/s). The absolute value of Vr should not exceed the maximum unambiguous velocity of the radar. The maximum unambiguous velocity is equal to the radar wavelength (in m) times the pulse repetition frequency (in Hz) divided by 4.

Data Types: double

SW — Spectrum width

nonnegative scalar | nonnegative length-N vector

Spectrum width, a measure of the velocity dispersion, specified as a scalar or length-N vector, where N is the number of resolution cells. Units are in meters per second (m/s).

Data Types: double

ZDR — Differential reflectivity

scalar | length-N vector

Differential reflectivity between horizontally and vertically polarized signals, specified as a scalar or length-N vector, where N is the number of resolution cells. Units are in decibel (dB).

Data Types: double

Rho — Correlation coefficient

positive scalar | 0 < Rho ≤ 1 | positive length-N vector

Correlation coefficient between horizontally and vertically polarized signals specified as a scalar or length-N vector, where N is the number of resolution cells. The value of Rho, which is commonly denoted as ρ_hv, must be greater than 0 and less than or equal to 1. Units are dimensionless.

Data Types: double

Phi — Differential phase

scalar | 0 ≤ Phi < 360 | length-N vector

Differential phase between horizontally and vertically polarized signals, specified as a scalar or length-N vector, where N is the number of resolution cells. The value of Phi, which is commonly denoted as $${\varphi }_{DP}$$, must be greater than or equal to 0 and less than 360. Units are in degrees (deg).

Data Types: double

Name-Value Arguments

Specify optional pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is the argument name and Value is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Example: PRF=1000,NumPulses=64,NumTrials=100

PRF — Pulse repetition frequency

1e3 (default) | positive scalar

Pulse repetition frequency (PRF), specified as a positive scalar. Units are in hertz (Hz).

Data Types: double

NumPulses — Number of pulses within the radar dwell time

64 (default) | positive integer scalar

Number of pulses per trial, specified as a positive integer scalar. Units are dimensionless.

Data Types: double

NumTrials — Number of trials

1 (default) | positive integer scalar

Number of trials, specified as a positive integer scalar. Units are dimensionless.

Data Types: double

Output Arguments

Vh — Complex voltage for horizontal polarization

N-by-M-by-L array

Complex voltage for horizontal polarization returned as a complex-valued N-by-M-by-L array where N is the number of resolution cells, M is the number of pulses (NumPulses), and L is the number of trials (NumTrials). By default, the number of pulses is set to 64, and the number of trials is set to 1. The resolution cells are characterized by polarimetric variables as input arguments Pr, Vr, SW, ZDR, Rho, and Phi. Each trial produces an independent radar return from a particular resolution volume, or cell, that is modeled as a complex stationary Gaussian random process using Monte Carlo method.

Data Types: double

Vv — Complex voltage for vertical polarization

N-by-M-by-L array

Complex voltage for vertical polarization returned as a complex-valued N-by-M-by-L array where N is the number of resolution cells, M is the number of pulses (NumPulses), and L is the number of trials (NumTrials). By default, the number of pulses is set to 64, and the number of trials is set to 1. The resolution cells are characterized by polarimetric variables as input arguments Pr, Vr, SW, ZDR, Rho, and Phi. Each trial produces an independent radar return from a particular resolution volume, or cell, that is modeled as a complex stationary Gaussian random process using Monte Carlo method.

Data Types: double

Algorithms

weatherTimeSeries simulates I/Q signals using a numerical Monte Carlo approach combined with a statistical model that is based on the expected behavior of radar returns from weather phenomena [2] [3]. weatherTimeSeries calculates scattering amplitudes for every time step using relevant scattering parameters specified as input arguments to generate a time series.

Statistical model assumptions for radar returns from weather phenomena:

• I/Q signals follow Gaussian distribution

• Signal amplitude follows Rayleigh distribution

• Signal phase follows uniform distribution

The first step in the Monte Carlo approach is to specify scattering amplitudes and generate random numbers for particle position and motion. Scattering amplitudes for a given resolution volume are modeled as combinations from random multiple scattering centers. Assume that scatterers within the resolution volume have the same size and a canting angle of zero (which implies that off-diagonal elements of the backscattering matrices have a value of zero). Then we have

$$P_r=N{\left|V_{hh}\right|}^2,$$

where P_r is the received power for horizontal polarization, N is the number of scatterers in the resolution volume, and V_hh is the complex voltage for an individual scatterer for horizontal polarization. At least 20 scatterers (N ≥ 20) are necessary to adequately model Gaussian random signals.

The amplitude of V_hh can be obtained as

$$\left|V_{hh}\right|=\sqrt{P_r/N}.$$

Accordingly, the amplitude of the complex voltage for an individual scatterer for vertical polarization can be calculated as

$$\left|V_{vv}\right|=\left|V_{hh}\right|/\sqrt{Z_{dr}},$$

where Z_dr is the differential reflectivity on a linear scale.

In addition, the correlation coefficient ρ_hv is assumed to be reduced by a factor of $${\text{e}}^{{}^{-{\sigma }_{\delta }^2/2}}$$ due to the random scattering phase difference, where σ_δ is the standard deviation of random scattering phase difference, that is, $${\rho }_{hv}={\text{e}}^{{}^{-{\sigma }_{\delta }^2/2}}$$. Then we have

$${\sigma }_{\delta }=\sqrt{-2\ln \left({\rho }_{hv}\right)}=\sqrt{{\sigma }^2{}_{\delta h}+{\sigma }^2{}_{\delta v}}.$$

For simplicity, it is assumed that $${\sigma }_{\delta h}={\sigma }_{\delta v}={\sigma }_{\delta }/\sqrt{2}$$, where σ_δh and σ_δv are standard deviations of the backscattering phase δ_h for horizontal polarization and the backscattering phase δ_v for vertical polarization, respectively.

Next, the scattered wave field for each particle is calculated. The amplitudes of complex voltage for an individual scatterer can be expressed as

$$V_{hh}=\left|V_{hh}\right|\cdot {\text{e}}^{-j{\delta }_h}\cdot {\text{e}}^{{}^{-j{\varphi }_{DP}}}$$

$$V_{vv}=\left|V_{vv}\right|\cdot {\text{e}}^{-j{\delta }_v},$$

where $${\varphi }_{DP}$$ is the differential phase.

The final step is to calculate the total scattered wave fields. The total complex voltage of the resolution volume is

$$V_h={\displaystyle {\sum }_{l=1}^N{V}_{hh}}\cdot {\text{e}}^{-j2k_i\cdot r_l}$$

$$V_v={\displaystyle {\sum }_{l=1}^N{V}_{vv}}\cdot {\text{e}}^{-j2k_i\cdot r_l},$$

where $$k_i$$ is the incident wave vector and $$r_l$$ is the random position of each scatterer.

References

[1] Doviak, R. and D. S. Zrnic. Doppler Radar and Weather Observations, 2nd Ed. New York: Dover Publications, 2006.

[2] Zhang, G. Weather Radar Polarimetry. Boca Raton: CRC Press, 2016.

[3] Li, Z, et al. Polarimetric phased array weather radar data quality evaluation through combined analysis, simulation, and measurements. IEEE Geosci. Remote Sens. Lett., vol. 18, no. 6, pp. 1029–1033, Jun. 2021.

Version History

Introduced in R2024a