基于准线性近似方法的SAR图像海浪参数反演研究❋

袁超文, 张彦敏❋❋, 姜文正, 王运华,3

(1.中国海洋大学信息科学与工程学部, 山东 青岛 266100; 2.自然资源部第一海洋研究所, 山东 青岛 266061; 3.青岛海洋科学与技术试点国家实验室 区域海洋动力学与数值模拟功能实验室, 山东 青岛 266237)

海浪作为海洋表面普遍存在的波动现象,与人类的海上活动和气候研究息息相关。合成孔径雷达(Synthetic Aperture Radar, SAR)具备全天时、全天候获取海面高分辨率微波后向散射图像的能力,是遥感观测海浪的有力工具。

海浪的SAR微波成像机制主要包括三种调制机制,即:流体力学调制、倾斜调制和速度聚束调制[1-3]。基于以上调制理论,Hasselmann等推导给出了SAR图像谱和海浪谱之间的非线性映射关系,建立了基于SAR图像的海浪谱迭代反演算法,即MPI反演法[4]。然而,由于MPI算法中的非线性映射关系无法直接求逆,因此,使用MPI方法反演海浪谱时依赖于初猜谱的输入。此外,基于MPI反演算法所得海浪谱存在原点对称的两个海浪谱谱峰,仅依靠MPI方法无法取舍两个海浪谱谱峰,因此,MPI SAR海浪谱反演算法存在海浪谱180°模糊问题。为了解决MPI方法中存在的海浪谱180°模糊问题,基于SAR子孔径图像的交叉谱方法被进一步应用于SAR海浪谱反演方法之中[5-7]。由于速度聚束效应存在显著非线性特征,这导致偏离SAR距离向传播的风浪在SAR图像中难以被成像,因此,为了进一步提升SAR海浪谱反演精度,研究者通过引入风场信息,提出了半参数化法海浪反演算法[8-10]。近年随着加拿大Radarsat-2和我国高分3号卫星在轨运行,为海浪谱和海浪参数提取提供了全极化数据。进而,诸多学者提出了基于极化SAR图像的海浪谱反演方法[11-16]。基于极化SAR图像的海浪谱反演算法可以消除流体力学调制效应的影响。近年来,有学者也提出了基于SAR图像参数的海浪有效波高经验反演算法[17-20],然而,这些海浪有效波高经验反演算法不能获得二维海浪谱。本质上基于单极化或多极化SAR图像的海浪谱反演算法都是在SAR图像纹理特征的基础上开展的,因此,这些反演算法的精确性必然会受到SAR图像中非海浪纹理的干扰。Li等[21]指出:SAR图像中的大尺度大气或海洋现象会导致SAR图像谱低波数域谱密度偏大。由SAR图像谱和海浪谱间的映射关系可知,低波数域的非海浪纹理信号谱密度在所反演的海浪谱中更容易被放大,从而影响海浪反演精度。此外,具有随机特征的相干斑噪声是海面SAR图像所固有的物理现象,然而这些随机相干斑噪声也会显著增大反演所得海浪谱的高波数域谱密度。

根据以上分析可见,当前基于SAR图像的海浪谱高精度反演算法依然存在挑战,其中,复杂的非线性映射关系依赖于初猜谱,而且计算过程较为困难;此外,SAR图像中的非海浪纹理和相干斑噪声也会影响反演所得海浪谱的精度。相对而言,准线性映射关系较为精简,本文分析了直接基于准线性映射关系从Sentinel-1 SAR波模式数据中反演海浪参数的可行性;同时,为了抑制SAR图像中大尺度非海浪纹理以及相干斑噪声的影响,根据大尺度非浪海洋现象的时间慢变特征和相干斑噪声的时间快变特征,基于SAR图像交叉谱提出了相应的抑制方案。

1 数据介绍

1.1 Sentinel-1卫星SAR数据

Sentinel-1由A星和B星构成,在欧洲航天局哥白尼计划的支持下两颗星分别发射于2014和2016年。Sentinel-1搭载了可多模式成像的C波段SAR,其中,该SAR的波模式是专门针对海浪观测的成像模式。波模式图像的中心入射角为23°(WV1)或36°(WV2),WV1和WV2间隔100 km交替对海成像。Sentinel-1波模式图像的空间分辨率和图像大小分别为5 m×5 m和20 km×20 km,其极化方式为单极化(HH或VV)。本文海浪反演时所使用的SAR数据为2020年Sentinel-1获取的全球VV极化WV2波模式L2级产品,该产品中提供了SAR图像交叉谱、方位向截断波长等信息。 其中,1月的数据用于误差校正模型的拟合,其余月份数据用于模型验证。需要说明的是:为了筛选出海浪纹理较清晰的SAR图像,同时避免海冰影响,在反演过程中,只保留了位于纬度(-60°N—+60°N)之间且图像归一化方差(nv)介于1～1.4之间的SAR数据[19]。在此,归一化方差nv定义为:

(1)

其中:I为SAR图像强度;var表示方差。

1.2 ECMWF再分析数据

本文使用了ECMWF提供的第五代再分析风场数据和有效波高数据用于海浪有效波高校正模型的拟合和结果验证,其中,风场数据的时空分辨率为1 h和0.25°×0.25°;海浪有效波高的时空分辨率为1 h和0.5°×0.5°。本文在与Sentinel-1 SAR数据进行时空匹配时,选取了ECMWF数据中与SAR数据的时空 间隔最近点。

1.3 NDBC浮标数据

为了进一步检验本文海浪反演结果的可靠性, NDBC提供的浮标数据也被用于最终的海浪有效波高结果检验。NDBC浮标每间隔1 h提供一个海浪有效波高数据,因此,本文选取浮标数据与SAR图像成像之间的时间间隔不超过30 min。同时,被选取的浮标位置与SAR图像中心位置不超过100 km。图1中的蓝点标注了匹配成功的浮标位置。

2 反演理论与方法

2.1 准线性近似法

SAR图像谱和海浪谱之间的非线性映射关系为[4]:

{1+fR(r)+ikyβ[fRv(r)-fRv(-r)]+

(kyβ)2[fRv(r)-fRv(0)][fRv(-r)-fRv(0)]}dr。

(2)

式中:PS(k)为SAR图像谱;i为虚数单位;r为空间位置矢量;k为海浪波数向量;ky为方位向波数;ξ2=β2fv(0),ξ为目标成像位置沿方位向偏移量的均方根;v为散射单元的雷达视向运动速度;β=R/V,R为平台到目标的斜距,V为平台运动速度。

(3)

(4)

(5)

(6)

(7)

(8)

直接应用式(2)反演海浪谱是困难的,然而,如果忽略式(2)中的高阶非线性项后,则可得到SAR图像谱和海浪谱之间的准线性映射关系[4]:

(9)

其中:

(10)

(11)

(12)

其中,方位向截断波长λc可通过高斯函数对SAR图像沿方位向的自相关函数进行拟合得到[22]。然而,在实际反演过程中,为了更好地抑制方位向高波数域噪声,本文将指数因子e-(kyξ)2的值设置为:

(13)

2.2 基于交叉谱的噪声抑制

SAR图像交叉谱和海浪谱之间的准线性映射关系为[7]:

(14)

(15)

由于海面上不同散射点的随机运动导致SAR不同时刻获取的子视图像上的相干斑噪声不相关,这一特征使得交叉谱受相干斑噪声的影响很小[6]。对比式(14)、式(15)和式(9)可见:

(16)

其中,Psc(k)为实际SAR图像的交叉谱。与SAR强度图像谱PS(k)相比,由于交叉谱Psc(k)能很好地抑制相干斑噪声,因此,本文应用交叉谱Psc(k)的模值代替图像谱PS(k)。

