CN109283505B - 一种修正雷达回波外推图像发散现象的方法 - Google Patents
一种修正雷达回波外推图像发散现象的方法 Download PDFInfo
- Publication number
- CN109283505B CN109283505B CN201811019917.0A CN201811019917A CN109283505B CN 109283505 B CN109283505 B CN 109283505B CN 201811019917 A CN201811019917 A CN 201811019917A CN 109283505 B CN109283505 B CN 109283505B
- Authority
- CN
- China
- Prior art keywords
- rimg
- image
- radar
- echo
- radar echo
- 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.)
- Active
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
- G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
- G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
- G01S7/41—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00 using analysis of echo signal for target characterisation; Target signature; Target cross-section
-
- 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)
- Computer Networks & Wireless Communication (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明提出一种修正雷达回波外推图像发散现象的方法，以改善雷达回波外推图像因发散而影响降水量等预报准确性的问题，由于数值模式及其产品能够较好地反映大气运动的客观规律，特别是对风暴发生、发展的变化趋势具有较好地整体性和连续性，因此本发明对数值模式产品这些特性加以利用。为了克服雷达外推预报出的回波强度“过度增加”或“过度减小”的问题，本发明通过融合数值模式预报产品，对回波图像在起报时刻与预报时刻的变化幅度进行修订，进而提高了雷达降水预报的准确性。
Description
技术领域
本发明属于地球科学领域，针对雷达回波图像外推预报随预报时效的延长逐渐发散进而影响降水等预报准确率的问题提出改进的方法。
背景技术
利用雷达回波图像进行外推预报是短时临近天气预报预警的关键技术之一，在气象领域得到广泛、深入的研究与应用。该技术起步于上世纪60年代，用于监测强对流天气的发生、发展，以保障机场空管和飞行的安全。该技术发展至今，得益于雷达探测技术的发展和计算机运算性能的急速提升，雷达探测的空间分辨率越来越高，反演出的产品越来越多，进而推动了雷达外推预测准确性的不断提升，并且预报的时效性也越来越强。
雷达回波外推预报的基本思想是利用雷达探测到的各仰角面上的基本反射率信息，对回波所表征的风暴体进行识别，然后根据相邻时刻雷达资料识别出的风暴体形态变化和移动趋势，对风暴体未来发展的走向、强度做出预测。在此基础上，利用Z-R关系等气象学方法将预测出的基本反射率数值换算成降水量等信息，进而达到预报的目的。
一般情况下，外推预报是建立在线性力学和运动关系的基础上，利用已知的实况雷达回波图像和计算出的风暴体移动速度矢量，根据预报的时长进行外推。对于风暴主体的移动趋势，这类方法具有良好的预报效果，但随着预报时效的延长，外推出的回波图像会越发显现出发散的现象，如图2和图3所示。其中，图2是一幅外推预报的回波图像，图3是与图2同一时刻、同一高度层的实况雷达回波图像。这两幅图像对比可以看出，图2所显示的回波区域中存在若干不规则的白点噪声，即所谓的发散现象。这主要是因为如下一些原因造成：
1)外推预报的回波图像，其每一个像素的值都取决于实况雷达回波图像中的某一或某些像素值以及计算出的风暴体移动速度矢量VECTOR，不同坐标位置的像素与VECTOR经计算后可能会重叠到预报回波图像的同一像素上。很显然，这种位置重叠发生的次数越多，预报回波图像上缺失的像素值也就越多，进而在图像上表现出一定程度地发散。
2)由于外推预报是建立在线性力学和运动关系的基础上，而实际大气和大气中粒子的运动有着复杂的动力学和热力学关系，其运动轨迹并不符合单纯线性运动的规律。因此，对于仅利用实况雷达资料进行外推预报的计算方法而言，很难克服外推预报出的回波强度“过度增加”或“过度减小”的问题。为此，通常的方法是设定一些阈值参数，来限制这种“过度”现象的发生。举例来说，将外推结果为负值的回波强度值修正为0，并将外推结果过大的回波强度值根据阈值参数设定做适当修正，如此也会造成外推预报出的回波图像呈现发散和回波强度“失真”的现象。
尽管这种图像发散的问题，可通过插值计算等数学方法在一定程度上解决，但很显然，这种“失真”的回波图像与实际大气状况并不相符，并且，利用这种回波图像换算出的降水量与实际也存在较大的偏差。
需要补充说明的是，雷达回波图像上的像素值与雷达探测到的基本反射率之间存在数学上的对应关系，并且这种对应关系对于本领域的专业人员而言是公知的。上述发散现象是从雷达回波图像的视觉层面进行描述，对于雷达基本反射率数值而言，就是存在很多本应有数值的点或小区域，当前的数值却为0或缺测。不论采用基本反射率数值或雷达回波图像上的像素值进行外推计算，没有本质区别，为了便于表述，本发明以下内容统一采用后者，即回波图像中的像素值作为计算的对象。
发明内容
为了改善雷达回波外推图像因发散而影响降水量等预报准确性的问题，本发明提出一种利用数值模式产品修正雷达回波外推图像的方法。
本发明所采用的技术方案如下：
一种修正雷达回波外推图像发散现象的方法，包括如下步骤：
步骤1：取某起报时刻雷达探测资料并采用雷达回波外推预报算法，计算未来某一时刻 t0+Δt各高度层的雷达基本反射率数据，再将雷达基本反射率数据绘制成对应各高度层的雷达回波外推图像，记为RImg(t0+Δt,h,x,y)，简称RImg；其中，t0表示起报时间；Δt表示外推预报的时效，0<Δt<3小时；h表示高度层的高度，0<h<20km；x和y分别表示雷达回波外推图像上任一像素的横坐标和纵坐标；
步骤2：取与步骤1雷达探测资料起报时刻t0时间上最相近的数值预报产品，从数值预报产品中提取与t0+Δt时刻最相近的降水要素并计算出与步骤1雷达回波外推图像各个相同高度层对应的基本反射率数据；
步骤3：将步骤2计算出的各个高度层的基本反射率数据绘制成回波图像记为MImg(t0+Δt,h,x’,y’)，简称MImg；其中，x’和y’分别表示该回波图像上任一像素的横坐标和纵坐标；再利用空间插值算法，将MImg(t0+Δt,h,x’,y’)转为与RImg(t0+Δt,h,x,y)具有相同地理范围、相同图像大小和相同坐标系的回波图像MImg(t0+Δt,h,x,y)；
步骤4：对雷达回波外推图像RImg(t0+Δt,h,x,y)中的发散点进行识别，判断是否存在下式情况：
其中，RImg＝255表示RImg中某一坐标(x,y)处的像素值为255；Count(RImg,r)表示RImg 中以r为半径，(x,y)为圆心的区域内像素值小于255的像素的数量；Sum(RImg,r)表示RImg中以r为半径，(x,y)为圆心的区域内像素的总数量；r取[2,5]的正整数；Ta为经验阈值；
步骤5：如果式1情况存在，即雷达回波外推图像RImg中(x,y)处为发散点，则对该点的像素值进行修补，修补过程如下：
首先，定义一个数据集DSR(x,y,r)，用于记录RImg中以r为半径，以坐标(x,y)为圆心的区域内所有像素值小于255的坐标；
然后，对数据集DSR(x,y,r)中的每一个坐标，从MImg图像中提取相应坐标位置的像素值，并对这些像素值按其数值大小进行排序，形成有序序列MList(i,x,y)；其中，i为该序列中各项的序号；提取MImg图像中坐标(x,y)处的像素值，再从MList(i,x,y)中找到与该像素值差值绝对值最小的一项的序号，并将该序号记为Pos(x,y)；
接着，对数据集DSR(x,y,r)中的每一个坐标，从RImg图像中提取相应坐标位置的像素值，并对这些像素值按其数值大小进行排序，形成有序序列RList(i,x,y)，i为有序序列各项的序号；
最后，按下式计算RImg中(x,y)处的像素值：
其中，RList(i,x,y)[Pos(x,y)-1]表示有序序列RList(i,x,y)中Pos(x,y)-1项的像素值；类似的， RList(i,x,y)[Pos(x,y)+1]表示有序序列RList(i,x,y)中Pos(x,y)+1项的像素值；
步骤6：对雷达回波外推图像RImg中的每个像素点，逐层根据步骤4的方法识别是否为发散点；如果为发散点，则按照步骤5的方法对发散点的像素值逐个进行修补计算，当所有高度层所有发散点全部修补后，计算结果对应存储到RImg中替换修补前的像素值，如果不是发散点，不进行替换；由此得到修补后的雷达回波外推图像RImg’(t0+Δt,h,x,y)，简称RImg’；
步骤7：将步骤1起报时刻雷达探测资料绘制成雷达回波图像，记为RImg(t0,h,x,y)；将上述与起报时刻t0时间上最相近的数值预报产品按步骤2和步骤3的方法绘制成回波图像，记为MImg(t0,h,x,y)；计算雷达回波图像RImg(t0,h,x,y)与RImg’(t0+Δt,h,x,y)各像素的变化值，以及MImg(t0,h,x,y)与MImg(t0+Δt,h,x,y)各像素的变化值，计算方法为：
RDif(x,y)＝RImg(t0,h,x,y)-RImg’(t0+Δt,h,x,y)
MDif(x,y)＝MImg(t0,h,x,y)-MImg(t0+Δt,h,x,y)
其中，RDif(x,y)是一个数据集，记录了雷达回波图像上的每一个像素从t0时刻到t0+Δt 时刻像素值的差值；类似的，MDif(x,y)也是一个数据集，记录了数值预报产品回波图像上的每一个像素从t0时刻到t0+Δt时刻像素值的差值；
步骤8：计算数据集RDif(x,y)的中值，记为RMV，计算方法是将RDif(x,y)中各项按其数值大小排序后，提取处于该序列中间位置的值；类似地，计算数据集MDif(x,y)的中值，记为 MMV；
步骤9：按下式计算MMV与RMV的倍率MP：MP＝MMV/RMV；
步骤10：利用倍率MP对修补后的雷达回波外推图像RImg’的像素值按下式逐层逐一进行修正：
RImg_New(t0+Δt,h,x,y)＝RImg’(t0+Δt,h,x,y)+Tc×MP×RDif(x,y)
其中，修正后的雷达回波外推图像RImg_New(t0+Δt,h,x,y)的物理意义表示：以t0时刻为起报时间，外推预测Δt时刻，h高度上，任一坐标(x,y)处的像素值，简称RImg_New；Tc为经验系数。
本发明进一步设计在于：
步骤1中雷达回波外推预报算法采用最大相关法或交叉相关法。
步骤2中数值预报产品中降水要素包括湿度、降水量、气压、温度、水凝物和冰粒子气象信息。
步骤3中空间插值算法采用反距离权重算法。
步骤5中当有序序列MList中存在数值大小相同的像素值时，这些像素值之间的序列顺序可随机排序；当MList(i,x,y)中与MImg(x,y)像素差值绝对值最小存在多项时，多项中可任意选择其中一项的序号记为Pos(x,y)。
步骤5计算过程中，当Pos(x,y)＝0时，则RImg＝RList(i,x,y)[Pos(x,y)+1]；当Pos(x,y)＝Count(MImg,r)时，则RImg＝RList(i,x,y)[Pos(x,y)-1]。
对步骤10得到的修正后的雷达回波外推图像RImg_New(t0+Δt,h,x,y)，采用与基本反射率数据转化为回波图像像素值相反的计算过程，将RImg_New(t0+Δt,h,x,y)中各个像素的像素值转化为基本反射率数据，再通过雷达气象学中的Z-R关系理论方法，将基本反射率计算转化为降水方面的预报信息。
本发明相比现有技术具有如下有益效果：
1、本发明改善了雷达回波图像随外推时效的增加越发呈现发散的现象，为后续进行降水量等预报提供更加准确的基本反射率数据。
2、不同于一般的空间插值方法，本发明考虑到数值模式预报产品能够较好地反映大气运动的客观规律，特别是对风暴发生、发展的变化趋势具有较好地整体性和连续性，但监测和预报的实时性不如雷达外推预报等特点，融合数值模式预报产品对雷达回波外推图像中的发散点和发散小区域进行修补，降低了传统插值方法对回波图像造成的“失真”。
3、为了克服雷达外推预报出的回波强度“过度增加”或“过度减小”的问题，本发明通过融合数值模式预报产品，对回波图像在起报时刻与预报时刻的变化幅度进行修订，进而提高了雷达降水预报的准确性。
附图说明
图1是本发明主要算法过程的流程图。
图2是一幅未经本发明技术方案改进的外推预报的回波图像。
图3是与图2同一时刻、同一高度层的实况雷达回波图像，即由实际探测到的雷达基本反射率数据绘制出的回波图像。
图4是应用实例一修正后的雷达回波外推图像。
具体实施方式
下面结合附图对本发明的技术方案进一步说明。
实施例一：
如图1所示，本发明的一种修正雷达回波外推图像发散现象的方法，其具体步骤如下：
步骤1：取某起报时刻雷达探测资料并采用本领域已知的雷达回波外推预报算法，如最大相关法或交叉相关法，计算出未来某一时刻的各高度层的雷达基本反射率数据，然后采用现有的算法，将基本反射率数据制成对应各高度层的雷达回波外推图像，记为RImg(t0+Δt,h,x,y)，简称RImg。其中，t0表示起报时间，Δt表示外推预报的时效，0<Δt<3小时，h表示高度层的高度，0<h<20km，x和y分别表示雷达回波外推图像上任一像素的横坐标和纵坐标，每个(x,y) 对应一个地理位置坐标。RImg(t0+Δt,h,x,y)的物理意义为：以t0时刻为起报时间，外推预测的预测时效为Δt，h高度上，(x,y)坐标处的像素值。对于灰阶图像而言，一般情况下，RImg的值域为[0,255]，且值越大表示回波强度越弱，当RImg＝255时，表示当前位置的基本反射率数据小于或等于0dBZ。
步骤2：取与步骤1雷达探测资料起报时刻t0时间上最相近的数值预报产品，从数值预报产品中提取与t0+Δt时刻最相近的降水要素并计算出与步骤1雷达回波外推图像各个相同高度层对应的基本反射率数据；
步骤2：取与步骤1雷达探测资料起报时刻t0时间上最相近(与t0最相近时刻，视为t0) 的数值预报产品，从该数值预报产品中提取与t0+Δt时刻最相近(视为t0+Δt)降水要素(包括湿度、降水量、气压、温度、水凝物以及冰粒子等信息)并换算成与步骤1雷达回波外推图像各个相同高度层对应的基本反射率数据。(本发明中数值预报产品涉及2个与时间有关的概念，一个是数值预报产品的起报时间T1，另一个是数值预报产品的预报时效T2；每个数值预报产品只有1个起报时间，并且有多个预报时效。该步骤中数值预报产品的起报时间T1＝t0，预报时效T2＝Δt；该步骤中t0+Δt对应的是起报时间为t0的数值模式产品中，预报时效为Δt的气象数据)
步骤3：将步骤2计算出的各个高度层的基本反射率数据采用现有算法绘制成回波图像，记为MImg(t0+Δt,h,x’,y’)，简称MImg。其中，x’和y’分别表示由数值预报产品绘制出的回波图像上任一像素的横坐标和纵坐标。补充说明：一般情况下，数值预报产品覆盖的空间地理范围远大于单部雷达的空间地理范围，但数值模式的空间分辨率往往低于雷达的空间分辨率。再利用空间插值算法，将上述两种资料的空间分辨率处理成同一标准，即将MImg(t0+Δt,h,x’,y’) 转化为MImg(t0+Δt,h,x,y)。例如，假设RImg的空间分辨率为1km/像素，而MImg的空间分辨率为2.5km/像素，那么可采用反距离权重算法，将MImg的空间分辨率转换为1km/像素。作为本发明的约束，假定数值预报产品覆盖的空间地理范围大于单部雷达的空间地理范围，那么 RImg中的任一像素(x,y)都能从MImg中找到与(x,y)地理位置相一致的像素(x’,y’)，为了便于表述，本发明中RImg与MImg具有相同的地理范围、相同的图像大小和相同的坐标系，下述步骤中，RImg中的(x,y)与MImg中的(x,y)表示同一地理位置。
步骤4：对雷达回波外推图像RImg(t0+Δt,h,x,y)中的发散点进行识别，识别的方法为：
其中，RImg＝255表示RImg中某一坐标(x,y)处的像素值为255(像素值和基本反射率数据存在对应关系，可进行换算)；Count(RImg,r)表示RImg中以r为半径，(x,y)为圆心的区域内像素值小于255的像素的数量；Sum(RImg,r)表示RImg中以r为半径，(x,y)为圆心的区域内像素的总数量；r为大于1的正整数，一般情况下，r的建议值域为[2,5]；Ta是一个经验阈值， Ta∈[0,1]，Ta的物理意义为：在以r为半径，(x,y)为圆心的区域内，像素值小于255的像素数量占该区域内全部像素数量的比例。式1的物理意义为：当某一坐标(x,y)处的像素值为255，且在以r为半径，(x,y)为圆心的区域内，像素值小于255的像素数量与该区域内全部像素数量的比大于Ta时，则判定当前坐标点(x,y)为发散点。
步骤5：如果式1情况存在，即雷达回波外推图像RImg中(x,y)处为发散点，则对该点的像素值进行修补，修补的方法为：
首先，定义一个数据集DSR(x,y,r)，用于记录RImg中以r为半径，以坐标(x,y)为圆心的区域内所有像素值小于255的坐标。
然后，对数据集DSR(x,y,r)的每一个坐标，从MImg图像中提取相应坐标位置的像素值，并对这些像素值按其数值大小进行排序，形成有序序列MList(i，x,y)；其中，i为该序列中各项的序号。提取MImg图像中对应坐标(x,y)处的像素值，再从MList(i，x,y)中找到与该像素值差值绝对值最小的一项的序号，并将该最小项的序号记为Pos(x,y)。当有序序列MList中存在数值大小相同的像素值时，这些像素值之间的序列顺序可随机排序；当MList(i，x,y)中与 MImg(x,y)像素差值绝对值最小存在多项时，多项中可任意选择其中一项的序号记为Pos(x,y)。
接着，对数据集DSR(x,y,r)中的每一个坐标，从RImg图像中提取相应坐标位置的像素值，并对这些像素值按其数值大小进行排序，与MList(i，x,y)排序方式相同，形成有序序列RList(i， x,y)，i为有序序列各项的序号。
最后，计算RImg中(x,y)处的像素值，计算方法为：
其中，RList(i,x,y)[Pos(x,y)-1]表示有序序列RList(i,x,y)中Pos(x,y)-1项的像素值；类似的， RList(i,x,y)[Pos(x,y)+1]表示有序序列RList(i,x,y)中Pos(x,y)+1项的像素值。计算过程中，存在一些特殊情况，可另外处理，如：当Pos(x,y)＝0时，RImg＝RList(i,x,y)[Pos(x,y)+1]；当 Pos(x,y)＝Count(MImg,r)时，RImg＝RList(i,x,y)[Pos(x,y)-1]。
步骤6：对雷达回波外推图像RImg中的每个像素点，逐层根据步骤4进行判别是否为发散点；如果为发散点，则按照步骤5的方法对发散点的像素值进行逐个修补计算，当所有高度层所有发散点全部修补后，计算结果对应存储到RImg中替换修补前的像素值，如果不是发散点，则不进行替换；由此得到修补后的雷达回波外推图像，记为RImg’(t0+Δt,h,x,y)，简称RImg’。
由于大气中风暴体的生消、演变过程十分复杂，而仅依靠雷达探测资料进行外推预测的算法难以表征这一复杂的过程，因此，再利用数值预报产品换算出的基本反射率及其回波图像，对由步骤6修补后的回波图像RImg’进行调整，计算过程按步骤7至步骤10所述。
步骤7：将起报时刻的雷达探测资料采用现有算法绘制成雷达回波图像，记为RImg(t0,h,x,y)，再将上述与起报时刻t0时间上最相近的数值预报产品，按步骤2和步骤3的方法绘制回波图像，记为MImg(t0,h,x,y)。计算雷达回波图像RImg(t0,h,x,y)与RImg’(t0+Δt,h,x,y) 各像素的变化值，以及MImg(t0,h,x,y)与MImg(t0+Δt,h,x,y)各像素的变化值，计算方法为：
RDif(x,y)＝RImg(t0,h,x,y)-RImg’(t0+Δt,h,x,y)
MDif(x,y)＝MImg(t0,h,x,y)-MImg(t0+Δt,h,x,y)
其中，RDif(x,y)是一个数据集，记录了雷达回波图像上的每一个像素从t0时刻到t0+Δt 时刻，像素值的数学差值；类似的，MDif(x,y)也是一个数据集，记录了数值预报产品输出的回波图像上的每一个像素从t0时刻到t0+Δt时刻，像素值的数学差值。
步骤8：计算数据集RDif(x,y)的中值，即中位数，记为RMV，计算方法是将记录于RDif(x,y) 中的各项按值大小排序后，提取处于该序列中间位置的值。该值反映了雷达回波图像从t0时刻到t0+Δt时刻演变的总体水平。类似地，计算数据集MDif(x,y)的中值，记为MMV。
步骤9：计算MMV与RMV的倍率MP，计算方法为：
MP＝MMV/RMV
步骤10：利用MP对修补后的雷达回波外推图像称RImg’的像素值逐层逐一进行修正，修正方法为：
RImg_New(t0+Δt,h,x,y)＝RImg’(t0+Δt,h,x,y)+Tc×MP×RDif(x,y)
其中，修正后的雷达回波外推图像称RImg_New(t0+Δt,h,x,y)的其物理意义表示：以t0时刻为起报时间，外推预测Δt时刻，h高度上，任一坐标(x,y)处的像素值，该值是通过本发明技术方案修订后的像素值，简称RImg_New。Tc是一个经验系数，用于调整修正的幅度，Tc∈[-1,1]，具体取值需根据季节、降水类型并结合以往外推预报的准确率做出适度调整。
实施例二：
本例在实施例一基础上，进一步对步骤10修正后的各层雷达回波外推图像RImg_New，采用与基本反射率数据转化为回波图像像素值相反的计算过程，将RImg_New中各个像素的像素值转化为基本反射率数据，再通过雷达气象学中的Z-R关系理论方法，将基本反射率计算转化为可降水量预报信息，为相关领域的专业人员提供辅助决策。
应用实例一：
本应用实例选用2018年5月18日16至18时南京一部天气雷达部分时刻的探测资料作为实施对象。
步骤1：采用本领域公知的交叉相关法，对当日08:33、08:39和08:45(均为世界时，与北京时间相差8小时，起报时刻t0为08:45)时刻的雷达探测资料进行预报时效为1小时(Δt 为1小时)的外推计算，得到9:45(世界时)高度层为1250米的雷达基本反射率数据，再将雷达基本反射率数据绘制成该高度层的雷达回波外推图像，记为RImg(t0+Δt,h,x,y)，简称RImg，该图像如图2所示。其中，t0表示起报时间，值为“2018-05-18 08:45(世界时)”；Δt表示外推预报的时效，值为“60(单位：分钟)”；h表示高度层，值为“1250(单位：米)”，x和y 分别表示回波图像上任一像素的横坐标和纵坐标。RImg(t0+Δt,h,x,y)的物理意义为：以t0 时刻为起报时间，外推预测Δt时刻，h高度上，(x,y)坐标处的像素值。本实施例处理的所有图像均为256阶的灰阶图像，因此，RImg的值域为[0,255]，且值越大表示回波强度越弱，当 RImg＝255时，表示当前位置的基本反射率数据小于或等于0dBZ。
步骤2：当前，气象相关领域的数值模式有很多种，本实施例选用GFS(GlobalForecast System全球预报***)数值预报产品。选择时间上与t0最相近的一次预报结果，即9:00时刻的预测结果，从该数据中提取10:00时刻(即Δt取1小时)850hPa高度场(经变换可换算为 1250米)上与降水有关的要素(湿度、降水量、气压、温度、水凝物以及冰粒子等信息)，换算成雷达反射率数据。该换算采用本领域公知的算法，如参考文献[1]。
步骤3：将步骤2计算出的基本反射率数据绘制成雷达回波图像，绘制的算法与步骤1相同，绘制出的图像记为MImg(t0+Δt,h,x’,y’)，简称MImg；其中，x’和y’分别表示由数值预报产品绘制出的回波图像上任一像素的横坐标和纵坐标；再利用空间插值算法中的反距离权重法，将MImg(t0+Δt,h,x’,y’)转为与RImg(t0+Δt,h,x,y)具有相同地理范围、相同图像大小和相同坐标系的回波图像MImg(t0+Δt,h,x,y)；
步骤4：对雷达回波外推图像RImg中的发散点和发散小区域进行识别，识别的方法为：
其中，RImg＝255表示RImg中某一坐标(x,y)处的像素值为255；Count(RImg,r)表示RImg中以r为半径，(x,y)为圆心的区域内像素值小于255的像素的数量；Sum(RImg,r)表示RImg中以r为半径，(x,y)为圆心的区域内像素的总数量；r为大于1的正整数，一般情况下，r的建议值域为[2,5]；Ta是一个经验阈值，Ta∈[0,1]，Ta的物理意义为：在以r为半径，(x,y)为圆心的区域内，像素值小于255的像素数量占该区域内全部像素数量的比例。式1的物理意义为：当某一坐标(x,y)处的像素值为255，且在以r为半径，(x,y)为圆心的区域内，像素值小于255的像素数量与该区域内全部像素数量的比大于Ta时，则判定当前坐标(x,y)处存在发散点。本应用实例中，r＝3，Ta＝0.88。
步骤5：如果式1成立，即雷达回波外推图像RImg中(x,y)处为发散点，则对该点的像素值进行修补，修补的方法为：
首先，定义一个数据集DSR(x,y,r)，用于记录RImg中以r为半径，以坐标(x,y)为圆心的区域内所有像素值小于255的坐标。
然后，对数据集DSR(x,y,r)的每一个坐标，从MImg图像中提取相应坐标位置的像素值，并对这些像素值按其数值大小进行排序，形成有序序列MList(i，x,y)；其中，i为该序列中各项的序号。提取MImg图像中对应坐标(x,y)处的像素值，再从MList(i，x,y)中找到与该像素值差值绝对值最小的一项的序号，并将该最小项的序号记为Pos(x,y)。当有序序列MList中存在数值大小相同的像素值时，这些像素值之间的序列顺序可随机排序；当MList(i，x,y)中与 MImg(x,y)像素差值绝对值最小存在多项时，多项中可任意选择其中一项的序号记为Pos(x,y)。
接着，对数据集DSR(x,y,r)中的每一个坐标，从RImg图像中提取相应坐标位置的像素值，并对这些像素值按其数值大小进行排序，与MList(i，x,y)排序方式相同，形成有序序列RList(i， x,y)，i为有序序列各项的序号。
最后，计算RImg中(x,y)处的像素值，计算方法为：
其中，RList(i,x,y)[Pos(x,y)-1]表示有序序列RList(i,x,y)中Pos(x,y)-1项的像素值；类似的， RList(i,x,y)[Pos(x,y)+1]表示有序序列RList(i,x,y)中Pos(x,y)+1项的像素值。计算过程中，存在一些特殊情况，可另外处理，如：当Pos(x,y)＝0时，RImg＝RList(i,x,y)[Pos(x,y)+1]；当 Pos(x,y)＝Count(MImg,r)时，RImg＝RList(i,x,y)[Pos(x,y)-1]。
步骤6：对雷达回波外推图像RImg中的每个像素点，根据步骤4进行判别是否为发散点；如果为发散点，则按照步骤5，对发散点的像素值逐个进行修补计算。当该高度层中所有发散点全部修补后，计算结果对应存储到RImg中替换修补前的像素值，如果不是发散点，则该像素点跳过步骤5的计算，不进行替换。由此得到修补后的外推回波图像，记为RImg’(t0+Δt,h, x,y)，简称RImg’。
由于大气中风暴体的生消、演变过程十分复杂，而仅依靠雷达探测资料进行外推预测的算法难以表征这一复杂的过程，因此，再利用数值预报产品换算出的基本反射率及其回波图像，对由步骤6修补后的回波图像RImg’进行调整，计算过程按步骤7至步骤11所述。
步骤7：将步骤1起报时刻的雷达探测资料采用现有方法绘制成雷达回波图像，记为 RImg(t0,h,x,y)；将步骤2中与起报时刻时间上最相近的9：00的GFS数值预报产品，按步骤 2和步骤3的方法绘制成回波图像，记为MImg(t0,h,x,y)。计算雷达回波图像RImg(t0,h,x,y) 与RImg’(t0+Δt,h,x,y)各像素的变化值，以及MImg(t0,h,x,y)与MImg(t0+Δt,h,x,y)各像素的变化值，计算方法为：
RDif(x,y)＝RImg(t0,h,x,y)-RImg’(t0+Δt,h,x,y)
MDif(x,y)＝MImg(t0,h,x,y)-MImg(t0+Δt,h,x,y)
其中，RDif(x,y)是一个数据集，记录了雷达回波图像上的每一个像素从t0时刻到t0+Δt 时刻，像素值的数学差值；类似的，MDif(x,y)也是一个数据集，记录了数值预报产品输出的回波图像上的每一个像素从t0时刻到t0+Δt时刻，像素值的数学差值。
步骤8：计算数据集RDif(x,y)的中值，即中位数，记为RMV，计算方法是将记录于RDif(x, y)中的各项按值大小排序后，提取处于该序列中间位置的值。该值反映了雷达回波图像从t0 时刻到t0+Δt时刻演变的总体水平。类似地，计算数据集MDif(x,y)的中值，记为MMV。
步骤9：计算MMV与RMV的倍率，计算方法为：
MP＝MMV/RMV
步骤10：利用MP对修补后的雷达回波外推图像称RImg’的像素值逐一进行修正，修正方法为：
RImg_New(t0+Δt,h,x,y)＝RImg’(t0+Δt,h,x,y)+Tc×MP×RDif(x,y)
其中，RImg_New(t0+Δt,h,x,y)为计算结果，其物理意义表示：以t0时刻为起报时间，外推预测Δt时刻，h高度上，任一坐标(x,y)处的像素值，该值是通过本发明技术方案修正后的像素值，简称RImg_New。Tc是一个经验系数，用于调整修正的幅度，Tc∈[-1,1]，具体取值需根据季节、降水类型并结合以往外推预报的准确率做出适度调整,本应用实例中，Tc＝0.2。
步骤11：通过上述步骤的计算，可得到修正后的雷达回波外推图像，如图4所示。为了便于比对分析，图3显示了与图4和图2相同时刻，即2018-05-18 09:46的实况雷达回波图像。(由于雷达探测的周期理论值为6分钟，实际上存在数秒的误差，因此出现基于8:45外推1小时即9:45的预测，雷达实际探测的时间为9:46，相差了约1分钟。)可以看出，图2所显示的回波区域中存在若干不规则的白点噪声，即所谓的发散现象，而图4中的发散现象相较图2 有明显的改善。并且，相较图2，图4所示的回波强度及其空间分布与图3所示的实况具有更佳的一致性。
进一步地，采用与基本反射率数据转化为回波图像像素值相反的计算过程，将RImg_New 中各个像素的像素值转化为基本反射率数据，再通过雷达气象学中的Z-R关系理论方法，将基本反射率计算转化为可降水量预报信息，为相关领域的专业人员提供辅助决策。
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。
参考文献
[1]Kober K,Craig G C,Keil C,et al.Blending a probabilistic nowcastingmethod with a high-resolution numerical weather prediction ensemble forconvective precipitation forecasts[J]. Quarterly Journal of the RoyalMeteorological Society,2012,138(664):755–768。
Claims (7)
1.一种修正雷达回波外推图像发散现象的方法，其特征在于：包括如下步骤：
步骤1：取某起报时刻雷达探测资料并采用雷达回波外推预报算法，计算未来某一时刻t0+Δt各高度层的雷达基本反射率数据，再将雷达基本反射率数据绘制成对应各高度层的雷达回波外推图像，记为RImg(t0+Δt,h,x,y)，简称RImg；其中，t0表示起报时间；Δt表示外推预报的时效，0<Δt<3小时；h表示高度层的高度，0<h<20km；x和y分别表示雷达回波外推图像上任一像素的横坐标和纵坐标；
步骤2：取与步骤1雷达探测资料起报时刻t0时间上最相近的数值预报产品，从数值预报产品中提取与t0+Δt时刻最相近的降水要素并计算出与步骤1雷达回波外推图像各个相同高度层对应的基本反射率数据；
步骤3：将步骤2计算出的各个高度层的基本反射率数据绘制成回波图像记为MImg(t0+Δt,h,x’,y’)，简称MImg；其中，x’和y’分别表示该回波图像上任一像素的横坐标和纵坐标；再利用空间插值算法，将MImg(t0+Δt,h,x’,y’)转为与RImg(t0+Δt,h,x,y)具有相同地理范围、相同图像大小和相同坐标系的回波图像MImg(t0+Δt,h,x,y)；
步骤4：对雷达回波外推图像RImg(t0+Δt,h,x,y)中的发散点进行识别，判断是否存在下式情况：
RImg＝255，且
式1
其中，RImg＝255表示RImg中某一坐标(x,y)处的像素值为255；Count(RImg,r)表示RImg中以r为半径，(x,y)为圆心的区域内像素值小于255的像素的数量；Sum(RImg,r)表示RImg中以r为半径，(x,y)为圆心的区域内像素的总数量；r取[2,5]的正整数；Ta为经验阈值；
步骤5：如果式1情况存在，即雷达回波外推图像RImg中(x,y)处为发散点，则对该点的像素值进行修补，修补过程如下：
首先，定义一个数据集DSR(x,y,r)，用于记录RImg中以r为半径，以坐标(x,y)为圆心的区域内所有像素值小于255的坐标；
然后，对数据集DSR(x,y,r)中的每一个坐标，从MImg图像中提取相应坐标位置的像素值，并对这些像素值按其数值大小进行排序，形成有序序列MList(i,x,y)；其中，i为该序列中各项的序号；提取MImg图像中坐标(x,y)处的像素值，再从MList(i,x,y)中找到与该像素值差值绝对值最小的一项的序号，并将该序号记为Pos(x,y)；
接着，对数据集DSR(x,y,r)中的每一个坐标，从RImg图像中提取相应坐标位置的像素值，并对这些像素值按其数值大小进行排序，形成有序序列RList(i,x,y)，i为有序序列各项的序号；
最后，按下式计算RImg中(x,y)处的像素值：
其中，RList(i,x,y)[Pos(x,y)-1]表示有序序列RList(i,x,y)中Pos(x,y)-1项的像素值；RList(i,x,y)[Pos(x,y)+1]表示有序序列RList(i,x,y)中Pos(x,y)+1项的像素值；
步骤6：对雷达回波外推图像RImg中的每个像素点，逐层根据步骤4的方法识别是否为发散点；如果为发散点，则按照步骤5的方法对发散点的像素值逐个进行修补计算，当所有高度层所有发散点全部修补后，计算结果对应存储到RImg中替换修补前的像素值，如果不是发散点，不进行替换；由此得到修补后的雷达回波外推图像RImg’(t0+Δt,h,x,y)，简称RImg’；
步骤7：将步骤1起报时刻雷达探测资料绘制成雷达回波图像，记为RImg(t0,h,x,y)；将上述与起报时刻t0时间上最相近的数值预报产品按步骤2和步骤3的方法绘制成回波图像，记为MImg(t0,h,x,y)；计算雷达回波图像RImg(t0,h,x,y)与RImg’(t0+Δt,h,x,y)各像素的变化值，以及MImg(t0,h,x,y)与MImg(t0+Δt,h,x,y)各像素的变化值，计算方法为：
RDif(x,y)＝RImg(t0,h,x,y)-RImg’(t0+Δt,h,x,y)
MDif(x,y)＝MImg(t0,h,x,y)-MImg(t0+Δt,h,x,y)
其中，RDif(x,y)是一个数据集，记录了雷达回波图像上的每一个像素从t0时刻到t0+Δt时刻像素值的差值；MDif(x,y)也是一个数据集，记录了数值预报产品回波图像上的每一个像素从t0时刻到t0+Δt时刻像素值的差值；
步骤8：计算数据集RDif(x,y)的中值，记为RMV，计算方法是将RDif(x,y)中各项按其数值大小排序后，提取处于该序列中间位置的值；计算数据集MDif(x,y)的中值，记为MMV；
步骤9：按下式计算MMV与RMV的倍率MP：MP＝MMV/RMV；
步骤10：利用倍率MP对修补后的雷达回波外推图像RImg’的像素值按下式逐层逐一进行修正：
RImg_New(t0+Δt,h,x,y)＝RImg’(t0+Δt,h,x,y)+Tc×MP×RDif(x,y)
其中，修正后的雷达回波外推图像RImg_New(t0+Δt,h,x,y)的物理意义表示：以t0时刻为起报时间，外推预测Δt时刻，h高度上，任一坐标(x,y)处的像素值，简称RImg_New；Tc为经验系数。
2.根据权利要求1所述修正雷达回波外推图像发散现象的方法，其特征在于：步骤1中雷达回波外推预报算法采用最大相关法或交叉相关法。
3.根据权利要求2所述修正雷达回波外推图像发散现象的方法，其特征在于：步骤2中数值预报产品中降水要素包括湿度、降水量、气压、温度、水凝物和冰粒子气象信息。
4.根据权利要求3所述修正雷达回波外推图像发散现象的方法，其特征在于：步骤3中空间插值算法采用反距离权重算法。
5.根据权利要求4所述修正雷达回波外推图像发散现象的方法，其特征在于：步骤5中当有序序列MList中存在数值大小相同的像素值时，这些像素值之间的序列顺序可随机排序；当MList(i,x,y)中与MImg(x,y)像素差值绝对值最小存在多项时，多项中可任意选择其中一项的序号记为Pos(x,y)。
6.根据权利要求5所述修正雷达回波外推图像发散现象的方法，其特征在于：步骤5计算过程中，当Pos(x,y)＝0时，则RImg＝RList(i,x,y)[Pos(x,y)+1]；当Pos(x,y)＝Count(MImg,r)时，则RImg＝RList(i,x,y)[Pos(x,y)-1]。
7.根据权利要求1-6任一所述修正雷达回波外推图像发散现象的方法，其特征在于：对于步骤10得到的修正后的雷达回波外推图像RImg_New(t0+Δt,h,x,y)，采用与基本反射率数据转化为回波图像像素值相反的计算过程，将RImg_New(t0+Δt,h,x,y)中各个像素的像素值转化为基本反射率数据，再通过雷达气象学中的Z-R关系理论方法，将基本反射率计算转化为降水方面的预报信息。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811019917.0A CN109283505B (zh) | 2018-09-03 | 2018-09-03 | 一种修正雷达回波外推图像发散现象的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201811019917.0A CN109283505B (zh) | 2018-09-03 | 2018-09-03 | 一种修正雷达回波外推图像发散现象的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109283505A CN109283505A (zh) | 2019-01-29 |
| CN109283505B true CN109283505B (zh) | 2022-06-07 |
Family
ID=65183843
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201811019917.0A Active CN109283505B (zh) | 2018-09-03 | 2018-09-03 | 一种修正雷达回波外推图像发散现象的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109283505B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN111352113B (zh) * | 2020-04-01 | 2022-08-26 | 易天气(北京)科技有限公司 | 一种强对流天气短临预报方法及***、存储介质和终端 |
| CN112363168B (zh) * | 2021-01-13 | 2021-03-26 | 南京满星数据科技有限公司 | 一种基于雷达外推和模式预报的同化融合方法 |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JP5115912B2 (ja) * | 2001-02-23 | 2013-01-09 | 独立行政法人日本原子力研究開発機構 | 高速ゲート掃引型３次元レーザーレーダー装置 |
| CN100487522C (zh) * | 2005-01-13 | 2009-05-13 | 奥林巴斯映像株式会社 | 模糊校正方法及摄像装置 |
| CN105717491B (zh) * | 2016-02-04 | 2018-09-21 | 象辑知源（武汉）科技有限公司 | 天气雷达回波图像的预测方法及预测装置 |
| CN106886023B (zh) * | 2017-02-27 | 2019-04-02 | 中国人民解放军理工大学 | 一种基于动态卷积神经网络的雷达回波外推方法 |
| CN107632295A (zh) * | 2017-09-15 | 2018-01-26 | 广东工业大学 | 一种基于时序卷积神经网络的雷达回波外推方法 |
-
2018
- 2018-09-03 CN CN201811019917.0A patent/CN109283505B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN109283505A (zh) | 2019-01-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN113496104A (zh) | 基于深度学习的降水预报订正方法及*** | |
| Atencia et al. | A comparison of two techniques for generating nowcasting ensembles. Part II: Analogs selection and comparison of techniques | |
| CN109917394B (zh) | 一种基于天气雷达的短临智能外推方法 | |
| CN109283505B (zh) | 一种修正雷达回波外推图像发散现象的方法 | |
| CN104237890B (zh) | 一种由“列车效应”引起的暴雨识别及预报方法 | |
| WO2019225064A1 (ja) | 気象予測装置、気象予測方法、並びに風力発電出力推定装置 | |
| CN110942111A (zh) | 一种识别强对流云团的方法及装置 | |
| CN115660137B (zh) | 一种船舶风浪航行能耗精准估算方法 | |
| CN116863274A (zh) | 一种基于半监督学习的钢板表面缺陷检测方法及*** | |
| Squitieri et al. | On the forecast sensitivity of MCS cold pools and related features to horizontal grid spacing in convection-allowing WRF simulations | |
| CN115691049A (zh) | 一种基于深度学习的对流初生预警方法 | |
| CN108732648B (zh) | 一种面向山地暴雨预报的渐进决策方法 | |
| Brisson et al. | Lagrangian evaluation of convective shower characteristics in a convection-permitting model | |
| Phan-Van et al. | Evaluation of the NCEP climate forecast system and its downscaling for seasonal rainfall prediction over Vietnam | |
| Mandal | Geo-information based spatio-temporal modeling of urban land use and land cover change in Butwal municipality, Nepal | |
| Case et al. | An objective technique for verifying sea breezes in high-resolution numerical weather prediction models | |
| Li et al. | Improving multi-model ensemble probabilistic prediction of Yangtze River valley summer rainfall | |
| CN116451554A (zh) | 考虑多种气象因子的电网气象风险预测方法 | |
| JP2004317173A (ja) | 雷観測システム | |
| CN114325879A (zh) | 一种基于分级概率的定量降水订正方法 | |
| JPH0972965A (ja) | 気象擾乱の予測方法 | |
| Christodoulou et al. | Prediction of rainfall rate based on weather radar measurements | |
| Singh et al. | An analogue method for simultaneous prediction of surface weather parameters at a specific location in the Western Himalaya in India | |
| CN114004426B (zh) | 一种短时暴雨预报释用模型的动态调整方法 | |
| CN116955964B (zh) | 一种连续路面气象状况分析及推演方法 |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |