CN103424783A - 一种基于极轨卫星微波遥感的去除模糊风场的方法 - Google Patents

Abstract

本发明涉及一种基于极轨卫星微波遥感的去除模糊风场的方法，包括以下步骤：一、对测量得到的非气旋风场初始化；二、确定滤波器参数；三、分区圆中数滤波；四、迭代计算：重复步骤三，直到当前风场不再改变或迭代次数达到最大迭代次数为止，从而得到每个风元上的去除了模糊风场的真实风矢量。本发明能够有效去除块状模糊，取得良好效果；且定义简单明确，计算量小，收敛速度较快，具有广阔的通用性和普适性。

Description

一种基于极轨卫星微波遥感的去除模糊风场的方法

技术领域

本发明涉及一种基于极轨卫星微波遥感的去除模糊风场的方法，属于大气产品反演算法技术领域。

背景技术

极轨卫星微波遥感主要采用星载微波散射计探测海面风场，它是目前唯一能够同时探测海面风速和风向的主动微波遥感仪器，国内外主要的星载微波散射计是欧洲METOP卫星搭载的ASCAT，日本QuickSCAT卫星搭载的Seawinds和中国FY-3卫星搭载的WindRAD。工作原理是顺轨方向上的同一风元(Wind VectorCell，WVC)一般会有多个雷达后向散射截面积NRCS(normalized radar crosssection)实测值，利用这些不同观测条件下获取的NRCS值和相应的观测参数，即可反演出海面风场。

由于NRCS与风速存在着对数正相关，与相对风向χ存在双调和关系，所以在理论上只要有2个以上实测NRCS，就可反演出一个确定的风矢量。但由于大气、云水和雨衰，以及卫星运行状态的不稳定性，得到的NRCS测量值会受到污染，其误差满足均值为零的高斯分布。这些受污染的NRCS，使得在近似风速条件下，出现风向180°变化的解，通过最大似然法进行风场反演，也往往只能得到一组模糊解。模糊解一般为2～4个，其中只有一个接近真实风矢量，其他的称为伪解。所以，要得到每个风元上的真实风矢量，必须对初步反演的结果进行模糊去除。

模糊去除通常采用圆中数滤波法，但是常规的圆中数滤波法无法去除块状模糊，同时传统的圆中数滤波法还会由于第一风场的缺陷使得滤波迭代毁坏相邻的风矢量值。

因此，在模糊风场去除的应用方面，需提出一种更好的模糊风场去除方法。

发明内容

有鉴于此，本发明的目的是基于极轨卫星微波遥感多角度观测特点，建立一种能够有效去除风场模糊特性的反演方法。

本发明为实现上述目的，采用如下技术方案，包括以下步骤：

一、对测量得到的非气旋风场初始化，该非气旋风场包括M行N列个风元，具体方法为：

按照风向角度将360°分为8个区，0～45°，45°～90°，90°～135°，135°～180°，180°～225°，225°～270°，270°～315°，315°～360°，每个区的代表角度为θ₁＝22.5，θ₂＝67.5，θ₃＝112.5，θ₄＝157.5，θ₅＝202.5，θ₆＝247.5，θ₇＝292.5，θ₈＝337.5，它们分别代表8个角度范围的特征值；

对通过最大似然估计得到的具有模糊特性的第一风场风向进行统计，统计出各个角度范围中风向的个数，即某个代表角度θ_i区间里有n_i个角，将这8个区间中具有最多角度个数的区间提取出来，并且提取其左右相邻两边区间的风向数量，计算三个区间的风向数量之和，以每个区间风向数量与该和值的比值作为权值，对三个区间的特征值进行加权求和，作为风向均值；

将最大似然估计得到的每个风元的4个模糊解的角度

i∈M，j∈N，k∈4依次和该风向均值进行比较，取差值Δθ_min(i，j)最小的风向