除相干斑噪声外,如图2(a)所示,SAR海浪图像还常常呈现出大尺度大气或海洋现象的纹理,这些纹理会造成SAR图像谱密度在低波数域增大。由于SAR海浪成像调制函数在低波数域的幅值非常小,因此,在通过对式(9)求逆反演海浪谱的过程中,低波数域的非海浪特征谱密度更容易被放大,从而影响海浪反演结果。另一方面,由于大尺度大气或海洋现象具有慢时变特征,因此,基于实测SAR数据获得的交叉谱Psc(k)中并不能有效消除这些低频纹理的影响。由此可见,为了进一步提升海浪反演结果的精度,需要针对大尺度大气或海洋现象的纹理提出相应抑制方法。由式(15)可见,海浪SAR图像交叉谱虚部是由于海浪在两个子视间隔时间内的规则运动产生的,因此,一些变化相对缓慢的大气或海洋现象,由于其角频率ω较小,在相同子视图像时间间隔τ的情况下,传递到交叉谱虚部的能量比重更低。如图2(b)所示,在图像谱的低波数域有一对噪声峰,而在对应的交叉谱虚部图2(c)中看不到这个现象。根据这一特征,论文提出了针对大尺度大气或海洋现象的抑制方法。

图2 Sentinel-1 VV极化图像(a),SAR图像谱(b),交叉谱虚部(c)和抑制低频噪声后的图像谱(d)Fig.2 Sentinel-1 VV polarization image(a), SAR image spectrum(b), Imaginary part of the SAR cross spectrum(c) and SAR image spectrum after suppressing the low-frequency noise(d)

首先定义交叉谱虚部的模与SAR图像谱的比值作为阈值,来判断是否去除该区域的能量,由于交叉谱虚部上存在噪声,在求比值前先减去一个与虚部能量成正比的平均噪声平面。定义该阈值参数Rk为:

(17)

在此,mean(·)为求平均值。并通过下式对SAR图像谱进行滤波:

(18)

其中λ为海浪波长。经过对大量数据的观察验证,本文中α设置为4。β1设置为0.12,β2设置为0.05。

图2(d)为低波数域噪声去除后的图像谱,可以看到,噪声峰被去除而海浪图像谱的边缘被较好地保留。

3 仿真分析

为了检验准线性近似法的反演能力,本节通过仿真实验对不同海况下的反演情况进行分析。仿真中使用的雷达和平台参数与Sentinel-1 WV2数据相同,入射角为36°,β为116 s,极化方式为VV。输入的海浪谱由风浪和涌浪两部分组成,其中风浪使用Elfouhaily 谱,涌浪使用窄带谱。生成二维海浪谱后,使用SAR图像谱和海浪谱之间的非线性映射关系式(2)前向映射得到仿真图像谱。进而将仿真图像谱与准线近似法相结合,反演海浪谱,并将反演所得海浪谱有效波高参数与输入海浪参数进行对比,检验该方法的准确性。

图3展示了在海面10 m高处风速(U)为5 m/s,风向和涌浪主波波向与距离向夹角为60°,涌浪有效波高(Hswell)为3.00 m,主波波长(λswell)为250 m时的反演情况。此时,风浪和涌浪总有效波高为3.06 m,由于风速较低,风浪谱的值远小于涌浪谱的值,且风浪主波部分位于高波数域,故图3(a)中未显示出风浪谱。对比图3(b)和图3(c)可以发现,使用非线性映射关系仿真得到的图像谱和使用准线性映射关系仿真得到的涌浪区域图像谱相似度很高,说明在该海况下,在进行准线性近似时舍弃的高阶项所占比重并不大,但从图3(e)中可以看到,由于高阶项的存在,导致涌浪谱附近一些区域谱值偏大,在去除谱值小于输入海浪谱峰值千分之一的区域的能量之后得到图3(f)所示的海浪谱,计算其有效波高为3.02 m,与输入的有效波高3.06 m差别不大,说明该方法成功地反演出了输入海浪谱的涌浪谱部分。图4展示了其它参数不变,风速增加到15 m/s时的情况,此时,风浪谱的谱值增大,导致了方位向截断效应增强,SAR图像谱上峰值处的波向比海浪谱上的更靠近距离向,准线性映射图像谱(见图4(c))与SAR图像谱图(见图4(b))仍有一定的相似性,但相对于低风速时,相似度更低。从图4(f)可以看到,最终反演的海浪谱发生了畸变,由于方位向高波数域海浪谱的缺失,反演有效波高也显著小于输入海浪有效波高。

