CN111175709B - 一种基于误差抑制的面向大范围气象雷达的拼图方法 - Google Patents

Abstract

本申请公开了一种基于误差抑制的面向大范围气象雷达的拼图方法，包括：步骤1，根据雷达站点之间的地表距离和距离阈值，确定雷达对和雷达环，获取各雷达站的雷达回波信号，并生成雷达站无向图。步骤2，在雷达站无向图中提取相应的雷达回波信号，依次计算符合雷达对线性约束条件、环路约束条件和对称性约束条件的第一矩阵和第二矩阵，以及各雷达站雷达回波信号的配准系数；步骤3，根据配准系数和获取的雷达回波信号，更新各雷达站的雷达回波信号，并根据各雷达站的经纬度，进行雷达回波信号拼接，记作雷达拼接图。通过本申请中的技术方案，利用雷达回波的时空关联特性，消除、抑制不同区域与时刻的误差因素，提高围雷达拼图的准确性与精度。

Description

一种基于误差抑制的面向大范围气象雷达的拼图方法

技术领域

本申请涉及电流检测装置的技术领域，具体而言，涉及一种基于误差抑制的面向大范围气象雷达的拼图方法。

背景技术

随着全球气候环境变化日益剧烈，极端天气频发，大范围气象监测预报对提高灾害天气预报精度、加强抗灾救灾能力具有重要意义。气象雷达覆盖范围广、测量精度高，可直接获取云、雨等天气现象的完整信息，是目前实现全天候大范围气象观测的重要手段。大范围的气象雷达拼图决定了其监测的准确性，但雷达信号中会引入误差，如杂波、噪声等，降低了雷达回波图像准确性，这种降低气象雷达回波质量的主要因素包括：地杂波、异常数据、多站配准误差等。

而现有技术中，对杂波、异常数据的处理，多采用频域的数字滤波器或者固定阈值法，这种处理方法，不能有效利用杂波、云的时空分布特性，不利于提高大范围气象雷达拼图效果。

并且，对于多雷达站拼图的配准误差，通常采用插值法过渡，如双线性插值等，未有效利用雷达站间的关联信息。

发明内容

本申请的目的在于：通过雷达回波的时空关联特性，消除与抑制不同区域与时刻的误差因素，提高大范围雷达拼图的准确性与精度。

本申请的技术方案是：提供了一种基于误差抑制的面向大范围气象雷达的拼图方法，该拼图方法包括：步骤1，获取指定区域内各雷达站点的经纬度，根据雷达站点之间的地表距离和距离阈值，确定雷达对和雷达环，获取各雷达站的雷达回波信号，并根据雷达站组成图，生成雷达站无向图，其中，雷达站组成图由各雷达站的雷达回波信号确定；步骤2，根据雷达对和雷达环，在雷达站无向图中提取相应的雷达回波信号，并依次计算符合雷达对线性约束条件、环路约束条件和对称性约束条件的第一矩阵和第二矩阵，并根据第一矩阵和第二矩阵计算各雷达站雷达回波信号的配准系数；步骤3，根据配准系数和获取的雷达回波信号，进行乘法计算，利用计算结果，更新各雷达站的雷达回波信号，并根据更新后的雷达回波信号和各雷达站的经纬度，进行雷达回波信号拼接，将拼接结果记作指定区域的雷达拼接图。

上述任一项技术方案中，进一步地，步骤2中，计算符合雷达对线性约束条件的第一矩阵和第二矩阵，具体包括：进行零值初始化，得到初始化第一矩阵A₀和第二矩阵C₀，采用遍历算法，计算任一个雷达对的符合雷达对线性约束条件的第一增量矩阵ΔA₁和ΔC₁，根据第一增量矩阵ΔA₁和ΔC₁，引入强度绝对值偏差，更新得到符合雷达对线性约束条件的第一矩阵A₁和第二矩阵C₁，其中，第一增量矩阵ΔA₁和ΔC₁的计算公式为：

ΔA₁＝{a_xy|a_xy＝0,xy≠m；a_mm＝1}_M×M

ΔC₁＝{c_x|c_x＝0,i≠m；c_mm＝Σlog[g_j(r_ij)/g_i(r_ij)],f(i,j)＝m}_M×1

