CN113325423A - 一种多普勒气象雷达数据采集及三维拼图方法 - Google Patents

一种多普勒气象雷达数据采集及三维拼图方法 Download PDF

Info

Publication number
CN113325423A
CN113325423A CN202110493990.7A CN202110493990A CN113325423A CN 113325423 A CN113325423 A CN 113325423A CN 202110493990 A CN202110493990 A CN 202110493990A CN 113325423 A CN113325423 A CN 113325423A
Authority
CN
China
Prior art keywords
radar
coordinate system
grid point
grid
cartesian coordinate
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN202110493990.7A
Other languages
English (en)
Inventor
卢明
郭志民
李哲
梁允
姚德贵
刘善峰
王超
王津宇
王磊
李帅
杨磊
苑司坤
高阳
崔晶晶
马建伟
耿俊成
吕卓
陈岑
刘莘昱
张小斐
刘昊
杨文�
刘枫琪
马彩霞
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
State Grid Corp of China SGCC
Electric Power Research Institute of State Grid Henan Electric Power Co Ltd
Henan Jiuyu Enpai Power Technology Co Ltd
Original Assignee
State Grid Corp of China SGCC
Electric Power Research Institute of State Grid Henan Electric Power Co Ltd
Henan Jiuyu Enpai Power Technology Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by State Grid Corp of China SGCC, Electric Power Research Institute of State Grid Henan Electric Power Co Ltd, Henan Jiuyu Enpai Power Technology Co Ltd filed Critical State Grid Corp of China SGCC
Priority to CN202110493990.7A priority Critical patent/CN113325423A/zh
Publication of CN113325423A publication Critical patent/CN113325423A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/95—Radar or analogous systems specially adapted for specific applications for meteorological use
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02A—TECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00—Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10—Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本申请公开了一种多普勒气象雷达数据采集及三维拼图方法，包括：采集目标区域范围内球坐标系下多部多普勒气象雷达数据；在笛卡尔坐标系上划分若干网格点；检索出覆盖每一笛卡尔坐标系网格点的雷达及其雷达数据；根据笛卡尔坐标系网格点坐标计算对应每一雷达的球坐标系网格点坐标；插值计算每一雷达球坐标系网格点的反射率因子值；笛卡尔坐标系网格点赋值，完成本次三维拼图。本发明综合应用把多个雷达资料进行拼图处理，实现了区域范围内多部多普勒雷达拼接，扩大了区域范围内的天气雷达网覆盖范围，解决了地面固定的单部雷达的探测范围有限和强对流预警范围及准确度受限的问题。

Description

一种多普勒气象雷达数据采集及三维拼图方法
技术领域
本发明属于电力气象数据处理技术领域，涉及一种多普勒气象雷达数据采集及三维拼图方法。
背景技术
电网的快速发展及规模的不断扩大和延伸，并随着全球气候变暖、环流异常，灾害性天气风险逐年增大风险，灾害性天气对电网安全运行的影响越来越大。电力气象是电力与气象领域相结合形成一个交叉学科，并不仅仅是气象部门针对电力行业的气象服务，而是面向电力生产所需开展的专门的气象预报、气象监测、电网灾害预测以及电网灾害风险预警工作。
目前，中国气象局已经在全国范围内建成了新一代天气雷达网，这些雷达的基数据通过气象系统专网及时传送到省气象局和中国气象局，且雷达绝大部分范围覆盖情况很好，是监测并预警强对流天气的重要工具，这就为开展地面监测和预警提供了条件。地面固定的单部雷达的探测范围有限,不足以覆盖更大尺度的天气系统,例如梅雨锋、台风、飑线等。这一缺陷使天气系统的发生发展机理研究无法做到全面、客观、细致，也造成了强对流预警范围和准确度受限。因此,为了提高对中尺度灾害性天气的研究以及预警能力，必须把来自多部雷达的资料进行组网拼图。
由多个雷达重复取样获得的天气信息要比只由单个雷达取样获得的更加精确,还可以很大程度解决因单部雷达观测的波束几何学原因(例如,静锥区、波束展宽、波束高度、波束阻挡等)引起的很多问题，扩大天气雷达网的范围是推动风暴尺度天气预报研究的主要动力。
发明内容
为解决现有技术中的不足，本申请提供一种多普勒气象雷达数据采集及三维拼图方法。
为了实现上述目标，本发明采用如下技术方案：
一种多普勒气象雷达数据采集及三维拼图方法，包括以下步骤：
步骤1，采集目标区域范围内球坐标系下多部多普勒气象雷达数据；
步骤2，根据目标区域定义笛卡尔坐标系网格范围和分辨率，在笛卡尔坐标系上划分若干网格点；
步骤3，根据笛卡尔坐标系网格点坐标，逐一检索出覆盖每一笛卡尔坐标系网格点的雷达及其雷达数据；
步骤4，根据笛卡尔坐标系网格点坐标计算对应每一雷达的球坐标系网格点坐标；
步骤5，插值计算每一雷达球坐标系网格点的反射率因子值，作为雷达在笛卡尔坐标系网格点的反射率分析值；
步骤6，将覆盖笛卡尔坐标系网格点的雷达反射率分析值赋值给该笛卡尔坐标系网格点，得到该笛卡尔坐标系网格点的雷达反射率因子值，所有笛卡尔坐标系网格点完成赋值后，即完成本次三维拼图。
本发明进一步包括以下优选方案：
优选地，步骤2中，所述笛卡尔坐标系网格在垂直方向选取的范围为0-17km，共21层；
5km以下的垂直分辨率为0.5km，5-17km之间的垂直分辨率为1km。
在水平方向的经纬度分辨率为0.01°×0.01°，范围是覆盖目标区域。
优选地，步骤4中，利用笛卡尔坐标系下网格点的经度、纬度和高度计算其对应雷达在球坐标系中的仰角e、方位角a和斜距r；
设笛卡尔坐标系三维网格中任意网格点的坐标为(αg，βg，hg)，其中αg为网格点纬度，βg为网格点经度，hg为网格点高度；
雷达天线所在点的坐标为(αr，βr，hr)，其中αr为雷达天线纬度，βr为雷达天线经度，hr为雷达天线高度；
则使用雷达波束传播和大圆几何学理论确定笛卡尔坐标系三维网格中网格点对应雷达天线所在点的球坐标位置(r，a，e)，其中r为斜距，a为方位角，e为仰角，具体的：
由球面三角公式得出：
sina＝cos(αg)sin(βg-βr)/sin(s/R) (公式1)
其中R为地球半径，s为大圆距离，其表达式为：
s＝Rcos-1(sin(αr)sin(αg)+cos(αr)cos(αg)cos(βg-βr)) (公式2)
设C＝sina，则：
Figure BDA0003053548560000031
Figure BDA0003053548560000032
其中Rm为等效地球半径，
Figure BDA0003053548560000033
r＝sin(s/Rm)(Rm+hg-hr)/cos(e) (公式5)。
优选地，步骤5中，采用最近邻居法计算球坐标系网格点的反射率因子值，具体为：
计算球坐标系中的网格点的中心与雷达距离库中心之间的距离，用最靠近网格点的雷达距离库的观测值，即球坐标系下的雷达反射率因子值去填充球坐标系网格点的雷达反射率因子值，作为雷达在笛卡尔坐标系网格点的反射率分析值；
所述雷达距离库中心即雷达有效照射体积的中心。
优选地，步骤5中，采用径向和方位上的最近邻居和垂直线性内插法球坐标系网格点的反射率因子值，具体为：
假设(r,a,e)是某一网格点在雷达球坐标系中的位置，r为斜距，a为方位角，e为仰角；e位于其上下相邻仰角e2和e1之间；
(r,a,e2)和(r,a,e1)分别是经过该网格点的垂线与其上下仰角波束轴线的交点，那么：
该球坐标系网格点的反射率因子值fa(r,a,e)用(r,a,e2)点反射率因子值fa(r,a,e2)和(r,a,e1)点反射率因子值fa(r,a,e1)进行垂直线性内插得到，即：
fa(r,a,e)＝(we1fa(r,a,e1)+we2fa(r,a,e2))/(we1+we2) (公式6)
其中，
we1、we2分别为fa(r,a,e1)和fa(r,a,e2)的内插权重：
we1＝(e2-e)/(e2-e1) (公式7)
we2＝(e-e1)/(e2-e1) (公式8)
fa(r,a,e2)和fa(r,a,e1)分别等于最靠近点(r,a,e2)和(r,a,e1)的雷达距离库的观测值，即球坐标系下的雷达反射率因子值，采用径向和方位上的最近邻居法获取。
优选地，假设球坐标系中，ri-1、ri、ri+1为相邻径向距离库，aj-1、aj、aj+1为相邻方位角，由半功率线和半距离库所围成的梯形区是距离库ri的影响区；
在径向、方位方向上落在所述梯形区的点(r,a)的反射率因子值fa(r,a)用距离库ri的观测值fo(ri,aj)来赋值，即fa(r,a)＝fo(ri,aj)。
优选地，步骤5中，采用垂直水平线性内插法球坐标系网格点的反射率因子值，具体为：
假设(r,a,e2)和(r,a,e1)分别是经过网格点(r,a,e)的垂线与其上下仰角波束轴线的交点，(r1,a,e2)、(r2,a,e1)分别是经过该网格点的水平线与其相邻上下仰角波束轴线的交点，那么：
该网格点的反射率因子值fa(r,a,e)用(r,a,e2)、(r,a,e1)、(r1,a,e2)、(r2,a,e1)点反射率因子值fa(r,a,e2)、fa(r,a,e1)、fa(r1,a,e2)、fa(r2,a,e1)通过垂直和水平内插得到；
其中，fa(r,a,e2)、fa(r,a,e1)、fa(r1,a,e2)、fa(r2,a,e1)通过径向和方位的最近邻居法得到，则：
Figure BDA0003053548560000041
其中，wr1、wr2分别是fa(r1,a,e2)、fa(r2,a,e1)的内插权重：
wr1＝(r2-r)/(r2-r1) (公式10)
wr2＝(r-r1)/(r2-r1) (公式11)
且有r1＝rsine/sine2，r2＝rsine/sine1。
优选地，步骤5中，采用8点插值法球坐标系网格点的反射率因子值，具体为：
假设(r,a,e2)和(r,a,e1)分别是经过网格点(r,a,e)的垂线与其上下仰角波束轴线的交点，(r1,a,e2)、(r2,a,e1)分别是经过该网格点的水平线与其相邻上下仰角波束轴线的交点，那么：
若某一球坐标系网格点(r,a,e)落在由(r1,a1,e1)、(r2,a1,e1)、(r1,a2,e1)、(r2,a2,e1)、(r1,a1,e2)、(r2,a1,e2)、(r1,a2,e2)、(r2,a2,e2)点围成的锥体内，则该网格点的反射率因子值fa(r,a,e)由这8个点的观测值
Figure BDA0003053548560000051
Figure BDA0003053548560000052
Figure BDA0003053548560000053
进行双线性内插获得：
Figure BDA0003053548560000054
其中，wa1、wa2为方位内插权重：
wa1＝(a2-a)/(a2-a1) (公式13)
wa2＝(a-a1)/(a2-a1) (公式14)；
wr1是
Figure BDA0003053548560000055
的内插权重，wr2是
Figure BDA0003053548560000056
的内插权重；
we1是
Figure BDA0003053548560000057
的内插权重；
we2是
Figure BDA0003053548560000058
的内插权重。
优选地，步骤6中，根据步骤5的反射率因子值计算笛卡尔坐标系每个网格点i的反射率因子值，实现目标区域范围内多部多普勒气象雷达数据拼接。
笛卡尔坐标系网格点i的雷达反射率因子值通过下面公式得到：
Figure BDA0003053548560000061
其中，
fm(i)是网格点i的合成反射率因子值；
Figure BDA0003053548560000062
是在网格点i处来自第n个雷达的反射率分析值；
wn是反射率分析值
Figure BDA0003053548560000063
的权重；
Nrad是在网格点i处有反射率分析值的总的雷达个数。
优选地，所述权重采用指数权重函数：
Figure BDA0003053548560000064
其中D为适当的长度比例，取D＝150，d为网格点到雷达的距离。
优选地，所述权重采用Cressman权重函数：
Figure BDA0003053548560000065
其中D为影响半径，取D＝300，d为网格点到雷达的距离。
优选地，步骤6中，采用最近邻居法计算笛卡尔坐标系网格点的雷达反射率因子值，具体的：
将来自最靠近笛卡尔坐标系网格点的那个雷达的反射率分析值的权重赋为1，其它的权重全赋为0，即把最靠近笛卡尔坐标系网格点的那个雷达的反射率分析值赋给网格点。
优选地，步骤6中，采用最大值法计算笛卡尔坐标系网格点的雷达反射率因子值，具体的：
将覆盖同一笛卡尔坐标系网格点的多个雷达反射率分析值中的最大值的权重赋为1，其它的权重全赋为0，即把覆盖同一网格点的多个雷达反射率分析值中的最大值赋给网格点。
本申请所达到的有益效果：
本发明综合应用把多个雷达资料进行拼图处理，实现了区域范围内多部多普勒雷达拼接，扩大了区域范围内的天气雷达网覆盖范围，解决了地面固定的单部雷达的探测范围有限和强对流预警范围及准确度受限的问题，为开展大尺度的天气系统,例如梅雨锋、台风、飑线的监测预警提供了条件。
由多个雷达重复取样获得的天气信息要比只由单个雷达取样获得的更加精确,提升了基于雷达监测开展天气预测预警的精确性和准确性。很大程度解决因单部雷达观测的波束几何学原因(例如,静锥区、波束展宽、波束高度、波束阻挡等)引起的雷达监测问题，扩大天气雷达网的范围，在一定程度上促进风暴尺度天气预报研究，可有效提高社会防灾减灾水平。
附图说明
图1为本发明一种多普勒气象雷达数据采集及三维拼图方法流程图。
具体实施方式
下面结合附图对本申请作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本申请的保护范围。
单部多普勒雷达覆盖范围有限，因此对于面积较大的区域需要部署多普勒雷达，以便覆盖和监测更大范围的区域。而每部雷达扫描数据都是独立的，因此要对整块区域进行统一监测，需要将此区域内多部多普勒雷达的雷达监测资料进行拼接后使用。如图1所示，本发明的一种多普勒气象雷达数据采集及三维拼图方法，包括：
单个雷达监测资料是球坐标系(即方位角、俯仰角和距离)，不同雷达监测资料之间的同一方位角、俯仰角和距离所标识反射率因子值并不相同也代表不是同一个位置，为了进行数据规范化统一和方便应用，需要进行开展雷达拼图工作，将多部雷达进行坐标转换后开展三维拼图工作，雷达拼图后的数据则是采用的笛卡尔坐标系(经度、纬度、高度)。
其次，雷达拼图首先要根据情况定义范围，例如定义某一省的范围定义为经度为东经111至115度，纬度为北纬33至36度范围，然后定义经纬度分辨率为0.01度，这样就在平面笛卡尔坐标系上划分了若干网格点。具体为步骤1-2。
步骤1，采集目标区域范围内球坐标系下多部多普勒气象雷达数据；
步骤2，根据目标区域定义笛卡尔坐标系网格范围和分辨率，在笛卡尔坐标系上划分若干网格点；
笛卡尔坐标系提供了一个统一框架，多部多普勒气象雷达资料能够在该框架下相互融合，这有利于多部多普勒气象雷达资料的综合应用，以便提供比单个观测系统对气象现象更加真实和科学合理的描述。
具体实施时，笛卡尔坐标系网格的分辨率和范围根据不同的需求进行不同的选择或根据球坐标系下的资料分辨率来决定。
笛卡尔坐标系网格在垂直方向选取的范围为0-17km，共21层；
5km以下的垂直分辨率为0.5km，5-17km之间的垂直分辨率为1km。这样垂直分层的原因是VCP21扫描方式的最低仰角0.5°在斜距230km处的高度约为5km，在斜距460km处的高度约为17km。在5km以下高度(离雷达230km以内)，雷达资料的空间密度比较大，所以垂直分辨率取的高一些，而在5-17km高度(离雷达230-460km)，雷达资料的空间密度比较小，所以垂直分辨率取的低一些。
在水平方向的经纬度分辨率为0.01°×0.01°(约1×1km)，平面方向，一般定义为能完整覆盖目标区域的经纬度范围，例如某一省的范围定义为经度为东经111至115度，纬度为北纬33至36度范围，然后定义经纬度分辨率为0.01度。
因为雷达拼图是三维的，上述垂直分辨率用1km表示，平面(即经度和纬度构成的平面)分辨率用0.01度表示。
为了给单个笛卡尔坐标系网格点去赋值反射率因子值等参数，则需要根据逐一对笛卡尔坐标系某一网格点的坐标信息，检索覆盖该笛卡尔坐标系网格点的雷达及其雷达基数据，然后根据该笛卡尔坐标系某一网格点的坐标信息计算对应某一雷达的球坐标系(即方位角、俯仰角和距离)网格点信息，然后插值出球坐标系网格点的反射率因子值等参数，然后在赋值给该笛卡尔坐标系网格点。具体为步骤3-6.
步骤3，根据笛卡尔坐标系网格点坐标，逐一检索出覆盖每一笛卡尔坐标系网格点的雷达及其雷达数据；
步骤4，根据笛卡尔坐标系网格点坐标计算对应每一雷达的球坐标系网格点坐标，具体的：
利用笛卡尔坐标系下网格点的经度、纬度和高度计算其对应雷达在球坐标系中的仰角e、方位角a和斜距r；
设笛卡尔坐标系三维网格中任意网格点的坐标为(αg，βg，hg)，其中αg为网格点纬度，βg为网格点经度，hg为网格点高度；
雷达天线所在点的坐标为(αr，βr，hr)，其中αr为雷达天线纬度，βr为雷达天线经度，hr为雷达天线高度；
则使用雷达波束传播和大圆几何学理论确定笛卡尔坐标系三维网格中网格点对应雷达天线所在点的球坐标位置(r，a，e)，其中r为斜距，a为方位角，e为仰角，具体的：
由球面三角公式得出：
sina＝cos(αg)sin(βg-βr)/sin(s/R) (公式1)
其中R为地球半径，s为大圆距离，其表达式为：
s＝Rcos-1(sin(αr)sin(αg)+cos(αr)cos(αg)cos(βg-βr)) (公式2)
设C＝sina，则：
Figure BDA0003053548560000091
Figure BDA0003053548560000092
其中Rm为等效地球半径，
Figure BDA0003053548560000093
r＝sin(s/Rm)(Rm+hg-hr)/cos(e) (公式5)。
步骤5，插值计算每一雷达球坐标系网格点的反射率因子值，作为雷达在笛卡尔坐标系网格点的反射率分析值；
具体实施时，可使用如下四种方法得到球坐标系网格点的反射率因子值：
(1)采用最近邻居法计算球坐标系网格点的反射率因子值，具体为：
计算球坐标系中的网格点的中心与雷达距离库中心之间的距离，用最靠近网格点的雷达距离库的观测值(即球坐标系下的雷达反射率因子值)去填充球坐标网格点的雷达反射率因子值，作为雷达在笛卡尔坐标系网格点的反射率分析值；
所述雷达距离库中心即雷达有效照射体积的中心。
(2)采用径向和方位上的最近邻居和垂直线性内插法球坐标系网格点的反射率因子值，具体为：
假设(r,a,e)是某一网格点在雷达球坐标系中的位置，r为斜距，a为方位角，e为仰角；e位于其上下相邻仰角e2和e1之间；
(r,a,e2)和(r,a,e1)分别是经过该网格点的垂线(仰角低于20°时，垂直方向可用仰角方向近似)与其上下仰角波束轴线的交点，那么：
该球坐标系网格点的反射率因子值fa(r,a,e)用(r,a,e2)点反射率因子值fa(r,a,e2)和(r,a,e1)点反射率因子值fa(r,a,e1)进行垂直线性内插得到，即：
fa(r,a,e)＝(we1fa(r,a,e1)+we2fa(r,a,e2))/(we1+we2) (公式6)
其中，
we1、we2分别为fa(r,a,e1)和fa(r,a,e2)的内插权重：
we1＝(e2-e)/(e2-e1) (公式7)
we2＝(e-e1)/(e2-e1) (公式8)
fa(r,a,e2)和fa(r,a,e1)分别等于最靠近点(r,a,e2)和(r,a,e1)的雷达距离库的观测值，即球坐标系下的雷达反射率因子值，采用径向和方位上的最近邻居法获取。
假设球坐标系中，ri-1、ri、ri+1为相邻径向距离库，aj-1、aj、aj+1为相邻方位角，由半功率线和半距离库所围成的梯形区是距离库ri的影响区；
在径向、方位方向上落在所述梯形区的点(r,a)的反射率因子值fa(r,a)用距离库ri的观测值fo(ri,aj)来赋值，即fa(r,a)＝fo(ri,aj)。
(3)采用垂直水平线性内插法球坐标系网格点的反射率因子值，具体为：
假设(r,a,e2)和(r,a,e1)分别是经过网格点(r,a,e)的垂线(仰角低于20°时，垂直方向可用仰角方向近似)与其上下仰角波束轴线的交点，(r1,a,e2)、(r2,a,e1)分别是经过该网格点的水平线与其相邻上下仰角波束轴线的交点，那么：
该网格点的反射率因子值fa(r,a,e)用反射率因子值fa(r,a,e2)、fa(r,a,e1)、fa(r1,a,e2)、fa(r2,a,e1)通过垂直和水平内插得到；该网格点的反射率因子值fa(r,a,e)用(r,a,e2)、(r,a,e1)、(r1,a,e2)、(r2,a,e1)点反射率因子值fa(r,a,e2)、fa(r,a,e1)、fa(r1,a,e2)、fa(r2,a,e1)通过垂直和水平内插得到；
其中，fa(r,a,e2)、fa(r,a,e1)、fa(r1,a,e2)、fa(r2,a,e1)通过径向和方位的最近邻居法得到，则：
Figure BDA0003053548560000111
其中，wr1、wr2分别是fa(r1,a,e2)、fa(r2,a,e1)的内插权重：
wr1＝(r2-r)/(r2-r1) (公式10)
wr2＝(r-r1)/(r2-r1) (公式11)
且有r1＝rsine/sine2，r2＝rsine/sine1。
(4)采用8点插值法球坐标系网格点的反射率因子值，具体为：
假设(r,a,e2)和(r,a,e1)分别是经过网格点(r,a,e)的垂线与其上下仰角波束轴线的交点，(r1,a,e2)、(r2,a,e1)分别是经过该网格点的水平线与其相邻上下仰角波束轴线的交点，那么：
若某一球坐标系网格点(r,a,e)落在由(r1,a1,e1)、(r2,a1,e1)、(r1,a2,e1)、(r2,a2,e1)、(r1,a1,e2)、(r2,a1,e2)、(r1,a2,e2)、(r2,a2,e2)点围成的锥体内，则该网格点的反射率因子值fa(r,a,e)由这8个点的观测值
Figure BDA0003053548560000112
Figure BDA0003053548560000113
Figure BDA0003053548560000114
进行双线性内插获得：
Figure BDA0003053548560000121
其中，wa1、wa2为方位内插权重：
wa1＝(a2-a)/(a2-a1) (公式13)
wa2＝(a-a1)/(a2-a1) (公式14)；
wr1是
Figure BDA0003053548560000122
的内插权重，wr2是
Figure BDA0003053548560000123
的内插权重；
Figure BDA0003053548560000124
Figure BDA0003053548560000125
的内插权重；
we2是
Figure BDA0003053548560000126
的内插权重。
步骤6，将覆盖笛卡尔坐标系网格点的雷达反射率分析值赋值给该笛卡尔坐标系网格点，得到该笛卡尔坐标系网格点的雷达反射率因子值，所有笛卡尔坐标系网格点完成赋值后，即完成本次三维拼图。
具体实施时，可采用如下方法实现笛卡尔坐标系网格点幅值：
(1)采用权重函数法计算笛卡尔坐标系网格点的雷达反射率因子值：
根据步骤5的反射率因子值计算笛卡尔坐标系每个网格点i的反射率因子值，实现目标区域范围内多部多普勒气象雷达数据拼接。
笛卡尔坐标系网格点i的雷达反射率因子值通过下面公式得到：
Figure BDA0003053548560000127
其中，
fm(i)是网格点i的合成反射率因子值；
Figure BDA0003053548560000128
是在网格点i处来自第n个雷达的反射率分析值；
wn是反射率分析值
Figure BDA0003053548560000129
的权重；
Nrad是在网格点i处有反射率分析值的总的雷达个数。
如果Nrad＝0，那么网格点不被任何一个雷达覆盖，该网格点的fm(i)被赋一个缺值符号；
如果Nrad＝1，那么网格点的fm(i)等于那个雷达在该网格点的反射率分析值；
如果Nrad>1，那么网格点的fm(i)使用多个雷达反射率分析值的权重平均。
即笛卡尔坐标系网格点不被任何一个雷达覆盖时，则赋一个缺值符号；
被一个雷达覆盖时，则赋该雷达在该网格点的反射率分析值；
被至少两个雷达覆盖时，则覆盖雷达反射率分析值按权重幅值。
所述权重可采用指数权重函数：
Figure BDA0003053548560000131
其中D为适当的长度比例，取D＝150，d为网格点到雷达的距离。
所述权重还可采用Cressman权重函数：
Figure BDA0003053548560000132
其中D为影响半径，取D＝300，d为网格点到雷达的距离。
(2)采用最近邻居法计算笛卡尔坐标系网格点的雷达反射率因子值，具体的：
将来自最靠近笛卡尔坐标系网格点的那个雷达的反射率分析值的权重赋为1，其它的权重全赋为0，即把最靠近笛卡尔坐标系网格点的那个雷达的反射率分析值赋给网格点。
(3)采用最大值法计算笛卡尔坐标系网格点的雷达反射率因子值，具体的：
将覆盖同一笛卡尔坐标系网格点的多个雷达反射率分析值中的最大值的权重赋为1，其它的权重全赋为0，即把覆盖同一网格点的多个雷达反射率分析值中的最大值赋给网格点。
本发明申请人结合说明书附图对本发明的实施示例做了详细的说明与描述，但是本领域技术人员应该理解，以上实施示例仅为本发明的优选实施方案，详尽的说明只是为了帮助读者更好地理解本发明精神，而并非对本发明保护范围的限制，相反，任何基于本发明的发明精神所作的任何改进或修饰都应当落在本发明的保护范围之内。

