基于双多普勒雷达反演的区域三维风场拼图方法与流程

本申请属于大气遥感监测分析技术领域，具体涉及一种基于双多普勒雷达反演的区域三维风场拼图方法。

背景技术：

多普勒天气雷达除了提供回波强度外，还能够提供降水粒子的径向速度。鉴于每部多普勒天气雷达都有一定的观测范围，为了监测分析大范围的天气系统，需要将多部多普勒天气雷达的信息进行拼图分析。

目前，国内外多普勒雷达拼图应用主要集中在回波强度拼图，即组合反射率拼图，其原理是对覆盖某点上方的所有雷达的所有体扫层的回波强度进行比较，选取其中的最大值作为该点的组合反射率值。然而，由于多普勒雷达径向速度扫描半径远小于回波强度扫描半径，例如，SA雷达的径向速度扫描半径通常为230km，SA雷达的回波强度扫描半径通常为460km，SA雷达的径向速度扫描半径远小于回波强度扫描半径；以及雷达风场拼图技术发展还不够成熟，因此雷达风场拼图还未能得到较好的应用。

中国气象科学研究院灾害天气国家重点实验室在公益性行业专项“热带西太平洋观测试验与我国高影响天气可预报性研究”的支持下，在自主研究的雷达资料质量控制和三维组网技术以及兰州大学的3DVAR单部达风场反演方法的基础上，研发了基于3DVAR的区域多普勒天气雷达组网风场反演技术，实现了反演区域高时、空分辨率三维风场的目的。该项技术充分发挥了单雷达风场反演区域范围大的优势，能够较为快速的实现大范围的雷达风场拼图。但是，由于其基于单雷达风场反演方法，只用到单部雷达的径向速度，相对于双雷达风场反演来说，存在较大的误差。对于多普勒雷达密集区，这项技术未能充分利用多部雷达的径向速度信息。

技术实现要素：

为至少在一定程度上克服相关技术中存在的问题，本申请提供了一种基于双多普勒雷达反演的区域三维风场拼图方法。

根据本申请实施例的第一方面，本申请提供了一种基于双多普勒雷达反演的区域三维风场拼图方法，其包括以下步骤：

根据两多普勒雷达之间的距离对多部多普勒雷达进行分组；

对每组双多普勒雷达均进行双多普勒雷达三维风场反演；

对多组双多普勒雷达的反演结果进行权重计算，得到反演网格点处的径向速度，根据径向速度和质量连续方程得到反演网格点处的垂直速度；

根据所有反演网格点处的径向速度和垂直速度得到区域三维风场的拼图。

进一步地，所述根据两多普勒雷达之间的距离对多部多普勒雷达进行分组的具体过程为：

根据用户对雷达风场反演区域的设定，统计所述雷达风场反演区域内的多普勒雷达；

统计所述雷达风场反演区域四周预定范围内的多普勒雷达；

对统计的所有多普勒雷达中的任意两个多普勒雷达，计算二者之间的距离；

将距离小于或等于预设的距离阈值的两个多普勒雷达作为一组双多普勒雷达。

进一步地，所述对每组双多普勒雷达均进行双多普勒雷达三维风场反演的具体过程为：

确定动态地球坐标系；

采用与动态地球坐标系相应的等经纬度、等海拔高度网格作为反演网格；

读取每组双多普勒雷达中多普勒雷达的基数据，将双多普勒雷达的径向速度和回波强度分别插值到等经纬度、等海拔高度的网格反演网格点上；

利用插值后的双多普勒雷达探测到的目标点的径向速度和回波强度计算得到目标点处的三维风速分量。

更进一步地，所述动态地球坐标系的坐标原点O设置在地球球心，Z轴随反演目标点的变化而变化。

更进一步地，所述动态地球坐标系中，Z轴由原点O通过反演目标点指向天顶方向，Y轴由原点O指向反演目标点的正北方向，X轴由原点O指向反演目标点的正东方向。