((a)输入的海浪谱;(b)由式(2)仿真得到的SAR图像谱;(c)由式(9)仿真得到的准线性SAR图像谱;(d)准线性映射窗反演得到的海浪谱;(f)去除(e)中低谱值区能量后的海浪谱。(a) The input ocean wave spectrum; (b) Simulated SAR image spectrum by Eq.(2); (c) Simulated quasi-linear SAR image spectrum by Eq.(9); (d) Quasi-linear mapping (e) Retrieved ocean wave spectrum; (f) The ocean wave spectrum after removing the energy in the low spectral density region.)图3 基于海浪谱的SAR图像谱仿真结果Fig.3 The ocean wave spectrum and the image spectrum

图4 与图3相同,但风速为15 m/sFig.4 Same as Figure 3, but with a wind speed of 15 m/s

图5展示了准线性近似法在不同海况下的反演情况,在此,反演有效波高误差Herror定义为:去除低谱值区后反演的有效波高与输入的有效波高之差。 从图5(a)和图5(b)可看到,风速为5 m/s时,在涌浪波长较长,波高较小的情况下,有效波高反演精度高,在涌浪波高极高,波长又很短时,反演结果出现较大误差,而此类海况在现实中很少出现。而从图5(c)和图5(d)可看到,风速为15 m/s时,反演的有效波高通常偏低,浪向与距离向的夹角为90°时偏低得比30°时更多,这是因为海浪谱越靠近方位向,缺失的成分就越多。

((a) 风速5 m/s,风向和涌浪主波波向与SAR距离向夹角为30°时,反演误差随涌浪有效波高从和涌浪主波波长的变化;(b) 同(a),但是风向和涌浪主波波向与距离向夹角为90°;(c) 同(a),但是风速为15 m/s;(d) 同(c),但是风向和涌浪主波波向与距离向夹角为90°。 (a) When the wind speed is 5 m/s, the wind direction and the dominant wave direction of the swell are at an angle of 30° from the SAR range direction, the variation of the error with significant wave height and the wave wavelength of swell; (b) Same as (a), but the wind direction and the dominant wave direction of the swell are at an angle of 90° from the range direction; (c) Same as (a), but with a wind speed of 15 m/s; (d) Same as (c), but the wind direction and the dominant wave direction of the swell are at an angle of 90° from the range direction.)图5 不同风速条件下有效波高反演误差Fig.5 The error of the retrieved significant wave height at different wind speed

本文进一步分析了风场对准线性近似法反演效果的影响,图6展示了输入的涌浪有效波高为3.00 m,主波波长为200 m,主波波向与距离向的夹角为45°时反演误差与风速和风向与距离向的夹角(Wind direction)的关系。随着风速的增加,反演的有效波高逐渐偏低,而风向与距离向的夹角越大,偏低得越多。

图6 不同风向时有效波高反演误差随风速的变化Fig.6 The errors of the retrieved significant wave height versus wind speed for different wind directions

通过对仿真图像谱的反演可以发现,该方法在涌浪主导情况下反演结果较好,但是在高风速情况下,由于方位向高波数域海浪信息缺失和准线性映射关系与非线性映射关系相似度降低,导致反演结果不理想。

4 Sentinel-1波模式数据反演

4.1 准线性近似反演

在完成噪声抑制后,首先筛选出风速小于7 m/s的数据用准线性近似法进行反演,并将反演的准线性近似有效波高(Hs ql)与ECMWF提供的有效波高(Hs EC)进行对比验证。本文中使用的风场数据为ECMWF提供的第五代再分析数据,对于Radarsat-2和高分3号等卫星获取的全极化数据,风场可直接从SAR图像反演[24]。图7(a)和图7(b)分别为对1月数据使用理论倾斜调制函数和经验倾斜调制函数的有效波高反演结果。本文用于评价反演结果的参数包括偏差(BIAS)、均方根误差(RMSE)和相关系数(COR)。对比发现,使用经验倾斜调制函数的反演结果虽然比使用理论倾斜调制的反演结果具有更大的负偏差,这与使用该方法反演的海浪谱缺失高波数域分量有关,但仍然具有更低的均方根误差和更高的相关性。

