一种基于灰度共生矩阵的航海雷达反演海面风向的方法

1.本发明属于海面风向遥感技术领域，具体地说，是一种基于灰度共生矩阵的航海雷达反演海面风向的方法。


背景技术：

2.海面风场是海-汽界重要的动力学因素，对海洋军事作业、航海作业及海洋动力学研究有着至关重要的作用，故基于纹理分析的海面风场辨识极具实用价值与研究前景。传统海面风向测量主要基于船舶和浮标，但此方法受空间覆盖面积的限制，故风场研究热点主要集中在合成孔径雷达(sar)和航海雷达遥感图像上，sar可以全天时、全天候成像，但其测量范围很大。航海雷达可以对海表面的空间和时间变化进行实时监测，且不受光线影响，能够测量近船舶1千米的风场信息，为舰船舰载机提供起飞和降落时1km内的平均风场信息，具有重要的工程应用价值。
3.现阶段基于航海雷达图像反演海面风场主要有三种方法：第一类方法是基于雷达回波强度与海面风场关系建立模型获取海面风向信息。2012年lund等应用航海雷达图像回波强度与海面风场建立谐波函数，根据逆风时函数值最大这一特性提取出海面风向信息。2017年xilong liu等提出基于集合经验模态分解方法，建立雷达图像固有模态函数分量与海面风场函数关系获取海面风向信息。2019年开始x chen和weiminhuang等相继提出各类基于深度学习直接建立雷达图像与风场的模型，根据模型由雷达图像直接获取海面风向信息。通过分析发现这类方法的雷达设备适用性比较差，针对不同型号雷达数据有较大的差异，影响海面风向的反演精度。
4.第二类方法是通过雷达图像频谱获取海面风向信息，2016年yali wang提出利用波数谱区域值分别获取降雨与非降雨条件下的海面风向信息，由此可以得到波数能量谱与海面风场是存在一定函数关系的。2021年王慧等提出根据航海雷达海面风场条纹特征，利用风场能量谱直接获得海面风向信息，取得了较好的精度和数据适用性。但此方法对条纹尺度特征的依赖性很高，在特殊天气或者条纹尺度不在预估范围的精度有很大降低，部分数据不适用的情况。而且在进行频谱变换时需要对极坐标图像进行笛卡尔坐标变换，在插值变化过程中图像会存在一定的畸变，影响海面风向的反演精度。
5.第三类方法是由图像空间域特征获得海面风向信息，2004年dankert提出基于航海雷达图像中的阵风信号特征，基于光流法提取出海面风场信息，但风场属于静态特征信号，很难在单幅图像中提取出风场空间域特征，导致光流法反演精度无法达到工程要求。2007年dankert提出根据风条纹特征与海面风向平行特征，基于局部梯度法从海面静态特征图像中提取出海面风向信息。2015年卢志忠等提出一种自适应缩减算法改进的局部梯度算法，设能够连续对图像进行缩减，增强图像适应小尺度风条纹尺度变化的能力。局部梯度法依然存在依赖条纹尺度特征，同样需要极坐标图像向笛卡尔坐标变换，降低了算法精度和数据适用性。
6.纹理是人们对于事物表面不同颜色部分的方向性、周期性和精细度的认知。而在
这三个角度中，方向是纹理感知中的一个重要元素。纹理方向可被应用于多种实际场景，如海域风向、风速分析领域、陆地遥感图像辨识领域以及图像模式识别领域。灰度共生矩阵的方法可解决以上问题，但传统的灰度共生矩阵方法具有以下问题：
7.(1)辨识范围限制。传统的灰度共生矩阵计算方法能准确辨识的角度只有与水平方向成90
°
、45
°
、0
°
和-45
°
的角度。而在实际应用中，仅有部分纹理分布在这些特殊方向上，导致灰度共生矩阵通常不足以解析纹理方向。
8.(2)应用范围有限。传统的灰度共生矩阵计算仅局限于笛卡尔坐标系下，而在海面风向分析问题中，航海雷达以旋转式扫描成像，生成的图像多基于极坐标。而关于灰度共生矩阵在极坐标下的应用较少。
9.(3)方向精度不足。传统灰度共生矩阵的方向结果置信区间为10度，而基于孔径雷达的海面方向辨识则需要满足置信区间为1度的精度条件。传统的灰度共生矩阵，尚无法满足实际的精度需求。


