一种双多普勒雷达三维风场反演方法

摘要

本发明涉及一种双多普勒雷达三维风场反演方法，其包括以下步骤：确定动态地球坐标系；采用与动态地球坐标系相应的等经纬度、等海拔高度网格作为反演网格；分析大气折射对风场反演的影响，得到目标点处的雷达波束与该处水平面的夹角；利用目标点处的雷达波束与该处水平面的夹角，得到考虑了标准大气折射对垂直风速投影影响的情况下的径向速度；将同一网格点上两部雷达观测到的径向速度近似作为水平风矢量进行合成，生成合成风速；借助质量连续方程以及降水粒子的下落速度经验公式，进行迭代计算，直至两次迭代的误差小于预置值，最终获得降水回波区域的三维风场。本发明考虑标准大气折射对双雷达风场反演的影响，能够提高风场反演的准确性。

说明书

技术领域

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

背景技术

多普勒雷达应用多普勒频移原理，通过测量降水粒子反射回来的电磁波与发射电磁波间频率的变化，来测定降水粒子相对于雷达移近或移远的速度。这种速度称为径向速度。由于无法测量降水粒子的真实矢量速度，因此单雷达在观测风场方面受到很大的限制。为了突破单雷达的限制，国内外科研人员发展出双多普勒雷达风场反演技术，即将处于一定距离的两部多普勒雷达所观测的径向速度进行合成，并借助相关假定和方程反演出降水系统内部的三维风场。

目前，国内外普遍采用的是笛卡尔坐标系下的双多普勒雷达风场反演方法。该方法采用笛卡尔坐标系，坐标原点O设置在双雷达中的某一雷达处，或者两部雷达连线的中点上，X轴指向正东方，Y轴指向正北方，Z轴指向天顶，反演网格采用的是等X、Y和Z坐标点。这种方法主要存在以下缺点：1)由于地球是近似球体，地面近似为球面，因此在现有反演技术采用的笛卡尔坐标系下， X-O-Y平面与地球表面不平行，目标点距离原点O越远，X-O-Y平面距离地面就越高。采用该方法反演出的网格点与气象业务科研常用的等经纬度、等海拔高度的网格点不一致，虽然这对小范围天气系统三维风场结构分析影响不大，但风场反演产品难以与数值预报等其他气象信息进行综合分析，直接影响业务应用效果。2)随着我国多普勒天气雷达网的进一步完善，雷达风场拼图分析显得越来越迫切。对于不同的双雷达组合，利用笛卡尔坐标系下的反演技术就需要设置不同的坐标系，导致不同双雷达反演的三维风场产品无法组合到统一坐标系里进行组网分析应用，这在一定程度上影响了雷达观测网建设效益的发挥。

为克服笛卡尔坐标系下双多普勒雷达风场反演方法存在的缺点，有人提出了动态地球坐标系下双多普勒雷达风场反演方法，该方法采用动态地球坐标系，以及等经纬度和等海拔高度的网格。这种方法既便于将反演风场与其他产品进行综合分析，又便于将不同双雷达反演的风场进行组网分析，从而能够有效地克服笛卡尔坐标系下反演方法存在的缺点。然而，由于大气实际上是非均质流体，导致雷达发射的电磁波在大气中存在折射现象，电磁波的传播路径并不是直线，因此会对降水粒子速度在电磁波传播方向上的投影(即径向速度)产生影响。目前普遍采用的上述两种风场反演方法均没有考虑大气折射的影响，而是假定电磁波按直线传播。虽然在较小的空间范围内该假定的影响是微乎其微的，但在实际业务应用中，双雷达之间的距离大多超过100公里，大气折射的影响是有必要考虑的。

发明内容

为了解决现有技术存在的上述问题，本发明提供了一种双多普勒雷达三维风场反演方法。

本发明所采用的技术方案为：一种双多普勒雷达三维风场反演方法包括以下步骤：

确定动态地球坐标系；

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

分析大气折射对风场反演的影响，得到目标点处的雷达波束与该处水平面的夹角；

利用目标点处的雷达波束与该处水平面的夹角，得到考虑了标准大气折射对垂直风速投影影响的情况下的径向速度；

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

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

进一步地，所述动态地球坐标系的坐标原点O设置在地球球心，Y轴由原点O指向北极，Z轴由原点O指向雷达站天顶方向，X轴与Z轴成90°并指向目标点一侧。

进一步地，所述步骤分析大气折射对风场反演的影响的具体过程为：

推导雷达波束方程；

在已知雷达天线仰角、目标点与雷达相对于地球球心所成的夹角、目标点的海拔高度中的两个参量时，通过雷达波束方程求取第三个参量；

计算目标点处雷达波束与该处水平面的夹角。

更进一步地，所述步骤推导雷达波束方程的具体过程为：

设P(x,z)是雷达天线仰角为δ的波束上任一点，在ΔOCO₁中，由正弦公式得到即

又由β＝α-γ及余弦公式得到|OC|²＝|OO₁|²+|O₁C|²-2|OO₁||O₁C|cosβ，设>

即

式中，Δh表示标准大气中由于折射现象而引起的高度的变化调整；α表示目标点P与雷达P₀相对于地球球心O所成的夹角；β表示点A与点C相对于等效地球球心O₁所成的夹角；γ表示地球球心O与地球球心O₁相对于点C所成的夹角；>m表示标准大气地球等效半径；点A表示雷达P₀所在经纬度海平面点，点B表示目标点P所在经纬度海平面点，点C表示地球球心点O与目标点P的连线与标准大气地球等效表面的交点；h₀表示雷达P₀的海拔高度，h>

得到极坐标形式的雷达波束方程为：

更进一步地，雷达波束方程的简化形式为：

式中，H₀＝R+h₀，F＝2R_m(R_m-R)，

更进一步地，所述步骤计算目标点处雷达波束与该处水平面的夹角的具体过程为：

计算雷达波束上目标点P(x,z)处的斜率：

波束上目标点P(x,z)处的斜率为：

其中，

波束上目标点P处的径向速度方向与X轴的夹角为目标点P处的水平线与X轴的夹角为-α；设目标点P处雷达波束与该处水平面的夹角为θ，则有

更进一步地，所述考虑了标准大气折射对垂直风速投影影响的情况下的径向速度为：

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

更进一步地，所述步骤获得降水回波区域的三维风场的具体过程为：

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

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

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

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

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

由于采用以上技术方案，本发明的有益效果为：本发明在动态地球坐标系下，采用等经纬度、等海拔高度网格，考虑标准大气折射对双雷达风场反演的影响，能够提高风场反演的准确性；既便于将反演风场结果与其他产品进行综合分析，又便于将不同双雷达反演的风场进行组网分析。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1是本发明一实施例提供的一种双多普勒雷达三维风场反演方法的流程图；

图2是本发明一实施例提供的一种双多普勒雷达三维风场反演方法中雷达波束路径受大气折射影响示意图；

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将对本发明的技术方案进行详细的描述。显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所得到的所有其它实施方式，都属于本发明所保护的范围。

本发明双多普勒雷达三维风场反演方法将常用的笛卡尔坐标系变换为动态地球坐标系，同时，将风场反演网格点由等XYZ网格改进为等经纬度和等海拔高度网格。这样得到的反演结果与地理位置对应关系更为简单，便于与台站或数值产品进行对比分析，同时也便于将反演结果用于数值模式计算以及雷达组网拼图的实现。其中，动态地球坐标系的原点位于地球球心，Z轴随反演目标点的变化而变化。

在不考虑大气折射的情况下，雷达波束是直线；考虑大气折射后，雷达波束变为曲线。大气折射对双雷达风场反演至少有以下三点影响：1)波束上目标点径向速度方向并不是雷达与目标点的连线方向，而是波束的切线方向；2)雷达扫描图像上库点与雷达的距离，不再是直线距离，而是波束路径长度即弧长； 3)波束上目标点海拔高度发生变化。因此，本发明提出了一种动态地球坐标系下考虑标准大气折射影响的双多普勒雷达三维风场反演方法。

如图1所示，本发明提供了一种双多普勒雷达三维风场反演方法，其包括以下步骤：

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

动态地球坐标系的坐标原点O设置在地球球心，Z轴随反演目标点的变化而变化。具体地，Y轴由原点O指向北极，Z轴由原点O指向雷达站天顶方向，X轴与Z轴成90°并指向目标点一侧。

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

S3、分析大气折射对风场反演的影响，得到目标点处的雷达波束与该处水平面的夹角。

如图2所示，利用标准大气地球等效半径的近似，分析大气折射对风场反演的影响。地球球心为点O，地球半径为R。不考虑大气折射时标准大气等效地球半径所对应的球心为点O₁，标准大气地球等效半径为R_m。雷达P₀所在经纬度海平面点为A，雷达P₀的海拔高度为h₀。目标点P所在经纬度海平面点为B，目标点P的海拔高度为h。雷达P₀的波束经过目标点P，雷达P₀的仰角为δ。不考虑大气折射时的δ仰角波束与目标点P天顶方向的交点为Q。地球球心点O>

下面给出标准大气折射对风场反演的影响的具体计算过程：

首先，推导雷达波束方程。

设P(x,z)是雷达天线仰角为δ的波束上任一点，在ΔOCO₁中，由正弦公式得到即

又由β＝α-γ及余弦公式得到|OC|²＝|OO₁|²+|O₁C|²-2|OO₁||O₁C|cosβ，设>

即

式中，Δh表示标准大气中由于折射现象而引起的高度的变化调整，故有 |PQ|＝|BC|＝Δh；α表示目标点P与雷达P₀相对于地球球心O所成的夹角；β表示点A与点C相对于等效地球球心O₁所成的夹角；γ表示地球球心O与地球球心O₁相对于点C所成的夹角。

由此，得到极坐标形式的雷达波束方程为：

令H₀＝R+h₀，F＝2R_m(R_m-R)，则上述雷达波束方程可以简化为：

特殊地，当时，α＝0，由雷达波束方程得到h＝R，与实际情况不符。此时的波束方程可改用α＝0表示。

其次，在已知雷达天线仰角δ、目标点与雷达相对于地球球心所成的夹角α、目标点的海拔高度h三个参量中的两个时，通过雷达波束方程可以求取第三个参量。

最后，计算目标点处雷达波束与该处水平面的夹角。

在标准大气折射的假定下，雷达波束向下弯曲，导致所测径向速度方向发生变化，同时导致垂直速度在波束上的投影发生变化。下面求取雷达波束上目标点P(x,z)处的斜率，即切线方向，亦即径向速度方向。由图2可知：

波束上目标点P(x,z)处的斜率为：

其中，

因此，如图2所示，波束上目标点P处的径向速度方向与X轴的夹角为目标点P处的水平线与X轴的夹角为-α。设目标点P处雷达波束与该处水平面的夹角为θ，则有

S4、利用目标点处的雷达波束与该处水平面的夹角，得到考虑了标准大气折射对垂直风速投影影响的情况下的径向速度。

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

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

考虑大气折射后，由于波束向下弯曲，导致垂直速度在波束切线方向即雷达所观测到的径向速度方向上投影的变化。用sinθ代替和得到考虑了标准大气折射对垂直风速投影影响的情况下的径向速度：

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

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

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

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

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

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

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

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以所述权利要求的保护范围为准。