CN103839104A - 一种海浪有效波高反演模型建模方法 - Google Patents

Abstract

本发明涉及一种海浪有效波高反演模型建模方法，其特征在于：步骤1：设置粒子群粒子位置初值和速度初值；步骤2：初始化粒子群各参数；步骤3：设置粒子群速度位置更新公式；步骤4：确定判断宽度；步骤5：基于判断宽度judgewidth和粒子的位置值X_i两个参数，建立适应值函数，适应值函数对样本数据进行分段，得到分段模型，并计算出与分段模型对应的整体样本数据的残差平方和作为粒子的适应值；步骤6：利用粒子群算法寻找拐角值最优解；步骤7：将步骤6.3得到的P_g值和judgewidth代入适应值函数，得到最优分段线性模型的各段斜率和截距。

Description

一种海浪有效波高反演模型建模方法

技术领域

本发明涉及一种海浪有效波高反演模型建模方法。 

背景技术

海浪是一种与人类关系最直接、最密切的海洋现象，其波高、波向、波周期等因素对航运、港口、以及海上石油平台的安全都具有非常重要的意义。船用X波段导航雷达回波形成的海杂波图像中包含丰富的海浪信息，可以利用雷达的回波强度反演海浪谱和海浪参数。1985年，Young等人首次提出了根据“海杂波”雷达图像序列提取海浪信息的方法。该方法一经发现，就引起了人们的极大兴趣。在此后10年间，Zimer、Rosenthal和Günther等人也纷纷开展了基于X波段导航雷达的海浪信息反演研究工作。1995年，德国GKSS实验室研制出了基于导航雷达的海浪监测系统WaMoS(Wave Monitoring System)；1996年，挪威Miros公司也研制出了类似的产品WAVEX。除以上主流研究机构外，目前美国、日本、丹麦、中国也在积极从事该方面的研究。 

有效波高是海浪信息的一种。由于海浪成像机制的非线性，利用X波段雷达图像进行海浪参数反演时，只能获得海浪谱能量的相对值，难以直接获取海浪有效波高。1982年，Alpers和Hasselmann提出了利用合成孔径雷达(SAR)信息估计有效波高的方法，该方法认为有效波高与雷达图像信噪比的平方根存在线性关系，可建立线性模型，通过计算得到雷达图像信噪比的平方根，进而根据线性模型得到有效波高。1994年，Ziemer和Günther将该方法推广到X波段导航雷达图像，计算获取了有效波高。到目前为止，该方法一直被作为基于X波段导航雷达图像有效波高的标准反演方法而使用。其线性模型如以下公式所示： 

$\mathrm{Hs}=A+B*\sqrt{\mathrm{SNR}}$

SNR＝SIG/BGN 

其中，Hs是有效波高，A和B是待定系数，SNR是雷达图像的信噪比，SIG是海浪波谱的能量，BGN是背景噪声的能量。 

在实际使用中发现，由于信噪比的计算方法不同，雷达系统的差异，以及海域的环境差异等因素，虽然海浪有效波高是随着雷达图像信噪比的增大而增大，但在整个变化范围内，海浪有效波高与雷达图像信噪比的平方根之间并不是完全线性的。因此，采用线性模型表达雷达图像信噪比的平方根与海浪有效波高之间的关系是不准确的，即使用雷达图像信噪比平方根反演海浪有效波高的精度偏低。 

发明内容

本发明目的在于提供一种海浪有效波高反演模型建模方法，能够有效提高使用雷达图像信噪比平方根反演海浪有效波高的精度。 

实现本发明目的技术方案： 

一种海浪有效波高反演模型建模方法，其特征在于： 

步骤1：设置粒子群粒子位置初值和速度初值，粒子位置代表拐角值； 

步骤2：初始化粒子群各参数，设种群数为N≥2；自身学习因子为c₁≥0，全局学习因子为c₂≥0；惯性权重为ω，0≤ω≤1；最大迭代次数T_max满足T_max≥2； 

步骤3：设置粒子群速度位置更新公式，更新公式为以下公式（5）， 