作为真解即初始圆中数，并且记录对应的风速值，对各个风元的模糊解按照与风向均值

的差值进行重新排序，初始化整个M行N列风元后得到第二风场；

二、确定滤波器参数：定义窗口大小为H×H，H取值为5或7，定义最大迭代次数＜＝100；

三、分区圆中数滤波，具体方法为：

1)定义h＝(H-1)/2，把第二风场的矩形风场区域分为9个子区域，即分为三行，每一行又分为左区、中区、右区；每个子区域的范围是：

第一行左区：1＜i＜h+1，1＜j＜h+1；

第一行中区：1≤i≤h+1，h+1≤j≤N-h

第一行右区：1≤i≤h+1，N-h≤j≤N，

第二行左区：h+1≤i≤M-h，1≤j≤h+1，

第二行中区：h+1＜i＜M-h，h+1＜j＜N-h，

第二行右区：h+1≤i≤M-h，N-h≤j≤N

第三行左区：M-h≤i≤M，1≤j≤h+1，

第三行中区：M-h≤i≤M，h+1≤j≤N-h，

第三行右区：M-h≤i≤M，N-h≤j≤N

2)对于每个风元，根据风元所处的子区域，按照以下公式计算出该风元对应的四个模糊解的

的值，并按照该值大小进行重新排序，选取

值最小的模糊解做为真值，代替原窗口中心的风矢量真值，形成新的圆中数；

