CN106526558B - 基于多普勒天气雷达资料的阵风锋自动识别方法 - Google Patents

基于多普勒天气雷达资料的阵风锋自动识别方法 Download PDF

Info

Publication number
CN106526558B
CN106526558B CN201610858401.XA CN201610858401A CN106526558B CN 106526558 B CN106526558 B CN 106526558B CN 201610858401 A CN201610858401 A CN 201610858401A CN 106526558 B CN106526558 B CN 106526558B
Authority
CN
China
Prior art keywords
doubtful
curve
gust front
region
gust
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.)
Expired - Fee Related
Application number
CN201610858401.XA
Other languages
English (en)
Other versions
CN106526558A (zh
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.)
Tianjin University
Original Assignee
Tianjin University
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 Tianjin University filed Critical Tianjin University
Priority to CN201610858401.XA priority Critical patent/CN106526558B/zh
Publication of CN106526558A publication Critical patent/CN106526558A/zh
Application granted granted Critical
Publication of CN106526558B publication Critical patent/CN106526558B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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
    • G01S7/415—Identification of targets based on measurements of movement associated with the target
    • 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/28—Details of pulse systems
    • G01S7/285—Receivers
    • G01S7/292—Extracting wanted echo-signals
    • G01S7/2923—Extracting wanted echo-signals based on data belonging to a number of consecutive radar periods
    • 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

本发明公开了基于多普勒天气雷达资料的阵风锋自动识别方法，所述方法包括以下步骤：1)依据阵风锋的雷达表现特征利用局部二值化算法提取出弱窄带回波疑似区域；2)对弱窄带回波疑似区域进行分割、连接和筛选处理，得到弱脊状带对应的骨架图像；3)由当前时刻和前一时刻的两幅低仰角雷达图像得到光流场，将步骤2)得到的弱脊状带对应的骨架图像中前后时刻匹配的骨架拟定为疑似阵风锋，根据该疑似阵风锋的位置和速度与风暴单体的位置和速度的关系及该疑似阵风锋的走向与速度的关系识别出阵风锋。本方法实现了阵风锋的自动检测，对灾害进行及时的预警，减少了经济损失和人员伤亡；并通过实验验证了本方法的有效性。

Description

基于多普勒天气雷达资料的阵风锋自动识别方法
技术领域
本发明涉及气象学领域，特别涉及一种基于多普勒天气雷达资料的阵风锋自动识别方法。
背景技术
天气雷达是对强对流天气进行监测和预警的主要工具之一。天气雷达发射脉冲形式的电磁波，当电磁波遇到降水物质(雪花、雨滴和冰雹等)，大部分能量继续前进，而有一小部分能量被降水物质向四面八方散射，向后散射的能量被雷达接收。
在对流单体的成熟阶段，冷性下沉气流作为一股冷空气，在近地面的底层向外扩展，与单体运动前方的暖湿气流交汇形成锋面。这种现象反映在雷达反射率图上，通常可以观测到一条若隐若现、粗细不等、强度较弱且取值不定的弱脊带状区域，气象上称之为阵风锋或雷暴的出流边界。阵风锋常引起气压突变、风速风向剧烈变化以及气温骤降等强烈天气现象，会导农作物倒伏、大树或树枝折断、威胁飞机的起降，严重影响人民大众的生命财产安全。
从图像角度看，做阵风锋区域的横剖面，其反射率值的分布呈脊状，但由于区域自身反射率值普遍较低，致使由区域两侧过渡到区域外围的反射率值的梯度变化微弱，且时常出现部分区域段混于大面积的弱回波区域之中的现象，致使原本就不太强的梯度变化在这部分区域段基本消失。这就使常规的基于边缘和区域图像分割方法用于检测这种阵风锋区域时奏效甚微,另外，弱回波区域中通常含有脊状非阵风锋带的干扰，这些干扰甚至与真正的阵风锋相邻或相交。使识别阵风锋方法的空识率与击中率之间的矛盾愈发突出。
在基于多普勒雷达数据的阵风锋识别算法中，Delanoy和Troxel[1]较早提出了基于模糊理论的函数模板相关法,该方法将模板尺度固定为17×7，为使检测算法与阵风锋的方向无关，该模板需围绕其中心多次旋转。Osama Alkhouli和Victor DeBrunner[2]应用熵函数模板法自动识别边界层辐合线，且方法无需旋转模板。郑佳峰等[3]结合速度场与强度场数据检测阵风锋。在速度场数据中，检测辐合带；在反射率强度场中，该方法主要使用的是双向梯度法寻找弱窄带回波，两者结合确定阵风锋。
发明人在实现本发明的过程中，发现现有技术中至少存在以下缺点和不足：
首先，人为观测费时费力，会影响到预报的时效性。文献[1]只适合检测宽度不大于3与周边区域存在较大梯度的阵风锋。文献[2]应用熵函数模板法提高了检测速度，但仅对典型的阵风锋有效，对于部分阵风锋区域融入背景的情况，会使检测结果出现断裂。文献[3]的方法在速度场数据中出现无效速度的区域时无法发挥；在反射率强度场时，双向梯度法限定窄带回波与母体回波距离。至今为止，在对阵风锋检测的文献中尚未见到阵风锋与其他弱窄带回波碰撞产生交叠现象的解决方案，也未见到在阵风锋的检测中出现断裂现象的处理方法。
[参考文献]
[1]Delanoy R L,Troxel S W.A machine intelligent gust front algorithmfor doppler weather radars[C].Contributions to the American MeteorologicalSociety’s 26th International Conference on Radar Meteorology.1993:9.
[2]Alkhouli O,DeBrunner V.Gust front detection in weather radarimages by entropy matched functional template[C].Image Processing,2005.ICIP2005.IEEE International Conference on.IEEE,2005,1:I-645-8.
[3]郑佳锋,张杰,朱克云,等.阵风锋自动识别与预警[J].应用气象学报,2013,24(1):117.
[4]Ojala T,M,Harwood D.A comparative study of texturemeasures with classification based on featured distributions[J].Patternrecognition,1996,29(1):51-59.
[5]G.Two-frame motion estimation based on polynomialexpansion[M].Image Analysis.Springer Berlin Heidelberg,2003:363-370.
发明内容
本发明公开了一种基于多普勒天气雷达资料的阵风锋自动识别方法，可以解决的技术问题包括：自动识别不同宽度的阵风锋；能识别与周边区域的梯度较小或部分域融入背景的阵风锋；能分离与其他弱窄带回波碰撞产生交叠的阵风锋；对阵风锋检测中出现断裂现象进行有效连接的处理；实现准确、完整自动识别阵风锋的目的。
为了解决上述技术问题，本发明提出的一种基于多普勒天气雷达资料的阵风锋自动识别方法，包括以下步骤：
步骤一、依据阵风锋的雷达表现特征，利用局部二值化算法提取出弱窄带回波疑似区域；步骤如下：
1-1)设一幅低仰角雷达图像的大小为N×N，对中间[N-(2n+1)×(2n+1)]×[N-(2n+1)×(2n+1)]区域内的像素点pij进行其是否属于弱脊状区域的判断；
1-2)求出区域[N-(2n+1)×(2n+1)]×[N-(2n+1)×(2n+1)]内大于等于35dBZ的连通区域ωi，i＝1,2，…s，计算该连通区域ωi的面积si及该连通区域ωi外包矩形长度li，对于面积si小于S1或外包矩形长度li小于L1的连通区域进行标记，得到ωj′，j＝1,2，…m，m≤s；
1-3)若像素点pij的反射率值f(i,j)∈x1或f(i,j)∈x2且f(i,j)∈ωj′，其中x1＝[5,35)dBZ，x2＝[35,40)dBZ，执行1-4)，否则，执行1-5)
1-4)在以该像素点pij为中心的区域根据式(1)做卷积运算，得到卷积运算结果g1(i,j)和g2(i,j)，当f(i,j)≥g1(i,j)且f(i,j)>g2(i,j)，则认为像素点pij属于弱窄带回波疑似区域Ω，并将该像素点pij置为前景，否则，置为背景；
1-5)将像素点pij置为背景，至此将低仰角雷达图像转化为一张二值图；
1-6)计算通过步骤1-4)和步骤1-5)所形成的二值图中各个连通区域的面积和外接矩形长度，将连通区域面积小于S2或外接矩形长度小于L2的连通区域置为背景；从而提取出弱窄带回波疑似区域；
步骤二、对弱窄带回波疑似区域进行分割、连接和筛选处理，得到弱脊状带对应的骨架图像；步骤如下：
2-1)对步骤一提取出的弱窄带回波疑似区域的轮廓进行除毛刺处理，再进行细化得到区域的骨架图像A；
2-2)在上述骨架图像A中的骨架交叉点和折点处断开骨架，得到骨架图像B，包括：
依照折点的特点通过计算骨架上某点两侧切线夹角识别出折点，即沿骨架行进，若由某点到其前侧第n个点的向量与该点到其后侧第n个点的向量的夹角小于135度，则认为该点为折点；
利用局部二值模式算子检测出端点和交叉点，即针对步骤2-1)得到的骨架图像，骨架图像中取值为1的点p为可能的端点或交叉点，以该点p为中心，考察其3×3区域及5×5区域边界的取值分布；其中，5×5区域边界点若取值为1，但在5×5范围内与区域中心不连通，则将该点置为0，分别从上述3×3区域及5×5区域的左上角开始沿逆时针方向形成8位01链码和16位01链码来描述所述两个区域边界的取值分布；然后分别沿两个区域边界循环一周，记录依次取值变化的次数n3(p)和n5(p)，若n3(p)＝2，则点p为端点；若n3(p)≥6或n5(p)≥6，则点p为交叉点；
在上述折点和交叉点处断开骨架，从而得到骨架图像B；
2-3)通过端点匹配的方法将骨架图像B中任意两段不连通的曲线连接；
根据阵风锋走向平缓的特点，设：曲线li的端点A和曲线lj(j≠i)的端点B的匹配条件如下：
式(2)中，L＝30km，Φ1＝Φ2＝Φ3＝-0.7，li和lj表示曲线i和j的长度，指由端点A沿曲线i在端点A处切线的方向指向远离该曲线的向量，指由端点B沿曲线j在端点B处切线的方向指向远离该曲线的向量；
当端点A仅与一个端点符合匹配条件时，则该端点为端点B，连接端点A和端点B；
当端点A与多个端点符合匹配条件时，则将其中值最小的端点作为端点B，连接端点A和端点B，其中，端点C是与端点A符合匹配条件的任意一端点，指由端点C沿其所在曲线在端点C处切线的方向指向远离该曲线的向量；
从而形成骨架图像C；
2-4)判断低仰角雷达图像是否存在单侧脊区域，若存在，则剔除单侧脊区域的骨架；
对骨架图像C中的每条曲线分别进行PCA处理，至少得到第二主成分方向的单位向量e2；
当距离k分别为2,3,4,5km时，设变量αk的初值为0，对每条曲线上m个点中每一点pi，i＝1,2...m,求得所述曲线一侧的位置qi，qi＝pi+k e2对于每一点pi计算：
式(3)中，f(x)表示点x的反射率值；
设，所述曲线一侧的加权比例
当距离k分别为-2,-3,-4,-5km时，设，所述曲线另一侧的加权比例为β0，根据上述过程求得加权比例β0；
若α0<Th且β0<Th，则认为该曲线与弱窄带回波疑似区域中对应的区域是弱脊状带，否则，该曲线与弱窄带回波疑似区域中对应的区域是单侧脊区域，滤除该曲线；其中，Th＝0.75；形成骨架图像D；
2-5)判断低仰角雷达图像是否存在虚线回波，若存在，则剔除虚线回波的骨架；
用2-4)中对骨架图像C中的每条曲线的PCA处理结果，能得到骨架图像D中每条曲线的第一主成分特征值λ1、第一主成分方向的单位向量e1和第二主成分特征值λ2，当第二主成分特征值λ2/第一主成分特征值λ1<0.05，且雷达中心点到曲线所拟合的直线的距离小于dl，则认为该曲线与低仰角雷达图像中对应的区域是虚线回波，滤除该曲线；其中dl＝3km；最终得到的骨架图像即为弱脊状带对应的骨架图像；
步骤三、由当前时刻和前一时刻的两幅低仰角雷达图像得到光流场，将步骤二得到的弱脊状带对应的骨架图像中前后时刻匹配的骨架拟定为疑似阵风锋，根据该疑似阵风锋的位置和速度与风暴单体的位置和速度的关系及该疑似阵风锋的走向与速度的关系识别出阵风锋；步骤如下：
3-1)利用光流法由当前时刻和前一时刻的两幅低仰角雷达图像得到光流场；
3-2)利用光流场信息将前一时刻的弱脊状带对应的骨架图像的曲线移动到由步骤二得到的当前时刻骨架图像中的对应位置，同时满足下述条件一和条件二的曲线为当前时刻和前一时刻的同一条弱窄带回波区域对应，即为疑似阵风锋；
条件一：前后时刻弱窄带回波区域的交叠长度大于30％；
条件二：利用PCA处理后得到的前后时刻弱窄带回波区域的第一主成分方向夹角小于30度；
3-3)若前一时刻弱窄带回波某端比当前时刻弱窄带回波长10km以上时，则利用多出的此段曲线对步骤3-2)中认定的疑似阵风锋进行延长补全；
3-4)根据该疑似阵风锋的位置和速度与风暴单体的位置和速度的关系及该疑似阵风锋的走向与速度的关系识别是否具有阵风锋，包括
判断当前时刻疑似阵风锋走向与疑似阵风锋速度方向锐角大于45度；
同时，当前时刻疑似阵风锋与风暴单体的速度方向与位置关系满足下述条件(1)至条件(5)中的一条；则认定当前时刻疑似阵风锋为阵风锋；并针对下一时刻低仰角雷达图像，直接将步骤3-3)补全后的该时刻的疑似阵风锋直接认定为阵风锋；
(1)当当前时刻疑似阵风锋仅位于风暴单体运动的前侧40km内时，疑似阵风锋速度方向与风暴单体移动方向夹角小于30度；
(2)当当前时刻疑似阵风锋仅位于风暴单体运动的右侧40km内时，疑似阵风锋速度方向与风暴单体移动方向右转90度后的方向夹角小于30度；
(3)当当前时刻疑似阵风锋仅位于风暴单体运动的左侧40km内时，疑似阵风锋速度方向与风暴单体移动方向左转90度后的方向夹角小于30度；
(4)当当前时刻疑似阵风锋既位于风暴单体前侧40km内且位于风暴单体右侧40km内时，疑似阵风锋速度方向位于风暴单体移动方向与风暴单体移动方向右转90度后的方向之间；
(5)当当前时刻疑似阵风锋既位于风暴单体前侧40km内且位于风暴单体左侧40km内时，疑似阵风锋速度方向位于风暴单体移动方向与风暴单体移动方向左转90度后的方向之间；
上述步骤3-2)中，利用光流场信息将前一时刻的弱脊状带对应的骨架图像的曲线移动到由步骤二得到的当前时刻骨架图像中的对应位置后，将没有同时满足条件一和条件二的前一时刻的弱脊状带对应的骨架图像的曲线保留两个体扫叠加到下一循环的疑似阵风锋匹配过程中的前一时刻的骨架图像中；
上述步骤3-3)补全后的疑似阵风锋用于替换下一循环中前一时刻对应位置的曲线。
与现有技术相比，本发明的有益效果是：
根据阵风锋在雷达强度场中的回波特性和几何特征，首先通过局部区域二值化等方法提取弱窄带回波，再根据阵风锋与风暴单体之间的关系从弱窄带回波中筛选出阵风锋，实现了准确、完整自动检测阵风锋，对灾害进行及时的预警，减少了经济损失和人员伤亡；并通过实验验证了本方法的有效性。
附图说明
图1(a)至图1(e)为天气雷达反射率图图例与低仰角反射率图上的阵风锋，其中，图1(a)至图1(e)的左侧显示了图例(下同)，右侧显示了低仰角反射率图上的阵风锋；
图2为卷积模板，其中每个方格代表一个像素点，也就是雷达图中1km*1km的区域(下同)，图中的黑色点代表要输出的像素方格，大的加粗框表示一个卷积操作数矩阵，小的加粗框代表另一个卷积操作数矩阵。
图3(a)至图3(f)为局部二值化算法提取的前景区域示意图，其中，图3(a)显示了单侧脊，图3(b)显示了与其他弱脊状区域交叠的阵风锋，图3(c)显示了虚线回波，图3d)显示了非阵风锋的边界层辐合线，图3(e)显示了非边界层辐合线，图3(f)显示了真正的阵风锋；
图4(a)至图4(c)为骨架的局部区域，其中，图4(a)中3×3区域的中心点为骨架的端点，图4(b)中5×5区域的中心点为骨架的交叉点，图4(c)中5×5区域的中心点既不是骨架的端点也不是骨架的交叉点；
图5为本发明提供的一种用于气象学中寻找阵风锋方法的流程图；
图6(a)和图(b)为本发明提供的测试效果图。
具体实施方式
下面结合附图和具体实施例对本发明技术方案作进一步详细描述，所描述的具体实施例仅对本发明进行解释说明，并不用以限制本发明。
本发明提供了一种基于多普勒天气雷达资料的阵风锋自动识别方法，本方法能自动检测出阵风锋，对灾害进行及时的预警，减少经济损失和人员伤亡。
实施例：在多普勒雷达低仰角反射率图中观测到的阵风锋如图1(a)至图1(e)所示，其图像特征包括但不限于：
1)阵风锋通常表现为在强回波外围的一条弱窄带回波，一般位于雷暴母体运动方向的前端或侧面，它与雷暴母体脱离或部分粘连，见图1(a)；
2)宽度变化范围：2km-10km；
3)长度变化范围：40km-300km；
4)反射率取值不定：在同一时刻，一条阵风锋区域上的回波强度通常不同，见图1(b)，不同时刻，同一条阵风锋上的回波强度取值和分布也会变化(如图1(b)和图1(c)，日期为2006.06.14其中图1(b)时间为09:49和图1(c)时间为09:31)，但多集中在10dBz—30dBz这个范围；
5)脊特征：阵风锋上的反射率值整体上高于其两侧的反射率值，但有时差值不大，呈现很弱的脊特征，且不排除局部区域或局部点这种特征消失；
6)不连续性：当阵风锋与非降水回波混合在一起时，可能会被较高反射率值区域或者较低的反射率值区域“污染”而发生断裂，见图1(d)；
7)交叠性：一条阵风锋可能与其他弱窄带回波相交甚至发生碰撞，使两者相连或相交，见图1(e)。
本发明基于多普勒天气雷达资料的阵风锋自动识别方法，如图5所示，包括以下步骤：
101、依据阵风锋的雷达表现特征，利用局部二值化算法提取出弱窄带回波疑似区域；具体内容如下：
1-1)设一幅低仰角雷达图像的大小为N×N，对中间[N-(2n+1)×(2n+1)]×[N-(2n+1)×(2n+1)]区域内的像素点pij进行其是否属于弱脊状区域的判断；
1-2)求出区域[N-(2n+1)×(2n+1)]×[N-(2n+1)×(2n+1)]内大于等于35dBZ的连通区域ωi，i＝1,2，…s，计算该连通区域ωi的面积si及该连通区域ωi外包矩形长度li，对于面积si小于S1或外包矩形长度li小于L1的连通区域进行标记，得到ωj′，j＝1,2，…m，m≤s；
1-3)若像素点pij的反射率值f(i,j)∈x1或f(i,j)∈x2且f(i,j)∈ωj′，其中x1＝[5,35)dBZ，x2＝[35,40)dBZ，执行1-4)，否则，执行1-5)
1-4)则用图2所示的模板在在以该像素点pij为中心的区域根据式(1)做卷积运算，得到卷积运算结果g1(i,j)和g2(i,j)，当f(i,j)≥g1(i,j)且f(i,j)>g2(i,j)，则认为像素点pij属于弱窄带回波区域Ω，并将该像素点pij置为前景，否则，置为背景；
1-5)将像素点pij置为背景，至此将低仰角雷达图像转化为一张二值图；
1-6)计算通过步骤1-4)和步骤1-5)所形成的二值图中各个连通区域的面积和外接矩形长度，将连通区域面积小于S2或外接矩形长度小于L2的连通区域置为背景；其中参数k＝3，S1＝15，L1＝5，S2＝45，L2＝15；图3(a)、图3(b)、图3(c)、图3(d)、图3(e)和图3(f)显示了局部二值化算法提取的前景区域示意图，图3(a)显示了单侧脊，图3(b)显示了与其他弱脊状区域交叠的阵风锋，图3(c)显示了虚线回波，图3(d)显示了非阵风锋的边界层辐合线，图3(e)显示了非边界层辐合线，图3(f)显示了真正的阵风锋；从而提取出弱窄带回波疑似区域。
102、对弱窄带回波疑似区域进行分割、连接和筛选处理，得到弱脊状带对应的骨架图像；具体内容如下：
2-1)对步骤一提取出的弱窄带回波疑似区域的轮廓进行除毛刺处理，再进行细化得到区域的骨架图像A；
2-2)在上述骨架图像A中的骨架交叉点和折点处断开骨架，得到骨架图像B，包括：
依照折点的特点通过计算骨架上某点两侧切线夹角识别出折点，即沿骨架行进，若由某点到其前侧第n个点的向量与该点到其后侧第n个点的向量的夹角小于135度，则认为该点为折点，本例中n＝3；
利用局部二值模式(Local Binary Pattern，LBP)[4]算子检测出端点和交叉点，即针对步骤2-1)得到的骨架图像，骨架图像中取值为1的点p为可能的端点或交叉点，以该点p为中心，考察其3×3区域及5×5区域边界的取值分布；其中，5×5区域边界点若取值为1，但在5×5范围内与区域中心不连通，如图4(c)中覆盖较深(第4行最右边)的方格，则将该点置为0，分别从上述3×3区域及5×5区域的左上角开始沿逆时针方向形成8位01链码和16位01链码来描述所述两个区域边界的取值分布；然后分别沿两个区域边界循环一周，记录依次取值变化的次数n3(p)和n5(p)，若n3(p)＝2，如图4(a)，则点p为端点；若n3(p)≥6或n5(p)≥6，如图4(b)，则点p为交叉点；
在上述折点和交叉点处断开骨架，从而得到骨架图像B；
2-3)通过端点匹配的方法将骨架图像B中任意两段不连通的曲线连接；
根据阵风锋走向平缓的特点，设：曲线li的端点A和曲线lj(j≠i)的端点B的匹配条件如下：
式(2)中，L＝30km，Φ1＝Φ2＝Φ3＝-0.7，li和lj表示曲线i和j的长度，指由端点A沿曲线i在端点A处切线的方向指向远离该曲线的向量，指由端点B沿曲线j在端点B处切线的方向指向远离该曲线的向量；
当端点A仅与一个端点符合匹配条件时，则该端点为端点B，连接端点A和端点B；
当端点A与多个端点符合匹配条件时，则将其中值最小的端点作为端点B，连接端点A和端点B，其中，端点C是与端点A符合匹配条件的任意一端点，指由端点C沿其所在曲线在端点C处切线的方向指向远离该曲线的向量；
从而形成骨架图像C；
2-4)判断低仰角雷达图像是否存在单侧脊区域，若存在，则剔除单侧脊区域的骨架；
对骨架图像C中的每条曲线分别进行PCA处理，至少得到第二主成分方向的单位向量e2；
当距离k分别为2,3,4,5km时，设变量αk的初值为0，对每条曲线上m个点中每一点pi，i＝1,2...m,求得所述曲线一侧的位置qi，qi＝pi+k e2对于每一点pi计算：
式(3)中，f(x)表示点x的反射率值；
设，所述曲线一侧的加权比例
当距离k分别为-2,-3,-4,-5km时，设，所述曲线另一侧的加权比例为β0，根据上述过程求得加权比例β0；
若α0<Th且β0<Th，则认为该曲线与弱窄带回波疑似区域中对应的区域是弱脊状带，否则，该曲线与弱窄带回波疑似区域中对应的区域是单侧脊区域，滤除该曲线；其中，Th＝0.75；形成骨架图像D；
2-5)判断低仰角雷达图像是否存在虚线回波，若存在，则剔除虚线回波的骨架；
用2-4)中对骨架图像C中的每条曲线的PCA处理结果，能得到骨架图像D中每条曲线的第一主成分特征值λ1、第一主成分方向的单位向量e1和第二主成分特征值λ2，当第二主成分特征值λ2/第一主成分特征值λ1<0.05，且雷达中心点到曲线所拟合的直线的距离小于dl，则认为该曲线与低仰角雷达图像中对应的区域是虚线回波，滤除该曲线；其中dl＝3km；最终得到的骨架图像即为弱脊状带对应的骨架图像。
103、由当前时刻和前一时刻的两幅低仰角雷达图像得到光流场，将弱脊状带对应的骨架图像中前后时刻匹配的骨架拟定为疑似阵风锋，根据该疑似阵风锋的位置和速度与风暴单体的位置和速度的关系及该疑似阵风锋的走向与速度的关系识别出阵风锋，具体内容如下：
3-1)利用光流法[5]由当前时刻和前一时刻的两幅低仰角雷达图像得到光流场，
3-2)利用光流场信息将前一时刻的弱脊状带对应的骨架图像的曲线移动到由步骤二得到的当前时刻骨架图像中的对应位置，同时满足下述条件一和条件二的曲线为当前时刻和前一时刻的同一条弱窄带回波区域对应，即为疑似阵风锋；
条件一：前后时刻弱窄带回波区域的交叠长度大于30％；
条件二：利用PCA处理后得到的前后时刻弱窄带回波区域的第一主成分方向夹角小于30度；
3-3)若前一时刻弱窄带回波某端比当前时刻弱窄带回波长10km以上时，则利用多出的此段曲线对步骤3-2)中认定的疑似阵风锋进行延长补全；
3-4)根据该疑似阵风锋的位置和速度与风暴单体的位置和速度的关系及该疑似阵风锋的走向与速度的关系识别是否具有阵风锋，包括
判断当前时刻疑似阵风锋走向与疑似阵风锋速度方向所夹锐角大于45度；
同时，当前时刻疑似阵风锋与风暴单体的速度方向与位置关系满足下述条件(1)至条件(5)中的一条；则认定当前时刻疑似阵风锋为阵风锋；并针对下一时刻低仰角雷达图像，直接将步骤3-3)补全后的该时刻的疑似阵风锋直接认定为阵风锋；
(1)当当前时刻疑似阵风锋仅位于风暴单体运动的前侧40km内时，疑似阵风锋速度方向与风暴单体移动方向夹角小于30度；
(2)当当前时刻疑似阵风锋仅位于风暴单体运动的右侧40km内时，疑似阵风锋速度方向与风暴单体移动方向右转90度后的方向夹角小于30度；
(3)当当前时刻疑似阵风锋仅位于风暴单体运动的左侧40km内时，疑似阵风锋速度方向与风暴单体移动方向左转90度后的方向夹角小于30度；
(4)当当前时刻疑似阵风锋既位于风暴单体前侧40km内且位于风暴单体右侧40km内时，疑似阵风锋速度方向位于风暴单体移动方向与风暴单体移动方向右转90度后的方向之间；
(5)当当前时刻疑似阵风锋既位于风暴单体前侧40km内且位于风暴单体左侧40km内时，疑似阵风锋速度方向位于风暴单体移动方向与风暴单体移动方向左转90度后的方向之间；
上述步骤3-2)中，利用光流场信息将前一时刻的弱脊状带对应的骨架图像的曲线移动到由步骤二得到的当前时刻骨架图像中的对应位置后，将没有同时满足条件一和条件二的前一时刻的弱脊状带对应的骨架图像的曲线保留两个体扫叠加到下一循环的疑似阵风锋匹配过程中的前一时刻的骨架图像中；
上述步骤3-3)补全后的疑似阵风锋用于替换下一循环中前一时刻对应位置的曲线。
下面以具体的实验来验证本发明实施例提供的一种基于多普勒天气雷达资料的阵风锋自动识别方法的可行性，测试样本由中国天津气象台提供，分两部分验证，详见下文描述：
图6(a)和图6(b)为本发明提供的测试效果图，其中，图6(b)标注的白色区域即最后自动检测的真实阵风锋，第一部分是天津地区含有阵风锋的6个过程115组雷达基数据。第二部分是天津2012年6月含有强对流天气的7天1680个基数据，其中包含5个阵风锋过程的88个基数据。利用击中率POD、虚警率FAR和临界成功指数CSI对检验结果进行评价(见表1和表2)。
表1第一部分样本数和识别情况
表1分别描述了6个阵风锋过程的识别情况，对其进行统计后得到样本总数为115个，成功识别样本数93个，未识别样本数22，误识别样本数3，得出击中率、虚警率和临界成功指数分别为80.87％,3.13％,78.81％。
表2第二部分2012年6月天津样本数和识别情况
表2分别描述了7天强对流天气的识别情况，击中了88个阵风锋中的51个，在1592个数据中空报出了15个阵风锋，得出击中率、虚警率和临界成功指数分别为57.95％,22.73％,49.51％。
本领域技术人员可以理解附图只是一个优选实施例的示意图，上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种基于多普勒天气雷达资料的阵风锋自动识别方法，其特征在于，包括以下步骤：
步骤一、依据阵风锋的雷达表现特征，利用局部二值化算法提取出弱窄带回波疑似区域；步骤如下：
1-1)设一幅低仰角雷达图像的大小为N×N，对中间[N-(2n+1)×(2n+1)]×[N-(2n+1)×(2n+1)]区域内的像素点pij进行其是否属于弱脊状区域的判断；
1-2)求出区域[N-(2n+1)×(2n+1)]×[N-(2n+1)×(2n+1)]内大于等于35dBZ的连通区域ωi，i＝1,2，…s，计算该连通区域ωi的面积si及该连通区域ωi外包矩形长度li，对于面积si小于S1或外包矩形长度li小于L1的连通区域进行标记，得到ωj′，j＝1,2，…m，m≤s；
1-3)若像素点pij的反射率值f(i,j)∈x1或f(i,j)∈x2且f(i,j)∈ωj′，其中x1＝[5,35)dBZ，x2＝[35,40)dBZ，执行1-4)，否则，执行1-5)
1-4)在以该像素点pij为中心的区域根据式(1)做卷积运算，得到卷积运算结果g1(i,j)和g2(i,j)，当f(i,j)≥g1(i,j)且f(i,j)>g2(i,j)，则认为像素点pij属于弱窄带回波疑似区域Ω，并将该像素点pij置为前景，否则，置为背景；
1-5)将像素点pij置为背景，至此将低仰角雷达图像转化为一张二值图；
1-6)计算通过步骤1-4)和步骤1-5)所形成的二值图中各个连通区域的面积和外接矩形长度，将连通区域面积小于S2或外接矩形长度小于L2的连通区域置为背景；从而提取出弱窄带回波疑似区域；
步骤二、对弱窄带回波疑似区域进行分割、连接和筛选处理，得到弱脊状带对应的骨架图像；步骤如下：
2-1)对步骤一提取出的弱窄带回波疑似区域的轮廓进行除毛刺处理，再进行细化得到区域的骨架图像A；
2-2)在上述骨架图像A中的骨架交叉点和折点处断开骨架，得到骨架图像B，包括：
依照折点的特点通过计算骨架上某点两侧切线夹角识别出折点，即沿骨架行进，若由某点到其前侧第n个点的向量与该点到其后侧第n个点的向量的夹角小于135度，则认为该点为折点；
利用局部二值模式算子检测出端点和交叉点，即针对步骤2-1)得到的骨架图像，骨架图像中取值为1的点p为可能的端点或交叉点，以该点p为中心，考察其3×3区域及5×5区域边界的取值分布；其中，5×5区域边界点若取值为1，但在5×5范围内与区域中心不连通，则将该点置为0，分别从上述3×3区域及5×5区域的左上角开始沿逆时针方向形成8 位01链码和16位01链码来描述所述两个区域边界的取值分布；然后分别沿两个区域边界循环一周，记录依次取值变化的次数n3(p)和n5(p)，若n3(p)＝2，则点p为端点；若n3(p)≥6或n5(p)≥6，则点p为交叉点；
在上述折点和交叉点处断开骨架，从而得到骨架图像B；
2-3)通过端点匹配的方法将骨架图像B中任意两段不连通的曲线连接；
根据阵风锋走向平缓的特点，设：曲线li的端点A和曲线lj， j≠i的端点B的匹配条件如下：
式(2)中，L＝30km，Φ1＝Φ2＝Φ3＝-0.7，li和lj表示曲线i和j的长度，指由端点A沿曲线i在端点A处切线的方向指向远离该曲线的向量，指由端点B沿曲线j在端点B处切线的方向指向远离该曲线的向量；
当端点A仅与一个端点符合匹配条件时，则该端点为端点B，连接端点A和端点B；
当端点A与多个端点符合匹配条件时，则将其中值最小的端点作为端点B，连接端点A和端点B，其中，端点C是与端点A符合匹配条件的任意一端点，指由端点C沿其所在曲线在端点C处切线的方向指向远离该曲线的向量；
从而形成骨架图像C；
2-4)判断低仰角雷达图像是否存在单侧脊区域，若存在，则剔除单侧脊区域的骨架；
对骨架图像C中的每条曲线分别进行PCA处理，至少得到第二主成分方向的单位向量e2；
当距离k分别为2,3,4,5km时，设变量αk的初值为0，对每条曲线上m个点中每一点pi，i＝1,2...m,求得所述曲线一侧的位置qi，qi＝pi+k e2对于每一点pi计算：
式(3)中，f(x)表示点x的反射率值；
设，所述曲线一侧的加权比例
当距离k分别为-2,-3,-4,-5km时，设，所述曲线另一侧的加权比例为β0，根据上述过程求得加权比例β0；
若α0<Th且β0<Th，则认为该曲线与弱窄带回波疑似区域中对应的区域是弱脊状带，否则，该曲线与弱窄带回波疑似区域中对应的区域是单侧脊区域，滤除该曲线；其中，Th＝0.75；形成骨架图像D；
2-5)判断低仰角雷达图像是否存在虚线回波，若存在，则剔除虚线回波的骨架；
用2-4)中对骨架图像C中的每条曲线的PCA处理结果，能得到骨架图像D中每条曲线的第一主成分特征值λ1、第一主成分方向的单位向量e1和第二主成分特征值λ2，当第二主成分特征值λ2/第一主成分特征值λ1<0.05，且雷达中心点到曲线所拟合的直线的距离小于dl，则认为该曲线与低仰角雷达图像中对应的区域是虚线回波，滤除该曲线；其中dl＝3km；最终得到的骨架图像即为弱脊状带对应的骨架图像；
步骤三、由当前时刻和前一时刻的两幅低仰角雷达图像得到光流场，将步骤二得到的弱脊状带对应的骨架图像中前后时刻匹配的骨架拟定为疑似阵风锋，根据该疑似阵风锋的位置和速度与风暴单体的位置和速度的关系及该疑似阵风锋的走向与速度的关系识别出阵风锋；步骤如下：
3-1)利用光流法由当前时刻和前一时刻的两幅低仰角雷达图像得到光流场；
3-2)利用光流场信息将前一时刻的弱脊状带对应的骨架图像的曲线移动到由步骤二得到的当前时刻骨架图像中的对应位置，同时满足下述条件一和条件二的曲线为当前时刻和前一时刻的同一条弱窄带回波区域对应，即为疑似阵风锋；
条件一：前后时刻弱窄带回波区域的交叠长度大于30％；
条件二：利用PCA处理后得到的前后时刻弱窄带回波区域的第一主成分方向夹角小于30度；
3-3)若前一时刻弱窄带回波某端比当前时刻弱窄带回波长10km以上时，则利用多出的此段曲线对步骤3-2)中认定的疑似阵风锋进行延长补全；
3-4)根据该疑似阵风锋的位置和速度与风暴单体的位置和速度的关系及该疑似阵风锋的走向与速度的关系识别是否具有阵风锋，包括
判断当前时刻疑似阵风锋走向与疑似阵风锋速度方向锐角大于45度；
同时，当前时刻疑似阵风锋与风暴单体的速度方向与位置关系满足下述条件(1)至条件(5)中的一条；则认定当前时刻疑似阵风锋为阵风锋；并针对下一时刻低仰角雷达图像，直接将步骤3-3)补全后的该时刻的疑似阵风锋直接认定为阵风锋；
(1)当当前时刻疑似阵风锋仅位于风暴单体运动的前侧40km内时，疑似阵风锋速度方向与风暴单体移动方向夹角小于30度；
(2)当当前时刻疑似阵风锋仅位于风暴单体运动的右侧40km内时，疑似阵风锋速度方向与风暴单体移动方向右转90度后的方向夹角小于30度；
(3)当当前时刻疑似阵风锋仅位于风暴单体运动的左侧40km内时，疑似阵风锋速度方向与风暴单体移动方向左转90度后的方向夹角小于30度；
(4)当当前时刻疑似阵风锋既位于风暴单体前侧40km内且位于风暴单体右侧40km内时，疑似阵风锋速度方向位于风暴单体移动方向与风暴单体移动方向右转90度后的方向之间；
(5)当当前时刻疑似阵风锋既位于风暴单体前侧40km内且位于风暴单体左侧40km内时，疑似阵风锋速度方向位于风暴单体移动方向与风暴单体移动方向左转90度后的方向之间；
上述步骤3-2)中，利用光流场信息将前一时刻的弱脊状带对应的骨架图像的曲线移动到由步骤二得到的当前时刻骨架图像中的对应位置后，将没有同时满足条件一和条件二的前一时刻的弱脊状带对应的骨架图像的曲线保留两个体扫叠加到下一循环的疑似阵风锋匹配过程中的前一时刻的骨架图像中；
上述步骤3-3)补全后的疑似阵风锋用于替换下一循环中前一时刻对应位置的曲线。
CN201610858401.XA 2016-09-27 2016-09-27 基于多普勒天气雷达资料的阵风锋自动识别方法 Expired - Fee Related CN106526558B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610858401.XA CN106526558B (zh) 2016-09-27 2016-09-27 基于多普勒天气雷达资料的阵风锋自动识别方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610858401.XA CN106526558B (zh) 2016-09-27 2016-09-27 基于多普勒天气雷达资料的阵风锋自动识别方法

Publications (2)

Publication Number Publication Date
CN106526558A CN106526558A (zh) 2017-03-22
CN106526558B true CN106526558B (zh) 2019-03-08

Family

ID=58344287

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610858401.XA Expired - Fee Related CN106526558B (zh) 2016-09-27 2016-09-27 基于多普勒天气雷达资料的阵风锋自动识别方法

Country Status (1)

Country Link
CN (1) CN106526558B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109061649A (zh) * 2018-08-10 2018-12-21 中国气象局气象探测中心 一种冰雹监测方法及系统
CN112347956B (zh) * 2020-11-12 2022-05-06 上海交通大学 基于多无人机和机器视觉的云观测系统及方法
CN114518213A (zh) * 2020-11-19 2022-05-20 成都晟甲科技有限公司 基于骨架线约束的流场测量方法、系统、装置及存储介质
CN115542326A (zh) * 2022-09-14 2022-12-30 江苏省气象台 一种孤立式下击暴流识别方法

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105260591A (zh) * 2015-09-18 2016-01-20 江苏省气象科学研究所 一种多仰角的migfa阵风锋识别改进算法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JPH06273524A (ja) * 1993-03-24 1994-09-30 Mitsubishi Electric Corp ガストフロント検出方法
CN102645679B (zh) * 2012-03-13 2014-07-02 天津大学 一种基于多普勒天气雷达回波图像的中气旋识别方法
CN102662172A (zh) * 2012-03-29 2012-09-12 天津大学 一种基于多普勒雷达反射率图像的风暴云团的外推方法

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105260591A (zh) * 2015-09-18 2016-01-20 江苏省气象科学研究所 一种多仰角的migfa阵风锋识别改进算法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Gust front detection in weather radar images by entropy matched functional template;O. Alkhouli et al.;《IEEE International Conference on Image Processing, 2005》;20050914;全文
多普勒天气雷达对阵风锋的识别算法研究;薛震刚 等;《中国体视学与图像分析》;20140930;第19卷(第3期);全文
阵风锋自动识别与预警;郑佳锋 等;《应用气象学报》;20130228;第24卷(第1期);全文

Also Published As

Publication number Publication date
CN106526558A (zh) 2017-03-22

Similar Documents

Publication Publication Date Title
CN108510467B (zh) 基于深度可变形卷积神经网络的sar图像目标识别方法
CN106526558B (zh) 基于多普勒天气雷达资料的阵风锋自动识别方法
CN103442209B (zh) 一种输电线路的视频监控方法
CN101975940B (zh) 基于分割组合的sar图像自适应恒虚警率目标检测方法
CN110414411A (zh) 基于视觉显著性的海面船只候选区域检测方法
CN102645679B (zh) 一种基于多普勒天气雷达回波图像的中气旋识别方法
KR101258668B1 (ko) 한반도 통합형 기상 레이더 품질 관리 시스템 및 그 방법
CN100520362C (zh) 基于彩色ccd图像分析的森林火情烟雾检测方法
CN106650812B (zh) 一种卫星遥感影像的城市水体提取方法
CN107238821A (zh) 一种基于特征谱特征的机场跑道异物检测方法及装置
CN109117802A (zh) 面向大场景高分遥感影像的舰船检测方法
CN101339613B (zh) 一种遥感图像背景噪声减弱方法
CN109191432A (zh) 基于域变换滤波多尺度分解的遥感图像云检测方法
CN104408482A (zh) 一种高分辨率sar图像目标检测方法
CN109214439A (zh) 一种基于多特征融合的红外图像结冰河流检测方法
CN103927758B (zh) 一种基于对比度与角点最小凸包的显著性检测方法
CN105512622B (zh) 一种基于图分割和监督学习的可见光遥感图像海陆分割方法
CN109614936A (zh) 遥感图像飞机目标的分层识别方法
CN108229342A (zh) 一种海面舰船目标自动检测方法
CN102855485A (zh) 一种小麦抽穗的自动检测方法
CN108805050A (zh) 基于局部二值模式的电线检测方法
CN105139034B (zh) 一种结合光谱滤除的船舶检测方法
CN107067039A (zh) 基于超像素的sar图像舰船目标快速检测方法
CN109117776A (zh) 基于航迹信息的飞机与气象杂波分类识别方法
CN100463001C (zh) 一种类球形果蔬的识别方法

Legal Events

Date Code Title Description
C06 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
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20190308

Termination date: 20210927

CF01 Termination of patent right due to non-payment of annual fee