$\left\{\begin{array}{c}{V_i}^{t+1}=\mathrm{\ω}\·{V_i}^t+c_1{r}_1\·(P_i-X_i^t)+c_2{r}_2\·(P_g-X_i^t)\\ X_i^{t+1}=X_i^t+{V_i}^{t+1}\end{array}\right.---\left(5\right)$

其中V_i ^t表示第i个粒子的速度；

表示第i个粒子的位置，i＝1,2,...,N，N≥2为粒子个数；t表示粒子群的迭代次数，t≤T_max；ω为惯性权重；c₁≥0为自身学习因子；c₂≥0为全局学习因子；r₁,r₂为[0,1]区间的随机数；P_i为粒子自身寻到的最优值；P_g为种群寻到的最优值； 

步骤4：确定判断宽度；分段时是以判断宽度为步长进行数据域的判断并分段，根据以下公式（6）判断宽度取值范围， 

judgewidth＝n·minwidth，1≤n≤3       （6） 

judgewidth为判断宽度，minwidth为判断宽度的最小值，是样本数据所有相邻点的横轴距离的最大值； 

步骤5：求解适应值；基于判断宽度judgewidth和粒子的位置值X_i两个参数，建立适应值函数，适应值函数对样本数据进行分段，得到分段模型，并计算出与分段模型对应的整体样本数据的残差平方和作为粒子的适应值； 

步骤6：利用粒子群算法寻找拐角值最优解，通过如下步骤实现， 

步骤6.1:将每个粒子的位置值X_i和judgewidth代入适应值函数，求出粒子的适应值F(X_i,jugdewidth)，并赋值给P_i，然后将P_i(i＝1,2,…N,N≥2)的最小值赋给P_g； 

步骤6.2:进入循环，按照以上公式（5）进行更新，直到满足迭代结束条件； 

步骤6.3:当迭代次数t=T_max时，循环结束，此时P_g值为拐角turnangle的最优解； 

步骤7：将步骤6.3得到的P_g值和judgewidth代入适应值函数，得到最优分段线性模型的各段斜率和截距。 

步骤5中，适应值函数通过如下方法获得， 

步骤5.1:第一段的初步划分； 

从左端开始，以判断宽度judgewidth为步长进行判断；将起始判断宽度域称为当前拟合区，将紧邻当前拟合区的判断宽度域称为判断区，判断区的判断宽度始终是judgewidth，当前拟合区则会根据分段情况判断宽度变宽或是不变； 

划分第一段时，应用最小二乘法，根据以下公式（7）求解当前拟合区的拟合线段的直线和斜率， 

$k=\frac{\underset{j=1}{\overset{M}{\mathrm{\Σ}}}(x_j-\stackrel{\‾}{x})(y_j-\stackrel{\‾}{y})}{\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{(x_j-\stackrel{\‾}{x})}^2},b=\stackrel{\‾}{y}-k\stackrel{\‾}{x}---\left(\text{7}\right)$

其中， $\stackrel{\‾}{x}=\frac{1}{M}\left(\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{x}_j\right),\stackrel{\‾}{y}=\frac{1}{M}\left(\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{y}_j\right).$ M表示参加拟合的数据点个数，M≥2；（x_j，y_j）指的是参加拟合的数据点，i＝1,2,...,M，x_j表示雷达图像信 噪比平方根

y_j表示海浪波高（Hs）；k、b表示拟合线段的斜率和截距； 

应用最小二乘法，根据以下公式（8）求解判断区的拟合线段的斜率和截距，求解判断区的拟合线段的斜率和截距如公式（8）所示。 

$k=\frac{\underset{j=1}{\overset{M}{\mathrm{\Σ}}}(x_j-x_0)(y_j-y_0)}{\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{(x_j-x_0)}^2},b=y_0-k{x}_0---\left(8\right)$

其中，M表示参加拟合的数据点个数，M≥2；（x_j,y_j）指的是参加拟合的数据点，i＝1,2,...,M，x_j表示雷达图像信噪比平方根

y_j表示海浪波高（Hs）；k、b表示拟合线段的斜率和截距；（x₀,y₀）代表拐点坐标；把当前拟合区的最右端的数据点横轴坐标当做拐点横轴坐标x₀，拐点在当前拟合区的线段上，根据其斜率和截距求出拐点的纵轴坐标y₀； 

若判断区线段倾斜角和当前拟合区线段倾斜角之差大于或等于X_i值，则认为发生争议，此时将判断区称为争议区；若判断区线段倾斜角和当前拟合区线段倾斜角之差小于X_i值，将这两个区合并为新的当前拟合区，然后以judgewidth为步长继续判断，直到发生争议；发生争议后，对争议区进行处理，将争议区内的数据点从最右端开始逐一的减少，每减少一个点，根据以上公式（8）求争议区的直线斜率，直到争议区的直线倾斜角和当前拟合区的直线倾斜角之差小于X_i值；然后将当前拟合区和处理后的争议区合并为新的当前拟合区，第一段初步划分结束，当前拟合区成为第一段区域，应用以上公式（7）求出第一段的斜率和截距并求出末端拐点； 

步骤5.2：下一段的划分； 

以判断宽度judgewidth为步长进行判断；划分下一段时，当前拟合区的线段应过上一段的末端拐点，根据以上公式（8）求解当前拟合区和判断区的线段斜率和截距，其余划分过程同步骤5.1； 

步骤5.3：进一步判断已划分出来的段； 

已划分出来的段的斜率必须是正值,判断处理过程如下： 

若已划分出两段；如果第一段斜率为负值，则认为当前划分错误，从除去第 一段的数据最左侧开始将数据点逐个的加入到第一段，直到第一段的斜率为正值，此时第二段不再保留，重新求第一段的斜率、截距和末端拐点，并返回到步骤5.2；如果第一段斜率为正值，第二段和第一段的倾斜角差值大于或等于X_i值，则认为划分两段正确，返回到步骤5.2；若差值小于X_i值，则将此两段合并为一段，重新求斜率、截距和末端拐点，并返回到步骤5.2； 

若已划分三段或三段以上；先判断最后两段，若其倾斜角差值大于或等于X_i值，则认为分段正确，返回到步骤5.2；若小于X_i值，则将这两段合并成一段，并求出斜率、截距和末端拐角，然后重复继续判断，直到结束，然后返回到步骤5.2； 

步骤5.4：得到粒子X_i的分段模型； 

当剩余的数据区域宽度小于judgewidth，根据以上公式（8）拟合剩余的数据；若其线段倾斜角和前一段的倾斜角之差大于或等于X_i值，则剩余的区域为最后一段；若其线段倾斜角和前一段倾斜角之差小于X_i值，则将这两个区域合并为最后一段；各段划分过程结束，得到粒子X_i的分段线性模型的各段斜率和截距。 

步骤6.2中，拐角值搜索范围应为0°～10°，每次迭代过程中，对粒子进行位置更新后，判断其值是否在0°～10°范围，若超出，则对其在0°～10°范围内进行随机赋值，未超出则不变。 

步骤1中，根据以下公式（1）设置各粒子位置的初值： 

X_i=10·rand      （1） 

其中，i＝1,2,...,N，N≥2为粒子个数；rand为0到1之间的随机数； 

根据以下公式（2）设置各粒子速度的初值： 

V_i=10·rand      （2） 

其中，i＝1,2,...,N，N≥2为粒子个数；rand为0到1之间的随机数。 

所有相关系数的精度要求的指标为小于0.01的正数，标准差的精度要求的 指标为小于0.001的正数。 

本发明具有的有益效果： 

海浪有效波高是随着雷达图像信噪比平方根增加而增加的，但二者之间在不同的波高阶段并非完全呈线性关系。本发明提出的是一种基于拐角寻优的分段线性拟合的海浪有效波高反演模型建模方法，在进行分段时不同的拐角值可以划分出不同的分段线性模型，利用粒子群算法的全局智能搜索功能，可搜索出最优的拐角值，然后根据最优拐角值，自动的划分出各段，得到最优模型，并且解决分段位置确定和各段的连接等问题，进而有效提高使用雷达图像信噪比平方根反演海浪有效波高的精度。 

本发明利用粒子群算法智能搜索能力对波高进行分段反演，比传统线性反演方法的反演精度更高。本发明不需要人为确定分段的段数，并且可以精确的找出分段位置，不需人为调整，可靠性高，自主性强，提高了反演的精度。本发明只有拐点这一个寻优对象，所需设置的参数较少，寻优的解是一维的，计算复杂度低，算法负担轻，寻优过程简单。 

附图说明

图1为本发明建模方法的粒子群流程图； 

图2为本发明建模方法的适应值函数流程图； 

图3为采用本发明建模方法对样本数据点进行分段拟合所得的拟合效果图； 

图4为采用本发明建模方法对样本数据点进行分段拟合，并对数据点进行反演，得到的反演波高与实际波高的相关度曲线； 

图5为采用最小二乘法线性反演方法，对样本数据点进行拟合所得的拟合效果图； 

图6为采用最小二乘法线性反演方法，对样本数据点进行拟合，并对样本数据进行反演，得到的反演波高与实际波高的相关度曲线。 

具体实施方式

如图1所示， 

步骤1：设置粒子群粒子位置初值和速度初值，粒子位置代表拐角值； 

粒子群算法的寻优对象只有拐角（turnangle）这一个参数，所以解空间是一维的。粒子群算法的粒子位置X_i值即为拐角值解空间的一个解。根据海浪的物理信息，分段时各段斜率值都应为正值，各段斜率从左到右依次变大，故设定解空间的搜索范围为0°～10°。 

根据以下公式（1）设置各粒子位置的初值： 

X_i=10·rand      （1）其中，i＝1,2,...,N，N≥2为粒子个数；rand为0到1之间的随机数； 

根据以下公式（2）设置各粒子速度的初值： 

V_i=10·rand      （2）其中，i＝1,2,...,N，N≥2为粒子个数；rand为0到1之间的随机数。 

步骤2：初始化粒子群各参数，设种群数为N≥2；自身学习因子为c₁≥0，全局学习因子为c₂≥0；惯性权重为ω，0≤ω≤1；最大迭代次数T_max满足T_max≥2； 

步骤3：设置粒子群速度位置更新公式，更新公式为以下公式（5）， 

$\left\{\begin{array}{c}{V_i}^{t+1}=\mathrm{\ω}\·{V_i}^t+c_1{r}_1\·(P_i-X_i^t)+c_2{r}_2\·(P_g-X_i^t)\\ X_i^{t+1}=X_i^t+{V_i}^{t+1}\end{array}\right.---\left(5\right)$

对速度和位置进行更新时，采用粒子群的标准模型。其中V_i ^t表示第i个粒子的速度；

表示第i个粒子的位置，i＝1,2,...,N，N≥2为粒子个数；t表示粒子群的迭代次数，t≤T_max；ω为惯性权重；c₁≥0为自身学习因子；c₂≥0为全局学习因子；r₁,r₂为[0,1]区间的随机数；P_i为粒子自身寻到的最优值；P_g为种群寻到的最优值； 

步骤4：确定判断宽度； 

分段时是以判断宽度为步长进行数据域的判断并分段的，判断宽度judgewidth求值过程如下：为使各段拟合时有意义，保证每个判断宽度内至少有一个数据点，判断宽度的最小值minwidth应该是样本数据所有相邻点的横轴距离的最大值。 

根据以下公式（6）判断宽度取值范围， 

judgewidth＝n·minwidth，1≤n≤3（6） 

judgewidth为判断宽度，minwidth为判断宽度的最小值，是样本数据所有相邻点的横轴距离的最大值； 

步骤5：适应值函数求解适应值； 

适应值用来评价粒子的适应程度。适应值函数求值流程图，如图2所示。本发明适应值函数有两个参数：第一个是以上公式（6）中确定的判断宽度（judgewidth）；第二个参数则是粒子的位置值X_i，X_i属于0°～10°范围的解空间内的一个拐角值（turnangle），该参数为本发明的寻优对象。适应值函数根据这两个参数对样本数据进行分段，得到一个分段模型，并计算出与分段模型对应的整体样本数据的残差平方和作为粒子的适应值。具体求解过程如下： 

步骤5.1:第一段的初步划分； 

从左端开始，以判断宽度judgewidth为步长进行判断。本发明将起始判断宽度域称为当前拟合区，将紧邻当前拟合区的判断宽度域称为判断区，判断区的宽度一直是judgewidth，而当前拟合区则会根据分段情况变宽或是不变。 

划分第一段时，求解当前拟合区的拟合线段的直线和斜率应用最下二乘法，如以下公式（7）所示。 

$k=\frac{\underset{j=1}{\overset{M}{\mathrm{\Σ}}}(x_j-\stackrel{\‾}{x})(y_j-\stackrel{\‾}{y})}{\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{(x_j-\stackrel{\‾}{x})}^2},b=\stackrel{\‾}{y}-k\stackrel{\‾}{x}---\left(\text{7}\right)$

其中， $\stackrel{\‾}{x}=\frac{1}{M}\left(\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{x}_j\right),\stackrel{\‾}{y}=\frac{1}{M}\left(\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{y}_j\right).$ M表示参加拟合的数据点个数，M≥2；（x_j，y_j）指的是参加拟合的数据点，i＝1,2,...,M，x_j表示雷达图像信噪比平方根

y_j表示海浪波高（Hs）；k、b表示拟合线段的斜率和截距。 

求解判断区的拟合线段的斜率和截距，同样应用最小二乘法，但必须过当前拟合区的末端拐点以保证线段的连续性。拐点值是这样确定的：把当前拟合区的最右端的数据点横轴坐标当做拐点横轴坐标，拐点在当前拟合区的线段上，根据 其斜率和截距可以求出拐点的纵坐标，这样就得到拐点值。 

求解判断区的拟合线段的斜率和截距如以下公式（8）所示。 

$k=\frac{\underset{j=1}{\overset{M}{\mathrm{\Σ}}}(x_j-x_0)(y_j-y_0)}{\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{(x_j-x_0)}^2},b=y_0-k{x}_0---\left(8\right)$

其中，M表示参加拟合的数据点个数，M≥2；（x_j,y_j）指的是参加拟合的数据点，i＝1,2,...,M，x_j表示雷达图像信噪比平方根

y_j表示海浪波高（Hs）；（x₀,y₀）代表拐点坐标；k、b表示拟合线段的斜率和截距。 

第一段划分过程如下：若判断区线段倾斜角和当前拟合区线段倾斜角之差大于或等于X_i值，则认为发生争议，此时将判断区称为争议区；若判断区线段倾斜角和当前拟合区线段倾斜角之差小于X_i值，将这两个区合并为新的当前拟合区，然后以judgewidth为步长继续判断，直到发生争议。 

发生争议后，对争议区进行处理。将争议区内的数据点从最右端开始逐一的减少，每减少一个点，应用以上公式（8）求争议区的直线斜率，直到争议区的直线倾斜角和当前拟合区的直线倾斜角之差小于X_i值。然后将当前拟合区和处理后的争议区合并为新的当前拟合区，第一段初步划分结束，当前拟合区成为第一段区域。应用以上公式（7）求出第一段的斜率和截距并求出末端拐点。 

步骤5.2：下一段的划分； 

进行下一段的划分，同样以judgewidth为步长进行判断。为保证分段模型的连续性，划分下一段时，当前拟合区的线段应过上一段的末端拐点。划分过程和步骤5.1相同，不同之处是求解当前拟合区和判断区的线段斜率和截距都是应用以上公式（8）。 

步骤5.3：进一步判断已划分出来的段； 

已划分出来的段的斜率必须是正值,判断处理过程如下： 

若已划分出两段；如果第一段斜率为负值，则认为当前划分错误，从除去第一段的数据最左侧开始将数据点逐个的加入到第一段，直到第一段的斜率为正值，此时第二段不再保留，重新求第一段的斜率、截距和末端拐点，并返回到步骤5.2；如果第一段斜率为正值，第二段和第一段的倾斜角差值大于或等于X_i值， 则认为划分两段正确，返回到步骤5.2；若差值小于X_i值，则将此两段合并为一段，重新求斜率、截距和末端拐点，并返回到步骤5.2； 

若已划分三段或三段以上；先判断最后两段，若其倾斜角差值大于或等于X_i值，则认为分段正确，返回到步骤5.2；若小于X_i值，则将这两段合并成一段，并求出斜率、截距和末端拐角，然后重复继续判断，直到结束，然后返回到步骤5.2； 

步骤5.4：得到粒子X_i的分段模型； 

当剩余的数据区域宽度小于judgewidth，根据以上公式（8）拟合剩余的数据；若其线段倾斜角和前一段的倾斜角之差大于或等于X_i值，则剩余的区域为最后一段；若其线段倾斜角和前一段倾斜角之差小于X_i值，则将这两个区域合并为最后一段；各段划分过程结束，得到粒子X_i的分段线性模型的各段斜率和截距。 

各段划分过程结束，得到粒子X_i的分段线性模型的各段斜率和截距，同时求出所有数据点的残差平方和，作为粒子的适应值F(X_i,jugdewidth)。 

步骤6：利用粒子群算法寻找拐角值最优解； 

步骤6.1:计算每个粒子的适应值，并初始化P_i和P_g

将每个粒子的位置信息X_i值和judgewidth代入适应值函数，求出粒子的适应值F(X_i,jugdewidth)，并赋值给P_i。然后将P_i(i＝1,2,…N,N≥2)的最小值赋给P_g。 

步骤6.2:进入循环，按照以上公式（5）进行更新，直到满足迭代结束条件； 

每次迭代时，对所有粒子X_i(i＝1,2,…N,N≥2)按照以上公式（5）进行速度和位置的更新，判断位置值是否在0°～10°范围，若超出，则对其在0°～10°范围内进行随机赋值，未超出则不变。求出粒子的适应值F(X_i,jugdewidth)，若F(X_i,jugdewidth)＜P_i,则令P_i=F(X_i,jugdewidth)，否则，P_i不变。令 P_i(i＝1,2,…N,N≥2)的最小值minP_i，若minP_i＜P_g，则令P_g=minP_i，否则P_g不变。 

步骤6.3:当迭代次数t=T_max时，循环结束，此时P_g值为拐角turnangle的最优解； 

步骤7：将步骤6.3得到的P_g值和judgewidth代入适应值函数，得到最优分段线性模型的各段斜率和截距。因为适应值函数在设计时，不仅可以得到适应值，还能得到当前粒子的分段模型，所以把P_g和judgewidth带入适应值函数后，可以得到各段的斜率和截距。 

上述所有相关系数的精度要求的指标为小于0.01的正数，标准差的精度要求的指标为小于0.001的正数。 

下面结合具体实施例对本发明进一步详细说明。 

利用2009年10月在福建平潭进行科研实验时获得的雷达图像信噪比平方根与对应的的海浪有效波高真值数据（称为样本数据），采用本发明提出的基于拐角寻优的分段线性拟合的海浪有效波高反演模型建模方法与基于最小二乘法的线性方法分别进行建模，其结果如图3、图4、图5、图6所示。采用本发明提出的建模方法建立的模型为分段线性曲线，而采用线性建模方法得到一条直线。 

采用相关系数和标准差这两个评价指标来分别评价本发明提出的建模方法得到的反演模型的精度与线性建模方法得到的反演模型的精度。并将本发明提出的建模方法得到的反演模型的评价指标和线性建模方法得到的反演模型的评价指标、刘利强等人在2011年提出的一种基于径向基神经网络的海浪有效波高反演模型建模方法得到的反演模型的评价指标、刘利强等人在2012年提出的一种基于PSO自适应分段线性拟合的海浪有效波高反演模型建模方法得到的反演模型的评价指标做对比，如表1所示。 

表1评价指标 

本发明用LMS（Least Mean Square，最小均方算法）、PSO、RBF分别代替 基于最小二乘法的线性建模方法、基于PSO自适应分段线性拟合的海浪有效波高反演模型建模方法、基于径向基神经网络的海浪有效波高反演模型建模方法。 

由反演结果和表1可以看出，和线性反演方法相比，标准差和相关系数都较小，精度更高，说明本文的波高反演方法更有效，证明了利用粒子群算法寻找拐角值最优解方法在波高反演建模中的可行性。和RBF反演方法相比，相关系数相等，但标准差较小；和PSO反演方法相比，相关系数指标略差，但标准差同样比PSO反演方法小。说明本发明建立的海浪有效波高反演模型能够准确的表达雷达图像信噪比平方根与海浪有效波高之间的关系，相关系数表现良好，标准差较小，反演精度更高。 

Claims (5)

1.一种海浪有效波高反演模型建模方法，其特征在于：

步骤1：设置粒子群粒子位置初值和速度初值，粒子位置代表拐角值；

步骤2：初始化粒子群各参数，设种群数为N≥2；自身学习因子为c₁≥0，全局学习因子为c₂≥0；惯性权重为ω，0≤ω≤1；最大迭代次数T_max满足T_max≥2；

步骤3：设置粒子群速度位置更新公式，更新公式为以下公式（5），

$\left\{\begin{array}{c}{V_i}^{t+1}=\mathrm{\ω}\·{V_i}^t+c_1{r}_1\·(P_i-X_i^t)+c_2{r}_2\·(P_g-X_i^t)\\ X_i^{t+1}=X_i^t+{V_i}^{t+1}\end{array}\right.---\left(5\right)$

其中V_i ^t表示第i个粒子的速度；

表示第i个粒子的位置，i＝1,2,...,N，N≥2为粒子个数；t表示粒子群的迭代次数，t≤T_max；ω为惯性权重；c₁≥0为自身学习因子；c₂≥0为全局学习因子；r₁,r₂为[0,1]区间的随机数；P_i为粒子自身寻到的最优值；P_g为种群寻到的最优值；

步骤4：确定判断宽度；分段时是以判断宽度为步长进行数据域的判断并分段，根据以下公式（6）判断宽度取值范围，

judgewidth＝n·minwidth，1≤n≤3      （6）

judgewidth为判断宽度，minwidth为判断宽度的最小值，是样本数据所有相邻点的横轴距离的最大值；

步骤5：求解适应值；基于判断宽度judgewidth和粒子的位置值X_i两个参数，建立适应值函数，适应值函数对样本数据进行分段，得到分段模型，并计算出与分段模型对应的整体样本数据的残差平方和作为粒子的适应值；

步骤6：利用粒子群算法寻找拐角值最优解，通过如下步骤实现，

步骤6.1:将每个粒子的位置值X_i和judgewidth代入适应值函数，求出粒子的适应值F(X_i,jugdewidth)，并赋值给P_i，然后将P_i(i＝1,2,…N,N≥2)的最小值赋给P_g；

步骤6.2:进入循环，按照以上公式（5）进行更新，直到满足迭代结束条件；

步骤6.3:当迭代次数t=T_max时，循环结束，此时P_g值为拐角turnangle的最优解；

步骤7：将步骤6.3得到的P_g值和judgewidth代入适应值函数，得到最优分段线性模型的各段斜率和截距。

2.根据权利要求1所述的海浪有效波高反演模型建模方法，其特征在于：步骤5中，适应值函数通过如下方法获得，

步骤5.1:第一段的初步划分；

从左端开始，以判断宽度judgewidth为步长进行判断；将起始判断宽度域称为当前拟合区，将紧邻当前拟合区的判断宽度域称为判断区，判断区的判断宽度始终是judgewidth，当前拟合区则会根据分段情况判断宽度变宽或是不变；

划分第一段时，应用最小二乘法，根据以下公式（7）求解当前拟合区的拟合线段的直线和斜率，

$k=\frac{\underset{j=1}{\overset{M}{\mathrm{\Σ}}}(x_j-\stackrel{\‾}{x})(y_j-\stackrel{\‾}{y})}{\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{(x_j-\stackrel{\‾}{x})}^2},b=\stackrel{\‾}{y}-k\stackrel{\‾}{x}---\left(\text{7}\right)$

其中， $\stackrel{\‾}{x}=\frac{1}{M}\left(\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{x}_j\right),\stackrel{\‾}{y}=\frac{1}{M}\left(\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{y}_j\right).$ M表示参加拟合的数据点个数，M≥2；（x_j，y_j）指的是参加拟合的数据点，i＝1,2,...,M，x_j表示雷达图像信噪比平方根

y_j表示海浪波高（Hs）；k、b表示拟合线段的斜率和截距；

应用最小二乘法，根据以下公式（8）求解判断区的拟合线段的斜率和截距，求解判断区的拟合线段的斜率和截距如公式（8）所示。

$k=\frac{\underset{j=1}{\overset{M}{\mathrm{\Σ}}}(x_j-x_0)(y_j-y_0)}{\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{(x_j-x_0)}^2},b=y_0-k{x}_0---\left(8\right)$

其中，M表示参加拟合的数据点个数，M≥2；（x_j,y_j）指的是参加拟合的数据点，i＝1,2,...,M，x_j表示雷达图像信噪比平方根

y_j表示海浪波高（Hs）；k、b表示拟合线段的斜率和截距；（x₀,y₀）代表拐点坐标；把当前拟合区的最右端的数据点横轴坐标当做拐点横轴坐标x₀，拐点在当前拟合区的线段上，根据其斜率和截距求出拐点的纵轴坐标y₀；

若判断区线段倾斜角和当前拟合区线段倾斜角之差大于或等于X_i值，则认为发生争议，此时将判断区称为争议区；若判断区线段倾斜角和当前拟合区线段倾斜角之差小于X_i值，将这两个区合并为新的当前拟合区，然后以judgewidth为步长继续判断，直到发生争议；发生争议后，对争议区进行处理，将争议区内的数据点从最右端开始逐一的减少，每减少一个点，根据以上公式（8）求争议区的直线斜率，直到争议区的直线倾斜角和当前拟合区的直线倾斜角之差小于X_i值；然后将当前拟合区和处理后的争议区合并为新的当前拟合区，第一段初步划分结束，当前拟合区成为第一段区域，应用以上公式（7）求出第一段的斜率和截距并求出末端拐点；

步骤5.2：下一段的划分；

以判断宽度judgewidth为步长进行判断；划分下一段时，当前拟合区的线段应过上一段的末端拐点，根据以上公式（8）求解当前拟合区和判断区的线段斜率和截距，其余划分过程同步骤5.1；

步骤5.3：进一步判断已划分出来的段；

已划分出来的段的斜率必须是正值,判断处理过程如下：

若已划分出两段；如果第一段斜率为负值，则认为当前划分错误，从除去第一段的数据最左侧开始将数据点逐个的加入到第一段，直到第一段的斜率为正值，此时第二段不再保留，重新求第一段的斜率、截距和末端拐点，并返回到步骤5.2；如果第一段斜率为正值，第二段和第一段的倾斜角差值大于或等于X_i值，则认为划分两段正确，返回到步骤5.2；若差值小于X_i值，则将此两段合并为一段，重新求斜率、截距和末端拐点，并返回到步骤5.2；

若已划分三段或三段以上；先判断最后两段，若其倾斜角差值大于或等于X_i值，则认为分段正确，返回到步骤5.2；若小于X_i值，则将这两段合并成一段，并求出斜率、截距和末端拐角，然后重复继续判断，直到结束，然后返回到步骤5.2；

步骤5.4：得到粒子X_i的分段模型；

当剩余的数据区域宽度小于judgewidth，根据以上公式（8）拟合剩余的数据；若其线段倾斜角和前一段的倾斜角之差大于或等于X_i值，则剩余的区域为最后一段；若其线段倾斜角和前一段倾斜角之差小于X_i值，则将这两个区域合并为最后一段；各段划分过程结束，得到粒子X_i的分段线性模型的各段斜率和截距。

3.根据权利要求2所述的海浪有效波高反演模型建模方法，其特征在于：步骤6.2中，拐角值搜索范围应为0°～10°，每次迭代过程中，对粒子进行位置更新后，判断其值是否在0°～10°范围，若超出，则对其在0°～10°范围内进行随机赋值，未超出则不变。

4.根据权利要求3所述的海浪有效波高反演模型建模方法，其特征在于：步骤1中，根据以下公式（1）设置各粒子位置的初值：

X_i=10·rand      （1）

其中，i＝1,2,...,N，N≥2为粒子个数；rand为0到1之间的随机数；

根据以下公式（2）设置各粒子速度的初值：

V_i=10·rand      （2）

其中，i＝1,2,...,N，N≥2为粒子个数；rand为0到1之间的随机数。

5.根据权利要求4所述的海浪有效波高反演模型建模方法，其特征在于：所有相关系数的精度要求的指标为小于0.01的正数，标准差的精度要求的指标为小于0.001的正数。

Images