引言

自20世纪80年代末以来，大气探测技术迅猛发展，卫星、雷达等非常规探测手段逐渐成熟，很大程度地促进了基于云分析的多源观测资料同化系统的发展。多普勒雷达作为云探测的重要手段之一，能够较好地反映对流层中低层云的宏微观信息(韩成鸣等，2015)。为了将多种云观测信息客观、有效地用于天气分析，改善天气预报效果，美国大气海洋管理局预报系统实验室着手设计和开发了资料同化和预报系统LAPS(Local Analysis and Prediction System)(McGinley et al, 1991)。LAPS具有其独特的资料融合优势，能够融合地面、探空、雷达、卫星、飞机报文和浮标等多种观测资料，得到高分辨率分析场。Albers (1995)和Albers等(1996)将雷达反射率用于修正三维云场，还将雷达反射率用于雨水场的反演，有效地将雷达观测到的云降水空间结构及强度用于云分析并将其引入数值模式初始场中，从而更加准确地描述中尺度对流系统。国内最早是由武汉暴雨所引进并投入业务运行(李红莉等，2009)，随后在北京市气象局准业务应用，建立了LAPS-RUC并进行检验评估(高华等，2010；李红莉，2011)。此外，刘瑞霞等(2011)研究了LAPS融合FY-2C卫星、地面观测、雷达数据后三维云场的分布，结果表明LAPS云分析中，地面观测对云底结构起主要订正作用，雷达观测对云中、低部信息起主要修正作用，而卫星云图数据对云顶分布订正效果显著。彭菊香等(2011)对华中区域LAPS中尺度分析场做了检验与评估。白永清等(2013)将LAPS分析场用于WRF模式逐时气温精细化预报释用中，改善了预报效果。

LAPS若想实现本地化，需要使用我国新一代多普勒雷达观测的基数据。LAPS可以实现多部雷达同时融合，首先将单部雷达数据插值到笛卡尔坐标系上并进行质量控制，然后将单部雷达数据通过最近邻居法进行拼图处理。单部雷达基数据大多都是VCP21体扫数据(图 1)，由于扫描方式的局限性，在19.5°仰角以上有静锥区，最低扫描仰角0.5°以下有资料空白区，高仰角之间也有间隙，还受分辨率、探测范围和地形阻挡等影响(肖艳姣，2007)，导致资料不完整，而多部雷达拼图则可以在一定程度上弥补单部雷达的不足。但是，LAPS最近邻居的拼图方法无法使用其他雷达在其探测盲区的反射率资料，需要使用将多部雷达观测资料效能最大化的拼图方法。国外学者研究雷达拼图较早，Zhang等(2004；2011)基于美国NWS(National Weather Service)接收的130部WSR-88D雷达数据开发三维雷达反射率拼图系统，拼图方法采用距离指数权重法。后来，Langston等(2007)又继续开发了基于距离和时间权重的4D(Four-Dimensional)动态雷达反射率拼图技术。国内研究三维雷达拼图方法的较少，肖艳姣等(2006)针对新一代天气雷达网资料，对比了多种单部雷达插值方法和多部雷达拼图方法，最后试验得出基于NVI插值法和距离指数权重平均拼图法能够较好地提供空间连续的三维反射率拼图。王瑾等(2011)采用多项式局部插值方法和最大值的拼图方法，将得到的三维雷达反射率用于冰雹诊断研究中，识别效果较好。周海光(2007)、王红艳等(2009)在三维雷达组网拼图系统开发方面取得了一定的成果。汪玲等(2015)、吴汪毅等运用雷达组网拼图用于人工增雨作用，取得了一些成效。

多部雷达拼图虽能在一定程度上弥补单部雷达探测的不足，但由于雷达探测仰角、地形和布网密度的限制，经过三维雷达拼图之后，通常只能对周围雷达能够探测到的对流层中层进行填补，而流层中低层(3~4 km以下)和高层仍然存在空白区域，静锥区并未完全由周围雷达观测值填补。国内很少有针对雷达反射率静锥区填充方法的研究，本文基于LAPS雷达资料融合程序，依据前人对雷达拼图方法的研究成果，改进LAPS拼图方法，并且试图通过最小二乘拟合方法来改善雷达反射率静锥区的空白区域。

1 方法设计

本文用于改进LAPS反射率拼图的方法有两种，即最大值法和距离指数权重法。静锥区填充方法是运用最小二乘法拟合，建立所选区域内实测雷达反射率与三维坐标系的多元线性回归方程，将三维雷达反射率的空间分布形态用于填充满足一定条件下的静锥区。另外，本文试验选取LAPS自带的质量控制程序，仅在最大值和距离指数权重拼图模块中取350 km为单部雷达有效观测半径。

