基于SAR图像的海洋中尺度涡弱纹理特征检测方法研究

Research on Weak Texture Feature Detection Method for Mesoscale Eddies in the Ocean Based on SAR Images

摘要: 海洋中尺度涡是一种典型的海洋中尺度现象之一，对水下声纳探测、全球气候系统、海洋生态环境以及渔业等具有严重影响，一直是海洋界研究的热点。全天时、全天候、高分辨率、宽幅SAR在海洋中尺度涡动态监测和特征参数反演精度方面具有明显优势，但是面临复杂海洋背景下中尺度涡反差信号弱、边缘特征检测难的问题。本文提出了多参数阈值法海洋中尺度涡边缘检测方法，实验结果表明该方法能够有效抑制背景噪声，显著提升SAR图像信杂比，实现海洋中尺度涡弱纹理边缘特征检测，为海洋中尺度涡精细化观测提供了新思路。

Abstract: Oceanic mesoscale eddies, as one of the typical oceanic mesoscale phenomena, have a significant impact on underwater sonar detection, the global climate system, marine ecological environment, and fisheries, and have always been a hot topic in oceanography research. All-day, all-weather, high-resolution, wide-swath Synthetic Aperture Radar (SAR) possesses obvious advantages in the dynamic monitoring of oceanic mesoscale eddies and the inversion accuracy of characteristic parameters. However, it faces the challenges of weak contrast signals and difficult edge feature detection of mesoscale eddies in complex marine backgrounds. This paper proposes a multi-parameter threshold method for edge detection of oceanic mesoscale eddies. Experimental results show that this method can effectively suppress background noise, significantly improve the signal-to-noise ratio of SAR images, and achieve weak texture edge feature detection of oceanic mesoscale eddies, providing a new approach for refined observation of oceanic mesoscale eddies.

1. 引言

海洋中尺度涡作为重要的海洋现象广泛存在于世界大洋和边缘海中，其水平尺度从几十公里至几百公里，时间尺度为几天至数月，其流场、高度场和温盐分布等通常表现出与周围环境不同的特性[1]。表现为海水温盐场、高度场分布异常和海面粗糙度变化，可被红外辐射计、雷达高度计和SAR等卫星探测，国内外已有大量基于遥感海面温度、海面高度数据的海洋涡旋识别研究[2] [3]。但红外辐射计受光照条件和海上云雾覆盖影响，可成像窗口时间少；高度计数据时空分辨率低、时效性差，而SAR具有全天时、全天候、空间分辨率高等特点，在海洋中尺度现象动态监测和特征参数反演精度方面具有明显优势，对于SAR图像特征明显的螺旋形亚中尺度涡的检测识别方法研究比较多[4]，而对于复杂海洋背景下形态复杂、多样的弱纹理边缘中尺度涡[5]，目前还没有有效的检测方法，因此，本文构建海洋中尺度涡弱纹理边缘检测多阈值模型，实现了中尺度涡边缘检测和特征参数高精度反演，对海洋SAR数据综合应用效益提升、海洋中尺度涡精细化研究等具有重要参考价值。

2. 中尺度涡弱纹理边缘检测算法

海洋中涡旋区域表面流场的“辐聚辐散”作用，改变了海面粗糙度，导致SAR后向散射强度变化，但由于SAR侧视成像特点以及海面风场、海浪、海底地形及油污等复杂背景噪声的干扰，在高分辨率SAR图像上明或暗纹理结构特征存在不明显、边缘不连续等问题。本文针对海洋中尺度涡在高分辨率SAR图像上弱纹理特征，构建了一种包含SAR图像辐射分辨率增强、边缘检测、区域面积筛选、形态学膨胀腐蚀等多参数阈值自适应SAR图像边缘特征检测模型，表示为：

$eddy=f\left(M,T,S,B_1,B_2\right)$

式中，$eddy$ 表示中尺度涡边缘特征，$M$ 为SAR图像辐射分辨率增强因子，$T$ 为边缘检测梯度阈值，$S$ 为区域面积筛选阈值，可根据中尺度涡物理海洋学尺度定义设置，$B_1,B_2$ 为形态学膨胀和腐蚀系数。

2.1. 辐射分辨率增强