Claims (13)

1.一种多普勒气象雷达数据采集及三维拼图方法，其特征在于：
所述方法包括以下步骤：
步骤1，采集目标区域范围内球坐标系下多部多普勒气象雷达数据；
步骤2，根据目标区域定义笛卡尔坐标系网格范围和分辨率，在笛卡尔坐标系上划分若干网格点；
步骤3，根据笛卡尔坐标系网格点坐标，逐一检索出覆盖每一笛卡尔坐标系网格点的雷达及其雷达数据；
步骤4，根据笛卡尔坐标系网格点坐标计算对应每一雷达的球坐标系网格点坐标；
步骤5，插值计算每一雷达球坐标系网格点的反射率因子值，作为雷达在笛卡尔坐标系网格点的反射率分析值；
步骤6，将覆盖笛卡尔坐标系网格点的雷达反射率分析值赋值给该笛卡尔坐标系网格点，得到该笛卡尔坐标系网格点的雷达反射率因子值，所有笛卡尔坐标系网格点完成赋值后，即完成本次三维拼图。
2.根据权利要求1所述的一种多普勒气象雷达数据采集及三维拼图方法，其特征在于：
步骤2中，所述笛卡尔坐标系网格在垂直方向选取的范围为0-17km，共21层；
5km以下的垂直分辨率为0.5km，5-17km之间的垂直分辨率为1km；
在水平方向的经纬度分辨率为0.01°×0.01°，范围是覆盖目标区域。
3.根据权利要求1所述的一种多普勒气象雷达数据采集及三维拼图方法，其特征在于：
步骤4中，利用笛卡尔坐标系下网格点的经度、纬度和高度计算其对应雷达在球坐标系中的仰角e、方位角a和斜距r；
设笛卡尔坐标系三维网格中任意网格点的坐标为(αg，βg，hg)，其中αg为网格点纬度，βg为网格点经度，hg为网格点高度；
雷达天线所在点的坐标为(αr，βr，hr)，其中αr为雷达天线纬度，βr为雷达天线经度，hr为雷达天线高度；
则使用雷达波束传播和大圆几何学理论确定笛卡尔坐标系三维网格中网格点对应雷达天线所在点的球坐标位置(r，a，e)，其中r为斜距，a为方位角，e为仰角，具体的：
由球面三角公式得出：
sina＝cos(αg)sin(βg-βr)/sin(s/R) (公式1)
其中R为地球半径，s为大圆距离，其表达式为：
s＝Rcos-1(sin(αr)sin(αg)+cos(αr)cos(αg)cos(βg-βr)) (公式2)
设C＝sina，则：
Figure RE-FDA0003133486040000021
Figure RE-FDA0003133486040000022
r＝sin(s/Rm)(Rm+hg-hr)/cos(e) (公式5)
其中Rm为等效地球半径，
Figure RE-FDA0003133486040000023
4.根据权利要求1所述的一种多普勒气象雷达数据采集及三维拼图方法，其特征在于：
步骤5中，采用最近邻居法计算球坐标系网格点的反射率因子值，具体为：
计算球坐标系中的网格点的中心与雷达距离库中心之间的距离，用最靠近网格点的雷达距离库的观测值，即球坐标系下的雷达反射率因子值去填充球坐标系网格点的雷达反射率因子值，作为雷达在笛卡尔坐标系网格点的反射率分析值；
所述雷达距离库中心即雷达有效照射体积的中心。
5.根据权利要求1所述的一种多普勒气象雷达数据采集及三维拼图方法，其特征在于：
步骤5中，采用径向和方位上的最近邻居和垂直线性内插法计算球坐标系网格点的反射率因子值，具体为：
假设(r,a,e)是某一网格点在雷达球坐标系中的位置，r为斜距，a为方位角，e为仰角；e位于其上下相邻仰角e2和e1之间；
(r,a,e2)和(r,a,e1)分别是经过该网格点的垂线与其上下仰角波束轴线的交点，那么：
该球坐标系网格点的反射率因子值fa(r,a,e)用(r,a,e2)点反射率因子值fa(r,a,e2)和(r,a,e1)点反射率因子值fa(r,a,e1)进行垂直线性内插得到，即：
fa(r,a,e)＝(we1fa(r,a,e1)+we2fa(r,a,e2))/(we1+we2) (公式6)
其中，
we1、we2分别为fa(r,a,e1)和fa(r,a,e2)的内插权重：
we1＝(e2-e)/(e2-e1) (公式7)
we2＝(e-e1)/(e2-e1) (公式8)
fa(r,a,e2)和fa(r,a,e1)分别等于最靠近点(r,a,e2)和(r,a,e1)的雷达距离库的观测值，即球坐标系下的雷达反射率因子值，采用径向和方位上的最近邻居法获取。
6.根据权利要求5所述的一种多普勒气象雷达数据采集及三维拼图方法，其特征在于：
假设球坐标系中，ri-1、ri、ri+1为相邻径向距离库，aj-1、aj、aj+1为相邻方位角，由半功率线和半距离库所围成的梯形区是距离库ri的影响区；
在径向、方位方向上落在所述梯形区的点(r,a)的反射率因子值fa(r,a)用距离库ri的观测值fo(ri,aj)来赋值，即fa(r,a)＝fo(ri,aj)。
7.根据权利要求1所述的一种多普勒气象雷达数据采集及三维拼图方法，其特征在于：
步骤5中，采用垂直水平线性内插法球坐标系网格点的反射率因子值，具体为：
假设(r,a,e2)和(r,a,e1)分别是经过网格点(r,a,e)的垂线与其上下仰角波束轴线的交点，(r1,a,e2)、(r2,a,e1)分别是经过该网格点的水平线与其相邻上下仰角波束轴线的交点，那么：
该网格点的反射率因子值fa(r,a,e)用(r,a,e2)、(r,a,e1)、(r1,a,e2)、(r2,a,e1)点反射率因子值fa(r,a,e2)、fa(r,a,e1)、fa(r1,a,e2)、fa(r2,a,e1)通过垂直和水平内插得到；
其中，fa(r,a,e2)、fa(r,a,e1)、fa(r1,a,e2)、fa(r2,a,e1)通过径向和方位的最近邻居法得到，则：
Figure RE-FDA0003133486040000041
其中，wr1、wr2分别是fa(r1,a,e2)、fa(r2,a,e1)的内插权重：
wr1＝(r2-r)/(r2-r1) (公式10)
wr2＝(r-r1)/(r2-r1) (公式11)
且有r1＝rsine/sine2，r2＝rsine/sine1。
8.根据权利要求1所述的一种多普勒气象雷达数据采集及三维拼图方法，其特征在于：
步骤5中，采用8点插值法球坐标系网格点的反射率因子值，具体为：
假设(r,a,e2)和(r,a,e1)分别是经过网格点(r,a,e)的垂线与其上下仰角波束轴线的交点，(r1,a,e2)、(r2,a,e1)分别是经过该网格点的水平线与其相邻上下仰角波束轴线的交点，那么：
若某一球坐标系网格点(r,a,e)落在由(r1,a1,e1)、(r2,a1,e1)、(r1,a2,e1)、(r2,a2,e1)、(r1,a1,e2)、(r2,a1,e2)、(r1,a2,e2)、(r2,a2,e2)点围成的锥体内，则该网格点的反射率因子值fa(r,a,e)由这8个点的观测值f1 o(r1,a1,e1)、
Figure RE-FDA0003133486040000042
Figure RE-FDA0003133486040000051
进行双线性内插获得：
Figure RE-FDA0003133486040000052
其中，wa1、wa2为方位内插权重：
wa1＝(a2-a)/(a2-a1) (公式13)
wa2＝(a-a1)/(a2-a1) (公式14)；
wr1是f1 o、
Figure RE-FDA0003133486040000053
的内插权重，wr2是
Figure RE-FDA0003133486040000054
的内插权重；
we1是
Figure RE-FDA0003133486040000055
的内插权重；
we2是
Figure RE-FDA0003133486040000056
的内插权重。
9.根据权利要求1所述的一种多普勒气象雷达数据采集及三维拼图方法，其特征在于：
步骤6中，根据步骤5的反射率因子值计算笛卡尔坐标系每个网格点i的反射率因子值，实现目标区域范围内多部多普勒气象雷达数据拼接；
笛卡尔坐标系网格点i的雷达反射率因子值通过下面公式得到：
Figure RE-FDA0003133486040000057
其中，
fm(i)是网格点i的合成反射率因子值；
Figure RE-FDA0003133486040000058
是在网格点i处来自第n个雷达的反射率分析值；
wn是反射率分析值
Figure RE-FDA0003133486040000059
的权重；
Nrad是在网格点i处有反射率分析值的总的雷达个数。
10.根据权利要求9所述的一种多普勒气象雷达数据采集及三维拼图方法，其特征在于：
所述权重采用指数权重函数：
Figure RE-FDA0003133486040000061
其中D为适当的长度比例，取D＝150，d为网格点到雷达的距离。
11.根据权利要求9所述的一种多普勒气象雷达数据采集及三维拼图方法，其特征在于：
所述权重采用Cressman权重函数：
Figure RE-FDA0003133486040000062
其中D为影响半径，取D＝300，d为网格点到雷达的距离。
12.根据权利要求1所述的一种多普勒气象雷达数据采集及三维拼图方法，其特征在于：
步骤6中，采用最近邻居法计算笛卡尔坐标系网格点的雷达反射率因子值，具体的：
将来自最靠近笛卡尔坐标系网格点的那个雷达的反射率分析值的权重赋为1，其它的权重全赋为0，即把最靠近笛卡尔坐标系网格点的那个雷达的反射率分析值赋给网格点。
13.根据权利要求1所述的一种多普勒气象雷达数据采集及三维拼图方法，其特征在于：
步骤6中，采用最大值法计算笛卡尔坐标系网格点的雷达反射率因子值，具体的：
将覆盖同一笛卡尔坐标系网格点的多个雷达反射率分析值中的最大值的权重赋为1，其它的权重全赋为0，即把覆盖同一网格点的多个雷达反射率分析值中的最大值赋给网格点。
CN202110493990.7A 2021-05-07 2021-05-07 一种多普勒气象雷达数据采集及三维拼图方法 Pending CN113325423A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110493990.7A CN113325423A (zh) 2021-05-07 2021-05-07 一种多普勒气象雷达数据采集及三维拼图方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110493990.7A CN113325423A (zh) 2021-05-07 2021-05-07 一种多普勒气象雷达数据采集及三维拼图方法