式中，a_xy为编号为第x和y的雷达对的编号，m为雷达对编号，g_i为第i个雷达站的雷达回波信号，r_ij为雷达回波信号g_i和g_j之间的重叠区域；

用于计算配准系数的第一矩阵A和第二矩阵C的计算公式为：

A₃＝A₂+ΔA₃

式中，符合环路约束条件的第一矩阵A₂由符合环路约束条件的第二增量矩阵ΔA₂确定。

上述任一项技术方案中，进一步地，步骤2中，计算符合环路约束条件的第一矩阵A₂，具体包括：依据第k个雷达环中包括的雷达对，计算符合环路约束条件的第二增量矩阵ΔA₂，并根据符合雷达对线性约束条件的第一矩阵A₁，更新得到符合环路约束条件的第一矩阵A₂，其中，第二增量矩阵ΔA₂的计算公式为：

A₂＝A₁+ΔA₂

p_i＝f(k_i,k_i+1),1≤i≤m-1

p_m＝f(k_m,k₁)

式中，p_i为第k个雷达环中第i、i+1个雷达站之间的索引映射关系。

上述任一项技术方案中，进一步地，步骤1中，获取各雷达站的雷达回波信号后，拼图方法还包括：步骤11，提取雷达回波信号中的多个仰角数据，并根据仰角数据与雷达接收机之间的距离，统计各个仰角数据的雷达回波强度分布；

步骤12，采用差分算法，计算雷达回波强度分布的一阶差分值和二阶差分值，并根据一阶差分值和二阶差分值，采用搜索算法，确定最小距离r₀，其中，最小距离r₀的判断公式为：

式中，dp(r)为一阶差分值，d²p(r)为二阶差分值，ε和η分别为经验阈值；

步骤13，根据最小距离r₀，依次对提取出的多个仰角数据进行滤波处理，当判定仰角数据与雷达接收机之间的距离r小于最小距离r₀时，进行滤波处理，将相应仰角数据的雷达回波强度分布中小于强度阈值p₀的强度值置为0，其中，强度阈值p₀为特征回波强度。

上述任一项技术方案中，进一步地，步骤13之后，拼图方法还包括：步骤14，根据滤波处理后的仰角数据的雷达回波强度，构建2L×N_A的回波矩阵Re(m,n)，其中，L为滤除异常方位角采用的邻方位角个数，其取值根据经验确定，通常确定为5-11之间的奇数。N_A为第s个仰角角度e_s上方位角的个数，该值可在气象雷达参数中查到，通常为360或520，回波矩阵Re(m,n)的计算公式为：

式中，m、n为行列数，a_m+n-L为第m+n-L个方位角，e_s为第s个仰角，r_k为第k个仰角数据与雷达站之间的距离，t₀为仰角数据对应的时刻；

步骤15，统计各方向角对应的雷达回波强度，生成长度为N_A的分布特征向量Q，其中，分布特征向量Q的元素个数为N_A；

步骤16，根据分布特征向量Q的取值，确定回波矩阵Re(m,n)中相应列数的中位数或均值，根据中位数和均值，计算雷达回波信号的方位角阈值；

步骤17，根据方位角阈值，依次判断各个方位角对应的雷达回波信号是否大于方位角阈值，若是，将该雷达回波信号记作波位异常值，计算与该方位角相邻的方位角对应的雷达回波信号的均值，记作替换值，利用替换值代替波位异常值，进行波位异常滤除。

本申请的有益效果是：

通过本申请中的拼图方法，利用雷达回波的时空关联特性，消除与抑制不同区域与时刻的误差因素，在不损失云信号的前提下，有效抑制地杂波，去除异常干扰，实现多雷达站联合配准，有效降低雷达配准误差，提高大范围雷达拼图的信噪比。

本申请中的方法，针对杂波、异常数据与多雷达站匹配误差等影响雷达拼图的主要误差因素进行建模分析，通过雷达回波的时空关联特性，通过自适应方法消除与抑制不同区域与时刻的误差因素，提高大范围雷达拼图的准确性与精度。

对于地杂波与海杂波，通过获取特定仰角上杂波与回波信号的空间分布特性，自适应估计杂波空间与强度分布范围，确定有效滤除阈值，去除近场地杂波的干扰。对回波中的异常数据，则基于相邻方位角上回波信号的空间相关性，采用动态恒虚警阈值法，在滤除干扰同时尽量保留有效回波信息。