技术实现要素：

10.为了解决上述技术问题，本发明披露了一种基于灰度共生矩阵的航海雷达反演海面风向的方法，基于改进的灰度共生矩阵方法，更有效的直接从x-band航海雷达极坐标图像中提取出海面风向信息。
11.本发明采用的具体技术方案如下：
12.一种基于灰度共生矩阵的航海雷达反演海面风向的方法，包括以下几个步骤：
13.步骤1，极坐标海面静态特征提取。海面风场属于静态特征，故需对得到的各时段的雷达图像序列进行叠加，通过较长时间序列下的海面风向状态对整体的海面风向趋势进行判断，依据小尺度风条纹与海面风向平行这一特征，以实现包含海面风场信息的海面静态特征提取。由于极坐标系与图像极坐标系对中心坐标及坐标轴方向定义不同，故对图像进行处理时，需要重新建系。原始极坐标系的原点位于图像中心，偏移角正方向为逆时针，而图像极坐标系由于受图像矩阵坐标索引的限制，建立的极坐标系只能遵循偏移角坐标系为顺时针，且起始角度必须为水平方向。故基于图像坐标系应用灰度共生矩阵时，在进行后续角度运算过程中需要进行角度换算。后续的图像预处理、灰度共生矩阵运算以及趋势融合方法均基于本部分建立的极坐标系。
14.步骤2，图像预处理。现使用滤波方法对海面静态特征图像进行预处理，对风向特征所在像素集合进行提取，为后续运算构建基础。首先，本发明针对海面静态特征图像的特点，对像素值分布进行聚类，并获取优势像素类别。针对提取出的优势像素类别，进行高斯差分滤波处理，以获取主要趋势所在像素。本发明根据类别置信度将主要趋势像素集合进行再分配，并对各集合分别进行快速收敛的灰度共生矩阵运算，在得到各类别的角度后，根据各类别的聚类结果进行自适应角度融合，并最终得到辨识出的图像总体趋势。
15.步骤3，快速收敛灰度共生矩阵。基于原始灰度共生矩阵，本发明通过粗估计-细估计循环过渡的方式，实现灰度共生矩阵的快速收敛及精度提升。
16.步骤4，自适应趋势融合。对图像预处理中得到的各趋势类别分别进行基于灰度共生矩阵的方向估计后，将各趋势类别所包含像素的数量，作为该趋势方向的置信度，进行自适应趋势融合，得到最终的图像整体趋势。图像中各趋势类别所占的权重的表达式为：
[0017][0018][0019]
其中，i＝1,2,...,n
η
,k＝1,2,...,n
η
,ηk为第k个优势群的权重，nk为第k个优势群包含的像素数量，为全部优势群像素数量的总和，mk为第k个优势群。
[0020]
最终的图像整体趋势角度可表示为：
[0021][0022]
其中，θk为各优势群进行灰度共生矩阵方向估计后的角度，为最终的图像整体趋势角度，k是对应的优势群类别。
[0023]
基于灰度共生矩阵的航海雷达反演海面风向的方法中所述步骤2包括以下步骤：
[0024]
步骤2.1，图像特征提取；
[0025]
①
绘制输入图像i的像素直方图，可用如下公式表示：
[0026]
i(i,j)＝pkꢀꢀꢀꢀꢀꢀꢀꢀ
(4)
[0027]
其中，i(i,j)为输入的极坐标海面静态特征图像，pk为第k个像素的像素值。
[0028]
②
求取q个最显著的像素值，可用如下公式表示：
[0029]
f＝max(pk,q)
ꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
[0030]
其中，f为最终求得的显著像素值，q为显著像素值的个数。
[0031]
③
对第二步求取的q个像素值进行增强处理，为了方便后续的滤波处理，将像素增强为红色。
[0032]
④
对整个像素集，将除了第二步的q个显著像素值之外的其他信息去除，令其为白色。
[0033]
步骤2.2，高斯差分滤波处理；
[0034]
①
对输入图像i求取初始相角下的二维高斯谱函数g(i)0和相角为k的二维高斯谱函数g(i)k，高斯差分滤波后的图像i'的二维高斯谱函数g(i)为初始相角和相角为k的二维高斯谱函数之差，其公式分别可表示为：
[0035][0036][0037][0038]
其中，σk＝kσ0，且取σ0＝1.6，k＝0.8，g(i(i,j))0表示初始相角下的二维高斯谱函
数，g(i(i,j))k表示相角为k的二维高斯谱函数，g(i'(i,j))表示高斯差分滤波后的二维高斯谱函数。
[0039]
②
对高斯差分滤波后的图像i'，有其服从函数：
[0040][0041][0042]
其中，ε＝-0.1，i'(i,j)表示高斯差分滤波后的图像，表示经过幅值转换的图像。
[0043]
③
对经过幅值转换的图像进行二值化处理，即设置阈值对低于该界限的设置为负界值0，对高于该界限的设置为正界值1。上述二值化处理可表达为：
[0044][0045]
其中，为阈值，为滤波后图像的矩阵像素平均值，表示二值化处理后的图像。
[0046]
步骤2.3，四邻域连通域搜索算法；
[0047]
①
数据准备：
[0048]
图像的大小为nr*nc，则灰度矩阵也应为行数为nr，列数为nc的矩阵。令经过二值化后的图像矩阵为b。首先，生成一个标记矩阵c，其大小也为nr*nc。且此标记矩阵中各个点的像素值均为1(白色)，用于对已经寻找过的像素点进行标记。
[0049]
在此算法中，i为行索引值，j为列索引值，即c(i,j)表示标记矩阵c中，第i行，第j列的像素值。由先前对于矩阵大小的定义，有1≤i≤nr，1≤j≤nc，在初始状态下，对于任意取值的i和j，恒有：c(i,j)＝1。连通域坐标储存矩阵d用于储存某一连通域所属的所有像素点坐标，储存这些坐标的意义是为了计算该连通域的重心，为之后的提取操作做准备。
[0050]
②
连通域搜索：
[0051]
当二值化后图像矩阵b(i,j)＝0且标记矩阵c(i,j)＝1时，将行索引值i和列索引值j压入坐标储存堆栈s，并将坐标信息储存到坐标储存矩阵d中，最后将标记矩阵中对应的像素点置0，视为该点已经遍历过。
[0052]
当二值化后图像矩阵b(i,j)＝0且标记矩阵c(i,j)＝0时，该点已经标记过，虽然该点在二值化后的图像中为黑色像素点，但是已经遍历过，故不重复统计。
[0053]
③
深度寻优：
[0054]
从堆栈中出栈一组坐标u，v，以该组坐标为中心，遍历其四邻域中的像素点。四邻域为以(u,v)为中心来说的上(u+1,v)、下(u-1,v)、左(u,v-1)、右(u,v+1)四个点，若中心点(u,v)已经处于图像的边缘部分，则相应的舍弃超过矩阵索引范围的邻域点，以下的说明中将四邻域坐标以其中心坐标u，v来统一代替。
[0055]
将四邻域中满足b(u,v)＝0且标记矩阵c(u,v)＝1的像素点压入堆栈s，并将坐标信息储存到坐标储存矩阵d中，最后将标记矩阵中对应的像素点置0，视为该点已经标记过。
[0056]
当二值化后图像矩阵b(i,j)＝0且标记矩阵c(i,j)＝0时，该点已经标记过，虽然该点在二值化后的图像中为黑色像素点，但是已经遍历过，故不重复统计。
[0057]
当堆栈s中还有点时，重复上述步骤，直到堆栈为空。
[0058]
基于灰度共生矩阵的航海雷达反演海面风向的方法所述步骤3包括以下步骤：
[0059]
步骤3.1，改进型灰度共生矩阵方法具体实现方式；
[0060]
①
若图像的灰度范围为0n-1，则灰度共生矩阵为n
×
n的矩阵。对于与像素重合的相对位置(r,θ)的的灰度共生矩阵，其矩阵元素(m,n)可利用像素对进行计数，如下所示：
[0061][0062]
其中，(m,n)是灰度共生矩阵元素，i为行索引值，j为列索引值，card表示一个集合元素的数量；r表示相对位置两端的距离；θ表示相对位置相对于图像水平方向的角度，逆时针方向；ω表示矩形区域内的像素集；f(x,y)为直角坐标系中位置(x,y)处的灰度值，q(r,θ)为归一化因子。
[0063]
②
采用像素间距作为距离单位，则归一化因子q(r,θ)为
[0064]
q(r,θ)＝card{(i,j)|(i,j)∈ω,(i+rcosθ,j+rsinθ)∈ω}
ꢀꢀꢀꢀꢀ
(13)
[0065]
归一化因子q(r,θ)表示满足相对位置(r,θ)的像素对个数。
[0066]
③
与垂直于纹理方向计算的灰度共生矩阵相比，沿着纹理方向的灰度共生矩阵的对角线周围聚集了较大的矩阵元素。可通过测量灰度共生矩阵中较大元素相对于对角线的程度来定量参数化该特性，满足如下公式：
[0067][0068]
其中，ω(m,n)表示矩阵元素(m,n)与对角线间距的递增函数，本发明将ω(m,n)设置为(m-n)2。相关度越高，z(r,θ)越小。
[0069]
④
利用z(r,θ)对距离r的积分，满足如下公式：
[0070]
z'(θ)＝∫z(r,θ)dr
ꢀꢀꢀꢀꢀꢀ
(15)
[0071]
实际应用中，可以使用z'(θ)的离散形式，表示为如下公式：
[0072][0073]
其中，j为上式的积分范围。
[0074]
⑤
将对应于z'(θ)的最小值θ
min
的方向作为主导纹理方向进行搜索。
[0075][0076]
步骤3.2，快速收敛方法的具体实现方式；
[0077]
①
初始粗估计。选用较为粗略的方向池精度，一般为1度间隔步长进行初步灰度共生矩阵计算，并得到粗略的方向估计结果
[0078]
②
粗估计-细估计循环估计。以粗估计得到的结果为中心，生成相对于上一次步长更精细化的方向池。细估计池的间隔步长为上一次粗估计间隔步长的十分之一。对细估计池进行灰度共生矩阵运算，得到方向估计结果。最后，将此次估计结果作为粗估计，算法再次进入到粗估计-细估计循环，直到粗估计-细估计得到的结果差分值满足收敛性要求ξ。
[0079][0080]
其中，θc表示为上一次迭代值，θr表示为下一次迭代值，ξ的值为1e-3
。当两次迭代值差距不大时，则结果收敛。
[0081]
有益效果：与现有技术相比，本发明具有以下优点：
[0082]
(1)本发明建立了图像-极坐标系，将灰度共生矩阵应用到极坐标系下，从而能够直接利用雷达图像采集到的图像，无需进行插值处理，减少了图像转化带来的误差。
[0083]
(2)本发明对传统的灰度共生矩阵计算方法进行改进，使得灰度共生矩阵在各个角度都能够进行准确辨识。
[0084]
(3)本发明将精度可细化的改进灰度共生矩阵算法与快速收敛方法相结合，通过粗估计-细估计循环过渡的方式，实现了灰度共生矩阵的快速收敛及精度提升
附图说明
[0085]
图1为本发明具体实施方式流程图。
[0086]
图2为本发明航海雷达叠加图像。
[0087]
图3为本发明坐标系转换图。
[0088]
图4为本发明局部增强灰度化前后雷达图。
[0089]
图5为本发明差分滤波后雷达图。
[0090]
图6为本发明四邻域搜索结果图。
[0091]
表1为本发明结果估计与误差分析表。
具体实施方式
[0092]
为了使本技术领域的人员更好地理解本发明方案，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0093]
如图1所示，一种基于灰度共生矩阵的航海雷达反演海面风向的方法，分为极坐标海面静态特征提取、图像预处理、快速收敛灰度共生矩阵方法、自适应趋势融合这四大块。具体实施步骤共分为十四步，第一步到第三步为极坐标海面静态特征提取；第四步到第七步为图像预处理；第八步到第十二步是快速收敛的灰度共生矩阵方法；第十三步为自适应趋势融合；第十四步为海面风向信息提取及分析。具体步骤如下：
[0094]
第一步，由自制海浪监测设备采集2010年10月22日，共10组航海雷达图像序列，每组图像序列中包含32幅航海雷达图像。雷达图像幅向以0.1
°
接收第一列海面回波信号，这条线信号角度与船艏向重合，以此开始雷达图像幅向以0.1
°
为间隔成像，最终形成3600条
幅向线航海雷达图像，并且每条线拥有不同强弱回波信号。在3600个角度中寻找与1800个角度最相近的角度，并形成新的图像矩阵。
[0095]
第二步，由于海面风场属于静态特征，故需对得到的各时段的雷达图像序列进行叠加，通过较长时间序列下的海面风向状态对整体的海面风向趋势进行判断，以实现包含海面风场信息的海面静态特征提取。再根据岸基遮挡实际情况，研究区域选取半径为600m-2100m，幅向106
°‑
291
°
。具体叠加过程见图2，如图2所示，将32个短时段的雷达图像进行叠加后，所得图像风向趋势明显，可判断性有了较大提升。
[0096]
第三步，由于极坐标系与图像极坐标系对中心坐标及坐标轴方向定义不同，故对图像进行处理时，需要重新建系。如图3所示，原始极坐标系的原点位于图像中心，偏移角正方向为逆时针。而图像极坐标系由于受图像矩阵坐标索引的限制，建立的极坐标系只能遵循偏移角坐标系为顺时针，且起始角度必须为水平方向。故基于图像极坐标系应用灰度共生矩阵时，在进行后续角度运算过程中需要进行角度换算。后续的图像预处理、灰度共生矩阵运算以及趋势融合方法均基于本部分建立的坐标系。
[0097]
第四步，灰度化图像，并对局部特征增强；
[0098]
①
绘制输入图像i的像素直方图，可用如下公式表示：
[0099]
i(i,j)＝pkꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0100]
其中，i(i,j)为输入的极坐标海面静态特征图像，pk为第k个像素的像素值。
[0101]
②
求取q个最显著的像素值，可用如下公式表示：
[0102]
f＝max(pk,q)
ꢀꢀꢀꢀꢀꢀ
(2)
[0103]
其中，f为最终求得的显著像素值，q为显著像素值的个数。
[0104]
③
对第二步求取的q个像素值进行增强处理，为了方便后续的滤波处理，将像素增强为红色。
[0105]
④
对整个像素集，将除了第二步的q个显著像素值之外的其他信息去除，令其为白色。局部增强后的雷达图见图4。
[0106]
第五步，基于局部增强的高斯差分滤波；
[0107]
①
对输入图像i求取初始相角下的二维高斯谱函数g(i)0和相角为k的二维高斯谱函数g(i)k，高斯差分滤波后的图像i'的二维高斯谱函数g(i)为初始相角和相角为k的二维高斯谱函数之差，其公式分别可表示为：
[0108][0109][0110][0111]
其中，σk＝kσ0，且取σ0＝1.6，k＝0.8，g(i(i,j))0表示初始相角下的二维高斯谱函数，g(i(i,j))k表示相角为k的二维高斯谱函数，g(i'(i,j))表示高斯差分滤波后的二维高斯谱函数。
[0112]
②
对高斯差分滤波后的图像i'，有其服从函数：
[0113][0114][0115]
其中，ε＝-0.1，i'(i,j)表示高斯差分滤波后的图像，表示经过幅值转换的图像，高斯差分滤波后的图像i'见图5。
[0116]
③
对经过幅值转换的图像进行二值化处理，即设置阈值对低于该界限的设置为负界值0，对高于该界限的设置为正界值1。上述二值化处理可表达为：
[0117][0118]
其中，为阈值，为滤波后图像的矩阵像素平均值，表示二值化处理后的图像。
[0119]
第六步，四邻域连通搜索；
[0120]
①
数据准备；
[0121]
图像的大小为nr*nc，则灰度矩阵也应为行数为nr，列数为nc的矩阵。令经过二值化后的图像矩阵为b。首先，生成一个标记矩阵c，其大小也为nr*nc。且此标记矩阵中各个点的像素值均为1(白色)，用于对已经寻找过的像素点进行标记。
[0122]
在此算法中，i为行索引值，j为列索引值，即c(i,j)表示标记矩阵c中，第i行，第j列的像素值。由先前对于矩阵大小的定义，有1≤i≤nr，1≤j≤nc，在初始状态下，对于任意取值的i和j，恒有：c(i,j)＝1。连通域坐标储存矩阵d用于储存某一连通域所属的所有像素点坐标，储存这些坐标的意义是为了计算该连通域的重心，为之后的提取操作做准备。
[0123]
②
连通域搜索；
[0124]
当二值化后图像矩阵b(i,j)＝0且标记矩阵c(i,j)＝1时，将行索引值i和列索引值j压入坐标储存堆栈s，并将坐标信息储存到坐标储存矩阵d中，最后将标记矩阵中对应的像素点置0，视为该点已经遍历过。
[0125]
当二值化后图像矩阵b(i,j)＝0且标记矩阵c(i,j)＝0时，该点已经标记过，虽然该点在二值化后的图像中为黑色像素点，但是已经遍历过，故不重复统计。四邻域搜索后的结果图见图6左侧的框图。
[0126]
③
深度寻优；
[0127]
从堆栈中出栈一组坐标u，v，以该组坐标为中心，遍历其四邻域中的像素点。四邻域为以(u,v)为中心来说的上(u+1,v)、下(u-1,v)、左(u,v-1)、右(u,v+1)四个点，若中心点(u,v)已经处于图像的边缘部分，则相应的舍弃超过矩阵索引范围的邻域点，以下的说明中将四邻域坐标以其中心坐标u，v来统一代替。
[0128]
将四邻域中满足b(u,v)＝0且标记矩阵c(u,v)＝1的像素点压入堆栈s，并将坐标信息储存到坐标储存矩阵d中，最后将标记矩阵中对应的像素点置0，视为该点已经标记过。
[0129]
当二值化后图像矩阵b(i,j)＝0且标记矩阵c(i,j)＝0时，该点已经标记过，虽然
该点在二值化后的图像中为黑色像素点，但是已经遍历过，故不重复统计。当堆栈s中还有点时，重复上述步骤，直到堆栈为空。对四邻域搜索后的范围进行深度寻优，最终提取出两组特征组如图6所示，红色和橙色的椭圆形区域为寻找出的两组特征组。
[0130]
第七步，对像素值分布进行聚类，并获取优势像素类别。针对提取出的优势像素类别，进行高斯差分滤波处理，以获取主要趋势所在像素。再根据类别置信度将主要趋势像素集合进行再分配，并对各集合分别进行快速收敛的灰度共生矩阵运算，在得到各类别的角度后，根据各类别的聚类结果进行自适应角度融合，并最终得到辨识出的图像总体趋势。
[0131]
第八步，若图像的灰度范围为0n-1，则灰度共生矩阵为n
×
n的矩阵。对于与像素重合的相对位置(r,θ)的的灰度共生矩阵，其矩阵元素(m,n)可利用像素对进行计数，如下所示：
[0132][0133]
其中，(m,n)是灰度共生矩阵元素，i为行索引值，j为列索引值，card表示一个集合元素的数量；r表示相对位置两端的距离；θ表示相对位置相对于图像水平方向的角度，逆时针方向；ω表示矩形区域内的像素集；f(x,y)为直角坐标系中位置(x,y)处的灰度值，q(r,θ)为归一化因子。
[0134]
采用像素间距作为距离单位，则归一化因子q(r,θ)为
[0135]
q(r,θ)＝card{(i,j)|(i,j)∈ω,(i+rcosθ,j+rsinθ)∈ω}
ꢀꢀꢀꢀꢀꢀ
(10)
[0136]
其中，归一化因子q(r,θ)表示满足相对位置(r,θ)的像素对个数。
[0137]
第九步，与垂直于纹理方向计算的灰度共生矩阵相比，沿着纹理方向的灰度共生矩阵的对角线周围聚集了较大的矩阵元素。可通过测量灰度共生矩阵中较大元素相对于对角线的程度来定量参数化该特性，满足如下公式：
[0138][0139]
其中，ω(m,n)表示矩阵元素(m,n)与对角线间距的递增函数，本发明将ω(m,n)设置为(m-n)2。相关度越高，z(r,θ)越小。
[0140]
第十步，利用z(r,θ)对距离r的积分，满足如下公式：
[0141]
z'(θ)＝∫z(r,θ)dr
ꢀꢀꢀꢀꢀꢀꢀ
(12)
[0142]
实际应用中，可以使用z'(θ)的离散形式，表示为如下公式：
[0143][0144]
其中，j为上式的积分范围。
[0145]
将对应于z'(θ)的最小值θ
min
的方向作为主导纹理方向进行搜索。
[0146]
[0147]
第十一步，获得粗纹理方向
[0148]
选用较为粗略的方向池精度，一般为1度间隔步长进行初步灰度共生矩阵计算，并得到粗略的方向估计结果
[0149]
第十二步，迭代索引，细化角度；
[0150]
①
以粗估计得到的结果为中心，生成相对于上一次步长更精细化的方向池。细估计池的间隔步长为上一次粗估计间隔步长的十分之一。对细估计池进行灰度共生矩阵运算，得到方向估计结果。最后，将此次估计结果作为粗估计，算法再次进入到粗估计-细估计循环，直到粗估计-细估计得到的结果差分值满足收敛性要求ξ。
[0151][0152]
其中，θc表示为上一次迭代值，θr表示为下一次迭代值，ξ的值为1e-3
。当两次迭代值差距不大时，则结果收敛。
[0153]
第十三步，对图像预处理中得到的各趋势类别分别进行基于灰度共生矩阵的方向估计后，将各趋势类别所包含像素的数量，作为该趋势方向的置信度，进行自适应趋势融合，得到最终的图像整体趋势。图像中各趋势类别所占的权重的表达式为：
[0154][0155][0156]
其中，i＝1,2,...,n
η
,k＝1,2,...,n
η
,ηk为第k个优势群的权重，nk为第k个优势群包含的像素数量，mk为第k个优势群，为全部优势群像素数量的总和。
[0157]
最终的图像整体趋势角度可表示为：
[0158][0159]
其中，θk为各优势群进行灰度共生矩阵方向估计后的角度，为最终的图像整体趋势角度，k是对应的优势群类别。
[0160]
第十四步，对最终获取的图像趋势数据选取皮尔逊系数pccs和kl散度作为误差检验指标，验证本发明提出方法的准确性。计算公式如下：
[0161][0162][0163]
其中，s为海面参考风向，为海面估计风向，σ为协方差，nr×
nc为输入雷达图像的
矩阵。pccs越接近1，越接近于0，则参考海面风向和实测海面风向的分布越接近，模型越精准，海面风向反演精度越高。kld误差分析结果如表1所示。本发明海面风向与参考风向的相关系数达到了0.99，标准差为0.49
°
，偏差为-0.05
°
，足以达到工程和环境监测要求。
[0164]
表1结果估计与误差分析表
[0165][0166]
以上显示和描述了本发明的基本原理、主要特征及优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。