引言

中小尺度对流系统空间尺度小、时间演变快、致灾严重，而且很多天气现象出现在中下层，加上复杂地形的影响，因而利用观测周期为6 min、空间分辨率相对低的S波段双偏振天气雷达(S-band dual-polarization weather radar, S-POL)观测和预警灾害天气，可能存在观测盲区和快速演化信息的缺失。为此，我国多个地区布设了X波段双偏振相控阵天气雷达(X-band dual-polarization phased array weather radar, X-PAR)，以弥补S波段雷达网的不足。从探测特点分析，X-PAR的时空分辨率明显高于业务组网的S-POL，但单部X-PAR探测范围小，对大范围、强降水的天气存在观测缺失；S-POL受衰减影响小、探测范围大，数据质量较好，但时空分辨率比较低。如何通过不同波长雷达的组网拼图，获得大范围精细准确的三维雷达观测场，是组网雷达数据应用的关键问题。

对于S或C波段天气雷达网，目前已建立了比较成熟的组网拼图方法。常见的雷达观测网格化插值方法包括最近邻插值、线性内插、Cressman插值(Weygandt et al, 2002)和Barnes插值(Askelson et al, 2000)等。已有的S波段业务组网雷达拼图研究表明，对单部雷达进行水平方向最近邻插值，垂直方向线性插值的网格化方案效果较好，当各雷达观测偏差较小时，利用雷达距离指数加权平均能够较好实现多雷达拼图融合(肖艳姣, 2007; Zhang et al, 2005)。Sun et al(2019)针对三部S波段天气雷达反射率因子，发现与最近邻插值、线性内插和8点插值相比，自适应的Barnes插值方法能够获得更加连续且分辨率较高的网格化插值结果。Langston et al(2007)针对组网雷达观测存在时间差的问题，将插值距离、雷达径向距离和观测时间差三个方面的权重乘积作为拼图插值权重，获得了时空连续且能够保留细节特征的拼图结果。针对雷达观测点分布不均匀的问题，Scovell and Al-Sakka(2016)基于二叉空间分割树原理，结合Barnes连续校正，实现多雷达反射率因子三维拼图，在保留雷达观测数据密集区细节特征的同时避免观测稀疏区的欠平滑问题。针对华南组网S-POL的对比研究显示，相比于先对单部雷达网格化插值再进行多雷达拼图融合的方法，将多部雷达观测直接加权插值能够获得更加合理的拼图结果(李哲, 2022)。除了对雷达观测的加权插值方法外，利用变分方法实现对笛卡尔坐标下雷达观测的求解能够兼顾雷达拼图的连续性和细节特征的保留(Roca-Sancho et al, 2014; Brook et al, 2022)，在雷达观测稀疏区获得较好的拼图结果，但对于密集组网观测地区，加权插值方式更加高效。

针对S(或C)波段的雷达组网拼图研究表明，由于雷达的布设普遍距离较远，观测样本在雷达远端相对稀疏，雷达数据质量问题主要源于远端分辨率下降和波束非均匀(Ryzhkov, 2007)，因此雷达拼图融合主要利用雷达观测距离来设置权重。而对于密集组网布设的X波段天气雷达，其时空分辨率较高，共同覆盖区的观测密度较大，单部雷达的观测范围较小，数据质量的问题较S波段雷达更为复杂。以往针对S(或C)波段雷达网的拼图方法对X波段雷达已不再完全适用。

尽管利用最近邻插值法对多部X波段雷达观测反射率因子进行组网拼图，能够有效补充业务雷达探测的低空盲区(吴翀等, 2016)。但根据数据质量调整置信度阈值，利用雷达观测距离与衰减影响计算拼图融合权重，可以更好发挥X波段雷达观测优势，获得准确精细的组网拼图结果(吴翀等, 2021)。此外，也有研究根据X波段雷达观测距离、积分衰减、零滞后相关系数等反映雷达观测数据质量的参量，确定拼图插值权重进行组网拼图，获得了与附近S波段雷达接近的观测结果(Liu et al, 2007)。针对组网雷达观测时间差的问题，可使用归一化互相关方法估计移动向量，对时间和空间采用线性插值实现单部雷达观测的网格化处理，再进行多雷达拼图(Kim et al, 2019)。对长沙机场阵列天气雷达回波强度融合的试验表明，计算雷达每个方位向和仰角方向分辨率扩展系数，并依次进行插值填充，然后利用均值法对多个收发子阵反射率因子进行融合，能够得到更为准确、精细且细节突出的回波结构(叶开等, 2020)。相比于以往针对S(或C)波段雷达的组网拼图方法，针对X波段天气雷达网的插值拼图方法应避免低质量观测数据对拼图结果的“污染”，同时最大化保留回波细节特征。

已有的拼图方法普遍针对雷达反射率因子，对于双偏振物理量的拼图方法研究较少。此外，X波段天气雷达的优势在于获得更加精细的低层降水观测信息，与S波段天气雷达网形成优势互补，但目前对于S波段与X波段天气雷达网的拼图融合方法研究仍然不足。