((a)基于理论倾斜调制函数的反演结果; (b)基于经验倾斜调制函数的反演结果。(a) The retrieved result based on the theoretical tilt modulation function; (b) The retrieved result based on the empirical tilt modulation function.)图7 1月风速小于7 m/s时基于准线性近似反演所得有效波高与ECMWF结果的比较Fig.7 The comparisons between the retrieved significant wave height by the quasi-linear approximation and the results of ECMWF at wind speed less than 7 m/s in January

为了探索适用范围更广的有效波高反演方法,我们将反演范围扩展到风速小于20 m/s的数据,反演结果展示在图8中,此时负偏差和均方根误差相对于低风速情况下都有明显的增大。

((a)基于理论倾斜调制函数的反演结果; (b)基于经验倾斜调制函数的反演结果。(a) The retrieved result based on the theoretical tilt modulation function; (b) The retrieved result based on the empirical tilt modulation function.)图8 1月风速小于20 m/s时基于准线性近似反演所得有效波高与ECMWF结果的比较Fig.8 The comparisons between the retrieved significant wave height by the quasi-linear approximation and the results of ECMWF at wind speed less than 20 m/s in January

4.2 有效波高校正

为了提升有效波高反演精度,分析了反演误差与环境参数以及图像参数之间的联系,图9中展示了使用经验倾斜调制函数的有效波高反演误差与风速和截断波长的关系,图中颜色表示截断波长,随着风速的增大,反演的有效波高逐渐偏小,这与仿真实验中看到的现象相似。而在风速较高的情况下,截断波长较长的案例更多地位于散点图的上层。因此,我们将有效波高反演误差(Herror ql)与风速和截断波长进行多项式拟合,得到有效波高反演误差的经验关系,将反演的有效波高减去通过经验关系拟合的误差,从而提升有效波高反演精度。

(红线为风速与有效波高误差的二次多项式拟合曲线。The red line is the quadratic polynomial fit line of the wind speed and the error of the retrieved significant wave height.)图9 有效波高反演误差与风速和截断波长之间的关系Fig.9 The relationship between the error of the retrieved significant wave height and wind speed as well as cut off wavelength

拟合的公式和具体参数为:

(19)

其中,拟合系数为:

p00=-0.296 3;p10=-0.258 1;p01=0.022 08;

p20=0.051 62;p11=-1.905×10-5;

p02=-0.000 243;p30=-0.002 31;

p21=-0.000 430 9;p12=1.727×10-5;

p03=9.042×10-7;p40=0.000 561;

p31=-4.655×10-5;p22=2.865×10-6;

p13=-9.022×10-8;p04=-1.015×10-9;

p50=-6.942×10-6;p41=-4.343×10-7;

p32=7.412×10-8;p23=-3.819×10-9;

p14=1.099×10-10;

图10中给出了校正后的有效波高结果(Hs correct)与ECMWF再分析数据的比较,由图可见:不论是高风速还是低风速的情况,校正后的有效波高精度明显提升。

(色标代表风速, The color bar represents wind speed。(a)风速小于20 m/s时的结果; (b)风速小于7 m/s时的结果。(a) The results when wind speed is lower than 20 m/s; (b) The results when wind speed is lower than 7 m/s.)图10 校正后的1月份SAR数据反演所得有效波高与ECMWF有效波高对比Fig.10 The comparison between the corrected significant wave height retrieved based on the SAR data in January and the results of ECMWF