辐射分辨率增强是通过入射角效应校正、噪声抑制、辐射分辨率提升等处理，去除SAR成像机制、入射角、风场非均匀性、海浪背景杂波等多种因素的干扰，提升SAR图像质量和信杂比。其中，辐射分辨率提升是以牺牲空间分辨率为代价的，因此，在保持较高空间分辨率的情况下，合理去除噪声、提升辐射分辨率，可以有效解决SAR图像上海洋中尺度涡“看不见、看不清”的问题。

2.1.1. 入射角效应校正

辐射校正的目的是消除入射角引起的距离向上亮度的不均匀性，校正公式为：

${\sigma }^{0^{\prime }}\left(i,j\right)=\frac{{\sigma }^0\left(i,j\right)}{\overline{{\displaystyle \sum _{i=1}^N{\sigma }^0\left(i,j\right)}}}$

其中，${\sigma }^{0^{\prime }}\left(i,j\right)$ 为$\left(i,j\right)$ 点距离向校正之后的后向散射系数，${\sigma }^0\left(i,j\right)$ 为$\left(i,j\right)$ 点校正前的后向散射系数，为$\left(i,j\right)$ 点所在$i$ 列$N$ 个点的${\sigma }^0\left(i,j\right)$ 均值。

2.1.2. 噪声抑制

SAR图像噪声主要是由于SAR电荷耦合器探测元响应不一致产生的斑点噪声和周期性条带噪声，斑点噪声在辐射分辨率提升的同时得到抑制，周期性条带噪声去除方法是利用傅里叶变换和形态学处理精确定位频域图像中条带噪声对应的异常亮点，然后采用一致性扩散进行异常亮点修复，实现有效去除条带噪声的目的。

一致性扩散能够根据图像已知区域周围几何结构和纹理方向智能地、连贯地修复，有效保持了边缘信息。假设集合D表示需要进行插值修复的区域，其边界为$\Gamma$ 。将所有区域$D$ 内的像素按照距离边界$\Gamma$ 的欧氏距离$T\left(x\right)$ 的大小进行排序。$B_{\epsilon }(x)=\left\{y\in \Omega :\left|y-x\right|\le \epsilon \right\}$ 表示点$x$ 为$\epsilon$ 的邻域集合，其中$\Omega$ 表示整个图像域，即最终复原图像的像素的有限阵列。$|·|$ 为欧几里得范数。将集合$D$ 内的点进行排序后，得到$D=\left\{x_1,\cdots ,x_N\right\}$ 。令$B_{\epsilon ,h}^{<}\left(x_k\right)=B_{\epsilon ,h}\left(x_k\right)\left\{x_1,\cdots ,x_N\right\}$ ，$k=1,\cdots ,N$ 表示邻域内已经修复的点[6]。其核心思想是图像二维傅里叶变换后找到条带噪声对应亮点的位置并以合适的值替换亮点。一般处理步骤是对存在条带噪声的图像进行二维快速傅里叶变换；引用形态学处理方法，通过设置膨胀系数，利用膨胀运算增强亮点信息，随后利用边缘提取及闭运算获取亮点大区域范围，最后经一定次数的腐蚀处理得到需替换的亮点区域；针对条带信息的亮点采用一致性扩散法进行修补，从而消除亮点抑制条带噪声干扰；再经傅里叶逆变化后即可得到条带噪声处理后的图像。

2.1.3. 辐射分辨率提升

提高SAR图像辐射分辨率通过多视处理实现，即使用滑动窗口直接平均相邻像素来取代直接测量值，平均使用像素的个数称为视数[7]。视数的增加可以提升图像的辐射分辨率，抑制SAR图像的斑点噪声[8]，增强弱特征信号的强度。

2.2. 涡弱边缘特征检测

针对涡边缘具有明显的后向散射异常，采用梯度阈值法是一种有效的边缘检测方法，但阈值的选取受主观因素影响严重，本文采用二维最大类间方差法(OTSU)实现了边缘检测阈值的自适应，检测得到面积不等的系列“涡”边缘，通过设定区域面积阈值标准剔除小面积区域，获取符合物理海洋学定义的中尺度涡。另外，对于检测出的边缘线不连续、空洞、毛刺较多等现象，结合形态学算子调整膨胀和腐蚀，获取连续、光滑的边缘曲线。