1.1 反射率拼图方法

通过一个或多个客观分析方法把来自各个雷达的反射率场插值到统一的笛卡尔坐标系之后，需要把来自多个雷达的格点反射率场拼接起来形成三维拼图网格。在拼图网格的很多区域，特别是在对流层中高层，有来自多个雷达的资料重叠区域，在拼图网格中的每个网格单元i的反射率值可以通过下面公式得到

$ {f^m}\left(i \right) = \sum\limits_{n = 1}^N {{w_n}} f_n^a\left(i \right)/\sum\limits_{n = 1}^N {{w_n}} $

(1)

式中，f^m(i)是网格单元i的合成反射率值，$ f_{n}^{a}\left(i \right) $是在网格单元i出来自第n个雷达的分析值，w_n是分析值$ f_{n}^{a}\left(i \right) $的权重，N是在网格单元i处有分析值的总雷达数。

(a)最大值法

最大值法就是把覆盖同一网格单元的多个雷达反射率中的最大值赋给网格单元，即把分析值中的最大值权重赋为1，其他权重赋为0。

(b)距离指数权重法

指数权重函数

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

(2)

式中，R为影响半径，本文取R＝150，r为网格点到雷达的距离。

1.2 静锥区填充方法

最小二乘法是提供“观测组合”的主要工具之一，它依据对某事件的大量观测而获得“最佳”结果或“最可能”表现形式。如已知两变量为线性关系y＝ax+b，对其进行n(n＞2) 次观测而获得n对数据。若将这n对数据带入方程求解a、b的值则无确定解。最小二乘法提供了一个求解方法，其基本思想就是寻找“最接近”这n个观测点的直线(贾小勇等，2006)。

本文将建立雷达反射率f(x, y, z)与三维坐标系(x, y, z)之间的多元线性方程，试图通过最小二乘法求取多元线性回归方程的回归系数，然后将回归方程用于填充静锥区。

对于给定的数据点(x_i, y_i, z_i)i=1, 2, …, N，三维反射率与坐标系的关系可以用一个多元线性方程来表示

$ f\left({{x_i}, {y_i}, {z_i}} \right) = a{x_i} + b{y_i} + c{z_i} + d $

(3)

式中a、b、c、d为多元回归系数。运用最小二乘法的求解，即使得总误差$ Q\left({a, b, c, d} \right) = \sum\limits_{i = 1}^N {{{\left[ {{f_i} - \left({a{x_i} + b{y_i} + c{z_i} + d} \right)} \right]}^2}} $为最小。由微积分知识可知，这一问题的求解可归结为四元函数Q(a, b, c, d)的极值问题，即a、b、c和d应满足：$ \partial Q/\partial a = 0, \partial Q/\partial b = 0, \partial Q/\partial c = 0, \partial Q/\partial d = 0 $；得到拟合方程的法相方程组：

