基于Sentinel-1A的长江口近岸风矢量场反演研究

戚纤云, 周云轩, 田 波, 于 鹏

(华东师范大学河口海岸学国家重点实验室,上海 200062)

基于Sentinel-1A的长江口近岸风矢量场反演研究

戚纤云, 周云轩, 田 波, 于 鹏

(华东师范大学河口海岸学国家重点实验室,上海 200062)

选取长江口外近岸海域,以欧空局发射的Sentinel-1A为SAR图像数据源,利用多级小波变换和傅里叶变换获得高精度风向信息,再利用基于GMF的3种CMOD系列模型反演风速,最后将反演结果与同时相的ECMWF模式风场数据以及五组实测数据进行对比.结果表明,除个别包含复杂纹理特征的SAR子图像以外,SAR风场反演精度整体优于ECMWF模式风场.经过多级小波变换后的风向反演结果更优,均方根误差分别为31.6°,29.7°,23.5°.经过二级小波变换之后的整景SAR图像结合风向信息代入到CMOD系列模型中反演得到的风速精度最优,均方根误差控制在0.8 m/s.比较适合长江口外近海海域的是MOD-IFR2和CMOD4,均方根误差分别为1.08 m/s和1.05 m/s.

合成孔径雷达SAR; 风场反演; 小波变换; 傅里叶变换;CMOD系列模型

0 引 言

作为海洋环境的关键要素之一,海表风场数据的获取方法有很多.相比于浮标零星的数据覆盖度,散射计较低的分辨率,合成孔径雷达(Synthetic Aperture Radar,SAR)因其分辨率高、穿透能力强、重访周期短,具有十分广阔的前景.近年来,利用SAR影像对海气交互作用、海表流、内波、油迹监测和浅海水下地形反演等方面的研究快速发展[1].

利用SAR反演风场信息包括风向提取和风速反演两部分.SAR影像上的风致条纹最初是由Gerling等[2]发现的,之后Fetterer[3]和Muller[4]以及Wackerman等[5]进一步验证了风致条纹的风向提取能力.利用风致条纹提取风向主要有三种方法:①是Gerling等提出的根据傅里叶变换(Fast Fourier transform,FTT)分析频谱域的波谱能量峰提取风向[5];②是Horstmann等提出的利用局部梯度(Local Gradient,LG)的方法分析空间域上的梯度风[6-7];③是利用小波变换或与傅里叶变换相结合提取风向[8-11].利用SAR风速反演这种方法主要基于两种方式:一种是结合散射计实测数据的半经验模型(EP,Empirical Model),另一种是基于电磁计算,即建立电磁模型(EM,Electromagnetic Model)[12].其中,基于地球物理模式函数(Geophysical Model Function,GMF)的CMOD系列模型包括CMOD4、CMODIFR2、CMOD5、CMOD5.n、CMOD7等.各种研究表明CMOD系列模型获取中性风场信息是最简便、使用率最高的方法[13].国内研究也集中在应用CMOD系列模型并作相关改进的研究[14-15].随着SAR卫星技术的飞速发展,针对不同的研究区域和不同类型的SAR数据源,选择合适的反演方法,对于提高SAR反演海表风场的精度十分重要.同时,相应的精度验证和适用性研究也是非常必要的.在这些研究中,大部分研究者对于反演精度验证,采用的是将模式风场作为验证数据,但是两种风场数据的分辨率并不相同,需要将反演得到的风场插值到模式风场中,并且作为模式风场资料同化来源的散射计数据受陆地影响较大,并不能很好地反映出近岸海表风场的真实情况.所以,为了提高验证的合理性和准确性,本文结合出海采集到的实测点数据与模式风场数据共同作为验证数据进行分析.

1 研究区及数据