对于大范围多雷达拼图中的配准误差，采用递归图搜索方法，遍历多雷达站组成的环约束，提取站点之间的回波映射关系，提高雷达站间的配准精度，优化大范围雷达拼图准确度。

附图说明

本申请的上述和/或附加方面的优点在结合下面附图对实施例的描述中将变得明显和容易理解，其中：

图1是根据本申请的一个实施例的基于误差抑制的面向大范围气象雷达的拼图方法的示意流程图；

图2是根据本申请的一个实施例的各个仰角数据的雷达回波强度分布的示意图；

图3是根据本申请的一个实施例的滤波处理后的雷达回波的示意图；

图4是根据本申请的一个实施例的异常帧处理后的雷达回波的示意图；

图5是根据本申请的一个实施例的大范围气象雷达拼图方法的示意图；

图6是根据本申请的一个实施例的拼图结果示意图。

具体实施方式

为了能够更清楚地理解本申请的上述目的、特征和优点，下面结合附图和具体实施方式对本申请进行进一步的详细描述。需要说明的是，在不冲突的情况下，本申请的实施例及实施例中的特征可以相互结合。

在下面的描述中，阐述了很多具体细节以便于充分理解本申请，但是，本申请还可以采用其他不同于在此描述的其他方式来实施，因此，本申请的保护范围并不受下面公开的具体实施例的限制。

如图1所示，本实施例提供了一种基于误差抑制的面向大范围气象雷达的拼图方法，该方法包括：

步骤1，获取指定区域内各雷达站点的经纬度，根据雷达站点之间的地表距离和距离阈值，确定雷达对和雷达环，获取各雷达站的雷达回波信号，并根据雷达站组成图V，生成雷达站无向图(G，V)。

具体的，本实施例基于麦卡托投影，将某个经纬度范围内的方形区域D，记作指定区域，设定该方形区域D内雷达站的总数为N，N＞1，且N个雷达站中存在若干个覆盖区域有交集的雷达站，将第i个雷达站的雷达回波信号g_i组成的雷达站回波信号集合记作G＝{g_i，1≤i≤N}。

在该方形区域D中，根据各雷达站的经纬度，可以计算出各雷达站间的地表距离，当判定两个雷达站见得地表距离d_ij小于距离阈值d₀时，将这两个雷达站记作一个雷达对。

雷达环的确定方法为，以第i个雷达站为起点，通过深度优先搜索遍历该方形区域D内的所有雷达站，在任一次遍历过程中，若重新回到第i个雷达站，则此次遍历时经过的雷达站记作一个雷达环。

对于该方形区域D内的雷达站，可以生成雷达站组成图V＝{v_ij，1≤i，j≤N}，其中，元素v_ij的取值由第i、j个雷达站的雷达回波信号确定。

若第i个雷达站的雷达回波信号g_i与第j个雷达站的雷达回波信号g_j之间存在重叠区域，则元素v_ij＝1，并将该重叠区域记作r_ij，否则，元素v_ij＝0，其中，雷达站组成图V中元素v_ij为非零元素的个数为M，0≤M≤N(N-1)/2，即该方形区域D内的雷达站间存在M个重叠区域，重叠区域集合R＝{r_ij，1≤i,j≤N}。也就是说，第i、j个雷达站组成的雷达对之间，存在重叠区域r_ij，即M为雷达对的总对数。

本实施例中，在确定重叠区域后，对重叠区域进行筛选，通过对该重叠区域进行二维相关分析，计算相关系数，当判定相关系数低于0.8时，去除对应的重叠区域，当判定相关系数大于0.8时，进行平面位置校准。

因此，结合雷达站组成图V和雷达站回波信号集合G，可以生成雷达站无向图(G，V)，并且，结合雷达站的标号，可以建立第i、j个雷达站与雷达站组成图V中M个非零元素之间的索引映射关系f(i，j)。

进一步的，步骤1中，获取各雷达站的雷达回波信号后，该拼图方法还包括对任一个雷达站的雷达回波信号进行滤波处理，滤波处理的方法具体包括：