$\left\{\begin{array}{cc}{E}_{\mathrm{ij}}^k=\underset{m=i-h}{\overset{i+h}{\mathrm{\&Sigma;}}}\underset{n=j-h}{\overset{j+h}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& h+1<i<M-h,h+1<j<N-h,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=1}{\overset{i+h}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{j+h}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& 1<i<h+1,1<j<h+1,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=i-h}{\overset{i+h}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{j+h}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& h+1\&le;i\&le;M-h,1\&le;j\&le;h+1,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=i-h}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{j+h}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& M-h\&le;i\&le;M,1\&le;j\&le;h+1,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=1}{\overset{i+h}{\mathrm{\&Sigma;}}}\underset{n=j-h}{\overset{j+h}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& 1\&le;i\&le;h+1,h+1\&le;j\&le;N-h,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=i-h}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=j-h}{\overset{j+h}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& M-h\&le;i\&le;M,h+1\&le;j\&le;N-h,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=1}{\overset{i+h}{\mathrm{\&Sigma;}}}\underset{n=j-h}{\overset{N}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& 1\&le;i\&le;h+1,N-h\&le;j\&le;N,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=i-h}{\overset{i+h}{\mathrm{\&Sigma;}}}\underset{n=j-h}{\overset{N}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& h+1\&le;i\&le;M-h,N-h\&le;j\&le;N,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=i-h}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=j-h}{\overset{N}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& M-h\&le;i\&le;M,N-h\&le;j\&le;N,h=(H-1)/2\end{array}\right.$

其中，θ_mn为所计算风元周边的每个风元的圆中数；

四、迭代计算：重复步骤三，直到当前风场不再改变或迭代次数达到最大迭代次数为止，从而得到每个风元上的去除了模糊风场的真实风矢量。

有益效果：

本发明采用一种加强的圆中数滤波法对非气旋风场进行模糊去除，能够克服传统圆中数滤波法的苛刻条件(如各风元模糊风场必须随机分布，不可有块状模糊)，从该风元的几个模糊解中提取出和真实矢量最近的风矢量解，有效去除块状模糊，取得良好效果。

该方法定义简单明确，计算量小，且收敛速度较快，具有广阔的通用性和普适性。

附图说明

图1为将矩形风场区域分为9个区的示意图；

图2为未经过去模糊的第一风场；

图3为使用传统圆中数滤波直接解模糊处理过后的风场；

图4为使用本发明所述加强型圆中数滤波方案去除模糊后得到的风场；

表1采用欧空局的ASCAT数据，并基于加强型园中数滤波的结果与欧洲Eumetsat的同一区域同时段的风场产品相对/绝对偏差统计。

具体实施方式

下面结合附图，具体说明本发明的优选实施方式。

针对测量得到的第一风场，本实施例采用加强型圆中数滤波算法去除风场模糊矢量，该方法的具体步骤为：

一、非气旋风场初始化：

非气旋风场不同于气旋风场，在该范围内风向无较大突变。

按照风向角度对360°分为8个区，0～45°，45°～90°，90°～135°，135°～180°，180°～225°，225°～270°，270°～315°，315°～360°，每个区的代表角度为θ₁＝22.5，θ₂＝67.5，θ₃＝112.5，θ₄＝157.5，θ₅＝202.5，θ₆＝247.5，θ₇＝292.5，θ₈＝337.5，它们分别代表8个角度范围的特征值。对第一风场风向进行统计，统计出各个范围中风向的个数，即某个代表角度θ_i区间里有n_i个角，将这8个区间中具有最多角度个数的区间提取出来，并且提取其左右相邻两边区间的风向数量，按照以下公式进行加权求风向均值；

$\left\{\begin{array}{cc}\stackrel{\&OverBar;}{\mathrm{\&theta;}}={\mathrm{\&theta;}}_{i-1}\&times;\frac{n_{i-1}}{n_{i-1}+n_i+n_{i+1}}+{\mathrm{\&theta;}}_i\&times;\frac{n_i}{n_{i-1}+n_i+n_{i+1}}+{\mathrm{\&theta;}}_{i+1}\&times;\frac{n_{i+1}}{n_{i-1}+n_i+n_{i+1}}& i=\mathrm{2,3},..,7\\ \stackrel{\&OverBar;}{\mathrm{\&theta;}}={\mathrm{\&theta;}}_8\&times;\frac{n_8}{n_8+n_i+n_{i+1}}+{\mathrm{\&theta;}}_i\&times;\frac{n_i}{n_8+n_i+n_{i+1}}+{\mathrm{\&theta;}}_{i+1}\&times;\frac{n_{i+1}}{n_8+n_i+n_{i+1}}& i=1\\ \stackrel{\&OverBar;}{\mathrm{\&theta;}}={\mathrm{\&theta;}}_{i-1}\&times;\frac{n_{i-1}}{n_{i-1}+n_i+n_1}+{\mathrm{\&theta;}}_i\&times;\frac{n_i}{n_{i-1}+n_i+n_1}+{\mathrm{\&theta;}}_1\&times;\frac{n_1}{n_{i-1}+n_i+n_1}& i=8\end{array}\right.---\left(\text{1}\right)$

对第一风场进行初始化，将每个风元4个模糊解的角度

i∈M，j∈N，k∈4按照从最大可能解到最小可能解的顺序依次和该风向均值

进行比较，取差值Δθ_min(i，j)最小的风向

作为真解，

$\mathrm{\&Delta;}{\mathrm{\&theta;}}_{\min (i,j)}={\mathrm{\&theta;}}_{i,j}^k-\stackrel{\&OverBar;}{\mathrm{\&theta;}}---\left(2\right)$

并且记录对应的风速值，对各个风元的模糊解进行重新排序，初始化整个M行N列风元后得到第二风场；

二、确定滤波器参数：

定义窗口大小为5×5或7×7、最大迭代次数＜100；

三、分区圆中数滤波；

把矩形风场区域按照(3)式分为9个子区域(第一行：左区、中区、右区、第二行：左区、中区、右区、第三行：左区、中区、右区)，如图1。对于每个风元，根据风元所处的子区域，按照以下公式计算出该风元对应的四个模糊解的的值，并按照该值大小进行重新排序，选取

值最小的模糊解做为真值，代替原窗口中心的风矢量真值，形成新的圆中数；

$\left\{\begin{array}{cc}{E}_{\mathrm{ij}}^k=\underset{m=i-h}{\overset{i+h}{\mathrm{\&Sigma;}}}\underset{n=j-h}{\overset{j+h}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& h+1<i<M-h,h+1<j<N-h,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=1}{\overset{i+h}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{j+h}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& 1<i<h+1,1<j<h+1,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=i-h}{\overset{i+h}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{j+h}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& h+1\&le;i\&le;M-h,1\&le;j\&le;h+1,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=i-h}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{j+h}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& M-h\&le;i\&le;M,1\&le;j\&le;h+1,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=1}{\overset{i+h}{\mathrm{\&Sigma;}}}\underset{n=j-h}{\overset{j+h}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& 1\&le;i\&le;h+1,h+1\&le;j\&le;N-h,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=i-h}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=j-h}{\overset{j+h}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& M-h\&le;i\&le;M,h+1\&le;j\&le;N-h,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=1}{\overset{i+h}{\mathrm{\&Sigma;}}}\underset{n=j-h}{\overset{N}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& 1\&le;i\&le;h+1,N-h\&le;j\&le;N,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=i-h}{\overset{i+h}{\mathrm{\&Sigma;}}}\underset{n=j-h}{\overset{N}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& h+1\&le;i\&le;M-h,N-h\&le;j\&le;N,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=i-h}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=j-h}{\overset{N}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& M-h\&le;i\&le;M,N-h\&le;j\&le;N,h=(H-1)/2\end{array}\right.$

四、迭代计算：

重复上一个步骤，直到当前风场不再改变或迭代次数达到最大迭代次数为止。

在传统的圆中数算法中，圆中数θ满足分布在两半圆中的风向数目相等。在该定义下，会出现多个圆中数，需从中选择与圆平均数最接近的一个作为唯一的圆中数。因此，该定义下的圆中数计算较复杂。对于本发明所采用的方法，由于第二风场中每个风元都有几个模糊解可供选择，所以不必要在整个风矢量窗口中进行搜索确定圆中数，再在模糊解中选择一个与圆中数最接近的模糊解做为窗口中心的风矢量；可以直接从中选出一个使(2)式达到最小的模糊解代替窗口中心的风矢量。若选取(M×N)的二维风矢量窗口进行滤波运算，这里引入一种分区的圆中数定义，将一个矩形区域风场分为9个区，每个区的领域范围不同，导致了迭代取值范围的不同。

下面举一个实例来进一步说明本发明所提供的方法。

数据选取与参数设定：

本实例所使用的散射计是一种采用推帚型扫描方式，扇形波束探测的仪器，其工作频率为5.5GHz，对于雨衰不敏感，通过排列于系统两侧45°、90°和135°的各三个辐射波导天线，产生仰角宽度为40°、方位角宽度为0.4°的扇形波束。

该散射计每条轨道包括1629行、42列个风元数据，每个风元大小为50km×50km。行号顺序与卫星飞行方向一致，列号从右到左分别为1到42。本实例选取了散射计2009年9月15日轨道号为15092的部分数据(地理范围是介于N9～29，E142.8～153，跨度范围长1100*宽550km，该区域无雨)对本发明的模糊去除方法进行验证。由于散射计工作在C波段，该波段对云雨导致回波信号衰减并不灵敏，因此没有区分该区域内晴空、云水和雨区的单独影响，而是进行统一的反演考虑。滤波窗口大小选7×7较为适宜。故在本实验中，圆中数滤波器参数设定如下：窗口大小为7×7，最大迭代次数为100次。

结果分析

为了进行对比，先对所选数据全部用第一模糊解初始化风场进行圆中数滤波，然后采用加强型圆中数滤波方案进行滤波。实验结果如图2到4所示。各图中的箭头方向表示风向，箭头长度表示风速，长度越大表示风速越大，反之风速越小。

图2为经最大似然法算出的由第一模糊解构成的初始风场。从该图可以看出，中间区域出现块状风向反向和不连续现象(正确风向应当指向右边，图中大部分风向正确)，右上角区域出现少量风向不连续现象，整个区域有半数以上的风向与真实风向一致，满足传统圆中数滤波条件，即该区域中的伪解可以由圆中数滤波器去除。但是块状模糊问题由传统圆中数滤波算法无法解决，只能解决右上角少量不连续现象。

图3为图2的第一风场直接经过传统圆中数滤波过后的结果。该图与图2相比，右上角区域的少量风向模糊得到有效的去除。但是，中间区域大量连续模糊没有得到改善，出现了所谓的块状模糊，该区域内的风向与周围风向相差160°至180°之间。故用传统圆中数滤波方法在该区域失效。

图4为采用本文所提出的加强型圆中数滤波方案进行滤波后的结果。滤波迭代次数为9次，每次风矢量改变数目依次为662，221，160，151，129，96，65，26，0。和图3该区域进行比较，可以发现，图3中间区域的块状模糊得到有效去除，右上角少量不连续也得到纠正，整个风场风向分布较为连续，达到了加强型圆中数滤波去除风矢量模糊的目的。

数据偏差分析

为了对滤波结果精度进行定量分析，采用欧空局的ASCAT数据，并基于加强型园中数滤波的结果与欧洲Eumetsat的同一区域同时段的风场产品相对/绝对偏差统计。(本文反演数据与欧洲产品数据相减)，统计数据(1638个点)如下：

表1

表1中的各项指标基于相对/绝对偏差进行统计。对于相对偏差：风速偏差满足正态分布，风速误差均值-1.5858m/s，测量风速比真实风速普遍偏小，主要是由于对主动探测的回波功率没有做水汽、云水和氧气衰减修正，使得实际收到回波功率比真实值偏小，由于回波功率与风速存在着对数正相关，因此测得的风速偏小，本文重点是解决风向块状模糊问题，没有做大气、云水和雨衰订正，因此风速偏小。风向偏差满足正态分布，风向误差均值1.0424°，风向偏差较小，因为块状模糊得到良好解决。对于绝对偏差：平均偏差远小于风场反演精度要求的最大误差值(风速误差＜2m/s，风向误差＜20°)。

实验表明，模糊解之间风速差异很小，而风向相差较大。加强圆中数滤波去除模糊解的目的是从几个模糊解序列中选出一个最可能接近真实解的风矢量，并不改变模糊解的值，因此滤波精度主要体现在风向偏差上，而表中平均风向偏差确实很小(只有8.426°)。以上分析从定量的角度证明了本文采用的滤波算法是比较可靠的。

本发明的方法具有普适性，但在业务运行该方法前，应当使用其它星载仪器确定台风云图覆盖区域，在去除台风覆盖区域后，在非气旋风场区域使用本方法。因此该方法为业务化非气旋洋面风场数据提取提供了新方案。这就大大降低了风场反演过程中反演风向模糊的特性。

综上所述，以上仅为本发明的较佳实施例而已，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (1)

1.一种基于极轨卫星微波遥感的去除模糊风场的方法，其特征在于，包括以下步骤：

一、对测量得到的非气旋风场初始化，该非气旋风场包括M行N列个风元，具体方法为：

按照风向角度将360°分为8个区，0～45°，45°～90°，90°～135°，135°～180°，180°～225°，225°～270°，270°～315°，315°～360°，每个区的代表角度为θ₁＝22.5，θ₂＝67.5，θ₃＝112.5，θ₄＝157.5，θ₅＝202.5，θ₆＝247.5，θ₇＝292.5，θ₈＝337.5，它们分别代表8个角度范围的特征值；

对通过最大似然估计得到的具有模糊特性的第一风场风向进行统计，统计出各个角度范围中风向的个数，即某个代表角度θ_i区间里有n_i个角，将这8个区间中具有最多角度个数的区间提取出来，并且提取其左右相邻两边区间的风向数量，计算三个区间的风向数量之和，以每个区间风向数量与该和值的比值作为权值，对三个区间的特征值进行加权求和，作为风向均值；

将最大似然估计得到的每个风元的4个模糊解的角度

i∈M，j∈N，k∈4依次和该风向均值

进行比较，取差值Δθ_min(i，j)最小的风向

作为真解即初始圆中数，并且记录对应的风速值，对各个风元的模糊解按照与风向均值

的差值进行重新排序，初始化整个M行N列风元后得到第二风场；

二、确定滤波器参数：定义窗口大小为H×H，H取值为5或7，定义最大迭代次数＜＝100；

三、分区圆中数滤波，具体方法为：

1)定义h＝(H-1)/2，把第二风场的矩形风场区域分为9个子区域，即分为三行，每一行又分为左区、中区、右区；每个子区域的范围是：