更进一步地，所述利用插值后的双多普勒雷达探测到的目标点的径向速度和回波强度计算得到目标点处的三维风速分量的具体过程为：

设目标点为P(x,y,z)，由于Z轴通过目标点P(x,y,z)，因此，x＝y＝0；在动态地球坐标系下，如果不考虑大气折射，则两部多普勒雷达探测到的目标点P(x,y,z)的两个方向的径向速度分别为：

其中，R₁和R₂分别表示两部多普勒雷达到地球球心点O的距离，

w_t表示降水粒子的下落速度，降水粒子的下落速度经验公式为：w_t＝3.8I^0.072，I表示回波强度；(x₀₁,y₀₁,z₀₁)、(x₀₂,y₀₂,z₀₂)分别表示两部多普勒雷达天线的位置；

将同一网格点上两部雷达观测到的径向速度近似作为水平风矢量进行合成，生成合成风速；将合成风速作为水平风速初始值，垂直速度初始值设为0；利用质量连续方程以及降水粒子的下落速度经验公式，进行迭代计算，直至两次迭代的误差小于预置值，最终获得雷达风场反演区域的三维风场。

更进一步地，所述雷达风场反演区域的三维风场的计算过程为：

1)假设w(z＝z₀)＝0，将w(z＝z₀)＝0带入质量连续方程中，计算得到径向速度第一次迭代值u⁽¹⁾和v⁽¹⁾；

2)将径向速度第一次迭代值u⁽¹⁾和v⁽¹⁾带入不考虑标准大气折射对垂直风速投影影响的情况下的径向速度公式中，计算得到垂直速度的第一估计值w⁽¹⁾；

3)将垂直速度的第一估计值w⁽¹⁾带入质量连续方程，重新计算径向速度u和v，直到收敛为止。

更进一步地，所述步骤3)中收敛条件设置为：

max|w^(k)-w^(k-1)|<ε,ε＝0.0001m/s。

进一步地，所述对多组双多普勒雷达的反演结果进行权重计算得到反演网格点处的径向速度，根据径向速度和质量连续方程得到反演网格点处的垂直速度的具体过程为：

对于某反演网格点，如果有m组有效双多普勒雷达，则该反演网格点处的径向速度为：

式中，A_i表示夹角权重系数，其中，夹角为目标点与双多普勒雷达之间形成的夹角；B_i表示距离权重系数，B_i＝D_i/D_max，其中，D_i表示第i组双多普勒雷达之间的距离，D_max表示所有组双多普勒雷达中两多普勒雷达之间距离的最大值；

对于所有的反演网格点，利用径向速度和质量连续方程计算得到垂直速度。

更进一步地，当所述夹角大于20°且小于或等于90°时，夹角权重系数为：A_i＝100.-|90.-θ_i|；当所述夹角大于90°且小于160°时，夹角权重系数为：A_i＝100.-1.2*|90.-θ_i|；其中，θ_i表示夹角，100.表示浮点数，90.表示浮点数。

本申请的实施例提供的技术方案可以包括以下有益效果：本申请针对多普勒雷达密集区，采用动态地球坐标系下双多普勒雷达风场反演算法进行三维风场拼图。本申请采用动态地球坐标系，以及等经纬度、等海拔高度网格，便于将不同双多普勒雷达反演的风场进行组网分析。通过对区域内多普勒雷达进行分组，形成多组双多普勒雷达，对各组多普勒双雷达分别进行风场反演，最后对反演网格点上的多个反演值进行优选或根据权重系数进行计算确定，从而获得区域三维风场拼图。本申请充分利用多普勒雷达密集区的多部多普勒雷达的径向速度信息，能够获得比单雷达更为准确的风场反演产品。

应当理解的是，以上的一般描述和后文的细节描述仅是示例性和解释性的，并不能限制本申请。

附图说明

此处的附图被并入说明书中并构成本说明书的一部分，示出了符合本申请的实施例，并与说明书一起用于解释本申请的原理。