2014年,欧空局发射了Sentinel-1(C-band)系列雷达卫星,该系列包含Sentinel-1A和Sentinel-1B两颗卫星,接替和改进已经退役的ENVISAT-ASAR[16].Sentinel-1采用太阳同步轨道,轨道高度693 km,倾角98.18°;2016年Sentinel-1B发射后,Sentinel-1的重访周期变为6 d[17].本文用到的SAR数据类型为地距影像(Ground Range Detected,GRD)、干涉测量宽幅模式(IW)、幅宽为240 km、单视分辨率为5 m×20 m(距离向×方位向)、VV极化方式、轨道方向上行,图像成像时间是2016年6月8日09:54.两种验证数据分别是实测风场数据和欧洲中期天气预报中心(European Centre for Medium-Range Weather Forecasts,ECMWF)提供的模式风场数据.其中实测数据是Sentinel-1A过境时刻的佘山气象站实测数据以及出海采集的5个实测点风场数据(使用WindLog风向风速仪测量),如表1所示.覆盖了5个实测点的试验区域如图1所示,作为后续反演精度的验证区域.ECMWF模式风场是取2016年6月8日9:00和12:00两组数据后插值得到的同时相的距离海表10 m高的风场数据,分辨率为 0.125°×0.125°.

表1 实测数据Tab.1 In-situ data

图1 试验数据地理位置示意图Fig.1 Geographic locations of SAR images and in-situ data

2 方法及预处理

2.1 方法概述

本文首先对整景SAR图像进行一系列预处理,得到的后向散射系数波段经过多级小波变换后,提取各级包含主要信息的低频分量,将其进行二维傅里叶变换即利用风致条纹提取风向信息,再将风向信息以及相关的雷达参数一并代入到3种CMOD模型得到风速.最后将反演得到的风场结果与实测数据和ECMWF模式风场数据进行精度验证分析,流程图如图2所示.

2.2 图像预处理

ESA提供下载的Sentinel-1A图像是level1数据,需要进行的预处理包括辐射定标、重采样、滤波、地形校正、陆地掩膜和船只光斑滤除.其中的σ0波段是图像上的DN值,即遥感图像上的像元亮度值,需要根据辐射标定表(LUTs)中所对应的绝对校准常数进行辐射定标,转化为归一化雷达后向散射截面σ0.再将斜距分辨率约为8.5 m×9.9 m的SAR原始数据重采样为10 m×10 m.将SAR图像固有的椒盐斑点噪声通过滤波处理消除,采用了3×3像素的Lee滤波窗口对SAR图像上的高频噪点进行降噪处理.除此之外,SAR图像需要地形校正补偿因地形变化和卫星传感器的倾斜产生的扭曲.该景SAR图像覆盖了陆地和海面,因此需要进行陆地掩膜,消除陆地以及人造地物等对海表风场反演的干扰.长江口口外分布着众多从航道驶出的船只,因其不同的介电常数和特有的几何形状,在SAR图像上表现为光斑,本研究中根据其明显的后向散射系数特征进行船只光斑滤除.经过一系列预处理的SAR图像见图3,可以看到陆地已经被掩膜,仅剩部分近期促淤围垦的滩涂,因粗糙度远远高于海面,所以呈现亮度较高的状态,比较容易分辨.其次,长江口三级分汊、四口入海处的交界处海面线性条纹特征明显.除此之外,图像上仍包含一些光斑噪点,也是后续操作中需要去除的高频信息.

图2 Sentinel-1A图像反演海表风矢量场流程图Fig.2 Flowchart of sea surface wind vector field inversion based on Sentinel-1A

图3 预处理后的2016年6月7日长江口近岸Sentinel-1A图像Fig.3 Preprocessed Sentinel-1A image located in Yangtze Estuary offshore on June 7,2016

3 SAR海表风场反演

3.1 SAR海表风向提取

3.1.1 SAR提取海表风向物理机制

散射计具有多个天线,以不同入射角观测地面和海表,可以直接获取风矢量信息.单一天线的SAR只能间接获得风场信息[18].因此建立后向散射系数与风矢量场之间的关系成为了该研究方向的重点.研究表明SAR后向散射系数与海面粗糙度有关,这是由于海洋大气导致海表辐聚和辐散,从而引起粗糙度变化,显示出交错的明暗条纹,间接反映出风向信息[19].傅里叶变换和局部梯度法都是基于上述原理提取风向的.

3.1.2 小波变换

利用小波变换将SAR图像进行多级分解后提取具有主要信息的低频分量,选取合适的低频分量再进行傅里叶变换得到风向.经过信号重构后的风向信息再代入到CMOD模型中,缩小风向的反演误差,进而提高风速反演精度.

设f(t)、φ(t)是平方可积函数,L2(R)代表平方可积函数空间,且f(t)、φ(t)∈L2(R),则信号f(t)的一维离散小波变换为[20]