$ \left\{ \begin{array}{l} a\sum\limits_{i = 1}^N {x_i^2} + b\sum\limits_{i = 1}^N {{x_i}{y_i}} + c\sum\limits_{i = 1}^N {{x_i}{z_i}} + d\sum\limits_{i = 1}^N {{x_i} = \sum\limits_{i = 1}^N {{x_i}{f_i}} } \\ a\sum\limits_{i = 1}^N {{x_i}} {y_i} + b\sum\limits_{i = 1}^N {y_i^2} + c\sum\limits_{i = 1}^N {{y_i}} {z_i} + d\sum\limits_{i = 1}^N {{y_i}} = \sum\limits_{i = 1}^N {{y_i}{f_i}} \\ a\sum\limits_{i = 1}^N {{x_i}} {z_i} + b\sum\limits_{i = 1}^N {{y_i}{z_i}} + c\sum\limits_{i = 1}^N {z_i^2} + d\sum\limits_{i = 1}^N {{z_i}} = \sum\limits_{i = 1}^N {{z_i}{f_i}} \\ a\sum\limits_{i = 1}^N {{x_i}} + b\sum\limits_{i = 1}^N {{y_i}} + c\sum\limits_{i = 1}^N {{z_i}} + dN = \sum\limits_{i = 1}^N {{f_i}} \end{array} \right. $

(4)

若系数矩阵行列式不等于0，可解得a、b、c、d的值，进而可以得到拟合的多元线性回归方程(霍晓程等，2009)。

为了比较有效地填充静锥区，本文将采样区域设定为以雷达位置为中心的200 km范围内，仅当每个等压面上的缺测值数量小于该区域总格点数的1/3时才能填充静锥区。这样约束是为了使得每个高度上静锥区附近存在足够多的有效雷达反射率观测值，从而保证静锥区填充之后具有一定的连续性，不引入过多的噪音。

2 LAPS雷达反射率拼图试验 2.1 试验设计

采用美国大气海洋管理局(NOAA)开发的中尺度分析预报系统LAPS资料融合模块开展雷达反射率拼图试验。LAPS融合我国新一代雷达基数据，需要将其预处理为9个仰角的数据文件，然后通过插值程序将每部雷达数据插值到LAPS网格上(笛卡尔坐标)，最后将多部雷达数据进行拼图，得到LAPS网格上的三维雷达反射率。

该拼图试验区域中心点为40°N、116°E，水平分辨率为5 km，网格数为200×200，垂直分辨率为50 hPa，共43层。雷达基数据选取2012年7月21日06时(世界时)北京大兴、天津塘沽、河北石家庄和沧州的四部S波段多普勒雷达资料。该拼图试验将分别用LAPS原有的最近邻居法、改进的最大值法和距离指数权重法展开，对比这三种方法的拼图效果及优缺点。

2.2 结果分析

本文所选个例为2012年北京“7·21”特大暴雨初期的雷达回波进行拼图试验。最近邻居法、最大值法和距离指数权重法在200、500和850 hPa等压面上的拼图结果如图 2所示。从三种拼图结果来看，锋面云系呈东北—西南走向，在北京上空已经开始形成一个弓形回波，大于40 dBz的回波已经出现在对流层中层以上，对流发展较强，回波顶高已经超过了200 hPa。对流层底层除了大片对流云系，暖区一侧还出现一些层状云系。试验所选取的四部雷达基本上能够反映锋面云系的空间结构特征，三种雷达拼图方法结果大体相近，但是也存在很大的不同。

最近邻居法(图 2a~2c)在200、500和850 hPa都有观测空白区域，尤其在对流层高层和底层较为明显，出现了多个雷达静锥区空白区域，而对流层中层情况相对较好。最大值法(图 2d~2f)在三个等压面上的观测空白区域有明显的减少，大于40 dBz的回波区域面积较最近邻居法更大，将北京和石家庄雷达之间的一个主要回波区域显示出来，而且500 hPa(图 2e)上几乎看不出有观测空白区域，较为连续；200 hPa(图 2d)和850 hPa(图 2f)上有一定的改善，但是仍然存在静锥区空白区域。距离指数权重法(图 2g~2i)在空白区域填充上与最大值法相近，回波空间结构在对流层底层差异较大；就回波分布情况来看，与最近邻居法更为相近，这主要与拼图算法有关。另外，最近邻居法(图 2a~2c)中出现在北京雷达较远处的一圈杂波，在最大值法和距离指数权重法中予以剔除。

最近邻居法由于只将距离格点最近雷达的反射率赋予格点，同一个格点处没有考虑其他雷达的观测，没有最大程度地利用多部雷达观测资料，导致空白区域较多。最大值法和距离指数权重法都考虑了周围雷达的影响，不同的是，最大值法选取覆盖格点上所有雷达反射率中最大值，保留了原有的观测值，但是受雷达探测精度的影响较大；而距离指数权重法将雷达相对格点的距离作为权重，从而可以将所有覆盖格点的反射率都利用起来，从而消除单部雷达的观测误差，但是在静锥区就会出现由于最近雷达缺测值占的权重最大而导致分析值偏小的现象。

约40°N处三种拼图方法得到的东西向垂直剖面如图 3所示。图 3a是最近邻居法拼图的结果，从图上可以明显地看出，由于扫描仰角的影响而出现的静锥区和高仰角之间的资料空白区，静锥区东西两侧有两支较强回波，回波分布不连续。从图 3b和3c上可以看出，最大值法和距离权重法很大程度地改善了回波分布的连续性，高仰角之间的资料空白区几乎消除了，静锥区也填充了大部分空白区域，但是在对流层高层还是有一些空白区域因为多部雷达探测距离和仰角的限制，无法填充。从回波强度分布来看，最大值法(图 3b)较最近邻居法(图 3a)和距离指数权重法(图 3c)大于30 dBz的区域更多，而最近邻居法与距离指数权重法结果较为相近。

从整体来看，最大值法和距离权重法要优于最近邻居法，能够利用多部雷达资料一定程度上改善观测空白区域，最大程度地发挥多部雷达的效能。最大值法和距离指数权重法各有优缺点，可以根据不同的研究对象和观测条件进行选取。

3 雷达反射率静锥区填充试验 3.1 试验设计

基于经最大值法和距离指数权重法改进后的LAPS雷达反射率拼图结果，静锥区填充试验将建立雷达反射率与三维坐标系之间的多元线性方程，通过最小二乘法求取回归系数。本试验将会对每部雷达建立一个方程，在满足方法设计中的特定条件下，对每个雷达产生的静锥区进行填充。

3.2 结果分析

两种拼图方法经过填充之后在200、500和850 hPa上的结果如图 4所示。在200 hPa图上(图 4a和4d)，由于填充算法的设计，没有对其进行填充。在500 hPa图上(图 4b和4e)，仅对北京大兴雷达的静锥区进行填充，填充区域与周围回波的梯度较小，拟合效果较好，但是在静锥区东南侧由于没有观测值，填充后出现了一个圆形不连续边界。在850 hPa图上(图 4c和4f)，四部雷达产生的静锥区都得到了填充，其中北京大兴雷达处于对流云系区域，周围回波观测值较多，填充之后较为连续，看不出有很明显的填充痕迹；其余三部雷达都处于层状云系区域，周围回波观测值相对较少，回波强度较弱，填充之后存在一些不连续。最大值法(图 4a~4c)和距离指数权重法(图 4d~4f)填充效果大体相似，最大值法在对流层中层填充之后表现的更为连续，而距离指数权重法则在对流层低层填充之后效果更好。

为体现本填充试验对对流层低层云系的填充效果，本文选取约38°N的垂直剖面(图 5)作进一步对比分析。从图 5a和5c可以看出，此剖面上均为小于30 dBz的回波，没有很强的对流云回波，114°E以东大多为对流层低层的层状云回波，石家庄和沧州两部雷达产生了两个反射率静锥区分别位于114.7°E和117°E附近。从图 5b和5d可以看出，经过填充之后，两个雷达反射率静锥区有了一定的改善，填充反射率都小于20 dBz，这与静锥区周围回波都属于层状云弱降水回波，符合雷达反射率的空间分布，但是由于受114°E以西回波垂直分布的影响，拟合的结果相对于静锥区周围回波强度稍微偏大，而且梯度不是很明显。最大值法(图 5b)和距离指数权重法(图 5d)在静锥区填充的回波强度上稍有不同，这与两者本身的回波空间结构有很大的关系。

从总体上看，这次填充试验取得了一定的效果。由于算法的设计，只能在静锥区周围有足够观测值得情况下才能够进行填充，这符合雷达回波空间分布的连续性。试验结果表明，该填充方法对静锥区中对流层低层云系的填充效果较好，这也是这次研究的主要目的。因为对流层低层由于雷达仰角的限制，多部雷达数据重叠区域少，仅靠拼图算法的改进不能够很好地填充静锥区的资料空白区域，只能通过根据周围观测资料的空间分布形态来拟合静锥区的雷达反射率值。

4 结论与讨论

多普勒雷达作为云探测的重要手段之一，能够很好地反映云的三维结构。LAPS融合雷达资料，分析三维云场和雨水场，能够有效地将云信息引入数值预报初始场中，从而改善数值预报质量。本文基于LAPS雷达资料融合程序，利用最大值和距离指数权重的拼图方法改进LAPS原有的最近邻居法拼图模块。研究结果表明，相比最近邻居法而言，最大值法和距离指数权重法能够更多地利用周围雷达的观测信息，明显改善了高仰角之间的资料空白，一定程度低改善了静锥区空白，尤其在对流层中层。最大值法和距离指数权重法由于算法设计的不同，拼图结果有一定的差异，最大值法能够较好地保留观测的原始信息，但受雷达观测质量的影响大。距离指数权重法拼图结果相对较为连续，将所有覆盖到的雷达反射率都利用起来，但在雷达静锥区的格点分析值会因为静锥区缺测值的权重最大而偏小。两种方法各有优缺点，可以针对不同的研究对象和雷达资料进行选择。

最大值法和距离指数权重法改进的LAPS雷达反射率仍然存在一些资料空白区，尤其是对流层高层和低层。本文针对这个问题，尝试通过最小二乘法拟合来进行填充。研究结果表明，该填充方法对本文个例中对流层低层的静锥区有一定程度的改善，填充之后的低云系较为连续，但是梯度小，这个与建立的多元线性回归方程有很大的关系，还受取样范围的影响，不能较为精确地模拟出静锥区的真实值。

总体而言，本文针对LAPS雷达反射率拼图方法的改进较为成功，填充试验一定程度上改善了静锥区资料空白，但是填充之后会引入一些新的非真实的不连续回波，主要是因为填充方法存在一些不足。填充方法存在的局限和需要进一步改进的地方包括：(1) 多元线性回归方程不能够较为精确地模拟出雷达回波的空间分布形态，需要考虑其他方法；(2) 填充的条件也有待精细化，避免引入过多的非真实不连续回波。此外，想要得到更好的拼图效果，也可以尝试改进LAPS单部雷达插值程序。