Publications (1)

Publication Number Publication Date
CN113325423A true CN113325423A (zh) 2021-08-31

Family

ID=77414298

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110493990.7A Pending CN113325423A (zh) 2021-05-07 2021-05-07 一种多普勒气象雷达数据采集及三维拼图方法

Country Status (1)

Country Link
CN (1) CN113325423A (zh)

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102117227A (zh) * 2011-03-09 2011-07-06 南京恩瑞特实业有限公司 天气雷达数据的多核并行计算方法
CN108445464A (zh) * 2018-03-12 2018-08-24 南京恩瑞特实业有限公司 Nriet基于机器学习的卫星雷达反演融合方法
CN109254290A (zh) * 2018-08-17 2019-01-22 深圳市雅码科技有限公司 一种天气雷达并行拼图方法及系统
CN109633783A (zh) * 2018-11-14 2019-04-16 中国气象科学研究院 基于双偏振雷达网的参量拼图方法及装置
US20210063568A1 (en) * 2019-09-03 2021-03-04 Korea Meteorological Administration Apparatus and method for composition for dual-polarization weather radar observation data using earth spherical coordinate system
CN112731564A (zh) * 2020-12-26 2021-04-30 安徽省公共气象服务中心 一种基于多普勒天气雷达数据的雷电智能预报方法

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102117227A (zh) * 2011-03-09 2011-07-06 南京恩瑞特实业有限公司 天气雷达数据的多核并行计算方法
CN108445464A (zh) * 2018-03-12 2018-08-24 南京恩瑞特实业有限公司 Nriet基于机器学习的卫星雷达反演融合方法
CN109254290A (zh) * 2018-08-17 2019-01-22 深圳市雅码科技有限公司 一种天气雷达并行拼图方法及系统
CN109633783A (zh) * 2018-11-14 2019-04-16 中国气象科学研究院 基于双偏振雷达网的参量拼图方法及装置
US20210063568A1 (en) * 2019-09-03 2021-03-04 Korea Meteorological Administration Apparatus and method for composition for dual-polarization weather radar observation data using earth spherical coordinate system
CN112731564A (zh) * 2020-12-26 2021-04-30 安徽省公共气象服务中心 一种基于多普勒天气雷达数据的雷电智能预报方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
肖艳姣, 刘黎平: "新一代天气雷达网资料的三维格点化及拼图方法研究", 《气象学报》, vol. 64, no. 5, pages 647 - 657 *
肖艳姣;刘黎平;: "新一代天气雷达网资料的三维格点化及拼图方法研究", 气象学报, no. 05, pages 647 - 657 *