图1是根据一示例性实施例示出的一种基于双多普勒雷达反演的区域三维风场拼图方法的流程图。

图2是根据一示例性实施例示出的一种基于双多普勒雷达反演的区域三维风场拼图方法中动态地球坐标系下双多普勒雷达反演示意图。

具体实施方式

这里将详细地对示例性实施例进行说明，其示例表示在附图中。下面的描述涉及附图时，除非另有表示，不同附图中的相同数字表示相同或相似的要素。以下示例性实施例中所描述的实施方式并不代表与本申请相一致的所有实施方式。相反，它们仅是与如所附权利要求书中所详述的、本申请的一些方面相一致的装置和方法的例子。

图1是本申请一实施例提供的基于双多普勒雷达反演的区域三维风场拼图方法的流程图。如图1所示，本申请基于双多普勒雷达反演的区域三维风场拼图方法包括以下步骤：

S1、根据两多普勒雷达之间的距离对多部多普勒雷达进行分组。

首先，根据用户对雷达风场反演区域即经纬度变化范围的设定，统计该雷达风场反演区域内的多普勒雷达。

其次，为充分发挥多普勒雷达探测信息的作用，统计雷达风场反演区域四周预定范围内的多普勒雷达，对这些多普勒雷达也进行分组。其中，雷达风场反演区域四周预定范围可以为雷达风场反演区域四周100km。

再次，对统计的所有多普勒雷达中的任意两个多普勒雷达，计算二者之间的距离。

最后，将距离小于或等于预设的距离阈值的两个多普勒雷达作为一组双多普勒雷达。其中，预设的距离阈值可以为100km。

S2、对每组双多普勒雷达均进行双多普勒雷达三维风场反演，其具体过程为：

S21、确定动态地球坐标系。

动态地球坐标系的坐标原点O设置在地球球心，Z轴随反演目标点的变化而变化。具体地，Z轴由原点O通过反演目标点指向天顶方向，Y轴由原点O指向反演目标点的正北方向，X轴由原点O指向反演目标点的正东方向。

S22、采用与动态地球坐标系相应的等经纬度、等海拔高度网格作为反演网格。

S23、读取每组双多普勒雷达中多普勒雷达的基数据，将双多普勒雷达的径向速度和回波强度分别插值到等经纬度、等海拔高度的网格反演网格点上。

S24、利用插值后的双多普勒雷达探测到的目标点的径向速度和回波强度计算得到目标点处的三维风速分量，其具体过程为：

S241、如图2所示，设目标点为P(x,y,z)，由于Z轴通过目标点P(x,y,z)，因此，x＝y＝0。在动态地球坐标系下，如果不考虑大气折射，则两部多普勒雷达探测到的目标点P(x,y,z)的两个方向的径向速度分别为：

其中，R₁和R₂分别表示两部多普勒雷达到地球球心点O的距离，

w_t表示降水粒子的下落速度，降水粒子的下落速度经验公式为：w_t＝3.8I^0.072，I表示回波强度；(x₀₁,y₀₁,z₀₁)、(x₀₂,y₀₂,z₀₂)分别表示两部多普勒雷达天线的位置。

S242、将同一网格点上两部雷达观测到的径向速度近似作为水平风矢量进行合成，生成合成风速。将合成风速作为水平风速初始值，垂直速度初始值设为0。利用质量连续方程(大气不可压缩近似)以及降水粒子的下落速度经验公式，进行迭代计算，直至两次迭代的误差小于预置值，最终获得雷达风场反演区域的三维风场。

采用迭代算法计算三维风场，其具体过程为：

1)假设w(z＝z₀)＝0，将w(z＝z₀)＝0带入质量连续方程中，计算得到径向速度第一次迭代值u⁽¹⁾和v⁽¹⁾。

2)将径向速度第一次迭代值u⁽¹⁾和v⁽¹⁾带入不考虑标准大气折射对垂直风速投影影响的情况下的径向速度公式(1)和(2)中，计算得到垂直速度的第一估计值w⁽¹⁾。