步骤11，提取雷达回波信号中的多个仰角数据，并根据仰角数据与雷达接收机之间的距离，统计各个仰角数据的雷达回波强度分布，其中，仰角数据x(a,e_s,r,t0)包括设定的各个仰角上360°方位角与一定距离的雷达回波，其中，a为方位角，e_s为第s个仰角的角度，r为仰角数据与雷达接收机之间的距离，t0为仰角数据对应的时刻；

也就是说，统计雷达回波强度分布时，首先从第一个仰角开始，逐个统计360°方位角上的雷达回波，再统计第二个仰角中所有方位角对应的雷达回波，直至统计至最后一个仰角。

具体的，设定雷达站采用扫描VCP模式，仰角为0.5至19.5°，方位角为0至360°，回波距离为300至460km。雷达站附近的地杂波与海杂波(通常在20km内)可视作地面与海平面上大量连续物体的回波反射的融合。因此，雷达站在不同波位扫描的回波均包含了所有主瓣范围内的地物回波与少量旁瓣干扰，导致不同波位的杂波分布并不均匀。

在本实施例中，如图2所示，通过大量的数据研究和分析，将地物反射设定为漫反射，利用Gamma分布描述单个波位的杂波信息对应的回波强度x_clutter(r)，用高斯分布描述大气中包含水汽的云团的回波强度x_cloud(r)，因此，获取雷达回波信号中的多个仰角数据后，结合仰角数据与雷达接收机之间的距离，雷达回波强度x(r)的计算公式为：

x(r)＝x_clutter(r)+x_cloud(r)

x_cloud(r)～N(μ,σ)

式中，r为仰角数据与雷达接收机之间的距离，α、β为Gamma分布的分布参数，Γ(·)为伽马函数，e为自然常数，μ,σ为高斯分布中的参数，进而可以根据雷达回波强度x(r)，统计出各仰角数据的雷达回波强度分布用p(r)，具体过程此处不再赘述。

根据气象雷达实测结果，Gamma分布中的分布参数α的取值范围为0.1-2.5，即最大分布既有可能在零点，也有可能出现在均值处。

因此，实际雷达回波分布将出现以下三种情况：

(a)杂波与云信号显著分离，云信号分布离雷达较远，杂波分布在雷达附近。

(b)云信号与杂波信号重叠，通常云信号较杂波信号强，则杂波信号会被云信号覆盖。

(c)若未探测到云，则回波分布均为杂波。

可见对于特定仰角，杂波与有效云信号之间存在波谷。只要识别该波谷，即可显著去除杂波同时不影响云回波强度。

步骤12，采用差分算法，计算雷达回波强度分布的一阶差分值和二阶差分值，并根据一阶差分值和二阶差分值，采用搜索算法，确定最小距离r₀，其中，最小距离r₀的判断公式为：

式中，dp(r)为一阶差分值，d²p(r)为二阶差分值，ε和η分别为经验阈值，两者在本实施例中的取值均为0.1。

步骤13，根据最小距离r₀，依次对提取出的多个仰角数据进行滤波处理，当判定仰角数据与雷达接收机之间的距离r小于最小距离r₀时，进行滤波处理，将相应仰角数据的雷达回波强度分布中小于强度阈值p₀的强度值置为0，其中，强度阈值p₀为特征回波强度，对气象雷达而言，强度阈值p₀通常取20dBz。

具体的，如图3(a)所示，获取到的雷达回波信号，通常会受到雷达站附近的地杂波影响，因此，通过获取特定仰角上杂波与回波信号的空间分布特性，自适应估计杂波空间与强度分布范围，确定有效滤除阈值，进行滤波处理，去除近场地杂波的干扰，得到处理后的雷达回波信号如图3(b)所示。

进一步的，本实施例中对任一个雷达站的雷达回波信号的滤波处理方法还包括：

步骤14，根据滤波处理后的仰角数据的雷达回波强度，构建2L×N_A的回波矩阵Re(m,n)，其中，L为滤除异常方位角采用的邻方位角个数，其取值根据经验确定，通常确定为5-11之间的奇数。N_A为第s个仰角角度e_s上方位角的个数，该值可在气象雷达参数中查到，通常为360或520，回波矩阵Re(m,n)的计算公式为：

式中，m、n为行列数，因此，回波矩阵中第n列的2L个元素分别为[m-L,m-1]与[m+1,m+L]方位角上雷达回波强度之和。a_m+n-L为第m+n-L个方位角，e_s为第s个仰角，r_k为第k个仰角数据与雷达站之间的距离，t₀为仰角数据对应的时刻。

步骤15，统计各方向角对应的雷达回波强度，生成长度为N_A的分布特征向量Q，其中，分布特征向量Q的元素个数为N_A。

分布特征向量Q中各元素的取值，由以下方法确定：

当判定回波矩阵中第n列元素里、存在连续10个强度大于5dBz的雷达回波强度时，将Q(n)置为1，否则，将Q(n)置为0。

步骤16，根据分布特征向量Q的取值，确定回波矩阵Re(m,n)中相应列数的中位数或均值，根据中位数和均值，计算雷达回波信号的方位角阈值，其中，方位角阈值的计算方法为：

当判定分布特征向量Q中第n个元素对应的取值Q(n)＝0时，计算回波矩阵Re(m,n)中第n列元素的均值，计算方位角阈值，方位角阈值的计算公式为：

当判定分布特征向量Q中第n个元素对应的取值Q(n)＝0时，提取回波矩阵Re(m,n)中第n列元素的中位数，计算方位角阈值，方位角阈值的计算公式为：

式中，M为阈值因子，取值为5～10，x_*为中位数，通常为连续(2L+1)去除中心点后回波强度升序排列的第m个，即第L个。

比如第10列，即n＝10，若Q(10)＝1，统计回波矩阵Re(m,n)中第10列2L个元素的中位数，若Q(10)＝0，统计回波矩阵Re(m,n)中第10列2L个元素的均值。

步骤17，根据雷达回波信号的方位角阈值，依次判断各个方位角对应的雷达回波信号是否大于方位角阈值，若是，将该雷达回波信号记作波位异常值，计算与该方位角相邻的方位角对应的雷达回波信号的均值，记作替换值，利用替换值代替波位异常值，进行波位异常滤除。

具体的，如图4(a)所示，雷达回波基数据中由于外部干扰、测量误差等因素，存在大量的显著的异常数据，即在特定波位所有距离均有回波，且与相邻波位有很大差异。因此，通过上述方法，进行波位异常滤除，得到处理后的雷达回波信号如图4(b)所示。

如图5所示，在对各雷达站进行基数据解码后，利用步骤11至步骤17进行杂波滤除和异常数据处理后，将进行球坐标变换、空间插值、麦卡托投影、地球圆角补偿等等操作，输出单站完整的雷达回波图。再通过下述过程，进行雷达拼接图。实现了在不损失云信号的前提下，有效抑制地杂波，去除异常干扰，有助于实现多雷达站联合配准，有效降低雷达配准误差，进而提高大范围雷达拼图的信噪比。

步骤2，根据雷达对和雷达环，在雷达站无向图(G，V)中提取相应的雷达回波信号，并依次计算符合雷达对线性约束条件、环路约束条件和对称性约束条件的第一矩阵A和第二矩阵C，并根据第一矩阵A和第二矩阵C计算各雷达站雷达回波信号g_i的配准系数，其中，符合对称性约束条件的第三增量矩阵依次为：

ΔC₃＝0

进一步的，计算符合雷达对线性约束条件的第一矩阵A和第二矩阵C，具体包括：

将第一矩阵A和第二矩阵C进行零值初始化，得到A₀和C₀，采用遍历算法，计算任一个雷达对的符合雷达对线性约束条件的第一增量矩阵ΔA₁和ΔC₁，根据第一增量矩阵ΔA₁和ΔC₁，引入强度绝对值偏差，更新得到第一矩阵A₁和第二矩阵C₁，其中，第一增量矩阵ΔA₁和ΔC₁的计算公式为：

ΔA₁＝{a_xy|a_xy＝0,xy≠m；a_mm＝1}_M×M

ΔC₁＝{c_x|c_x＝0,x≠m；c_mm＝Σlog[g_j(r_ij)/g_i(r_ij)],f(i,j)＝m}_M×1

式中，g_i为第i个雷达站的雷达回波信号，r_ij为雷达回波信号g_i和g_j之间的重叠区域，m表示本次迭代计算的雷达对的编号，M表示所有可能的雷达对的总对数，即1≤m≤M。第m号雷达对对应的两个雷达站编号记为i，j。雷达站号与雷达对号之间的关系记为f，即f(i,j)＝m。a_xy表征本次迭代计算第x号雷达对与y号雷达对之间的关联系数，当且仅当x＝y＝m时a_xy取1，其他均为0；

第一矩阵A和第二矩阵C的计算公式为：

A₁U＝C₁＝log[g_j(r_ij)/g_i(r_ij)]

U＝{u|u_p＝logw_ij,p＝f(i,j)}∈R^M

A₂U＝0

A₃U＝0

A₁＝A₀+ΔA₁

C₁＝C₀+ΔC₁

A₂＝A₁+ΔA₂

A₃＝A₂+ΔA₃

式中，ΔA₂为符合环路约束条件的第二增量矩阵。

进一步的，计算符合环路约束条件的第一矩阵A时，具体包括：

依据第k个雷达环中包括的雷达对，计算符合环路约束条件的第二增量矩阵ΔA₂，并根据符合雷达对线性约束条件的第一矩阵A₁，更新得到第一矩阵A₂，其中，第二增量矩阵ΔA₂的计算公式为：

p_i＝(k_i,k_i+1),1≤i≤m-1

p_m＝f(k_m,k₁)

式中，p_i为第k个雷达环中第i、i+1个雷达站之间的索引映射关系。

此处需要说明的是，由于单次计算的ΔA₁、ΔA₂、ΔA₃与ΔC₁均为稀疏矩阵，可进一步降低本实施例中拼图算法的内存占用率。

具体的，方形区域D内的N个雷达站生成雷达站无向图(G，V)，包括M个雷达对。对于第m个雷达对而言，通过重叠覆盖区域内的测量结果进行配准，校正各雷达测量回波强度的系统误差。

气象雷达回波强度与云中水粒子尺度的对数成正比，因此可将不同雷达站对同一对象的测量结果差异用线性映射近似。即雷达组网可表述为，对于第i、j个雷达站，若元素v_ij＝1，则存在：

g_i(r_ij)w_ij＝g_j(r_ij)

式中，g_i(r_ij)为重叠区域r_ij的第i个雷达站的雷达回波信号，w_ij为第i个雷达站与第j个雷达站进行雷达回波信号拼接时的配准系数，引入强度绝对值偏差b_ij，则可以将上式改写为：

W＝{w_ij，1≤i，j≤N}

式中，W为配准系数矩阵。

并可以生成雷达站有向图(G，W)，与雷达站无向图(G，V)拓扑同构，且易知配准系数矩阵W包含2M个非零元素，且w_ij×w_ji＝1。因此，求解配准系数矩阵W即可实现雷达站间的雷达回波信号拼接。

由于该方形区域D中还存在多个雷达环，因此，将雷达站无向图(G，V)中的雷达环记作L＝(l_k)，l_k为第k个环路的雷达站编号组成的序列。设定W(l_k)为按照l_k顺序的线性映射权重集合，则应满足：

ΠW(l_k)＝1

则雷达回波信号拼接可转化为求解下述最优解的过程：

s.t.ΠW(l_k)＝1

w_ij×w_ji＝1

由于雷达站回波信号均为正值，即w_ij>0，则上式可转换为：

s.t.ΠlogW(l_k)＝0

logw_ij×logw_ji＝0

本实施例中，进行如下标记，U＝{u|u_p＝logw_ij,p＝f(i,j)}∈R^M，雷达对的线性约束条件为：

A₁U＝C₁＝log[g_j(r_ij)/g_i(r_ij)]

式中，A₁∈R^M×M，C₁∈R^M。

设定第i、j个雷达站之间存在重叠区域r_ij，则元素v_ij＝1，w_ij>0，u_p＝logw_ij，矩阵A₁中存在行q，使得A₁(q,p)＝1，p＝f(i,j)，该行其他元素为0，相应的，C₁(q)＝log[g_j(r_ij)/g_i(r_ij)]。

对于雷达环的环路约束条件，则转换为：

A₂U＝0

此时，设定雷达站无向图(G，V)中雷达环的数量为L，则A₂＝R^2L×M。设定第k个雷达环包含的雷达站的编号为k₁，k₂，…，k_m，则

A₂(k,k₁)＝A₂(k,k₂)＝…＝A₂(k,k_m)＝1

且第k行其余元素为零。采用深度优先搜索，对每个雷达站逐次遍历，即可得到经过该雷达站的所有环路。

综上所述，雷达站雷达回波信号拼接过程中的限定条件为：

AU＝C

经过最小二乘法计算可得：

由此，可以计算出雷达站雷达站的配准系数，以便用于后续的雷达拼接图。

步骤3，根据配准系数和获取的雷达回波信号，进行乘法计算，利用计算结果，更新各雷达站的雷达回波信号，并根据更新后的雷达回波信号和各雷达站的经纬度，进行雷达回波信号拼接，将拼接结果记作指定区域的雷达拼接图。

具体的，如图6所示，通过本实施例中的拼图方法，对于大范围多雷达拼图中的配准误差，采用递归图搜索方法，遍历多雷达站组成的环约束，提取站点之间的回波映射关系，提高雷达站间的配准精度，优化大范围雷达拼图准确度。本实施例中的方法，能有效降低气象雷达回波拼图中的干扰，提高基于雷达回波的气象目标识别准确度。

以上结合附图详细说明了本申请的技术方案，本申请提出了一种基于误差抑制的面向大范围气象雷达的拼图方法，包括：步骤1，根据雷达站点之间的地表距离和距离阈值，确定雷达对和雷达环，获取各雷达站的雷达回波信号，并生成雷达站无向图。步骤2，在雷达站无向图中提取相应的雷达回波信号，依次计算符合雷达对线性约束条件、环路约束条件和对称性约束条件的第一矩阵和第二矩阵，以及各雷达站雷达回波信号的配准系数；步骤3，根据配准系数和获取的雷达回波信号，更新各雷达站的雷达回波信号，并根据各雷达站的经纬度，进行雷达回波信号拼接，记作雷达拼接图。通过本申请中的技术方案，利用雷达回波的时空关联特性，消除、抑制不同区域与时刻的误差因素，提高围雷达拼图的准确性与精度。

本申请中的步骤可根据实际需求进行顺序调整、合并和删减。

本申请装置中的单元可根据实际需求进行合并、划分和删减。

尽管参考附图详地公开了本申请，但应理解的是，这些描述仅仅是示例性的，并非用来限制本申请的应用。本申请的保护范围由附加权利要求限定，并可包括在不脱离本申请保护范围和精神的情况下针对发明所作的各种变型、改型及等效方案。

Claims (5)

1.一种基于误差抑制的面向大范围气象雷达的拼图方法，其特征在于，该拼图方法包括：

步骤1，获取指定区域内各雷达站点的经纬度，根据所述雷达站点之间的地表距离和距离阈值，确定雷达对和雷达环，获取各所述雷达站的雷达回波信号，并根据雷达站组成图，生成雷达站无向图，其中，所述雷达站组成图由各雷达站的雷达回波信号确定；

步骤2，根据所述雷达对和所述雷达环，在所述雷达站无向图中提取相应的雷达回波信号，并依次计算符合雷达对线性约束条件、环路约束条件和对称性约束条件的第一矩阵和第二矩阵，并根据所述第一矩阵和所述第二矩阵计算各雷达站雷达回波信号的配准系数；

步骤3，根据所述配准系数和获取的所述雷达回波信号，进行乘法计算，利用计算结果，更新各雷达站的雷达回波信号，并根据所述更新后的雷达回波信号和各雷达站的所述经纬度，进行雷达回波信号拼接，将拼接结果记作指定区域的雷达拼接图。

2.如权利要求1所述的基于误差抑制的面向大范围气象雷达的拼图方法，其特征在于，步骤2中，所述计算符合雷达对线性约束条件的第一矩阵和第二矩阵，具体包括：

进行零值初始化，得到初始化第一矩阵A₀和第二矩阵C₀，采用遍历算法，计算任一个雷达对的符合所述雷达对线性约束条件的第一增量矩阵ΔA₁和ΔC₁，根据所述第一增量矩阵ΔA₁和ΔC₁，引入强度绝对值偏差，更新得到符合雷达对线性约束条件的第一矩阵A₁和第二矩阵C₁，其中，所述第一增量矩阵ΔA₁和ΔC₁的计算公式为：