目前，国际上开展的综合外场试验初步验证了1分钟分辨率的X波段天气雷达(机械扫描雷达)网在地形复杂地区和重要气象保障区局地强天气监测和预警方面的明显作用，形成了一系列数据处理方法和应用成果(Mahale et al, 2014；Brewster et al, 2017；Diederich et al, 2015；Anagnostou et al, 2018；Chandrasekar et al, 2018)，但X波段相控阵天气雷达并没有进入到这些组网研究中，S和X波段雷达组网技术也只处于综合试验阶段，并没有进行大范围推广。我国已经进入了X波段相控阵天气雷达组网业务应用阶段(程元慧等, 2020；肖靖宇等, 2022；Zhao et al, 2024)，但现有的业务应用还主要依靠单一波长的雷达观测或简单的组网拼图产品，尚未能很好发挥S波段与X波段天气雷达密集组网优势，实现不同波长雷达探测能力的互补，亟需建立多波段雷达组网拼图融合方法，获得相对准确精细的三维雷达观测场，为云降水物理精细结构及其演化特征分析、雷达定量估测降水、强对流灾害天气预报预警等提供数据支撑。

1 数据概况 1.1 雷达参数

本研究选用广东省S波段双偏振天气雷达(S-POL)和X波段双偏振相控阵天气雷达(X-PAR) 资料为主要研究数据。其中S-POL的数据经过了严格的标定，数据质量相对稳定可靠。X-PAR是一款全固态相参、有源相控阵、脉冲多普勒、双偏振天气雷达。S-POL与X-PAR均包含反射率因子Z_H(单位：dBz)、零滞后相关系数ρ_HV、差分反射率因子Z_DR(单位：dB)、双程差分传播相移Φ_DP(单位：°)、差分传播相移率K_DP(单位：°·km^-1)的三维体扫数据。

X-PAR在垂直方向采用相控阵电子扫描的方式，完成一个方位角RHI(range height indicator)扫描后切换到下一个方位角，时间分辨率得到了明显提升。相比于传统机械扫描天气雷达的VCP21模式，相控阵雷达的VRHI观测模式对目标不同高度的观测时差更小，避免了逐仰角PPI(plan position indication)机械扫描方式造成的雷达回波在垂直方向上的不连续。在窄发窄收的观测模式下，X-PAR扫描波束从0.9°开始，以1.8°步进扫描到20.7°，共12层仰角，方位步进0.9°，一次体扫时间为92 s。该X-PAR水平半功率波束宽度≤3.6°，垂直波束宽度≤1.8°，与常规机械扫描天气雷达1°的波束宽度相比，X-PAR的波束宽度很宽，采用了0.9°方位步进的超分辨率I/Q信号处理算法对方位分辨率进行了改善，得到0.9°方位角分辨率的雷达数据。两种型号雷达的主要参数如表 1所示。需要说明的是，X-PAR切换不同观测模式后的最大探测距离、方位步进、仰角范围和时间分辨率均可作相应调整，表 1中只给出了窄发窄收观测模式下的主要参数。

1.2 雷达分布

广东省附近X-PAR的分布概况如图 1所示，广东省及周边区域布设的S-POL如图 2。可见该地区基本均可实现三部以上X-PAR和S-POL的共同覆盖，这为获得精细准确的三维降水雷达观测场提供了有利条件。

1.3 数据质量控制

在进行拼图之前，需要对雷达数据进行质量控制。本研究对S-POL资料利用零滞后相关系数ρ_HV与信噪比数据剔除非降水回波(吴翀, 2019)，利用线性规划方法(Giangrande et al, 2013)，采用9点Savitzky-Golay滤波器，计算得到S-POL的K_DP。基于Geng and Liu(2023a;2023b)对X-PAR的Z_H和K_DP进行质量控制和计算，以S-POL观测数据转换为X波段后为基准，X-PAR质量控制后的Z_H平均偏差在1 dBz以内。由于X-PAR的差分反射率因子Z_DR未保存原始观测数据，本研究只对厂家衰减订正后的Z_DR进行滤波去噪和系统偏差订正。

2 不同波长雷达组网拼图方法

基于S-POL与X-PAR的探测能力，S-POL数据质量较好，但远端空间分辨率下降明显，中高层仰角存在较大的观测空隙，拼图的关键在于利用多雷达组网改善空间分辨率。X-PAR时空分辨率较高，但水平波束较宽，受衰减影响较大，数据质量问题复杂，拼图的关键在于保留精细观测信息的同时，避免质量差的观测数据对拼图结果的污染。为此，本研究提出针对密集组网的S-POL与X-PAR拼图方法。该方法的整体思路为：首先建立时空分辨率低的S-POL拼图和时空分辨率高的X-PAR拼图，然后将两组拼图进行时空匹配(包括空间匹配和S-POL拼图的时间外推)，将可靠性更高的S-POL拼图作为基准，与X-PAR拼图作差得到一组与S-POL空间分辨率相同的粗糙偏差网格，将该粗糙偏差网格进行插值得到与X-PAR空间分辨率一致的精细偏差网格，利用X-PAR拼图与精细偏差网格相加，即可得到整体回波强度与S-POL拼图一致，同时保留X-PAR拼图细节的拼图融合结果。最终拼图融合结果的时间分辨率与X-PAR观测(1.5 min)相同，水平空间分辨率为50 m，包含Z_H、Z_DR及K_DP三组物理量，主要分为以下几个步骤实现：

(1) 利用研究区内地面雨滴谱观测资料，模拟计算降水粒子的雷达物理量，计算得到不同波段雷达物理量拟合关系。

(2) 对研究区建立两组拼图网格，一组精细网格(水平分辨率为50 m)用于X-PAR的插值拼图，一组粗糙网格(水平分辨率为500 m)用于S-POL的插值拼图。根据不同波长和不同参量的数据质量问题，构建数据质量评价参数作为数据质量权重，结合雷达探测点与目标网格点的距离权重，将同一波段的多部雷达观测数据同时向拼图网格进行加权插值，分别得到各雷达物理量500 m分辨率的S-POL拼图结果和50 m分辨率的X-PAR拼图结果。