第一行左区：1＜i＜h+1，1＜j＜h+1；

第一行中区：1≤i≤h+1，h+1≤j≤N-h

第一行右区：1≤i≤h+1，N-h≤j≤N，

第二行左区：h+1≤i≤M-h，1≤j≤h+1，

第二行中区：h+1＜i＜M-h，h+1＜j＜N-h，

第二行右区：h+1≤i≤M-h，N-h≤j≤N

第三行左区：M-h≤i≤M，1≤j≤h+1，

第三行中区：M-h≤i≤M，h+1≤j≤N-h，

第三行右区：M-h≤i≤M，N-h≤j≤N

2)对于每个风元，根据风元所处的子区域，按照以下公式计算出该风元对应的四个模糊解的的值，并按照该值大小进行重新排序，选取

值最小的模糊解做为真值，代替原窗口中心的风矢量真值，形成新的圆中数；

$\left\{\begin{array}{cc}{E}_{\mathrm{ij}}^k=\underset{m=i-h}{\overset{i+h}{\mathrm{\&Sigma;}}}\underset{n=j-h}{\overset{j+h}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& h+1<i<M-h,h+1<j<N-h,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=1}{\overset{i+h}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{j+h}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& 1<i<h+1,1<j<h+1,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=i-h}{\overset{i+h}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{j+h}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& h+1\&le;i\&le;M-h,1\&le;j\&le;h+1,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=i-h}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=1}{\overset{j+h}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& M-h\&le;i\&le;M,1\&le;j\&le;h+1,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=1}{\overset{i+h}{\mathrm{\&Sigma;}}}\underset{n=j-h}{\overset{j+h}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& 1\&le;i\&le;h+1,h+1\&le;j\&le;N-h,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=i-h}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=j-h}{\overset{j+h}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& M-h\&le;i\&le;M,h+1\&le;j\&le;N-h,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=1}{\overset{i+h}{\mathrm{\&Sigma;}}}\underset{n=j-h}{\overset{N}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& 1\&le;i\&le;h+1,N-h\&le;j\&le;N,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=i-h}{\overset{i+h}{\mathrm{\&Sigma;}}}\underset{n=j-h}{\overset{N}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& h+1\&le;i\&le;M-h,N-h\&le;j\&le;N,h=(H-1)/2\\ E_{\mathrm{ij}}^k=\underset{m=i-h}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=j-h}{\overset{N}{\mathrm{\&Sigma;}}}\left|({\mathrm{\&theta;}}_{\mathrm{ij}}^k-{\mathrm{\&theta;}}_{\mathrm{mn}})\right|,& M-h\&le;i\&le;M,N-h\&le;j\&le;N,h=(H-1)/2\end{array}\right.$

其中，θ_mn为所计算风元周边的每个风元的圆中数；

四、迭代计算：重复步骤三，直到当前风场不再改变或迭代次数达到最大迭代次数为止，从而得到每个风元上的去除了模糊风场的真实风矢量。

Images