其中，的单位是dB/km，的单位是，单位是dBZ。这个关系易受其他因素影响，常不稳定，如校准calibration、部分波束遮挡partial beam blockage。系数a和b是需要由雨滴谱观测数据进行回波模拟进行拟合而来，没有这个数据的时候先根据默认系数来计算就好：对于C波段，根据wradlib库，a=1.67e-4, b=0.7；对于X波段，根据Delrieu et al (1997)的图3，计算出来a=1.409e-4，b=0.76。
defcalc_attenuation_hb(gateset, a=1.67e-4, b=0.7, gate_length=1.0): """Gate-by-Gate forward correction as described in :cite:`Kraemer2008` Parameters ---------- gateset : :class:`numpy:numpy.ndarray` Multidimensional array, where the range gates (over which iteration has to be performed) are supposed to vary along the last array-dimension. Data has to be provided in decibel representation of reflectivity [dBZ]. a : float proportionality factor of the k-Z relation (:math:`k=a \\cdot Z^{b}`). Per default set to 1.67e-4. b : float exponent of the k-Z relation ( :math:`k=a \\cdot Z^{b}` ). Per default set to 0.7. gate_length : float length of a range gate [km]. Per default set to 1.0. Returns ------- pia : :class:`numpy:numpy.ndarray` Array with the same shape as ``gateset`` containing the calculated path integrated attenuation [dB] for each range gate. """ gateset = np.where(np.isnan(gateset), -999, gateset) pia = np.zeros(gateset.shape) for gate in range(gateset.shape[-1] - 1): k = ( a * idecibel(gateset[..., gate] + pia[..., gate]) ** b * 2.0 * gate_length ) pia[..., gate + 1] = pia[..., gate] + k pia = np.where(np.isnan(pia), 0, pia) return pia
2.2 MK方法
针对HB的不稳定问题，可设置一个最大的允许修正量对此进行约束。Jacobi and Heistermann (2016)比较了3种对HB方法进行约束的方案，发现MK（modified Kraemer）方法最好。其调整了幂律关系的参数，使整个数据订正后的反射率不超过59 dBZ，PIA上限设为20 dB。
MK方法（Jacobi and Heistermann 2016）是基于Hitschfeld and Bordan（1954）迭代方法的逐门衰减订正，有两个约束条件：允许的最大为59 dBZ，最大PIA为20 dB，以避免过订正。代码见wradlib库的函数atten.correct_attenuation_constrained。在不符合上述约束条件的情况下，用于计算衰减的相邻波束的数量sector_thr的默认值为10。
其中，单位deg/km，和单位dB/deg，单位dB/km。上述表达式常被认为是线性的，即，这在大部分气象雷达波段都是合理的。式中系数与雨滴谱分布DSD、温度和滴谱大小-形状关系有关。由05-07年的2D滴谱仪数据、T矩阵计算出的散射幅度（温度10°C、波长3.22 cm的情况下）、Brandesetal（2002）给出的滴谱的size–shape relation，和分别为0.313和0.0483。在Park et al（2005）进行的X波段雷达模拟中，随着雨滴形状的变化差异很大，其数值在0.173-0.315 dB/deg的范围内变化（本文参考CHUVA feild campaign，时，取0.219）。根据Overeem et al (2021)，C波段为$0.081。
DP方法中不准确的也会导致PIA的误差在径向求和过程中积累。
defpia_from_kdp(kdp, dr, gamma=0.313): """Retrieving path integrated attenuation from specific differential \ phase (Kdp). The default value of gamma is based on :cite:`Carey2000`. Parameters ---------- kdp : :class:`numpy:numpy.ndarray` array specific differential phase
Range dimension must be the last dimension. dr : float gate length (km) gamma : float linear coefficient (default value: 0.08 for C band) in the relation between Kdp phase and specific attenuation (alpha) 0.313 for X band Returns ------- output : :class:`numpy:numpy.ndarray` array of same shape as kdp containing the path integrated attenuation """ alpha = gamma * kdp pia = 2 * np.cumsum(alpha, axis=-1) * dr #进行累加求和，即积分 return pia
2.4 ZPHI方法
ZPHI方法通过雨区的总变化来约束，然后根据径向上的分布来分摊这个衰减（Bringi and Chandrasekar2001），因此会比HB方法稳定得多，其名称“ZPHI”也是因为既用到了Z又用到PHI来进行约束。ZPHI需要两个先验值：关系中的系数[，见式(4)]和关系中的指数[b，见式(3)]。Gorgucci and Chandrasekar (2005) 和Tabary et al (2011)表明ZPHI方法一般比DP表现更好。
对于，ZPHI算法通过归一化的截距参数建立了和的关系：。其中，根据Testud et al（2000），X波段, 。
defpia_ZPHI(gateset, kdp, dr, b = 0.76, gamma = 0.313, debug = False): """Retrieving path integrated attenuation using ZPHI methodwhich \ constrains PIA_H by the total change in φDP along a radial through a rain cell. The attenuation is then apportioned according to the distribution of Z_H along the radial. The default value of gamma is based on :cite:`FIGURE 6.21 in book by Zhang`. Parameters ---------- gateset : :py:class:`numpy:numpy.ndarray` multidimensional array. The range gates (over which iteration has to be performed) are supposed to vary along the *last* dimension. The azimuth are supposed to vary along the 2nd last dimension. so, e.g., for a set of `l` radar images stored in polar form with `m` azimuths and `n` range-bins the input array's shape should be either (l,m,n). data has to be provided in decibel representation of reflectivity [dBZ] kdp : :class:`numpy:numpy.ndarray` array specific differential phase Range dimension must be the last dimension. dr : float gate length (km) gamma : float linear coefficient (default value: 0.08 for C band, 0.313 for X band) in the relation between Kdp phase and specific attenuation (alpha) b : float exponent of the k-Z relation ( :math:`k=a \\cdot Z^{b}` ). Per default
set to 0.76. """
gateset = np.where(np.isnan(gateset), -999, gateset) alpha = gamma * kdp pia = 2 * np.cumsum(alpha, axis=-1) * dr pia = np.where(np.isnan(pia), 0, pia) # needed? rN = np.argmax(pia, axis = -1) # solving `a` assuming that the attenuation is apportioned according to the distribution of ZH along the radial. a, numerator, denominator = np.zeros(gateset.shape[:-1]), np.zeros(gateset.shape[:-1]), np.zeros(gateset.shape[:-1])
for iray in range(gateset.shape[-2] - 1):
if rN[iray] > 0: numerator[iray] = 1 - np.exp( -0.23 * b * pia[iray,rN[iray]]) tmp = 0.46 * b * np.cumsum(np.power( idecibel(gateset[iray,0:rN[iray]]), b ), axis = -1) * dr denominator[iray] = tmp[-1] if denominator[iray] > 100: a[iray] = numerator[iray] / denominator[iray] #a越大，衰减越厉害 if debug: print('iray is %i, range : %i, last echo : %.1f, coefficient a : %.5f ' \ %(iray, rN[iray], idecibel(gateset[iray,0:rN[iray]][-1]), a[iray]) )
pia = np.zeros(gateset.shape) for iray in range(gateset.shape[-2] - 1): for gate in range(gateset.shape[-1] - 1): k = ( a[iray] * idecibel(gateset[iray, gate] + pia[iray, gate]) ** b * 2.0 * dr ) pia[iray, gate + 1] = pia[iray, gate] + k pia = np.where(np.isnan(pia), 0, pia) return a, pia
2.5 自适应约束算法(self-consistent with constraints, SCWC, Bringi et al 2001）
defpia_scwc(gateset, kdp, phi_dp, dr, gamma_max = 0.35, gamma_min = 0.13 , delta_gamma = 0.01, b = 0.76): """ calculate the Path Integrated Attenuation using self-consistent with constraint method. Parameters ---------- gateset : :py:class:`numpy:numpy.ndarray` multidimensional array. The range gates (over which iteration has to be performed) are supposed to vary along the *last* dimension. The azimuth are supposed to vary along the 2nd last dimension. so, e.g., for a set of `l` radar images stored in polar form with `m` azimuths and `n` range-bins the input array's shape should be either (l,m,n). data has to be provided in decibel representation of reflectivity [dBZ] kdp : :class:`numpy:numpy.ndarray` array specific differential phase Range dimension must be the last dimension. phi_dp : :class:`numpy:numpy.ndarray` array differential_phase should be the reconstructed phi_dp (filetered), phi_dp[0]should be 0. dr : float gate length (km) gamma_max : float iterate from gamma_max to gamma_min with the step of delta_gamma, to find the optimal gamma, which is the linear coefficient in the relation between Kdp phase and specific attenuation b : float exponent of the A-Z relation ( :math:`A=a \\cdot Z^{b}` ). Per default set to 0.76. Returns ------- gamma_opt : :class:`numpy:numpy.ndarray` the coefficient of the A-K_DP retation (math: `A-\gamma*K_DP`) for all irays. a_opt : the coefficient a of the A-Z retation (math: `A=a \\cdot Z^{b}`) for all irays. pia : :class:`numpy:numpy.ndarray` array of same shape as kdp containing the path integrated attenuation """ gateset = np.where(np.isnan(gateset), -999, gateset) alpha = gamma_max * kdp pia = 2 * np.cumsum(alpha, axis=-1) * dr pia = np.where(np.isnan(pia), 0, pia) # needed? rN = np.argmax(pia, axis = -1)
Jacobi, Stephan, and Maik Heistermann. "Benchmarking attenuation correction procedures for six years of single-polarized C-band weather radar observations in South-West Germany." Geomatics, Natural Hazards and Risk 7.6 (2016): 1785-1799.
Snyder, Jeffrey C., et al. "Attenuation correction and hydrometeor classification of high-resolution, X-band, dual-polarized mobile radar measurements in severe convective storms." Journal of Atmospheric and Oceanic Technology 27.12 (2010): 1979-2001.
Overeem, Aart, et al. "Rainfall-induced attenuation correction for two operational dual-polarization c-band radars in the Netherlands." Journal of Atmospheric and Oceanic Technology 38.6 (2021): 1125-1142.
Bringi, V. N., and V. Chandrasekar, 2001: Polarimetric Doppler Weather Radar: Principle and Applications. Cambridge Uni- versity Press, 636 pp.
Gu, Ji-Young, et al. "Polarimetric attenuation correction in heavy rain at C band." Journal of Applied Meteorology and Climatology 50.1 (2011): 39-58.
Ryzhkov, Alexander, et al. "Potential utilization of specific attenuation for rainfall estimation, mitigation of partial beam blockage, and radar networking." Journal of Atmospheric and Oceanic Technology 31.3 (2014): 599-619.
Delrieu, G., S. Caoudal, and J. D. Creutin. "Feasibility of using mountain return for the correction of ground-based X-band weather radar data." Journal of Atmospheric and Oceanic Technology 14.3 (1997): 368-385.
Lengfeld, Katharina, M. Berenguer, and D. Sempere Torres. "Intercomparison of attenuation correction algorithms for single-polarized X-band radars." Atmospheric Research 201 (2018): 116-132.