2.2.1. 梯度计算

SAR图像后向散射系数梯度计算公式为[9]：

$\begin{array}{l}{\left(\frac{\partial \sigma }{\partial x}\right)}_{i,j}=\frac{{\sigma }_{i+1,j}-{\sigma }_{i-1,j}}{2\Delta x}\\ {\left(\frac{\partial \sigma }{\partial y}\right)}_{i,j}=\frac{{\sigma }_{i,j+1}-{\sigma }_{i,j-1}}{2\Delta y}\\ \left|\nabla \sigma \right|=\sqrt{{\left(\frac{\partial \sigma }{\partial x}\right)}^2+{\left(\frac{\partial \sigma }{\partial y}\right)}^2}\end{array}$

式中，$\sigma$ 为SAR图像后向散射系数，${\left(\frac{\partial \sigma }{\partial x}\right)}_{i,j}$ 、${\left(\frac{\partial \sigma }{\partial y}\right)}_{i,j}$ 分别为$x$ 、$y$ 方向$\left(i,j\right)$ 点的微分，$\left|\nabla \sigma \right|$ 为后向散射系数梯度。

2.2.2. 梯度阈值分割

OTSU是一种经典的自适应阈值分割算法，算法以类间方差作为分割效果的判定标准，当类间方差为最大值时分割效果最好[10]，可降低人工选择阈值的主观性和不确定性。分割阈值表示为：

$k=OTSU\left(I\left(x,y\right)\right)$

$I\left(x,y\right)$ 为去除背景信号之后的图像。假设图像大小为$M\times N$ ，初始化阈值为$T$ ，图像中梯度值小于初始化阈值$T$ 的个数记作$N_0$ ，大于初始化阈值$T$ 的个数记作$N_1$ ，则边缘像素占整幅图像的比例${\omega }_0$ 为：

${\omega }_0=N_0/\left(M\times N\right)$

背景像素占整幅图像的比例${\omega }_1$ 为：

${\omega }_1=N_1/\left(M\times N\right)$

其中，$N_0+N_1=M\times N$ ，${\omega }_0+{\omega }_1=1$ 。

图像的总平均梯度值为$\mu$ ，边缘平均梯度值为${\mu }_0$ ，背景平均梯度值为${\mu }_1$ ，则类间方差$g$ 表示为：

$g={\omega }_0\times {\left({\mu }_0-\mu \right)}^2+{\omega }_1\times {\left({\mu }_1-\mu \right)}^2$

$\mu ={\omega }_0\times {\mu }_0+{\omega }_1\times {\mu }_1$

联立两式可得：

$g={\omega }_0\times {\omega }_1+{\omega }_1\times {({\mu }_0-{\mu }_1)}^2$

采用遍历的方法得到使类间方差$g$ 最大的阈值锋值$k$ 作为最佳梯度阈值$T_{\sigma }$ ，依据$T_{\sigma }$ 将后向散射系数梯度大于$T_{\sigma }$ 判定为涡特征区域，小于$T_{\sigma }$ 判定为海洋背景区域[11]，即：