Similar Documents

Publication Publication Date Title
Mukupa et al. A review of the use of terrestrial laser scanning application for change detection and deformation monitoring of structures
US10915673B2 (en) Device, method, apparatus, and computer-readable medium for solar site assessment
Gál et al. Computing continuous sky view factors using 3D urban raster and vector databases: comparison and application to urban climate
CN112558076B (zh) 基于组网天气雷达覆盖域的体扫描模式计算方法与应用
CN107356926A (zh) 基于Hu矩的差值云团外推降雨预测算法
CN113325422B (zh) 天基测雨雷达目标定位及降雨信息三维处理方法和系统
Rodler et al. Local climate zone approach on local and micro scales: Dividing the urban open space
Hinks et al. Flight optimization algorithms for aerial LiDAR capture for urban infrastructure model generation
CN116069882B (zh) 空域网格图生成方法
CN114091756B (zh) 一种基于泰森多边形的乡镇级海啸风险评估方法
US20140266856A1 (en) System and method for filling gaps in radar coverage
Dorman et al. shadow: R Package for Geometric Shadow Calculations in an Urban Environment.
CN109001846A (zh) 一种复杂地形条件下s波段与x波段雷达组网测雨方法
CN114660591A (zh) 一种基于多部天气雷达三维组网生成方法
CN113325423A (zh) 一种多普勒气象雷达数据采集及三维拼图方法
CN107341568B (zh) 一种台风风暴增水预测方法及系统
CN113156441B (zh) 气象雷达探测的有效立体空域细分逼近计算方法
CN114912171A (zh) 一种基于建筑遮挡的日照辐射计算方法
CN114137637A (zh) 基于闪电和雷达数据的雷暴中心踪迹的集合概率预报方法
Sládek Applicability of Python-based geospatial analysis for operational assessment of weather station deployment
Viana-Fons et al. Methodology for the calculation of the shadow factor on roofs and facades of buildings in urban areas
Yu et al. Evaluation of detection ability of regional radar network in complex terrain based on digital elevation model: a sample study of Fujian province in China
Janowski et al. 3D modelling of liquid fuels base infrastructure for the purpose of visualization and geometrical analysis
Liao et al. Fast and accurate estimation of solar irradiation on building rooftops in Hong Kong: A machine learning-based parameterization approach
Zhang et al. Three-and four-dimensional high-resolution national radar mosaic

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination