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

技术特征：

1.一种基于多普勒天气雷达资料的阵风锋自动识别方法，其特征在于，包括以下步骤：

步骤一、依据阵风锋的雷达表现特征，利用局部二值化算法提取出弱窄带回波疑似区域；步骤如下：

1-1)设一幅低仰角雷达图像的大小为N×N，对中间[N-(2n+1)×(2n+1)]×[N-(2n+1)×(2n+1)]区域内的像素点p_ij进行其是否属于弱脊状区域的判断；

1-2)求出区域[N-(2n+1)×(2n+1)]×[N-(2n+1)×(2n+1)]内大于等于35dBZ的连通区域ω_i，i＝1,2，…s，计算该连通区域ω_i的面积s_i及该连通区域ω_i外包矩形长度l_i，对于面积s_i小于S₁或外包矩形长度l_i小于L₁的连通区域进行标记，得到ω_j′，j＝1,2，…m，m≤s；

1-3)若像素点p_ij的反射率值f(i,j)∈x₁或f(i,j)∈x₂且f(i,j)∈ω_j′，其中x₁＝[5,35)dBZ，x₂＝[35,40)dBZ，执行1-4)，否则，执行1-5)

1-4)在以该像素点p_ij为中心的区域根据式(1)做卷积运算，得到卷积运算结果g₁(i,j)和g₂(i,j)，当f(i,j)≥g₁(i,j)且f(i,j)>g₂(i,j)，则认为像素点p_ij属于弱窄带回波疑似区域Ω，并将该像素点p_ij置为前景，否则，置为背景；

$\left\{\begin{array}{cc}{g}_1(i,j)=\underset{s=-n}{\overset{s=n}{\Σ}}\underset{t=-n}{\overset{t=n}{\Σ}}w(s,t)f(x+s,y+t)& w(s,t)=\frac{1}{n^2}\\ g_2(i,j)=\underset{s=-(2n+1)}{\overset{s=2n+1}{\Σ}}\underset{t=-(2n+1)}{\overset{t=2n+1}{\Σ}}w(s,t)f(x+s,y+t)& w(s,t)=\frac{1}{{(2n+1)}^2}\end{array}\right.---\left(1\right)$

1-5)将像素点p_ij置为背景，至此将低仰角雷达图像转化为一张二值图；

1-6)计算步骤1-5)所形成的二值图中各个连通区域的面积和外接矩形长度，将连通区域面积小于S₂或外接矩形长度小于L₂的连通区域置为背景；从而提取出弱窄带回波疑似区域；

步骤二、对弱窄带回波疑似区域进行分割、连接和筛选处理，得到弱脊状带对应的骨架图像；步骤如下：

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链码来描述所述两个区域边界的取值分布；然后分别沿两个区域边界循环一周，记录依次取值变化的次数n₃(p)和n₅(p)，若n₃(p)＝2，则点p为端点；若n₃(p)≥6或n₅(p)≥6，则点p为交叉点；

在上述折点和交叉点处断开骨架，从而得到骨架图像B；

2-3)通过端点匹配的方法将骨架图像B中任意两段不连通的曲线连接；

根据阵风锋走向平缓的特点，设：曲线l_i的端点A和曲线l_j(j≠i)的端点B的匹配条件如下：

式(2)中，L＝30km，Φ₁＝Φ₂＝Φ₃＝-0.7，l_i和l_j表示曲线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处理，至少得到第二主成分方向的单位向量e₂；

当距离k分别为2,3,4,5km时，设变量α_k的初值为0，对每条曲线上m个点中每一点p_i，i＝1,2...m,求得所述曲线一侧的位置q_i，q_i＝p_i+k e₂对于每一点p_i计算：

$\left\{\begin{array}{cc}{\mathrm{\α}}_k={\mathrm{\α}}_k& f\left(p_i\right)\>f\left(q_i\right)\\ {\mathrm{\α}}_k={\mathrm{\α}}_k+1& f\left(p_i\right)=f\left(q_i\right)\\ {\mathrm{\α}}_k={\mathrm{\α}}_k+2& f\left(p_i\right)\<f\left(q_i\right)\end{array}\right.,i=1,2...m---\left(3\right)$

式(3)中，f(x)表示点x的反射率值；

设，所述曲线一侧的加权比例

当距离k分别为-2,-3,-4,-5km时，设，所述曲线另一侧的加权比例为β₀，根据上述过程求得加权比例β₀；

若α₀<Th且β₀<Th，则认为该曲线与弱窄带回波疑似区域中对应的区域是弱脊状带，否则，该曲线与弱窄带回波疑似区域中对应的区域是单侧脊区域，滤除该曲线；其中，Th＝0.75；形成骨架图像D；

2-5)判断低仰角雷达图像是否存在虚线回波，若存在，则剔除虚线回波的骨架；

用2-4)中对骨架图像C中的每条曲线的PCA处理结果，能得到骨架图像D中每条曲线的第一主成分特征值λ₁、第一主成分方向的单位向量e₁和第二主成分特征值λ₂，当第二主成分特征值λ₂/第一主成分特征值λ₁<0.05，且雷达中心点到曲线所拟合的直线的距离小于d_l，则认为该曲线与低仰角雷达图像中对应的区域是虚线回波，滤除该曲线；其中d_l＝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)补全后的疑似阵风锋用于替换下一循环中前一时刻对应位置的曲线。