(3) 首先根据计算得到的不同波段雷达物理量拟合关系，将X-PAR拼图物理量转换到S波段以便于与S-POL拼图的融合；针对两种雷达时间分辨率不同的问题，对S-POL拼图进行线性外推实现与X-PAR拼图的时间匹配；然后对两波长雷达拼图的匹配点进行差值即可得到一组粗糙偏差(水平分辨率为500 m)，将该粗糙偏差进行距离指数权重插值得到精细偏差(水平分辨率为50 m)，使用精细偏差对X-PAR拼图进行订正，得到的即是S-POL拼图与X-PAR拼图的初步融合结果。

以上为多雷达拼图方法概述，主要流程如图 3所示，具体细节在2.1~2.3节进行详述。

2.1 不同波长雷达物理量拟合关系计算

选用OTT Parsivel激光雨滴谱仪(Wu and Liu, 2017)观测数据模拟计算雷达物理量，计算得到不同波长雷达物理量拟合关系。该雨滴谱仪布设在广东省惠州市龙门观测站，观测时间分辨率为1 min，采样面积为54 cm²，可测量的粒子尺度和下落末速度分为32个非等距区间，分别为0~25 mm，0~20.8 m·s^-1。

对雨滴谱资料的质量控制包括以下几个方面，首先将降水强度小于0.5 mm·h^-1或单次采样内粒子总数小于50的数据作为噪声剔除(Tokay et al, 2013)，然后利用粒子直径计算得到的粒子下落末速度理论值(Atlas et al, 1973)，将该下落末速度的±60%作为阈值，剔除严重偏离理论值的数据(Jaffrain and Berne, 2011)。

利用扩展边界条件法，可计算得到雨滴的复数散射函数(Barber and Yeh, 1975)。其中，粒子轴线假设为铅直方向，雨滴椭率利用(Pruppacher and Beard, 1970)给出的关系式进行计算。进一步，即可利用雨滴谱数据模拟得到雷达物理量，计算不同波长雷达参量的拟合关系。

对于单个降水粒子，雷达电磁波的后向散射矩阵可以表示为：

$ \left[\begin{array}{l} E_{\mathrm{H}}^{\mathrm{r}} \\ E_{\mathrm{V}}^{\mathrm{r}} \end{array}\right]=\frac{\mathrm{e}^{-\mathrm{j} k r}}{r}\left[\begin{array}{ll} S_{\mathrm{HH}} & S_{\mathrm{HV}} \\ S_{\mathrm{VH}} & S_{\mathrm{VV}} \end{array}\right]\left[\begin{array}{c} E_{\mathrm{H}}^{\mathrm{i}} \\ E_{\mathrm{V}}^{\mathrm{i}} \end{array}\right] $

(1)

式中：E_Hⁱ、E_Vⁱ为水平和垂直偏振波的入射电场强度，E_H^r、E_V^r为水平和垂直粒子散射的电场强度，j为虚数单位，k为波数，r为传播距离；S_HH、S_HV、S_VH、S_VV为散射矩阵的四个复数散射函数。分别用0和180表示前向散射和后向散射。对于雨滴谱分布为N(D)的雨区，雷达参量可利用以下公式(Bringi and Chandrasekar, 2001)计算：

$ \begin{gathered} Z_{\mathrm{H}, \mathrm{~V}}= \\ \frac{\lambda^4}{\pi^5|U|^2} \int_0^{D_{\max }} 4 \pi\left|S_{\mathrm{HH}, \mathrm{VV}}(180, D)\right|^2 N(D) \mathrm{d} D \end{gathered} $

(2)

$ Z_{\mathrm{DR}}=10.0 \lg \left(\frac{Z_{\mathrm{H}}}{Z_{\mathrm{V}}}\right) $

(3)

$ \begin{gathered} K_{\mathrm{DP}}= \\ \frac{180 \lambda}{\pi} \int_0^{D_{\max }} \operatorname{Re}\left[S_{\mathrm{HH}}(0, D)-S_{\mathrm{VV}}(0, D)\right] N(D) \mathrm{d} D \end{gathered} $

(4)

式中：$U=\frac{\varepsilon-1}{\varepsilon+2}$，ε为复相对介电常数, λ为雷达电磁波长。

2.2 单波长雷达网拼图方法

针对不同波长雷达探测特性，首先对单一波长雷达网进行拼图。其中X-PAR拼图的水平空间分辨率为50 m，S-POL拼图的水平空间分辨率为500 m，两组拼图的垂直分辨率均设置为200 m。拼图的雷达物理量包括Z_H、Z_DR、K_DP，其中Z_H在拼图时的单位由dBz转换为mm⁶·m^-3后进行插值(Warren and Protat, 2019)。

为便于计算，首先建立目标经纬度网格坐标，并计算各网格点对应的雷达观测极坐标，实现坐标转换(肖艳姣, 2007)。在雷达极坐标下对该网格点进行水平方向(方位角、距离方向)的最近邻匹配，提取每部雷达在目标网格点上方和下方相邻两个仰角的观测，然后根据数据质量和雷达探测数据与拼图网格点的距离设置插值权重，将多部雷达观测数据同时进行加权插值，获得目标网格点的单波长雷达网插值拼图结果。