公式中,〈f(t),φj,k(t)〉是 f(t)与 φj,k(t)的内积.如果,{φj,k(t)}j,k∈Z是 L2(R)的标准正交基,则公式为正交离散小波变换,公式变化为

其中,djk〈f(t),φj,k(t)〉是第j尺度下k位置上的离散小波变换系数.如图4所示,经过一次分解后得到的是一个包含主要信息的低频分量和沿水平、垂直以及对角线方向的3个高频.如果去掉部分高频相当于原信号的去噪和压缩处理,对得到的低频分量再次进行连续分解,便可得到信号不同分辨率下的低频分量,即多分辨率分析.本文采用Daubechies小波中的db2.通过大量的试验比较表明,Daubechies小波在SAR海面风向的反演中取得了较优的精度[21].经过三级小波变换分解后的低频分量和高频分量结果如图5.

图4 二维离散小波变换阶层式架构示意图Fig.4 The 2D-discrete wavelet transform structure diagram

3.1.3 傅里叶变换

式中,X为定标后图像的灰度值,即NRCS(dB),Y为计算得到的低波数谱,并且l,m=1,2,···,N[22].二维傅里叶变换后的SAR图像由时域方式变为以频谱域的方式展现出图像包含的线性条纹特征,也就是将SAR图像的频谱按照不同的频率叠加到一幅图中,根据研究目标设置阈值,根据不同的频谱对应的特定频率去除掉噪音,选取低波数部分得到风向信息.但是使用这种方法得到的海面风向有180°方向模糊,需要外部信息确定海面风向[23].经过小波变换和傅里叶变换后,试验区域内SAR提取的风向信息和实测点风向数据以及ECMWF模式风场中的风向信息如表2所示.为了更好地进行精度对比,将前者的风场插值到0.25°×0.25°,与ECMWF风场数据保持一致.

图5 3级小波变化分解后的低频分量和高频分量Fig.5 Multiscale components of wavelet transform decomposition

表 2 试验区域SAR提取风向信息与验证数据结果Tab.2 S1-S5Wind direction data from SAR and validation data 风向/(°)

3.2 SAR海表风速反演

3.2.1 SAR反演海表风速物理机制

讨论SAR反演海表风速物理机制实际上是讨论海面对雷达入射波束的后向散射机制.在20°—70°入射角条件下,布拉格散射占主导地位,雷达后向散射的主要散射体是由海面风场产生的微尺度波[24].

3.2.2CMOD系列模型

地球物理模式函数是基于ECMWF模式风场和ERS散射计数据以及浮标数据等众多风场数据,建立的雷达后向散射系数与海面长波参数、雷达观测入射角、极化方式、雷达入射波的频率等参数之间的关系,继而间接建立了海面风矢量与后向散射系数之间的关系[25].通常地球物理模式函数中σ0可以表示为傅里叶级数展开的形式[26]:

其中A0,A1,A2是入射角、风速、相对方位角的函数,可以通过散射计的测量数据与实测数据经曲线拟合获得.有研究表明[27],在小于15 m/s的风速条件下CMOD4效果较优,风速大于15 m/s时,CMOD-IF2比CMOD4效果更好;CMOD4模型反演近岸风速性能最好,CMODIFR2模型次之,CMOD5模型最差.本文将整景SAR图像分别代入3个CMOD系列模型中,验证上述有关CMOD系列模型的适应范围,风矢量场反演结果如图6和表3所示.这里值得注意的是,A1、A2、A3代表的是分别经过一级小波变换、二级小波变换、三级小波变换后的包含主要信息的低频分量.基于CMOD系列模型的整景SAR反演结果和同时相ECMWF模式风场如图7所示.

图6 基于CMOD-ifr2/4/5模型反演结果与验证数据风矢量对比图Fig.6 The inversion results of wind vector based on CMOD-ifr2/4/5 models compared with the validation data

表 3 实验区域SAR反演风向信息与验证数据结果Tab.3 S1-S5 Wind Speed data from SAR and validation data

图7 长江口近岸风场反演结果图Fig.7 Retrieved sea surface wind field of fthe Yangtze Estuary

4 结果与分析