我们将误差校正模型应用到4、7、10月的数据中以检验其适用性。反演结果展示在图11中展示,不同风速限制下校正前后的结果比较统计在表1中。可以看到,由1月数据拟合出的经验校正模型在其它季度也是有效的,校正之后各统计参数都有明显的改善,尤其是在高风速情况下。对于低风速情况,四个月的数据均无需校正就可以得到比较准确的有效波高反演值,但是从图11中可以看到,当风速极低时(小于2 m/s),出现了较多的偏离验证数据的样本。从7月的反演结果图11(c)和图11(d)中都能更明显地看到一团低风速的离群值,它们的反演结果高于验证数据,使得该月的RMSE相对于其它几个月稍大。

表1 反演所得有效波高的统计结果Table1 Statistics result of the retrieved significant wave height

((a)4月,风速小于20 m/s时的结果;(b)4月,风速小于7 m/s时的结果;(c)7月,风速小于20 m/s时的结果;(d)7月,风速小于7 m/s时的结果;(e)10月,风速小于20 m/s时的结果;(f)10月,风速小于7 m/s时的结果。(a) Results when Wind speed less than 20 m/s in April; (b) Results when Wind speed less than 7 m/s in April; (c) Results when Wind speed less than 20 m/s in July ;(d) Results when wind speed less than 7 m/s in July; (e) Results when wind speed less than 20 m/s in October; (f)Results when wind speed less than 7 m/s in October.)图11 4、7、10月等三个月SAR数据反演有效波高与ECMWF有效波高对比Fig.11 The comparison of the corrected significant wave height retrieved based on the SAR image in April、July and October and the results of ECMWF

4.3 浮标测量数据验证

Ren等[21]提出了截断波长和有效波高间的线性模型:

(20)

并在Radarsat-2数据上验证了其有效性。本节使用ECMWF数据拟合了适用于Sentinel-1 WV2模式的经验参数:c0=1.511 5,c1=-0.356。

使用Ren模型(即式(20))反演得到的有效波高(Hs empri)与NDBC浮标测量得到的有效波高(Hs ND)对比结果展示在图12(a)中,使用本文方法得到的校正有效波高(Hs correct)与浮标测量结果的对比展示在图12(b)中,共有458景SAR数据与浮标数据满足匹配条件。可以看到,与浮标测量数据相比,两种方法均得到了较好的反演结果,且Hs correct精度高于Hs empri。

((a) Ren模型结果; (b)本文方法结果。(a)The results retrieved by Ren model; (b) The results retrieved by the method in this paper.)图12 反演的有效波高与浮标测量有效波高对比Fig.12 The comparison of the retrieved significant wave height with the buoy measurements

5 结语

本文先通过仿真实验分析了准线性近似法反演海浪参数的准确性和适用范围,仿真结果表明:排除一些实际不存在的一些极端情况,在海面风速较低时,不同波高波长波向的涌浪的海浪谱和有效波高都能得到较好的反演结果,但随着风速增加,海浪谱和有效波高的反演结果都不再可靠。这是由于在高风速下,经过准线性近似的映射关系与非线性映射关系产生了较大的偏差,以及方位向截断波长过长导致该方法反演的海浪谱大面积缺失造成的。然后对2020年的Sentinel-1卫星WV2数据进行反演,在反演过程中,利用交叉谱对相干斑噪声和图像谱中频繁出现的低波数域噪声进行了抑制,反演结果表明,在低风速下,反演的有效波高与ECMWF验证数据吻合较好,使用经验倾斜调制函数比使用理论倾斜调制函数反演结果与验证数据具有更高的相关性,而在高风速时,反演的有效波高逐渐偏低,这一现象与仿真结果相似。截断波长与有效波高反演误差也有一定的关联,但其相关性弱于风速,因此,将风速和截断波长同时用于误差校正将有更好的效果。本文通过1月数据拟合了有效波高反演误差与风速和截断波长之间的关系,建立了校正模型,并将4、7、10月数据的反演结果与ECMWF数据经行了比较,结果表明,在校正之后,有效波高反演精度得到明显的提升。最后将反演结果与NDBC浮标测量数据经行了对比,对于满足匹配条件的458个样本,反演有效波高的BIAS为0.075 m,RMSE为0.592 m,COR为0.878,统计参数优于直接使用截断波长拟合经验关系得到的结果。