$\sigma =\{\begin{array}{cc}\sigma & \sigma >T_{\sigma }\\ 0& \sigma \le T_{\sigma }\end{array}$

利用阈值分割结果将图像$I\left(x,y\right)$ 转为二值图，其中疑似涡区信号的像素点值为1，非涡区信号的像素点值为0。

2.2.3. 二值形态学处理

二值形态学最基本的一对运算膨胀和腐蚀[12]，结合开运算、闭运算对梯度分割后的二值图像红存在的边缘毛刺、空洞、不连续等现象进行处理[13]。

(1) 膨胀

图像$A$ 用结构元素$B$ 来膨胀记为$A\oplus B$ ，定义为：

$A\oplus B=\left\{x|{\left(\widehat{B}\right)}_x\cap A\ne \phi \right\}$

其中$\widehat{B}$ 是$B$ 的反射，$B$ 的原点平移至$x$ 点时，$B$ 与$A$ 的交集非空，则保留该点，结构元素$B$ 遍历图像$A$ 得到膨胀后的集合。

(2) 腐蚀

图像$A$ 用结构元素$B$ 来腐蚀记为$A\otimes B$ ，定义为：

$A\otimes B=\left\{x|{\left(\widehat{B}\right)}_x\subseteq A\right\}$

当$B$ 的原点平移至$x$ 点时，如$B$ 完全包含在$A$ 中，则保留该点，结构元素$B$ 遍历图像$A$ 得到腐蚀后的集合。

(3) 开运算与闭运算

开运算是先腐蚀，后膨胀，目的是消除小面积区域、平滑涡区边界，定义为：

$A\circ B=\left(A\otimes B\right)\oplus B$

闭运算是先膨胀，后腐蚀，目的是填补空洞和断裂缝隙、平滑涡区边界，定义为：

$A•B=\left(A\oplus B\right)\otimes B$

2.3. 涡特征参数提取

根据图像空间几何位置关系，提取中尺度涡边缘经纬度位置，利用椭圆拟合得到的椭圆中心经纬度位置和长短半径作为涡特征参数。

3. 实验结果与分析

3.1. 实验数据

实验数据主要包括不同时间、不同海域的ENVISAT、Radarsat-2、GF3等国内外SAR数据，包含不同强度、形态的弱纹理特征中尺度涡。

3.2. SAR数据预处理

Figure 1. Radarsat-2 SAR raw image and Incident angle correction results

图1. Radarsat-2 SAR原始图像及入射角校正结果

图1(a)为2021年3月18日Radarsat-2卫星SAR原始图像，幅宽500 km，分辨率50 m。由于入射角的变化引起SAR图像近端和远端明暗差异较大，经地理编码、入射角校正等预处理结果如图1(b)所示，图像质量得到有效改善，可视性和纹理特征得到明显增强。

图2(a)为一景Envisat卫星原始SAR VV极化图像，空间分辨率为16 m，辐射分辨率约为1.9 dB，入射角为25.79˚~31.66˚，显然原始图像中含有斑点噪声，纹理特征不明显，通过噪声抑制、入射角效应校正等预处理，辐射分辨率提升至0.45 dB，可以看出图像中纹理特征得到明显增强，如图中标注的A、B、C、D区域，SAR图像中海洋现象的可视性进一步加强。

Figure 2. Envisat SAR raw image and radiation resolution enhancement results

图2. Envisat SAR原始图像及辐射分辨率增强结果

图3为一景GF3卫星原始SAR图像，图像中含有倾斜条带噪声，如图4(a)为白色方框区域对应的细节放大图像，可见明显的倾斜条带噪声，经本文算法处理后结果如图4(b)，表明在去除条带噪声的同时没有过多抑制图像纹理细节。

Figure 3. GF3 SAR local raw image and suppression of oblique stripes

图3. GF3 SAR局部原始图像及倾斜条带抑制结果

为了验证图像预处理方法的有效性，本文采用信噪比(Signal to Noise Ratio, SNR) [14]定量评价图像预处理质量，图像的信噪比定义为图像中信号与噪声的功率谱之比，通常用图像中信号和噪声的方差之比来近似。一般来说，信噪比越大，图像的去条带效果越好。与基于低秩单图像分解模型法(the Low-Rank Single Image Decomposition, LRISD) [15]和适合于斜条带的剪切式低秩模型法(A Shearing Low-Rank Model, ASL) [16]噪声抑制效果比较，如表1所示，在运行时间、信噪比和图像纹理细节保留方面本文算法优势明显。

Table 1. Evaluation of noise suppression result indicators

表1. 噪声抑制结果指标评价

方法	运行时间(s)	SNR
LRISD	29.763	22.302
ASL	33.238	19.109
本文	8.462	24.695

本文针对多源SAR图像中存在的近端和远端亮度不均衡，倾斜条带噪声、信杂比低等典型干扰，通过入射角校正、条带噪声抑制、辐射分辨率增强等预处理算法，在保持纹理细节的情况下，实现了SAR图像弱反差信号增强、信杂比提升，有效改善了图像可视性。

3.3. 海洋中尺度涡边缘检测

在图1预处理后的图像中截取涡区图像进行涡边缘弱纹理特征检测，涡区为暗斑区域，边缘纹理特征较弱。本文处理结果如图4所示，其中，图4(a)为截取的涡区域SAR图像，图4(b)为涡区域SAR后向散射系数梯度分布，图4(c)为最大类间方差法梯度阈值分割及小面积区域剔除结果，图4(d)、图4(e)为形态学膨胀、腐蚀结果，图4(f)为涡边缘椭圆拟合结果在SAR图像上的叠加，同时标识了中尺度涡中心位置($x_0,y_0$ )和尺度(椭圆长短半径a、b)。

Figure 4. Results of the detection process of the weak edge of the medium-scale vortex in Radarsat-2 SAR

图4. Radarsat-2 SAR中尺度涡弱边缘检测过程结果

图5(a)为2012年6月5日Radarsat-2卫星原始SAR图像，标称分辨率50 m、幅宽300 km，包含一个螺旋形海洋中尺度涡，截取涡区图像并进行预处理如图5(b)所示，中尺度涡边缘检测及椭圆拟合结果，如图5(c)、图5(e)所示。

图6(a)为2021年3月17日的GF3卫星原始SAR图像，宽幅扫描模式成像模式，标称分辨率100 m、幅宽500 km，包含一个海洋中尺度涡，截取涡区图像及预处理结果如图6(b)所示，涡边缘检测结果如图6(c)所示，涡心位置如红色点所示，椭圆拟合结果如图6(d)所示。

海洋中尺度涡边缘检测精度采用时空匹配卫星遥感海面温度数据边缘检测结果比对方法，首先获取时空匹配的SAR数据与海面温度数据，分别检测得到涡的边缘，然后采用均方根误差公式统计中尺度涡中心位置和尺度偏差。其中，涡中心位置是指涡区形心的位置，尺度是指过涡旋形心任意直线与涡边缘两侧交点之间的距离，示意图如图7所示，图7(a)为SAR后向散射分布图，图7(b)为海面温度梯度分布图，偏差统计如表2、表3所示，可见SAR中尺度涡中心位置偏差小于2 km，尺度偏差小于4 km，能够满足海上应用保障需求。

Figure 5. Results of the measurement of the weak edge in the Radarsat-2 SAR

图5. Radarsat-2 SAR中尺度涡弱边缘检测过程结果

Figure 6. Weak edge detection results of mesoscale vortex in GF3 SAR

图6. GF3 SAR中尺度涡弱边缘检测结果

Figure 7. Graphical representation of statistical deviations

图7. 偏差统计示意图

Table 2. Statistical deviation of the center of the mid-scale axis

表2. 中尺度涡中心位置偏差统计

日期	SAR检测的涡中心位置		SST检测的涡中心位置		距离偏差(km)
	纬度(˚)	经度(˚)	纬度(˚)	经度(˚)
2012.06.05	18.5386	110.9905	18.5405	110.9923	0.788
2019.03.13	18.5876	18.5876	18.5786	111.6935	1.21
2021.01.25	20.6578	114.9995	20.6409	114.9998	1.875
2021.03.18	20.2965	117.5499	20.286	117.5499	1.166
涡中心位置均方根误差(km)：1.523

Table 3. Statistical deviation of the mean and median scales

表3. 中尺度涡尺度偏差统计

日期	SAR检测的涡尺度(km)	SST检测的涡尺度(km)	尺度偏差(km)
2012.06.05	183.5336	184.1227	0.5891
2019.03.13	119.0697	123.0548	5.6049
2021.01.25	136.4602	139.6712	3.211
2021.03.18	22.574	23.292	0.718
涡尺度均方根误差(km)：3.768

4. 结论

本文在分析海洋中尺度涡SAR图像特征的基础上，针对环形弱纹理结构特征不明显、不连续、背景噪声干扰等问题，首先通过辐射分辨率增强方法，校正了因SAR入射角、成像机制引起的图像两侧明暗差异、斑点噪声和条带噪声，提升了图像辐射分辨率和可视性，基于最大类间方差法梯度阈值、区域面积阈值、形态学膨胀腐蚀系数阈值等多阈值自适应，实现了中尺度涡SAR图像弱边缘特征的有效检测。实验结果表明本文方法对SAR图像中弱纹理特征海洋中尺度涡边缘检测的有效性，为SAR图像预处理、海洋中尺度涡特征参数高精度反演奠定了基础，对海洋中尺度涡的精细化观测具有重要意义。