对于SAR图像风向提取,如表2所示,ECMWF模式风场的风向和SAR反演风向分别与实测点风向的平均绝对误差为 28.9°、29.2°、24.6°、24.6°、15.7°,均方根误差分别为35.1°、41.5°、31.6°、29.7°、23.5°.5个实验区域中 S1、S3、S4、S5的 4个实测数据与SAR图像反演得到的风向基本吻合,而在实验区域S2处,实测风向与SAR反演风向相差较大.这是因为在S2处的子图像上不仅包含风致条纹还有不同方向上尺度相近的线性条纹,这导致了频谱域内无法区分真正的风致条纹,进而得到了不准确的频谱峰值拟合线.因此,反演得到的风向是不准确的.除去试验区域S2,其他验证区域内,根据平均绝对误差和均方根误差可知,经过三级小波变换(A3)处理后的风向反演精度最高,其次是二级小波变换(A2)和一级小波(A1).经过小波变换后的低频信息有效地滤掉了包含不同方向的高频信息,风向提取的精度明显提高,可以看出船舶光斑等奇异点造成的误差是较大的.无小波变换处理的风向是将图像的后向散射波段直接进行傅里叶变换得到风向,其提取精度误差较大,原因可能是SAR图像覆盖区域靠近陆地.海流导致的条纹,以及岛屿、船只等光斑使得SAR图像上的线性条纹特征复杂,容易与风致条纹混淆,导致风向提取精度较差,进一步影响后续风速反演.

对于SAR图像风速反演,整体来说,经过小波变换之后的A1、A2、A3均优于没有进行小波变换的原始数据,这其中很大一部分的原因是风向提取精度的提高.因为有证据表明风向的准确性对利用地球物理模式函数反演风速至关重要,也就是说很小的误差也可能会给风速的反演结果带来明显偏差.多级小波变换降低噪声的效果是十分明显的,但是在小波变换消除噪点的同时,降低了每景子图象的分辨率,这会间接影响风速反演.研究表明,经过二级小波变换之后的整景SAR图像的相应参数结合风向信息代入到CMOD模型中反演得到的风速精度是比较合理的,平均绝对误差为0.72 m/s,均方根误差为0.80 m/s.CMOD系列模型中,比较适合的是CMOD-IFR2和CMOD4,原因是整景SAR图像覆盖的长江口近海海域海表风速较低,在10m/s风速以内;经过CMOD5模型反演得到的风速存在高估的情况,适用于高风速海表风速反演.

在实验区域S5处,实测数据来自于佘山岛气象站数据.佘山岛气象站海拔高度为70 m,需要将气象站的数据根据高度与风速的关系调整到海表10 m高,以便精度验证[19].研究表明,S5处的ECMWF模式风场数据与实测数据最相近,但是整体来讲,其风向信息与实测点风向之间的平均绝对误差和均方根误差均小于仅经过傅里叶变换的SAR图像提取的风向,但是大于经过各级小波变换之后的SAR图像提取的风向.风速反演的精度对比中,ECMWF模式风场的风速与实测点风速之间的平均绝对误差和均方根误差均大于SAR图像反演得到的风速.

从图6中可以看到,每个实验区域都有丁坝、船只光斑、油迹、岛屿(佘山岛)等明显的噪声干扰,但是在这些影响因素下,基于CMOD系列模型利用Sentinel-1A反演长江口近岸风场是可行的.图7中,整景SAR图像风场反演结果图与利用率较高的ECMWF模式风场对比可知,前者的风向在分辨率高的情况下略显得凌乱,后者相邻子图像之间的风向延续性更强.这表明傅里叶变换得到的频谱图像更能捕捉到分辨率较高的SAR图像的风致条纹,从而获得分辨率更高的风向信息.但是,这种方法的弊端是反演的风向信息结果中会存在个别与真实风向相差很大的情况,而出现这种情况往往是因为受到其他因素影响图像上出现的其他线性特征(大气边界层涡旋、Langmuir环流等因素都有可能引起海面线性条纹),并非是风纹理特征.并且这种线性特征的尺度正好在风纹理尺度范围内(例如海表波浪),因此直接干扰SAR图像提取风向,且这种干扰难以用同以频率域为基础的小波变换滤掉.所以,在后续研究中应注意剔除非局地风引起的纹理特征,以获得正确的SAR反演风向.

5 结论与展望

微波遥感图像区别于光学图像,难以根据肉眼直观地目视判读,只能借助频谱域的差异,也就是条纹的尺度去判读,进而进行滤波降噪,根据风致条纹的尺度设置阈值来获取低频的风向信息.SAR图像缺乏风致条纹,或者某些线性条纹并不是由风引起时,在频率相似的情况下,很难分辨条纹的属性.研究表明,适宜长江口近岸风场反演的最优流程是经过二级小波变换后进行傅里叶变换,提取低频风向信息,利用基于GMF的CMOD-IFR2和CMOD4模型反演风速;SAR图像风场反演总体上优于ECMWF模式风场,但是当个别子图像的条纹特征复杂时,容易造成比较大的误差.

未来,针对SAR反演海表风场这一研究方向的改进有以下3个方面.①虽然SAR全球覆盖周期太长,但可以与其他散射计作为数据补充.②建立针对SAR开发的GMF;由于SAR与散射计成像尺度不同,会导致反演过程中无法避免的误差.③因数据获取的限制,仅比较了IW模式的Sentinel-1A图像风场反演精度,未来研究可以尝试不同模式下的Sentinel-1图像反演风场操作.中国已经开启微波遥感时代,随着中国自主研发卫星的技术越来越先进,覆盖我国陆地、海洋的微波遥感卫星会越来越多,模式也会越来越齐全,同时也会进一步提高我们对微波遥感在风场反演领域的认知.

[1] JOHANNESSEN J A.Coastal observing systems:The role of synthetic aperture radar[J].Johns Hopkins Apl Technical Digest,2000,21(1):41-48.

[2] GERLING T W.Structure of the surface wind field from the Seasat SAR[J].Journal of Geophysical Research,1986:2308-2320.

[3]FETTERER F M,GINERIS D,WACKERMAN C C,et al.Validating a scatterometer wind algorithm for ERS-1 SAR[J].IEEE Transactions on Geoscience and Remote Sensing,1998,36(2):479-492.

[4] MULLER G,BRUMMER B,ALPERS W,et al.Roll convection within an arctic cold-air outbreak:Interpretation of in situ aircraft measurements and spaceborne SAR imagery by a three-dimensional atmospheric model[J].Monthly Weather Review,1999,127(3):363-380.

[5] WACKERMAN C C.Estimating wind vectors from RADARSAT synthetic aperture radar imagery[R].Verdian ERIM International Report 10032100-1T,prepared for Office of Naval Research for Contract,2000(00014).

[6]HORSTMANN J,KOCH W,LEHNER S,et al.Ocean winds from RADARSAT-1 scan SAR[J].Canadian Journal of Remote Sensing,2002,28(3):524-533.

[7] HORSTMANN J,KOCH W,LEHNER S,et al.Ocean wind fields retrieved from the advanced synthetic aperture radar aboard ENVISAT[J].Ocean Dynamics,2004,54(6):570-576.

[8] FICHAUX N,RANCHIN T.Combined extraction of high spatial resolution wind speed and wind direction from SAR images:A new approach using wavelet transform[J].Canadian Journal of Remote Sensing,2011,28(3):510-516.

[9] DU Y,VACHON P W,WOLFE J,et al.Wind direction estimation from SAR images of the ocean using wavelet analysis[J].Canadian Journal of Remote Sensing,2002,28(3):498-509.

[10] LEITE G C,USHIZIMA D,MEDEIROS F N,et al.Wavelet analysis for wind fields estimation[J].Sensors,2010,10(6):5994-6016.

[11] 张毅,蒋兴伟,宋清涛,等.C波段地球物理模式函数对比分析研究[J].遥感信息,2012(1):31-36.

[12] LA T V,KHENCHAF A,COMBLET F,et al.Study of wind speed retrievals from Sentinel-1 images using physical models[C]//Antennas and Propagation(EuCAP),201610th European Conference on.IEEE,2016:1-4.

[13] 孔毅,赵现斌,艾未华,等.基于墨西哥帽小波变换的机载SAR海面风场反演[J].解放军理工大学学报(自然科学版),2011(3):301-306.

[14] HE Y,ZOU Q.Validation of wind vector retrieval from ENVISAT ASAR images[C]//Geoscience and Remote Sensing Symposium,2004,IGARSS’04.IEEE,2004:3184-3187 vol.5.

[15] 宋贵霆,侯一筠,何宜军.Comparison of two wind algorithms of ENVISAT ASAR at high wind[J].中国海洋湖沼学报(英文版),2006,24(1):92-96.

[16] MONALDO F,JACKSON C,PICHEL W,et al.Early validation of operational SAR wind retrievals from Sentinel-1A[C]//Geoscience and Remote Sensing Symposium.IEEE,2015:1223-1226.

[17]PICHEL W G,MONALDO F M,JACKSON C,et al.NOAA operational SAR winds:Current status and plans for Sentinel-1A[C]//Geoscience and Remote Sensing Symposium.IEEE,2015:4916-4919.

[18] 范开国,黄韦艮,常俊芳,等.NCEP/QSCAT混合风向用于SAR图像反演高分辨率海面风速[J].遥感技术与应用,2010,25(6):873-876.

[19]ALPERS W,BR¨UMMER B.Atmospheric boundary layer rolls observed by the synthetic aperture radar aboard the ERS-1 satellite[J].Journal of Geophysical Research Oceans,1994,99(C6):12613-12621.

[20] 张毅,蒋兴伟,林明森,等.基于小波分析的近岸海面风场反演研究[J].高技术通讯,2011,21(10):1056-1061.

[21] 张雷,石汉青,杜华栋,等.用小波分析反演星载SAR图像海面风向[J].遥感学报,2013,18(1):215-230.

[22] 杨劲松,黄韦艮,周长宝,等.合成孔径雷达图像的近岸海面风场反演[J].遥感学报,2001,5(1):13-16.

[23] 程玉鑫,艾未华,孔毅,等.基于影像纹理特征和外部风向的星载SAR海面风场反演研究[J].海洋科学,2015,39(12):157-164.

[24] LIN H,XU Q,ZHENG Q.An overview on SAR measurements of sea surface wind[J].Progress in Natural Science,2008,18(8):913-919.

[25]WACKERMAN C C,RUFENACH C L,SHUCHMAN R A,et al.Wind vector retrieval using ERS-1 synthetic aperture radar imagery[J].IEEE Transactions on Geoscience and Remote Sensing,1996,34(6):1343-1352.

[26]PORTABELLA M,STOFFELEN A,JOHANNESSEN J A.Toward an optimal inversion method for synthetic aperture radar wind retrieval[J].Journal of Geophysical Research:Oceans,2002,107(C8):1-13.

[27] 王珂,洪峻,张问一,等.SAR反演邻近岸海面风场方法[J].电子与信息学报,2013,35(8):1800-1805.

(责任编辑：李万会)

Sea surface wind vector retrieval of fthe Yangtze Estuary based on Sentinel-1A

QI Xian-yun,ZHOU Yun-xuan,TIAN Bo,YU Peng
(State Key Laboratory of Estuarine and Coastal Research,East China Normal University,Shanghai200062,China)

This paper selected SAR images in the Yangtze Estuary from Sentinel-1A which was launched by ESA.Firstly,it extracts high-precise wind direction using multilevel wavelet transform and fourier transform.Secondly,wind speed information are obtained through three kinds of CMOD models of GMF,and compare with five sets of in-situ data and as well as wind vector from ECMWF.The results indicats that except for some sub-images containing complex texture,the precision of SAR-derived wind field is better than that of ECMWF.The multiscale wavelet transform is conducive to wind field retrieval(RMSE of wind direction are 31.6°,29.7°,23.5°respectively),especially the wind speed retrieval after the second level wavelet(RMSE of wind speed is 0.8 m/s).Among CMOD models,CMOD-IFR2 and CMOD4 are suitable for coastal area of the Yangtze Estuary(RMSE of wind speedis 1.08 m/s and 1.05 m/s respectively).

synthetic aperture radar; wind vector retrieval; wavelet transform;fourier transform; CMOD models

TP732

A

10.3969/j.issn.1000-5641.2017.06.012

1000-5641(2017)06-0126-10

2016-12-08

国家自然科学基金(41476151)

戚纤云,女,硕士研究生,研究方向为河口海岸微波遥感.E-mail:qxybeijing0319@hotmail.com.

周云轩,男,教授,研究方向为河口海岸变化遥感研究.E-mail:xyzhou@sklec.ecnu.edu.cn.