CN108761461B - 基于天气雷达回波时序图像的降水预报方法 - Google Patents
基于天气雷达回波时序图像的降水预报方法 Download PDFInfo
- Publication number
- CN108761461B CN108761461B CN201810531177.2A CN201810531177A CN108761461B CN 108761461 B CN108761461 B CN 108761461B CN 201810531177 A CN201810531177 A CN 201810531177A CN 108761461 B CN108761461 B CN 108761461B
- Authority
- CN
- China
- Prior art keywords
- img
- gray
- time
- image
- value
- 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
- 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
-
- 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
-
- 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
- G01S7/411—Identification of targets based on measurements of radar reflectivity
-
- 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)
- Computer Networks & Wireless Communication (AREA)
- General Physics & Mathematics (AREA)
- Electromagnetism (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明公开了一种利用天气雷达回波时序图像进行降水预报的方法,该方法提出基本反射率数值与图像灰度值之间的直接关联,同时2条合理假设并建立相应的计算模型,依次对时间上两两相邻的基本反射率灰阶图像进行回波主体移动特征计算,再通过卷积函数,对时序图像的移动速度矢量场进行叠加计算,并遵循发生时间越久远的移动速度矢量对当前及未来的影响越小的原则,接着利用一个灰阶图像变化预测模型,对未来灰阶图像的形态、强度进行预测,最终实现基于雷达资料的降水预测。该方法对于较长预报时效且时空变化不连续的回波主体,其降水预报结果的准确性有着明显改善。
Description
技术领域
本发明属于大气科学领域,涉及一种利用天气雷达回波时序图像进行降水预报的方法。
背景技术
基于天气雷达探测资料的临近预报技术对短时临近降水预报起到十分积极有效的作用。该技术从上世纪60年代开始被研究,随着雷达硬件基础设施和探测技术的不断更新,基于天气雷达的降水预报方法也取得了长足的进步。这类方法的基本思想是利用雷达反射率因子对强回波进行识别,再根据回波区域的移动趋势,对回波未来的方位、强度做出预判,进而达到预报的目的。从数理统计角度看,这类方法有效解决了雷达回波在未来时间发展变化的预测难题,但实际天气的发展演变过程往往更加复杂。从时间连续的回波图像可以看出,回波图像所表征的可降水区域每时每刻、从整体到局部都在不断地生消、发展、移动和变化,传统预报方法采用的预测模型通常是建立在线性运动关系的基础上,而实际大气及大气中粒子的运动轨迹并不符合线性运动的规律。因此,基于线性运动关系的降水预报,在预测较长时间的回波图像时,往往与实际偏差很大,进而降低了降水预报的准确性。
发明内容
本发明的目的是为了解决现有技术中存在的缺陷,提供一种降水预报结果更加准确的方法。
为了达到上述目的,本发明提供了一种基于天气雷达回波时序图像的降水预报方法,包括以下步骤:
(1)读取当前时刻T及前N各时刻的雷达基数据;
(2)读取(1)中当前时刻T雷达基数据中的各个探测仰角的基本反射率数值;
(3)将步骤(2)中各个仰角的基本反射率数值的值域扩大到[0,255]之间;
(4)将步骤(3)中扩大值域后的数值转换为灰阶图像;
(5)重复步骤(2)至步骤(4),对步骤(1)中读取的所有雷达基数据进行灰阶图像转换;
(6)根据以下原理,依次对时间上两两相邻雷达基数据转换后的基本反射率灰阶图像进行回波主体移动特征计算,从而计算出灰阶图像各像素的移动速度矢量:图像中任意位置相邻的像素,其移动趋势近乎一致;以及图像中表征雷达回波主体的灰阶图像区域,其总体移动趋势基本一致;
(7)根据发生时间越久远,对当前及未来影响越小的原则,对步骤(6)计算出的移动速度矢量进行叠加计算;
(8)构建灰阶图像变化预测模型:
且
式中,Δt表示预测的时间长度,Δt的取值以时间上两两相邻的基数据时间间隔为基本单位;Δx和Δy分别表示未来Δt时间后,原来T时刻坐标(x,y)分别在x轴和y轴的偏移量;MV(x,y)是移动速度矢量,它的值由步骤(7)计算得到;MV(x,y)x和MV(x,y)y分别表示MV(x,y)在x轴和y轴的速度分量;IMG(T,a,x,y)表示基数据BD(T)计算出的a仰角灰阶图像在直角坐标(x,y)处的灰阶强度值;同理,IMG(T-t,a,x,y)表示基数据BD(T-t)计算出的a仰角灰阶图像在直角坐标(x,y)处的灰阶强度值;h(x,y,t)是一个时域卷积函数,该函数的具体形式没有限定,但始终须遵循一个原则:发生时间越久远的雷达回波图像的像素灰度值对当前及未来的影响越小;*是卷积运算符;Img(T-t,a,x,y)*h(x,y,t)的物理意义为:T-t时刻的灰阶图像中,任一像素(x,y)处的灰度值在[0,N-1]时序中的t时刻对预测值P(T+Δt,a,x+Δx,y+Δy)的作用程度,该作用程度是由h(x,y,t)函数控制,换言之,h(x,y,t)函数是Img(T-t,a,x,y)的权重系数。P(T+Δt,a,x+Δx,y+Δy)是预测的结果,表示在未来Δt时刻,a仰角面上,坐标(x+Δx,y+Δy)处的灰阶值。
(9)对步骤(8)计算得到的P(T+Δt,a,x+Δx,y+Δy)进行平滑处理;
(10)对步骤(9)平滑处理后数值的值域进行收缩调整,将各像素的取值还原至雷达基数据的正常数量级;
(11)将步骤(10)还原后数值所表征的回波强度值转换为可降水量。
更加具体的,本发明降水预报方法如下:
步骤1)读取当前时刻T及前N个时刻的雷达基数据,其中,N是大于0的整数。不妨设BD(T-t)表示T-t时刻的雷达基数据(以下简称基数据),t为正整数,且0≤t<N。
步骤2)读取BD(T)中各个探测仰角的基本反射率数值。探测仰角的数量取决于不同雷达及雷达体扫的方式,不妨设某种雷达的基本反射率有效探测仰角的数量为EN,某一仰角的基本反射率数值为BRV(T,a,b,c),其中,a表示仰角,且a<EN;b表示方位角,且b∈[0,360);c表示雷达探测的距离库数,且c∈[1,MaxDis],MaxDis的数值取决于雷达技术参数。
步骤3)将BD(T)中各个仰角的基本反射率数值BRV(T,a,b,c)的值域扩大到[0,255]之间,计算方法为:
且当BRVG(T,a,b,c)<0时,BRVG(T,a,b,c)=0。其中,MAX(BRV(T,a))表示BRV(T,a,b,c)某一仰角a中,所有方位角和所有距离库数中基本反射率数值最大的一个值。
步骤4)将BRVG(T,a,b,c)采用通用的图像处理算法转化为灰阶图像IMG(T,a,x,y),其中,x,y分别表示灰阶图像的横坐标和纵坐标。上述极坐标形式的(b,c)转化为直角坐标形式的(x,y)可由公知的方法计算得到。IMG(T,a,x,y)表示基数据BD(T)计算出的a仰角灰阶图像在直角坐标(x,y)处的灰阶强度值。
步骤5)按照步骤2至步骤4,对步骤1中所有基数据RBD(T-t)依次进行上述处理,则任一时刻任一仰角的基本反射率灰阶图像可表示为IMG(T-t,a,x,y)。
步骤6)依次对时间上两两相邻的基本反射率灰阶图像进行回波主体移动特征计算,计算的原理是基于如下2点假设:
假设1:图像中任意位置相邻的像素,其移动趋势是近乎一致的,用公式可表示为:
V(T-t,x,y)=V(T-t-1,x+Δx,y+Δy)(Δx→0,Δy→0) (式2)
其中,V表示移动速度矢量,V(T-t,x,y)表示基本反射率灰阶图像中直角坐标系下(x,y)处的移动速度矢量,按照假设1,该值与V(T-t-1,x,y)的值近似相等,即V相邻坐标位置的梯度变化近似为0。由此,可建立式3所示的方程:
假设2:图像中表征雷达回波主体的灰阶图像区域,其总体移动趋势是基本一致的。换言之,虽然无法断定“图像中雷达回波主体完全按统一的方向和速度移动”,但可以通过约束项的极小化来进行限定,用公式可表示为:
其中,min[]表示遍历x方向和y方向所有像素后,求取数值最小的一项。
将式3和式4合并,构成用于求取基本反射率灰阶图像移动矢量的约束条件,如式5所示:
其中,μ∈[0,1],它是假设2的调节因子,用于调节灰阶图像的“噪声”与假设2约束之间的权重,图像中的“噪声”越大,说明灰阶图像本身的置信度越低,亦即雷达回波中的杂波越严重,为了能够计算出合理的E值,则需要增加假设2约束的权重。反之,则μ可以取较小的值。
上述关于灰阶图像中的“噪声”,并非传统意义上的图像光学噪点,而是主要受到反映回波强度和区域面积的水汽粒子不断生消变化的影响,使得探测时间相邻的雷达回波灰阶图像之间具有明显的差异,这种差异对于式5所示的假设来说,是非常严重的“噪声”。
为了便于计算,有必要将积分形式的式5简化为离散化形式:
其中,
式7中,D(i,j,T-t)x、D(i,j,T-t)y和D(i,j,T-t)t分别表示图像中的任一像素点(i,j)的灰度值在x、y和时间T轴方向上的变化率。为了简化计算,可采用如下的差分公式进行计算:
D(i,j,T-t)x=(IMG(T-t,a,i+1,j)-IMG(T-t,a,i,j)
+IMG(T-t-1,a,i+1,j)-IMG(T-t-1,a,i,j)
+IMG(T-t,a,i+1,j+1)-IMG(T-t,a,i,j+1)
+IMG(T-t-1,a,i+1,j+1)-IMG(T-t-1,a,i,j+1))/4δx
D(i,j,T-t)y=(IMG(T-t,a,i+1,j)-IMG(T-t,a,i,j)
+IMG(T-t-1,a,i+1,j)-IMG(T-t-1,a,i,j)
+IMG(T-t,a,i+1,j+1)-IMG(T-t,a,i,j+1)
+IMG(T-t-1,a,i+1,j+1)-IMG(T-t-1,a,i,j+1))/4δy
D(i,j,T-t)t=(IMG(T-t,a,i+1,j)-IMG(T-t,a,i,j)
+IMG(T-t-1,a,i+1,j)-IMG(T-t-1,a,i,j)
+IMG(T-t,a,i+1,j+1)-IMG(T-t,a,i,j+1)
+IMG(T-t-1,a,i+1,j+1)-IMG(T-t-1,a,i,j+1))/4δt
(式8)
式6至式8表达了求解连续两幅灰阶图像中(x,y)处移动速度矢量的计算过程,(x,y)是灰阶图像中任一像素的坐标,大写的X和Y表示灰阶图像的长(即横坐标方向)和宽(即纵坐标方向),单位为像素数。i和j是计算中的变量,其取值取决于式6中的x和y。δx、δy和δt分别表示灰阶图像中的像素在x轴、y轴和时间轴方向上的宽度。和的数学含义是灰阶图像中的像素点(i,j)在x和y轴方向上的移动特征分量。实际计算中,可使用以下参考值:δx=δy=1,δt=6(雷达数据的时间间隔)。由于D(i,j,T-t)x、D(i,j,T-t)y和D(i,j,T-t)t可以由式8计算得到,因此,式6的计算转化为对自变量和的求解,该求解过程数学上有多种方法,典型的如,最小二乘法。为了便于描述,不妨将T-t与T-t-1两个时刻的灰阶图像在(x,y)处求得的最优解表示为V(T-t,x,y),该值即为(x,y)处的移动速度矢量。
需要补充说明的是,式6至式8是为了便于计算对式5做出的简化,式6至式8并不等价于式5,且对式5简化的方法也不限于式6至式8。步骤6的关键内容,是基于上述2条原理,对时间上两两相邻的基本反射率灰阶图像进行回波主体移动特征的计算,计算出灰阶图像各像素的移动速度矢量。
步骤7)根据发生时间越久远,对当前及未来影响越小的原则,对步骤6计算出的移动速度矢量V(T-t,x,y)进行叠加计算,计算方法为:
其中,x0和y0分别表示每次叠加计算时灰阶图像中参与计算的空间邻域的大小,取值为大于1的整数。*是卷积运算符,h(i,j,t)是时域卷积函数,该函数的具体形式没有限定,但始终须遵循一个原则:发生时间越久远的移动速度矢量对当前及未来的影响越小。式9中V(T-t,i,j)*h(i,j,t)的物理意义为:由IMG(T-t,i,j)和IMG(T-t-1,i,j)通过步骤6计算得到的移动速度矢量V(T-t,i,j),在[0,N]时序中的t时刻对MV(x,y)的作用程度,该作用程序是由h(i,j,t)函数控制,换言之,h(i,j,t)函数是V(T-t,i,j)的权重系数。实际计算过程中,h(i,j,t)函数可采用基于时间t的线性或非线性约束关系。
步骤8)根据大气中水凝物等粒子的运动规律及其在雷达回波图像上所表征出的特点,构建一个灰阶图像变化预测模型,该模型的计算方法为:
且
其中,Δt表示预测的时间长度,Δt的取值以时间上两两相邻的基数据时间间隔为基本单位。举例来说,雷达某种探测方式下,基数据时间间隔为6分钟,如果预测未来6分钟,则Δt=1,如果预测未来30分钟,则Δt=5。Δx和Δy分别表示未来Δt时间后,原来T时刻坐标(x,y)分别在x轴和y轴的偏移量。MV(x,y)x和MV(x,y)y分别表示MV(x,y)在x轴和y轴的速度分量。P(T+Δt,a,x+Δx,y+Δy)表示在未来Δt时刻,a仰角面上,坐标(x+Δx,y+Δy)处的灰阶值。h(x,y,t)与步骤7中h(i,j,t)相似,为时域卷积函数,该函数的具体形式没有限定,但始终须遵循一个原则:发生时间越久远的雷达回波图像的像素灰度值对当前及未来的影响越小,*是卷积运算符,Img(T-t,a,x,y)*h(x,y,t)的物理意义为:T-t时刻的灰阶图像中,像素(x,y)处的灰度值在[0,N-1]时序中的t时刻对P(T+Δt,a,x+Δx,y+Δy)的作用程度,该作用程序是由h(x,y,t)函数控制,换言之,h(x,y,t)函数是Img(T-t,a,x,y)的权重系数。
步骤9)为了改善步骤8计算中,T+Δt时刻的坐标(x+Δx,y+Δy)(x∈[0,X),y∈[0,Y))不能完全覆盖整个灰阶图像,进而造成图像发散的问题,再对步骤8计算得到的P(T+Δt,a,x+Δx,y+Δy)做适当的平滑处理,使得所预测的回波图像在各个局,部区域有着更好的一致性。平滑处理可采用公开的图像平滑或图像滤波算法实现,如高斯滤波、中值滤波等。
步骤10)按照与步骤3相反的过程,对P(T+Δt,a,x,y)的值域进行收缩调整,将各像素的取值还原到基数据的正常量级,计算方法为:
当前步骤中,P(T+Δt,a,x,y)与式10的P(T+Δt,a,x+Δx,y+Δy)是同一变量,只是坐标的表述方法有所不同。式10中Δx和Δy是用于表述当前任一坐标(x+Δx,y+Δy)相对于Img(T,a,x,y)的坐标(x,y)的偏移量。P(T+Δt,a,x,y)则表示T+Δt时刻,回波图像中任一坐标(x,y)处的灰阶强度值;类似地,P’(T+Δt,a,x,y)表示T+Δt时刻,任一坐标(x,y)处的雷达基本反射率数值。
步骤11)采用雷达气象学相关领域的Z-R关系,将P’(T+Δt,a,x,y)所表征的回波强度值转换为可降水量,计算方法为:
R(x,y)=α×P′(T+Δt,a,x,y)β (式12)
其中,α和β是Z-R关系式中的可变参量,通常根据降水类型的不同,取相应的经验值。
通过上述步骤1至步骤11的计算,实现基于天气雷达回波时序图像的降水预报。
本发明的优点如下:
本发明提出一种基于天气雷达回波时序图像的降水预报方法,首先针对雷达回波图像中RGB各颜色通道的值与实际回波强度的大小关系不对应的问题,提出基本反射率数值与图像灰度值之间的直接关联,为后续计算提供更能反映基本反射率强弱连续性关系的灰阶图像。其次,提出2条合理假设并建立相应的计算模型,依次对时间上两两相邻的基本反射率灰阶图像进行回波主体移动特征计算,再通过卷积函数,对时序图像的移动速度矢量场进行叠加计算,并遵循发生时间越久远的移动速度矢量对当前及未来的影响越小的原则。接着,提出一个灰阶图像变化预测模型,对未来灰阶图像的形态、强度进行预测,并提出改善预测图像发散问题的方法。最终实现基于雷达资料的降水预测。该方法对于较长预报时效且时空变化不连续的回波主体,其降水预报结果的准确性有着明显改善。
附图说明
图1是本发明的主要计算步骤流程图;
图2是本实施例中,6个雷达基数据按步骤1至步骤4计算得到的IMG(T-t,a,x,y)绘制出的基本反射率灰阶图像,其中,a均为0.5°仰角;
图2中,a为05:39的灰阶图像IMG(T-5,0.5°,x,y),b为05:45的灰阶图像IMG(T-4,0.5°,x,y),c为05:51的灰阶图像IMG(T-3,0.5°,x,y),d为05:57的灰阶图像IMG(T-2,0.5°,x,y),e为06:04的灰阶图像IMG(T-1,0.5°,x,y),f为06:10的灰阶图像IMG(T,0.5°,x,y);
图3是采用图2所示的6幅灰阶图像,依次两两按照步骤6所述计算过程,得到的移动速度矢量,绘制出的矢量场图像;
图3中,a为05:39~05:45移动速度矢量场图像,b为05:45~05:51移动速度矢量场图像,c为05:51~05:57移动速度矢量场图像,d为05:57~06:04移动速度矢量场图像,e为06:04~06:10移动速度矢量场图像;
图4以图形化方式展示了本实施例中,Δt=1、Δt=5和Δt=10,即未来6分钟(a)、30分钟(b)和60分钟(c)的移动速度矢量场图像;
图5以图形化方式展示了本实施例中,Δt=1、Δt=5和Δt=10,即未来6分钟、30分钟和60分钟的雷达回波灰阶图像,其中,每行左侧的图像为预测图像,右侧为相同时间的实况图像,所有图像都经过了步骤3所述的转换;
图5中,a为预测出的未来6分钟的回波灰阶图像,b为与a同时刻的实况图像(06:16:19);c为预测出的未来30分钟的回波灰阶图像,d为与c同时刻的实况图像(06:40:58);e为预测出的未来60分钟的回波灰阶图像,f为与e同时刻的实况图像(07:11:43)。
具体实施方式
本实施例选用2016年8月10日05:39至07:11时间段南京地区的雷达基数据作为实验对象。
步骤1)读取当日06:10及前5个时刻的雷达基数据,为便于表述,分别定义这些基数据名为BD(T)、BD(T-1)、BD(T-2)、BD(T-3)、BD(T-4)、BD(T-5)。该部雷达的基数据为CINRAD SA格式,体扫时间间隔为6分钟,基本反射率有效探测仰角数为9个。
步骤2)读取BD(T)中各个探测仰角的基本反射率数值。定义某一仰角的基本反射率数值为BRV(T,a,b,c),其中,a表示仰角,且a∈[0.5°,1.5°,2.4°,3.4°,4.3°,6.0°,9.9°,14.6°,19.5°];b表示方位角,且b∈[0,360);c表示雷达探测的距离库数,且c∈[1,MaxDis],MaxDis的数值取决于雷达技术参数,如当前雷达0.5°仰角的最大距离库数MaxDis=460。
步骤3)将BD(T)中各个仰角的基本反射率数值BRV(T,a,b,c)的值域扩大到[0,255]之间,计算方法为:
且当BRVG(T,a,b,c)<0时,BRVG(T,a,b,c)=0。其中,MAX(BRV(T,a))表示BRV(T,a,b,c)某一仰角a中,所有方位角和所有距离库数中基本反射率数值最大的一个值。
步骤4)将BRVG(T,a,b,c)采用通用的图像处理算法转化为灰阶图像IMG(T,a,x,y),其中,x,y分别表示灰阶图像的横坐标和纵坐标。上述极坐标形式的(b,c)转化为直角坐标形式的(x,y)可由公知的方法计算得到。IMG(T,a,x,y)表示基数据BD(T)计算出的a仰角灰阶图像在直角坐标(x,y)处的灰阶强度值。
步骤5)按照步骤2至步骤4,对步骤1中所有基数据RBD(T-t)依次进行上述处理,则任一时刻任一仰角的基本反射率灰阶图像可表示为IMG(T-t,a,x,y)。上述6个雷达基数据由步骤1至步骤4,绘制出的灰阶图像(均为0.5°仰角)如图2所示。
步骤6)依次对时间上两两相邻的基本反射率灰阶图像进行回波主体移动特征计算,计算的原理是基于如下2条假设:
假设1:图像中任意位置相邻的像素,其移动趋势是近乎一致的,用公式可表示为:
V(T-t,x,y)=V(T-t-1,x+Δx,y+Δy)(Δx→0,Δy→0) (式2)
其中,V表示移动速度矢量,V(T-t,x,y)表示基本反射率灰阶图像中直角坐标系下(x,y)处的移动速度矢量,按照假设1,该值与V(T-t-1,x,y)的值近似相等,即V相邻坐标位置的梯度变化近似为0。由此,可建立式3所示的方程:
假设2:图像中表征雷达回波主体的灰阶图像区域,其总体移动趋势是基本一致的。换言之,虽然无法做出“图像中雷达回波主体完全按统一的方向和速度移动”的假设,但可以通过约束项的极小化来进行限定,用公式可表示为:
min[]表示遍历x方向和y方向所有像素后,求取数值最小的一项。
将上述2条假设的式3和式4合并,构成用于求取基本反射率灰阶图像移动矢量的约束条件,如式5所示:
其中,μ是假设2的调节因子,用于调节灰阶图像的“噪声”与假设2约束之间的权重,图像中的“噪声”越大,说明灰阶图像本身的置信度越低,亦即雷达回波中的杂波越严重,为了能够计算出合理的E值,则需要增加假设2约束的权重。反之,则μ可以取较小的值。在本实施例中,μ取值为0.5。
上述关于灰阶图像中的“噪声”,并非传统意义上的图像光学噪点,而是主要受到反映回波强度和区域面积的水汽粒子不断生消变化的影响,使得探测时间相邻的雷达回波灰阶图像之间具有明显的差异,这种差异对于式5所示的假设来说,是非常严重的“噪声”。
为了便于计算,有必要将积分形式的式5简化为离散化形式:
其中,
式7中,D(i,j,T-t)x、D(i,j,T-t)y和D(i,j,T-t)t分别表示图像中的任一像素点(i,j)的灰度值在x、y和时间T轴方向上的变化率。为了简化计算,可采用如下的差分公式进行计算:
D(i,j,T-t)x=(IMG(T-t,a,i+1,j)-IMG(T-t,a,i,j)
+IMG(T-t-1,a,i+1,j)-IMG(T-t-1,a,i,j)
+IMG(T-t,a,i+1,j+1)-IMG(T-t,a,i,j+1)
+IMG(T-t-1,a,i+1,j+1)-IMG(T-t-1,a,i,j+1))/4δx
D(i,j,T-t)y=(IMG(T-t,a,i+1,j)-IMG(T-t,a,i,j)
+IMG(T-t-1,a,i+1,j)-IMG(T-t-1,a,i,j)
+IMG(T-t,a,i+1,j+1)-IMG(T-t,a,i,j+1)
+IMG(T-t-1,a,i+1,j+1)-IMG(T-t-1,a,i,j+1))/4δy
D(i,j,T-t)t=(IMG(T-t,a,i+1,j)-IMG(T-t,a,i,j)
+IMG(T-t-1,a,i+1,j)-IMG(T-t-1,a,i,j)
+IMG(T-t,a,i+1,j+1)-IMG(T-t,a,i,j+1)
+IMG(T-t-1,a,i+1,j+1)-IMG(T-t-1,a,i,j+1))/4δt
(式8)
式6至式8表达了求解连续两幅灰阶图像中(x,y)处移动速度矢量的计算过程,(x,y)是灰阶图像中任一像素的坐标,大写的X和Y最示灰阶图像的长(即横坐标方向)和宽(即纵坐标方向),单位为像素数。i和j是计算中的变量,其取值取决于式6中的x和y。δx、δy和δt分别表示灰阶图像中的像素在x轴、y轴和时间轴方向上的宽度。和的数学含义是灰阶图像中的像素点(i,j)在x和y轴方向上的移动特征分量。实际计算中,可使用以下参考值:δx=δy=1,δt=6(雷达数据的时间间隔)。由于D(i,j,T-t)x、D(i,j,T-t)y和D(i,j,T-t)t可以由式8计算得到,因此,式6的计算转化为对自变量和的求解,该求解过程数学上有多种方法,典型的如,最小二乘法。为了便于描述,不妨将T-t与T-t-1两个时刻的灰阶图像在(x,y)处求得的最优解表示为V(T-t,x,y),该值即为(x,y)处的移动速度矢量。
需要补充说明的是,式6至式8是为了便于计算对式5做出的简化,式6至式8并不等价于式5,且对式5简化的方法也不限于式6至式8。步骤6的关键内容,是基于上述2条假设,对时间上两两相邻的基本反射率灰阶图像进行回波主体移动特征的计算,计算出灰阶图像各像素的移动速度矢量。
上述如图2所示的6幅灰阶图像,依次两两按照步骤6所述计算过程,得到的移动速度矢量,绘制出的矢量场图像如图3所示。
步骤7)根据发生时间越久远,对当前及未来影响越小的原则,对步骤6计算出的移动速度矢量V(T-t,x,y)进行叠加计算,计算方法为:
其中,x0和y0分别表示每次叠加计算时灰阶图像中参与计算的空间邻域的大小,取值为大于1的整数。*是卷积运算符,h(i,j,t)是时域卷积函数,该函数的具体形式没有限定,但始终须遵循一个原则:发生时间越久远的移动速度矢量对当前及未来的影响越小。式9中V(T-t,i,j)*h(i,j,t)的物理意义为:由IMG(T-t,i,j)和IMG(T-t-1,i,j)通过步骤6计算得到的移动速度矢量V(T-t,i,j),在[0,N]时序中的t时刻对MV(x,y)的作用程度,该作用程度是由h(i,j,t)函数控制,换言之,h(i,j,t)函数是W(T-t,i,j)的权重系数。实际计算过程中,h(i,j,t)函数可采用基于时间t的线性或非线性约束关系。本实验例中,h(i,j,t)函数采用如下形式:
其中,h=5,t∈[0,4],即:
可以看出,本实验例h(i,j,t)函数符合发生时间越久远,对当前及未来影响越小的原则。
步骤8)根据大气中水凝物等粒子的运动规律及其在雷达回波图像上所表征出的特点,构建一个灰阶图像变化预测模型,该模型的计算方法为:
且
其中,Δt表示预测的时间长度,Δt的取值以时间上两两相邻的基数据时间间隔为基本单位。举例来说,雷达某种探测方式下,基数据时间间隔为6分钟,如果预测未来6分钟,则Δt=1,如果预测未来30分钟,则Δt=5。Δx和Δy分别表示未来Δt时间后,原来T时刻坐标(x,y)分别在x轴和y轴的偏移量。MV(x,y)x和MV(x,y)y分别表示MV(x,y)在x轴和y轴的速度分量。P(T+Δt,a,x+Δx,y+Δy)表示在未来Δt时刻,a仰角面上,坐标(x+Δx,y+Δy)处的灰阶值。
图4以图形化方式展示了本实施例中,Δt=1、Δt=5和Δt=10,即未来6分钟、30分钟和60分钟的移动速度矢量场图像。为了便于移动速度矢量的显示,所示图像采取了一定程度的抽稀。
步骤9)为了改善步骤8计算中,坐标(x+Δx,y+Δy)(x∈[0,X),y∈[0,Y))不能完全覆盖整个灰阶图像,进而造成图像发散的问题,再对步骤8计算得到的P(T+Δt,a,x+Δx,y+Δy)做适当的平滑处理,使得所预测的回波图像在各个局部区域有着更好的一致性。本实施例中,平滑处理采用了公开的高斯滤波算法,其高斯模板为:
图5以图形化方式展示了本实施例中,Δt=1、Δt=5和Δt=10,即未来6分钟、30分钟和60分钟的雷达回波灰阶图像,其中,每行左侧的图像为预测图像,右测为相同时间的实况图像,所有图像都经过了步骤3所述的转换。这些图像中各个像素的颜色深浅是由式10所计算的P(T+Δt,a,x+Δx,y+Δy)确定,且P(T+Δt,a,x+Δx,y+Δy)的值已经过步骤9高斯滤波计算。从图5可以看出,根据本发明进行降水预报与实际降水具有很好的一致性。
步骤10)按照与步骤3相反的过程,对P(T+Δt,a,x,y)的值域进行收缩调整,将各像素的取值还原到基数据的正常量级,计算方法为:
当前步骤中,P(T+Δt,a,x,y)与式10的P(T+Δt,a,x+Δx,y+Δy)是同一变量,只是坐标的表述方法有所不同。式10中Δx和Δy是用于表述当前任一坐标(x+Δx,y+Δy)相对于Img(T,a,x,y)的坐标(x,y)的偏移量。P(T+Δt,a,x,y)则表示T+Δt时刻,回波图像中任一坐标(x,y)处的灰阶强度值;类似地,P’(T+Δt,a,x,y)表示T+Δt时刻,任一坐标(x,y)处的雷达基本反射率数值。
步骤11)采用雷达气象学相关领域的Z-R关系,将P’(T+Δt,a,x,y)所表征的回波强度值转换为可降水量,计算方法为:
R(x,y)=α×P′(T+Δt,a,x,y)β (式12)
其中,α和β是Z-R关系式中的可变参量,通常根据降水类型的不同,取相应的经验值。
通过上述步骤1至步骤11的计算,实现基于天气雷达回波时序图像的降水预报。
以上所述,仅为本发明较佳的具体实施方式,但本发明的保护范围并不局限于此,任何熟悉本技术领域的技术人员在本发明揭露的技术范围内,可轻易想到的变化或替换,都应涵盖在本发明的保护范围之内。因此,本发明的保护范围应该以权利要求的保护范围为准。
Claims (3)
1.基于天气雷达回波时序图像的降水预报方法,其特征在于:包括以下步骤:
(1)读取当前时刻T及前N各时刻的雷达基数据;
(2)读取(1)中当前时刻T雷达基数据中的各个探测仰角的基本反射率数值;
(3)将步骤(2)中各个仰角的基本反射率数值的值域扩大到[0,255]之间;
(4)将步骤(3)中扩大值域后的数值转换为灰阶图像;
(5)重复步骤(2)至步骤(4),对步骤(1)中读取的所有雷达基数据进行灰阶图像转换;
(6)根据以下原理,依次对时间上两两相邻雷达基数据转换后的基本反射率灰阶图像进行回波主体移动特征计算,从而计算出灰阶图像各像素的移动速度矢量:图像中任意位置相邻的像素,其移动趋势近乎一致;以及图像中表征雷达回波主体的灰阶图像区域,其总体移动趋势基本一致;
(7)根据发生时间越久远,对当前及未来影响越小的原则,对步骤(6)计算出的移动速度矢量进行叠加计算;
(8)构建灰阶图像变化预测模型:
且
式中,Δt表示预测的时间长度,Δt的取值以时间上两两相邻的基数据时间间隔为基本单位;Δx和Δy分别表示未来Δt时间后,原来T时刻坐标(x,y)分别在x轴和y轴的偏移量;MV(x,y)为步骤(7)中叠加计算的结果;MV(x,y)x和MV(x,y)y分别表示MV(x,y)在x轴和y轴的速度分量;IMG(T,a,x,y)表示T时刻的雷达基数据计算出的a仰角灰阶图像在直角坐标(x,y)处的灰阶强度值;同理,IMG(T-t,a,x,y)表示基数据T-t时刻的雷达基数据计算出的a仰角灰阶图像在直角坐标(x,y)处的灰阶强度值;h(x,y,t)是为时域卷积函数;*是卷积运算符;IMG(T-t,a,x,y)*h(x,y,t)的物理意义为:T-t时刻的灰阶图像中,任一像素(x,y)处的灰度值在[0,N-1]时序中的t时刻对预测值P(T+Δt,a,x+Δx,y+Δy)的作用程度;
(9)对步骤(8)计算得到的P(T+Δt,a,x+Δx,y+Δy)进行平滑处理;
(10)对步骤(9)平滑处理后数值的值域进行收缩调整,将各像素的取值还原至雷达基数据的正常数量级;
(11)将步骤(10)还原后数值所表征的回波强度值转换为可降水量;
所述步骤(6)中图像中任意位置相邻的像素,其移动趋势是近乎一致的,用公式可表示为:
V(x,y,T-t)=V(x+Δx,y+Δy,T-t-1)(Δx→0,Δy→0) (式2)
式2中,V表示移动速度矢量,V(T-t,x,y)表示T-t时刻的基本反射率灰阶图像中直角坐标系下(x,y)处的移动速度矢量;根据所述步骤(6)中的原理,V(T-t,x,y)与V(T-t-1,x,y)的值近似相等,即V相邻坐标位置的梯度变化近似为0,由此建立以下公式:
根据图像中表征雷达回波主体的灰阶图像区域,其总体移动趋势基本一致的原理,由以下公式表示:
式4中,min[]表示遍历x方向和y方向所有像素后,求取数值最小的一项;
将式3和式4合并,构成用于求取基本反射率灰阶图像移动矢量的约束条件,如式5所示:
式5中,μ∈[0,1];
将式5简化为离散化形式:
其中,
式7中,D(i,j,T-t)x、D(i,j,T-t)y和D(i,j,T-t)t分别表示图像中的任一像素点(i,j)的灰度值在x、y和时间T轴方向上的变化率;采用如下的差分公式进行简化计算:
D(i,j,T-t)x=(IMG(T-t,a,i+1,j)-IMG(T-t,a,i,j)+IMG(T-t-1,a,i+1,j)-IMG(T-t-1,a,i,j)+IMG(T-t,a,i+1,j+1)-IMG(T-t,a,i,j+1)+IMG(T-t-1,a,i+1,j+1)-IMG(T-t-1,a,i,j+1))/4δx
D(i,j,T-t)y=(IMG(T-t,a,i+1,j)-IMG(T-t,a,i,j)+IMG(T-t-1,a,i+1,j)-IMG(T-t-1,a,i,j)+IMG(T-t,a,i+1,j+1)-IMG(T-t,a,i,j+1)+IMG(T-t-1,a,i+1,j+1)-IMG(T-t-1,a,i,j+1))/4δy
D(i,j,T-t)t=(IMG(T-t,a,i+1,j)-IMG(T-t,a,i,j)+IMG(T-t-1,a,i+1,j)-IMG(T-t-1,a,i,j)+IMG(T-t,a,i+1,j+1)-IMG(T-t,a,i,j+1)+IMG(T-t-1,a,i+1,j+1)-IMG(T-t-1,a,i,j+1))/4δt (式8)
式中,x和y分别表示灰阶图像中任一像素在横坐标和纵坐标的位置,(x,y)是灰阶图像中任一像素的坐标;大写的X和Y表示灰阶图像的长(横坐标方向),和宽(纵坐标方向),单位为像素数;δx、δy和δt分别表示灰阶图像中的像素在x轴、y轴和时间轴方向上的宽度;和的数学含义是灰阶图像中的像素点(i,j)在x和y轴方向上的移动特征分量;
所述步骤(7)中对移动速度矢量的叠加计算公式为:
其中,x0和y0分别表示每次叠加计算时灰阶图像中参与计算的空间邻域的大小,取值为大于1的整数;*是卷积运算符;h(i,j,t)是时域卷积函数,采用如下形式:
其中,h=5,t∈[0,4],即:
3.根据权利要求2所述的降水预报方法,其特征在于:所述步骤(11)进行可降水量转换的计算公式为:
R(x,y)=α×P′(T+Δt,a,x,y)β (式12)
式中,α和β是Z-R关系式中的可变参量。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810531177.2A CN108761461B (zh) | 2018-05-29 | 2018-05-29 | 基于天气雷达回波时序图像的降水预报方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201810531177.2A CN108761461B (zh) | 2018-05-29 | 2018-05-29 | 基于天气雷达回波时序图像的降水预报方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN108761461A CN108761461A (zh) | 2018-11-06 |
CN108761461B true CN108761461B (zh) | 2022-02-18 |
Family
ID=64003418
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201810531177.2A Active CN108761461B (zh) | 2018-05-29 | 2018-05-29 | 基于天气雷达回波时序图像的降水预报方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN108761461B (zh) |
Families Citing this family (9)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109298424A (zh) * | 2018-11-28 | 2019-02-01 | 广东电网有限责任公司 | 一种基于雷达拼图的杆塔电力线路天气监测方法和装置 |
CN110031823B (zh) * | 2019-04-22 | 2020-03-24 | 上海禾赛光电科技有限公司 | 可用于激光雷达的噪点识别方法以及激光雷达系统 |
CN111175744B (zh) * | 2019-09-20 | 2023-08-15 | 中国船舶工业系统工程研究院 | 一种雷达图像快速生成与缩放方法 |
CN110824479A (zh) * | 2019-10-28 | 2020-02-21 | 兰州大方电子有限责任公司 | 一种用于短临预报的雷达数据处理方法 |
CN110879828B (zh) * | 2019-11-20 | 2020-11-24 | 上海眼控科技股份有限公司 | 雷达回波图的处理方法、装置、计算机设备及存储介质 |
CN111337928B (zh) * | 2020-03-20 | 2021-09-28 | 厦门市气象台(厦门市海洋气象台、海峡气象开放实验室) | 一种雷达回波移动信息计算方法和装置 |
CN111458769B (zh) * | 2020-05-26 | 2021-05-28 | 南京大学 | 用于输电线路环境气象数据预测的方法及系统 |
CN111814961A (zh) * | 2020-07-03 | 2020-10-23 | 深圳市赑玄阁科技有限公司 | 基于生成式对抗网络的降水预报的方法 |
CN113640769B (zh) * | 2021-08-27 | 2023-06-13 | 南京信息工程大学 | 基于深度神经网络的天气雷达基本反射率预测方法 |
Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105388801A (zh) * | 2014-08-27 | 2016-03-09 | 瑞萨电子株式会社 | 控制系统、中继装置和控制方法 |
CN105738873A (zh) * | 2015-11-16 | 2016-07-06 | 象辑知源(武汉)科技有限公司 | 天气雷达回波图像的处理方法及处理装置 |
Family Cites Families (14)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
WO2007005328A2 (en) * | 2005-06-30 | 2007-01-11 | Massachusetts Institute Of Technology | Weather radar echo tops forecast generation |
JP6281460B2 (ja) * | 2014-09-24 | 2018-02-21 | 株式会社デンソー | 物体検出装置 |
CN104977584B (zh) * | 2015-06-29 | 2017-08-15 | 深圳市气象台 | 一种对流天气的临近预报方法及系统 |
CN105717491B (zh) * | 2016-02-04 | 2018-09-21 | 象辑知源(武汉)科技有限公司 | 天气雷达回波图像的预测方法及预测装置 |
EP3252501B1 (en) * | 2016-06-03 | 2022-06-08 | Veoneer Sweden AB | Enhanced object detection and motion state estimation for a vehicle environment detection system |
CN105891078B (zh) * | 2016-06-08 | 2017-05-31 | 中国气象局乌鲁木齐沙漠气象研究所 | 基于风廓线雷达的沙尘暴沙尘质量浓度定量反演估算方法 |
US10514770B2 (en) * | 2016-06-17 | 2019-12-24 | Texas Instruments Incorporated | Hidden Markov model-based gesture recognition with FMCW radar |
US10078155B2 (en) * | 2016-06-24 | 2018-09-18 | Climacell Inc. | Real-time precipitation forecasting system |
CN105974418B (zh) * | 2016-07-08 | 2018-05-08 | 南京信息工程大学 | 基于天气雷达反射率特征匹配的降水估测方法 |
CN106324709B (zh) * | 2016-10-21 | 2019-01-01 | 中国人民解放军理工大学 | 微波链路、雨滴谱仪、雨量计与天气雷达多源联合的降雨场重构方法 |
CN106872981A (zh) * | 2017-02-17 | 2017-06-20 | 水利部南京水利水文自动化研究所 | 雨量雷达的降水强中心跟踪与预报方法 |
CN106908782B (zh) * | 2017-02-23 | 2019-08-06 | 公安部第三研究所 | 基于水面状态连续成像系统的波浪传播方向的提取方法 |
CN107038507A (zh) * | 2017-06-06 | 2017-08-11 | 孝感市青谷信息科技有限公司 | 一种测雨卫星雨量智能测量和预报方法 |
CN107703564B (zh) * | 2017-10-13 | 2020-04-14 | 中国科学院深圳先进技术研究院 | 一种降雨预测方法、系统及电子设备 |
-
2018
- 2018-05-29 CN CN201810531177.2A patent/CN108761461B/zh active Active
Patent Citations (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105388801A (zh) * | 2014-08-27 | 2016-03-09 | 瑞萨电子株式会社 | 控制系统、中继装置和控制方法 |
CN105738873A (zh) * | 2015-11-16 | 2016-07-06 | 象辑知源(武汉)科技有限公司 | 天气雷达回波图像的处理方法及处理装置 |
Also Published As
Publication number | Publication date |
---|---|
CN108761461A (zh) | 2018-11-06 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN108761461B (zh) | 基于天气雷达回波时序图像的降水预报方法 | |
Wang et al. | Automatic laser profile recognition and fast tracking for structured light measurement using deep learning and template matching | |
CN108073864B (zh) | 目标对象检测方法、装置及系统和神经网络结构 | |
CN104299229B (zh) | 一种基于时空域背景抑制的红外弱小目标检测方法 | |
CN102609904B (zh) | 双变量非局部平均滤波x射线图像消噪方法 | |
CN103024349B (zh) | 一种基于稀疏限制的mle视频目标跟踪方法 | |
CN108222749B (zh) | 一种基于图像分析的智能自动门控制方法 | |
CN104820997B (zh) | 一种基于分块稀疏表达与hsv特征融合的目标跟踪方法 | |
CN103729854B (zh) | 一种基于张量模型的红外弱小目标检测方法 | |
CN106338733B (zh) | 基于蛙眼视觉特性的前视声呐目标跟踪方法 | |
CN105809715B (zh) | 一种基于帧间累计变化矩阵的视觉运动目标检测方法 | |
CN103020918B (zh) | 基于形状自适应邻域均值的非局部均值去噪方法 | |
CN109993775B (zh) | 基于特征补偿的单目标跟踪方法 | |
CN104574439A (zh) | 一种融合卡尔曼滤波与tld算法的目标跟踪方法 | |
Jing et al. | AENN: A generative adversarial neural network for weather radar echo extrapolation | |
CN107563397A (zh) | 一种卫星云图中云团运动矢量的确定方法 | |
CN103077537B (zh) | 基于l1正则化的实时运动目标跟踪的新方法 | |
CN110070565A (zh) | 一种基于图像叠加的船舶轨迹预测方法 | |
CN103473755B (zh) | 基于变化检测的sar图像稀疏去噪方法 | |
CN106651908A (zh) | 一种多运动目标跟踪方法 | |
CN105354863A (zh) | 基于特征滤波和快速运动检测模板预测的自适应尺度图像序列目标跟踪方法 | |
CN105631899A (zh) | 一种基于灰度纹理特征的超声图像运动目标跟踪方法 | |
CN101482969A (zh) | 基于同质点计算的sar图像去斑方法 | |
Zhang et al. | Infrared small target tracking based on sample constrained particle filtering and sparse representation | |
CN110706208A (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 |