在雷达仰角20°以内，可以近似将仰角方向作为垂直方向考虑。对目标网格点P(α_p, β_p, h_p)，组网雷达在该点的极坐标为(e_k, a_k, r_k)，其中k代表雷达编号。如图 4所示，以两部雷达为例，对于网格点P，可以获得第一部雷达在P点相邻上下两个仰角对应的极坐标点(e₁₂, a₁, r₁)和(e₁₁, a₁, r₁)以及第二部雷达在P点相邻上下两个仰角对应的极坐标点(e₂₂, a₂, r₂)和(e₂₁, a₂, r₂)，其中方位角a和距离r为网格点对应雷达极坐标的最近邻匹配点。对于两部雷达在P点上下观测的四个点，根据雷达探测能力，计算观测点的数据质量评价系数w_kq，进一步根据观测数据的质量评价系数和垂直方向的距离计算插值权重，对各观测点进行加权平均得到目标点P的插值结果f^a(e, a, r)。

$ \begin{gathered}f^a(e, a, r)= \\ \frac{\sum\limits_{k=1}^n\left[w_{k 1} f_k^a\left(e_{k 1}, a_k, r_k\right)+w_{k 2} f_k^a\left(e_{k 2}, a_k, r_k\right)\right]}{\sum\limits_{k=1}^n\left(w_{k 1}+w_{k 2}\right)}\end{gathered} $

(5)

$ w_{k 1}=w_{k \mathrm{q} 1}^2 w_{k \mathrm{d} 1} $

(6)

$ w_{k 2}=w_{k \mathrm{q} 2}^2 w_{k \mathrm{d} 2} $

(7)

$ w_{k \mathrm{~d} 1}=\exp \left[-\frac{r_k^2\left(e_k-e_{k 1}\right)^2}{R_{v 0}^2}\right] $

(8)

$ w_{k \mathrm{~d} 2}=\exp \left[-\frac{r_k^2\left(e_{k 2}-e_k\right)^2}{R_{v 0}^2}\right] $

(9)

式中：n为能够覆盖拼图目标点的雷达数量，$f_k^a\left(e_{k 1}\right.$，$\left.a_k, r_k\right)$和$f_k^a\left(e_{k 2}, a_k, r_k\right)$ 分别为雷达k在目标点下方和上方的观测匹配点。w_k1和w_k2为相应的插值权重，由数据质量系数和垂直距离权重两部分构成，w_kq1和w_kq2为相应观测数据的质量系数，w_kd1和w_kd2为雷达观测点向目标网格点插值的垂直距离指数权重，r_k为目标点距离雷达的径向距离，(e_k-e_k1)与(e_k2-e_k)在此处的单位为rad，R_v0设定为500 m。

对于观测数据的质量评价系数w_kq，需要根据雷达探测能力和观测物理量的特性进行设置，主要考虑以下几项评价指标。

$ w_{\mathrm{r}}=\exp \left(-\frac{r^2}{R_{\mathrm{w}}^2}\right) $

(10)

$ w_{\mathrm{n}}=\frac{1}{2 / \mathrm{SNR}+1} $

(11)

$ w_{\mathrm{a}}=\exp \left(-0.69 \frac{\varphi_{\mathrm{DP}}^2}{\varphi_{\mathrm{DPT}}^2}\right) $

(12)

$ w_{\mathrm{d}}=\exp \left(-\frac{d_{\mathrm{v}}^2}{R_{\mathrm{v} 0}^2}\right) $

(13)

式(10)~式(13)中：w_r为雷达距离权重，主要考虑随着径向距离增加，波束展宽、空间分辨率下降、充塞系数降低等因素造成的数据质量问题，是雷达组网拼图中各物理量都要考虑的重要参数。本研究的距离参数R_w对S-POL设置为300 km，对X-PAR设置为30 km。

w_n为信噪比权重，对S-POL而言，降水区的信噪比普遍较高，在剔除非降水回波后，能够获得数据质量较好的观测结果。但对X-PAR而言，由于雷达在20 km处最小可测回波已经达到10 dBz以上，加上衰减的影响，强降水区的信噪比往往偏低，该权重主要抑制低信噪比的雷达观测数据的影响。

w_a为衰减权重，对于S-POL而言，忽略降水的衰减影响。但X-PAR的Z_H和Z_DR在降水区受衰减影响很大，衰减订正后依然存在误差，该权重用来抑制严重衰减后的雷达观测数据的影响，由于雨区衰减率和差分衰减率均与K_DP存在准线性关系，利用前向传播相位差φ_DP反映路径积分衰减，参考吴翀等(2021)的研究，参数φ_DPT设定为80°。

w_d为观测点距离权重，对于S-POL而言，雷达径向远端相邻两个仰角的波束垂直距离很大，在进行垂直内插时可能会引入较大误差。为此，引入观测点距离权重，用来进一步抑制与网格点距离过远的观测数据的影响。d_v为雷达探测点与目标网格点的距离，距离参数R_v0设定为500 m。

对于不同波长雷达的不同观测物理量，质量系数评估方法有所不同。对于S-POL，主要数据质量问题为雷达波束展宽造成的径向远端数据质量下降，雷达观测点与网格目标点间距过大造成的匹配误差，以及低信噪比造成的观测数据质量较差。对于X-PAR，主要数据质量问题为雷达径向远端的波束展宽、雨区信号衰减及灵敏度较低。本研究对Z_H进行了严格的衰减订正，但由于X-PAR的基数据只存储了基于Φ_DP和固定参数的Z_DR衰减订正结果，难以做进一步的衰减订正，所以提高Z_DR插值中衰减所占的权重。由于X-PAR雷达灵敏度较低，经过雨区衰减后，信噪比偏低的数据权重也需要进行抑制。

最后，考虑波束遮挡问题引入w_o(默认取w_o=1)，对各雷达观测中受到严重遮挡、地物、避雷针影响的部分进行剔除(w_o=0)，对有遮挡但不严重的部分降低该数据的权重(w_o=0.1)，优先应用其他雷达观测结果进行插值拼图。w_o的判定主要基于地形波束遮挡率(张亚萍等, 2002)的计算结果，对地形波束遮挡率大于30%的部分取w_o=0.1，对地形波束遮挡率大于50%的部分取w_o=0；针对建筑物或地物影响，对雷达历史探测信息进行统计得到反射率因子均值，对与相邻方位反射率因子差异达到1 dB的库取w_o=0.1，差异达到3 dB的库取w_o=0。对避雷针所在方位角，取w_o=0。对不同波段雷达各物理量的质量评价方法如表 2所示。

这样，即可计算得到每部雷达在目标网格点上方和下方两个观测点的数据质量评价系数，结合雷达观测点和目标点的距离计算得到插值权重，进行加权平均后即可获得该目标网格点的单一波长雷达网观测的插值结果。

2.3 S-POL与X-PAR雷达拼图融合方法

在完成了单一波长雷达组网拼图后，能够获得时间分辨率为1.5 min、水平空间分辨率为50 m的X-PAR拼图和时间分辨率为6 min、水平空间分辨率为500 m的S-POL拼图。然后以S-POL拼图作为相对粗糙但准确可靠的基准场，将S-POL拼图与X-PAR拼图(根据拟合关系转换为S波段)中匹配的网格点的差值作为X-PAR拼图的粗糙偏差，对粗糙偏差进行插值获得相对平滑的高分辨率X-PAR拼图的精细偏差，然后对X-PAR拼图进行订正，即可得到初步融合结果。

由于两组拼图的时间分辨率存在差异，在计算粗糙偏差时，首先需要对S-POL的拼图结果进行外推，与X-PAR的拼图结果进行匹配。基于实际观测对比，当降水回波水平方向梯度较大时，两组拼图的空间错配可能造成反射率因子±30 dBz以上的偏差，并表现出高低错落的条状高值分布。考虑到两组拼图的时间差非常小(5 min以内)，并且粗糙偏差向精细网格的插值处理对匹配误差具有抑制作用，在实际计算中，将雷达回波短时间内的空间变化简化为不同高度层回波在水平方向的刚体运动。

定义平均四次误差MQE(采用四次方运算，使MQE对偏差的高值更为敏感)，表征S-POL拼图结果与X-PAR拼图结果的一致性；以500 m为步长设置移动矢量的东西向分量和南北向分量，对比8 km范围内的所有移动矢量下两组拼图的MQE，取最小MQE对应的移动矢量作为当前高度层S-POL拼图的外推矢量。需要说明的是，受制于探测能力，S-POL拼图在对流层底层，X-PAR拼图在对流层中上层均可能出现较大范围的缺测。外推矢量的计算要求具有充足的观测样本，当S-POL拼图或X-PAR拼图在某一层的观测样本量偏低时(本研究将这一样本量阈值定为S-POL拼图在2 km高度层样本量的1/2)，则利用相邻层次的外推矢量作为该层次的外推矢量。

$ \mathrm{MQE}=\frac{\sum\limits_{i=1}^n\left(Z_{\mathrm{S} i}-Z_{\mathrm{X} i}\right)^4}{n} $

(14)

式中：Z_Si(单位：dBz)为S-POL拼图在i点坐标的反射率因子，Z_Xi(单位：dBz)为X-PAR拼图(经过波长转换为S波段后)在i点坐标的反射率因子。

在对S-POL拼图逐层进行矢量外推后，即可根据经纬度网格坐标索引实现S-POL拼图和X-PAR拼图的空间匹配。将两组拼图在匹配点的差值作为X-PAR拼图的粗糙偏差。进一步计算距离指数权重，将粗糙偏差插值到X-PAR拼图的网格坐标点，获得X-PAR拼图的精细偏差分布。

$ {D_{{\rm{H}}j}} = \frac{{\sum\limits_{{d_{{{\rm{h}}_{ij}}}} = 0,{d_{{{\rm{v}}_{ij}}}} = 0}^{{d_{{{\rm{h}}_{ij}}}} = 2,{d_{{{\rm{v}}_{ij}}}} = 0.4} {{w_{ij}}{D_{{\rm{L}}i}}} }}{{\sum\limits_{{d_{{{\rm{h}}_{ij}}}} = 0,{d_{{{\rm{v}}_{ij}}}} = 0}^{{d_{{{\rm{h}}_{ij}}}} = 2,{d_{{{\rm{v}}_{ij}}}} = 0.4} {{w_{ij}}} }} $

(15)

$ D_{\mathrm{L} i}=M_{\mathrm{S} i}-M_{\mathrm{X} i} $

(16)

$ w_{i j}=\exp \left(\frac{-d_{i j}^2}{\text { roi }^2}\right) $

(17)

$ d_{i j}^2=d_{\mathrm{h}_{i j}}^2+\left(z_{\mathrm{f}} \times d_{\mathrm{v}_{i j}}\right)^2 $

(18)

式中: D_H为待求的50 m分辨率精细偏差，D_L为500 m分辨率的粗糙偏差，w_ij为距离指数权重，M_Si为S-POL拼图在i点的偏振物理量，M_Xi为X-PAR拼图在i点的偏振物理量，d_ij代表粗糙偏差网格点与精细偏差网格点的距离参数，其中d_hij代表水平距离，d_vij代表垂直距离，设置垂直插值系数z_f=5，用来降低垂直方向插值的权重。本研究中，对目标网格点，选择水平2 km垂直400 m范围内的粗糙偏差进行插值，roi取2 km。

计算得到X-PAR拼图结果的精细偏差网格后，即可得到X-PAR拼图订正后的结果，即S-POL与X-PAR的初步融合结果M_XC。

$ M_{\mathrm{XC}}=M_{\mathrm{X}}+D_{\mathrm{H}} $

(19)

受波束充塞不足、数据样本量偏少等原因的影响，初步融合的结果在回波边缘或底层的连续性较差。因此，针对目标网格点可能存在的以下四种情况，按照相应的方法确定该点最终的雷达拼图融合结果。

(1) 如果该点存在X-PAR拼图数据和精细偏差数据，并且该点计算精细偏差时使用的粗糙偏差数据量超过指定阈值(本研究设定为200个)。直接采用M_XC作为该点的最终拼图融合结果。在实际拼图融合中，绝大多数网格点均满足此类情况。

(2) 如果该点X-PAR拼图数据缺失，但存在S-POL拼图数据。将相应位置的S-POL拼图进行最近邻插值作为拼图融合结果。

(3) 如果该点存在X-PAR拼图数据和精细偏差数据，但该点计算精细偏差时使用的粗糙偏差数据量小于指定阈值。这种情况普遍出现在X-PAR拼图的边缘，该区域周围可以匹配的S-POL拼图结果相对较少。当粗糙偏差样本量不足时，计算的精细偏差可靠性较差。对这部分网格点，根据粗糙偏差的样本量设置权重，对M_XC与S-POL拼图结果进行加权平均，得到的结果作为该点的拼图融合结果。

(4) 如果该点处于大气低层S-POL探测盲区，S-POL拼图数据缺失，无法获得精细偏差数据，但存在X-PAR拼图数据。使用该点坐标底层到2 km高度内的精细偏差数据的平均值作为该点的精细偏差，对X-PAR拼图结果进行订正，作为该点最终拼图融合结果。四种情况的计算方案如式(20)所示，最终得到的M即为S-POL与X-PAR雷达网拼图融合结果。

$ M_j= \begin{cases}M_{\mathrm{XC} j} &\;\;\; N_j \geqslant N_{\mathrm{th}} \\ M_{\mathrm{S} j} &\;\;\; M_{\mathrm{X} j} \text { 缺失 } \\ w_{\mathrm{X} j} M_{\mathrm{XC} j}+w_{\mathrm{S} j} M_{\mathrm{S} j} & \;\;\;N_j <N_{\mathrm{th}} \\ M_{\mathrm{X} j}+D_{\mathrm{H} j \mathrm{~s}} \quad& M_{\mathrm{S} j} \text { 缺失, Alt }<1.5 \mathrm{~km}\end{cases} $

(20)

$ w_{x_j}=1 /\left\{1+\exp \left[-2\left(N_j / 40-4\right)\right]\right\} $

(21)

$ w_{\mathrm{S} j}=1-w_{\mathrm{X} j} $

(22)

$ {D_{{\rm{H}}j{\rm{s}}}} = \frac{{\sum\limits_{{d_{{{\rm{h}}_{jk}}}} = 0,{\rm{Alt}} = 0}^{{d_{{{\rm{h}}_{jk}}}} = 0,{\rm{Alt}} = 2} {{D_{{\rm{H}}k}}} }}{{{{{\mathop{\rm Num}\nolimits} }_k}}} $

(23)

式中：M_Sj为S-POL拼图在精细网格点j处的最近邻插值，N_j为计算点j的精细偏差时能够应用到的粗糙偏差的样本数，N_th作为阈值在本研究中取为200；当存在初步融合结果M_XCj，但插值样本量不足N_th时，则对j点的拼图结果进行加权计算，w_Xj为初步融合结果M_XCj的权重，w_Sj为M_Sj的权重；当该点为S-POL低层盲区且存在X-PAR拼图结果时，D_Hjs代表X-PAR拼图在j点的替代精细偏差，它是该位置从底层到2 km高度层之间的精细偏差(D_Hk)的平均值，Num_k代表j点位置0~2 km高度层精细偏差的样本量。

3 结果 3.1 不同波长雷达参量的拟合结果

由于降水粒子群对不同波长雷达的散射特性存在差异，在进行X-PAR与S-POL的对比和不同波长雷达拼图融合时，需要对雷达观测物理量进行波长转换。利用2019—2022年雨滴谱观测资料分别计算X波段和S波段的雷达物理量，并对两波长的Z_H、K_DP进行指数拟合，对Z_DR进行多项式拟合(图 5)。从Z_H的分布来看，弱降水粒子尺度较小，对X波段和S波段电磁波的散射均为瑞利散射，二者呈现很好的线性关系，随着粒子尺度增加，瑞利散射的假设不再对两个波长同时有效，二者的Z_H分布开始变得相对离散。两波长的Z_DR在降水粒子较小时同样表现出一致的线性关系，但随着粒子尺度增大，两波长的水平偏振波与垂直偏振波的散射差异对Z_DR的值产生了较大影响。K_DP的值与雷达波长负相关，两波长K_DP呈现准线性关系。两波长雷达物理量拟合结果如下。

$ Z_{\mathrm{HS}} =1.194 Z_{\mathrm{HX}}^{0.948} $

(24)

$ K_{\mathrm{DPS}} =0.2733 K_{\mathrm{DPX}}^{1.041} $

(25)

$ Z_{\mathrm{DRS}}=\frac{p_1 Z_{\mathrm{DRX}}^3+p_2 Z_{\mathrm{DRX}}^2+p_3 Z_{\mathrm{DRX}}+p_4}{Z_{\mathrm{DRX}}^2+q_1 Z_{\mathrm{DRX}}+q_2} $

(26)

式中：p₁ =1.125，p₂=-5.976，p₃=9.997，p₄=-0.1347，q₁ =-5.385，q₂ =9.834

3.2 单一波长雷达拼图结果

强对流天气系统是短时临近预警业务中天气雷达观测的主要目标，回波细节结构对于研判对流系统特征和短时临近预警具有重要意义。本研究以2022年5月11日07:12(世界时，下同)广东一次强对流天气为例，对提出的拼图方法进行检验。

单部S-POL的观测可能存在径向远端分辨率下降，回波细节特征缺失等问题，这些均不利于强对流的观测分析。通过多雷达组网观测，增加同一气象目标的观测样本，能够有效改善上述问题。如图 6所示，各雷达对该强对流系统的观测具有较好的一致性，但受波束展宽的影响，Z_H和Z_DR表现出不同程度的方位分辨率偏低的特征。多雷达拼图结果整体保留了各雷达观测的共性特征，同时使Z9200观测中的波束遮挡问题(图 6a)、Z9662观测中K_DP高值中心东北侧的噪声(图 6f)等问题得到抑制，获得了空间相对连续，共性结构特征得以保留的强对流观测结果。

对于X-PAR，由于衰减严重和灵敏度较差，单部X-PAR对强降水后侧存在严重缺测(图 7a~7l)。利用组网雷达拼图，能够有效实现对降水系统的完整观测，同时较好保留单部雷达观测回波的细节特征。对比可见，各雷达计算得到的K_DP虽然整体分布具有相似性，但强度和空间结构均存在一定差异，主要原因一方面是径向远端波束展宽，导致波束非均匀问题及回波结构沿方位向平滑，另一方面是对Φ_DP噪声信息进行抑制的过程中，低通滤波可能会引入误差。多雷达拼图获得了整体较为连续且具有细节特征的K_DP结构，并保留了各雷达观测的共性特征。对比各雷达观测的Z_DR，可见径向远端产生了非常明显的高估，这可能与过度衰减订正有关，在多雷达拼图过程中，由于严格约束了衰减严重部分的插值权重，在Z_DR的拼图结果中，雷达径向远端Z_DR偏高的问题得到了有效抑制。整体来看，各物理量获得了较为完整精细的拼图结果。

3.3 偏差计算

以3.2节所展示的对流降水观测为例，根据2.3节所述的外推矢量计算方法，计算得到S-POL逐层外推矢量，分别计算S-POL拼图外推前和外推后与X-PAR拼图的差值，结果如图 8所示。强对流降水回波水平梯度较大，回波移动导致的错配会在两拼图的差值上表现出条状正负极大值错落分布的形态。对比外推前和外推后的拼图差值，可见外推前的Z_H差值和K_DP差值均表现为明显的正、负极大值错落分布，经过外推后，这种条状高值错落分布的现象得到明显缓解，粗糙偏差的空间分布更加连续。由于S-POL拼图是多部雷达6 min内多个仰角连续观测数据的融合结果，拼图并不代表某一固定时间点，再加上雷达观测空间分辨率的差异，外推处理后，两组雷达拼图的偏差仍然存在空间较为离散的极大值和极小值分布。但这些极端高低值会在加权平均过程中得到抑制，最终得到相对连续平滑的精细偏差计算结果。

进一步利用2.3节所述方法，对X-PAR的粗糙偏差进行插值，得到与X-PAR拼图网格(50 m×50 m×200 m)匹配的精细偏差(图 9)。其中Z_H和K_DP的误差高值中心均位于该对流系统中心部分，这是因为该对流中心位置距离各部雷达的距离均较远，各雷达在对对流中心部位的探测时，都受到了严重衰减、波束展宽、信噪比较低的影响，导致该部位数据质量较差，在本个例中产生了一定的低估。而X-PAR的Z_DR拼图偏差整体呈现出南侧偏低，北侧偏高的空间分布特征，经对比发现，该现象主要源于回波南侧雷达观测与北侧雷达观测的系统偏差订正不足，这可能与Z_DR的衰减订正与标定误差有关。

总体而言，插值得到的精细偏差分布与图 8中粗糙偏差的空间分布特征基本一致，同时图 8中那些由于空间错配造成的极端噪声得到了很好的抑制。进一步，根据2.3节所述方法，利用精细偏差对X-PAR拼图结果进行订正，即可将S-POL观测作为背景值与X-PAR拼图的精细结构进行融合。

需要说明的是，精细偏差分布图中，回波边缘处(如西侧、北侧)存在偏差高值，这主要是因为在雷达回波边缘，网格点周围的粗糙偏差样本量偏低，不利于噪声的抑制，按照2.3节中所述方法，对于这些插值样本量偏低的网格点，会采用与S-POL加权融合的方式进行赋值，避免质量较差的偏差估计对最终拼图融合结果的污染。

3.4 融合效果

S-POL与X-PAR雷达拼图融合的目的在于充分利用S-POL拼图数据质量可靠的特性和X-PAR拼图分辨率较高的优势，得到准确精细的三维雷达观测场。为此，雷达拼图融合结果的回波结构和强度应与S-POL拼图基本一致，同时保留X-PAR拼图中的细节特征。S-POL与X-PAR雷达拼图融合结果如图 10所示，可以看出，S-POL拼图与X-PAR拼图(已根据2.1节中拟合关系转换到S波段)的空间结构特征基本一致，但受到雨区严重衰减和波束展宽对数据质量的不利影响，X-PAR各物理量拼图在高值中心部分存在低估，融合后的各物理量的强度得到较好订正，与S-POL拼图基本一致，同时保留了X-PAR拼图中的细节特征。