ΔA₁＝{a_xy|a_xy＝0,xy≠m；a_mm＝1}_M×M

ΔC₁＝{c_x|c_x＝0,i≠m；c_mm＝Σlog[g_j(r_ij)/g_i(r_ij)],f(i,j)＝m}_M×1

式中，a_xy为编号为第x和y的雷达对的编号，m为雷达对编号，g_i为第i个雷达站的雷达回波信号，r_ij为雷达回波信号g_i和g_j之间的重叠区域；

用于计算配准系数的第一矩阵A和第二矩阵C的计算公式为：

A₃＝A₂+ΔA₃

式中，符合环路约束条件的第一矩阵A₂由符合环路约束条件的第二增量矩阵ΔA₂确定。

3.如权利要求2所述的基于误差抑制的面向大范围气象雷达的拼图方法，其特征在于，步骤2中，计算符合环路约束条件的第一矩阵A₂，具体包括：

依据第k个雷达环中包括的雷达对，计算符合所述环路约束条件的第二增量矩阵ΔA₂，并根据所述符合雷达对线性约束条件的第一矩阵A₁，更新得到所述符合环路约束条件的第一矩阵A₂，其中，所述第二增量矩阵ΔA₂的计算公式为：

A₂＝A₁+ΔA₂

p_i＝f(k_i,k_i+1),1≤i≤m-1

p_m＝f(k_m,k₁)

式中，p_i为第k个雷达环中第i、i+1个雷达站之间的索引映射关系。

4.如权利要求1至3中任一项所述的基于误差抑制的面向大范围气象雷达的拼图方法，其特征在于，步骤1中，获取各雷达站的所述雷达回波信号后，所述拼图方法还包括：

步骤11，提取所述雷达回波信号中的多个仰角数据，并根据仰角数据与雷达接收机之间的距离，统计各个仰角数据的雷达回波强度分布；

步骤12，采用差分算法，计算雷所述达回波强度分布的一阶差分值和二阶差分值，并根据所述一阶差分值和所述二阶差分值，采用搜索算法，确定最小距离r₀，其中，所述最小距离r₀的判断公式为：

式中，dp(r)为一阶差分值，d²p(r)为二阶差分值，ε和η分别为经验阈值；

步骤13，根据所述最小距离r₀，依次对提取出的多个所述仰角数据进行滤波处理，当判定仰角数据与雷达接收机之间的距离r小于所述最小距离r₀时，进行滤波处理，将相应仰角数据的雷达回波强度分布中小于强度阈值p₀的强度值置为0，其中，强度阈值p₀为特征回波强度。

5.如权利要求4所述的基于误差抑制的面向大范围气象雷达的拼图方法，其特征在于，步骤13之后，所述拼图方法还包括：

步骤14，根据滤波处理后的仰角数据的雷达回波强度，构建2L×N_A的回波矩阵Re(m,n)，其中，L为滤除异常方位角采用的邻方位角个数，其取值根据经验确定，通常确定为5-11之间的奇数； N_A为第s个仰角角度e_s上方位角的个数，该值可在气象雷达参数中查到，通常为360或520，回波矩阵Re(m,n)的计算公式为：

式中，m、n为行列数，a_m+n-L为第m+n-L个方位角，e_s为第s个仰角，r_k为第k个仰角数据与雷达站之间的距离，t₀为仰角数据对应的时刻；

步骤15，统计各方向角对应的雷达回波强度，生成长度为N_A的分布特征向量Q，其中，所述分布特征向量Q的元素个数为N_A；

步骤16，根据所述分布特征向量Q的取值，确定所述回波矩阵Re(m,n)中相应列数的中位数或均值，根据所述中位数和所述均值，计算雷达回波信号的方位角阈值；

步骤17，根据所述方位角阈值，依次判断各个方位角对应的所述雷达回波信号是否大于所述方位角阈值，若是，将该雷达回波信号记作波位异常值，计算与该方位角相邻的方位角对应的雷达回波信号的均值，记作替换值，利用所述替换值代替 所述波位异常值，进行波位异常滤除。

Images