3)将垂直速度的第一估计值w⁽¹⁾带入质量连续方程，重新计算径向速度u和v，直到收敛为止，收敛条件为：

max|w^(k)-w^(k-1)|<ε,ε＝0.0001m/s。

S3、对多组双多普勒雷达的反演结果进行权重计算，得到反演网格点处的径向速度，根据径向速度和质量连续方程得到反演网格点处的垂直速度，其具体过程为：

由于同一反演网格点所在的反演区域往往有多组双多普勒雷达，因此需要按照一定的规则对反演结果进行取舍或者按一定权重系数对反演结果进行计算。

本申请主要考虑反演网格点与两部多普勒雷达之间的夹角以及两部多普勒雷达之间的距离这两个因素：

(1)根据已有分析结果，反演误差的大小与目标点和双多普勒雷达之间形成的夹角有关，形成的夹角越接近90°，反演误差越小。此外，考虑到双多普勒雷达连线附近的反演误差较大，对反演网格点来说，只有与目标点形成夹角大于20°且小于160°的双多普勒雷达才对其进行权重考虑。

当夹角θ_i大于20°且小于或等于90°时，夹角权重系数A_i＝100.-|90.-θ_i|；

当夹角θ_i大于90°且小于160°时，夹角权重系数A_i＝100.-1.2*|90.-θ_i|。

其中，100.表示浮点数，90.表示浮点数。

(2)由于多普勒雷达的每根扫描线的宽度为1°，因此，与雷达距离越近，空间分辨率越高，距离较近的双多普勒雷达附近区域的反演风场也就更为精确。为此，设定第i组双多普勒雷达的距离权重系数为：B_i＝D_i/D_max，其中D_i表示第i组双多普勒雷达之间的距离，D_max表示所有组双多普勒雷达中两多普勒雷达之间距离的最大值。

对于某反演网格点而言，如果有m组有效双多普勒雷达，则该反演网格点处的径向速度为：

对于所有的反演网格点，计算得到径向速度后，就可以通过质量连续方程进行积分计算得到垂直速度。

S4、根据所有反演网格点处的径向速度和垂直速度得到区域三维风场的拼图。

可以理解的是，上述各实施例中相同或相似部分可以相互参考，在一些实施例中未详细说明的内容可以参见其他实施例中相同或相似的内容。

需要说明的是，在本申请的描述中，术语“第一”、“第二”等仅用于描述目的，而不能理解为指示或暗示相对重要性。此外，在本申请的描述中，除非另有说明，“多个”的含义是指至少两个。

流程图中或在此以其他方式描述的任何过程或方法描述可以被理解为，表示包括一个或更多个用于实现特定逻辑功能或过程的步骤的可执行指令的代码的模块、片段或部分，并且本申请的优选实施方式的范围包括另外的实现，其中可以不按所示出或讨论的顺序，包括根据所涉及的功能按基本同时的方式或按相反的顺序，来执行功能，这应被本申请的实施例所属技术领域的技术人员所理解。

本技术领域的普通技术人员可以理解实现上述实施例方法携带的全部或部分步骤是可以通过程序来指令相关的硬件完成，所述的程序可以存储于一种计算机可读存储介质中，该程序在执行时，包括方法实施例的步骤之一或其组合。

此外，在本申请各个实施例中的各功能单元可以集成在一个处理模块中，也可以是各个单元单独物理存在，也可以两个或两个以上单元集成在一个模块中。上述集成的模块既可以采用硬件的形式实现，也可以采用软件功能模块的形式实现。所述集成的模块如果以软件功能模块的形式实现并作为独立的产品销售或使用时，也可以存储在一个计算机可读取存储介质中。

上述提到的存储介质可以是只读存储器，磁盘或光盘等。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本申请的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

尽管上面已经示出和描述了本申请的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本申请的限制，本领域的普通技术人员在本申请的范围内可以对上述实施例进行变化、修改、替换和变型。