沿图 10中虚线进行剖面(结果见图 11)，对比拼图融合结果，可见由于S-POL雷达间距较大，在S-POL拼图中，大气底层存在探测盲区，1 km以下存在较大范围缺测，X-PAR雷达布设更加密集，大气底层探测能力更强，由拼图融合结果可见，S-POL拼图中底层缺测部分得到了较好补盲。另外，由于X-PAR灵敏度较低，且探测范围较小，所以存在弱回波缺测、静锥区缺测的问题，拼图融合后，对X-PAR的这部分缺测也得到了较好补盲。融合后的拼图实现了X-PAR和S-POL雷达的相互补盲，获得了更加完整的三维观测场信息。

从强度信息来看，剖面结构与CAPPI表现出相似的结果，以东侧的回波悬垂为例，X-PAR拼图的Z_H和Z_DR高值中心存在强度偏低的问题，融合后的Z_H、Z_DR一方面保留了X-PAR拼图中的悬垂细节特征，同时高值中心的强度也与S-POL拼图相一致。相比于S-POL拼图的K_DP剖面特征，X-PAR拼图中K_DP西侧悬垂高值中心范围偏大，东侧悬垂高值中心范围偏小，这种偏差也在融合后得到校正。

沿图 10中虚线2 km高度层的三组拼图结果进行分析，对比融合前后雷达参量变化，结果如图 12所示，整体来看，拼图融合结果(蓝色曲线)的各物理量，既包含了X-PAR拼图(灰色曲线)的细节变化特征，同时整体强度分布也与S-POL拼图(红色曲线)比较一致。其中X-PAR的反射率因子数据质量相对较好，X-PAR拼图与S-POL拼图的结果在整体强度分布上表现基本一致，部分X-PAR拼图中强度偏低的区域，经过融合后得到明显改善。对东部强回波区，Z_DR拼图和K_DP拼图中，X-PAR均表现出明显低估，融合后结果的强度和S-POL基本一致。拼图融合整体实现了以S-POL拼图为背景值，同时保留X-PAR拼图细节特征的目的。

4 结论与讨论

针对X-PAR和S-POL的探测能力，本研究分别提出针对组网X-PAR和组网S-POL的雷达拼图方法，实现对单波长雷达网的组网拼图，并进一步对S-POL与X-PAR雷达网拼图进行融合，发挥多波段雷达组网观测优势，形成相对准确精细的三维雷达观测场。得到以下结论：

(1) S-POL的特性在于数据质量较高，观测较为准确，但S-POL布设间距相对较远，在拼图过程中，为抑制雷达远端空间分辨率下降造成的插值误差及波束展宽导致的数据质量下降，本研究综合考虑雷达观测点与目标网格点间距、雷达径向距离、波束遮挡和信噪比四个方面的影响设置插值权重，并相应调整垂直方向线性内插为指数权重插值，将多部S-POL观测进行同时插值实现S-POL组网拼图。结果显示该拼图方法能够保留各雷达观测的共性特征，得到细节特征明显且空间连续的拼图结果。

(2) X-PAR的特性在于空间分辨率较高，但观测范围较小，波束宽度较宽，雷达灵敏度较低且受衰减影响较大，在拼图中，数据质量的问题比S-POL更为复杂。因此在多雷达拼图过程中，需要对不同物理量计算相应数据质量评价系数，对数据质量差的观测进行抑制，然后结合观测点与目标点的距离，设置插值权重进行拼图。结果显示该方法能够使多雷达拼图更多应用到数据质量可靠的X-PAR观测资料，避免质量较差的观测数据对拼图的污染，克服单部X-PAR对强降水区观测不完整不连续的问题，获得了相对连续、完整且细节特征得以保留的X-PAR拼图结果。

(3) 根据两种雷达时空分辨率和观测数据质量的差异，将S-POL拼图作为分辨率较低，细节特征不足但相对准确的背景场，与X-PAR拼图进行时间匹配并计算得到粗糙偏差，对粗糙偏差进行距离指数权重插值得到细节偏差，以此对X-PAR拼图进行订正，即可得到S-POL拼图与X-PAR拼图的融合结果。结果表明，该方法能够实现对S-POL观测中的低空补盲和对X-PAR观测中的弱回波区补盲，并且整体实现了最终融合拼图强度分布特征与S-POL拼图相一致，同时保留X-PAR拼图细节结构的目的，获得了相对准确精细的三维雷达观测场。

需要说明的是，基于Warren and Protat(2019)的研究结果，本研究在进行拼图时将反射率因子的单位由dBz转换为mm⁶·m^-3，但使用mm⁶·m^-3插值的结果会比使用dBz插值结果偏大，在肖艳姣(2007)的对比研究中，使用dBz插值的结果总体上更接近观测值，反射率因子在插值拼图时的单位如何选取，仍需基于实际观测开展进一步的对比验证；本研究对X-PAR的数据质量控制及其与S-POL的雷达拼图融合方法均建立在地面雨滴谱数据模拟雷达物理量分析的基础之上，利用雨滴谱数据统计得到的雷达参量拟合关系对固态及混合相态降水则不再适用。此外，本研究提出的S-POL拼图与X-PAR拼图时空匹配方法比较简单，对降水回波短时间(5 min)内的变化采用了刚体假设进行外推，该方法忽略了单层CAPPI中各网格点可能存在的时间差及5 min内回波可能出现的涡旋、生消等复杂演变，可能会给X-PAR拼图偏差的计算引入新的误差，在实际应用中，可综合考虑对精度和运算效率的需求，利用更准确的外推匹配方法或不进行外推，结合外推匹配效果，灵活调整S-POL和X-PAR的时空匹配方式，改进雷达拼图融合效果。