袁策，田根，熊祖强

(1.河南理工大学 资源环境学院，河南 焦作 454003；2.河南理工大学 能源科学与工程学院，河南 焦作 454003)

0 引言

全球性、区域性和局部性的海表变化状态的在线数据的实时获取是现代海洋学亟待解决的问题[1]。海表变化状态的光谱信息对海-气交界层、海水污染、人类活动对海表状态影响的监测等有重要的作用，获取它们是现代海洋学基础研究和应用的重要工作[2]。

空间谱、频率谱和空间-频率谱是海表状态光谱信息表达的重要手段，既能描述海洋近表层现象的变化，又能表达海表状态的能量分布[3]。Cox等[4]分析了从高分辨率航拍影像中提取的一维或二维海浪空间谱的合理性，研究蕴含海浪空间结构信息的海浪亮度场的分布规律，指出从光学影像中提取海浪谱的可行性。随后，国内外学者针对海浪谱展开研究，有针对海表状态分布特征模拟、海浪谱重构和海表状态演化，有偏重海浪亮度场与海表状态的数学函数的构建。如Yurovskaya等[5]提出立体光学影像中海浪亮度场与高度场数学函数的半经验模型，根据哨兵-2获取影像的时间差，给出海表状态与海浪亮度场的数学模型，分析了近岸区海水折射对海表状态和海浪能量谱中局部能量的影响[6-7]。Maria 等[8]深入研究海浪亮度场与高度场的数学关系，指出海浪谱能量分布，与波峰点位置的变化相关，也与风海变化规律一致。这些研究偏重于光学影像中海表状态模型或海浪谱特征，未考虑海浪成因对海表状态的影响。

Boukhanovsky等[9]根据海浪能量谱中能量集中方向(即主方向)的个数判别海表状态复杂度，但未考虑数据中薄云等噪声引起的海浪能量谱中能量少许集中形成的伪主方向；Bitner-Gregersen等[10]研究了由涌浪和风浪组成的双波海况中交叉谱特征及成因；Park 等[11]从物理建模的角度，研究了由不同类型海浪组成的复杂海况(波数3个及以上的海况)的反演及其海浪谱特征；Luxmoore 等[12]提出由风浪、涌浪组成的极端海况判断方法。这些研究偏重海表状态成分，或已知海浪的类型和各自占比所构成海况的海浪谱模拟等，利用海浪谱研究海表变化特征的研究很少，如Mori等[13]研究了海浪谱的峰度与波向的关系等。

海表变化状态与海浪类型、成分、影响因素等密切相关，但综合分析海表状态的研究较少。本文从资源三号光学影像中海表状态模拟出发，基于光学成像原理，采用数字图像处理技术和统计学方法，通过对频率域中海浪谱特征的分析，提出应用海浪谱的能量谱、峰度和偏度参数，研究海况的复杂度、海表分布状态，结合海浪的空间表象特征，确定海浪类型及成因，实现海表状态分析的方法，为实时获取海表状态数据提供方法支持。

1 数据与方法

1.1 数据

海表状态变化不仅受风力的影响，在近岸区域还受海底地形、海水深度等因素的影响。因此，选择海况、海水深度、数据质量不同的典型近岸浪、涌浪和风浪区的资源三号全色和多光谱光学影像作实验数据(图1至图3，大小300像素×300像素)。海水深度数据来自美国国家海洋和大气管理局的高精度栅格导航海图(http：//maps.ngdc.noaa.gov/viewers/bathymetry)。验证数据是欧洲中期天气预报中心的海浪模式数据，应用ArcGIS软件处理，统计实验区内海浪的有效波高、波龄(表1)。

图1是美国西海岸加利福尼亚附近，海水深度1 000 m以下，分辨率2.1 m的ZY-3光学影像切块，拍摄于2014年5月20日，无云。该区域位于太平洋的东部，常年受太平洋面的风和沿岸峭壁影响，浪高比较大，是典型的风浪区。

图1 风浪实例

表1 基于欧洲中期天气预报中心海浪模式数据的波参数

图2是中国南海七连屿附近，含云量7%的5.8 m分辨率的ZY-3多光谱影像，拍摄于2013年4月26日，图2(a)和图2(b)区域内海水深度分别是400～600 m和50～200 m。七连屿常年受风暴袭击，灾害性天气频发，周围海域海浪比较复杂，常出现涌浪、风浪和涌浪的混合浪。

图2 混合浪实例

图3(a)是中国獐子岛附近，海水深度20～30 m，分辨率2.1 m的ZY-3全色影像切块，拍摄于2012年7月6日。獐子岛周围海水深度平均40 m，最深处约140 m，是典型的近岸浪海域。图3(b)是分辨率5.8 m的ZY-3多光谱影像，拍摄于2016年8月12日，海水深度浅于50 m，位于波罗的海，也是典型的近岸浪区域。

图3 近岸浪实例

表1数据是按月统计，与实验影像的时间误差较大，因此，仅从波龄的变化趋势和有效波高对应海况分析是否与研究的海表状态一致。

1.2 海表状态模拟

海表的起伏是由任意波高(z=ζ(x，y，t))的海浪组成的集合。幅宽52 km的ZY-3卫星光学影像拍摄时长7 s。实验数据中，分辨率2.1 m的全色影像实地距离630 m，分辨率5.8 m的多光谱影像实地距离1 740 m，实验区拍摄时长最长不足0.3 s，因此，忽略实验数据的拍摄时长，影像是某时刻点(t=t0)的快照数据，海表的起伏如式(1)所示。

z=ξ(x，y)=ζ(x，y，t)|t=t0

(1)

式中：(x，y，z)是笛卡尔坐标系下，海表任一点的三维坐标；(x，y)是静止海平面上任一点坐标；ξ(x，y)是海浪高度场模拟函数。

模拟的海表状态在频率域中规律更显著，因此，应用快速傅里叶变换(fastfourier transform,FFT)，将海表模拟模型变换到频率域，获取海浪的梯度，即海浪谱，在频率-波数坐标系中的表示如式(2)所示。

(2)

式中：(kx，ky)是波矢量k分别在x，y向的波数；s、l是实验影像的第s行和第l列；m、n是实验影像的总行数和总列数；e是自然对数的底数；i是虚数。

(3)

影像中海浪光谱亮度是海表接受光照辐射能量强弱的表征，与海浪坡度直接相关[15]，海浪坡度是在一定方向的光照下，由海浪高低起伏引起，因此，由海浪坡度ζ(x，y)与海浪亮度的关系，基于光学成像原理，可构建海浪高度场与亮度场的数学函数。

沿海浪波向φ(φ=arctan(kx，ky))的海浪坡度如式(4)所示。

ζφ(x，y)=cosφξx(x，y)+sinφξy(x，y)

(4)

式中：ξx(x，y)和ξy(x，y)是海浪高度场模拟函数在x、y轴的偏度。

在傅里叶空间，沿方向φ的海表坡度ζφ(k)表达如式(5)所示。

(5)

(6)

T(φ)=B×G×cos(φ-φG)

(7)

由海浪高度与坡度的关系(式(5))和坡度与海浪亮度场的函数(式(6))，构建海浪光谱场与高度场间的数学关系，如式(8)所示。

(8)

(9)

式中：C是常数，在海表状态分析过程中不考虑。

海浪能量谱ES是傅里叶空间海浪谱模的平方，如式(10)所示。

ES=|fft(I(s，l))|2

(10)

式中：fft(·)是交互式数据语言中快速傅里叶变换函数；I(s，l)是光学影像的亮度场。

由式(9)和式(10)可知，谱密度函数与海浪能量谱成正比，如式(11)所示。

(11)

所以

ES=|T(φ)(cosφkx+sinφky)|2

(12)

将式(7)和式(12)，代入式(8)得式(13)。

(13)

(14)

式(12)和式(14)显示海表状态与蕴含海浪波向φ、波数等分布特征的海浪谱能量相关，与传感器图像构成无关[16]，因此海表状态不受光学影像定标结果的影响；另一方面，越高分辨率影像记录海浪细节越详细，米级以上的影像的分辨率差异小，不影响海浪状态。

当海表状态遵循瑞利分布时，呈线性分布，用高斯函数模拟，海浪谱的峰度极值是3，偏度极值是0[17]；相反，海浪谱的峰度极值不等于3，且偏度极值不等于0时，呈非线性分布，需高斯函数、峰度和偏度曲线组合模拟[18]。因此，海浪谱的峰度和偏度的极值，是判断海浪分布状态的依据。

海浪谱峰度(kurtosis)和偏度(skewness)的计算如式(15)至式(16)所示。

kurtosis=fft(I(s，l))4

(15)

skewness=fft(I(s，l))3

(16)

1.3 海表状态分析过程

根据海表状态模型及参数计算方法，对资源三号光学影像进行裁剪和快速傅里叶变换，获取实验数据的海浪谱；计算海浪谱的能量谱、峰度和偏度，绘制它们统计图，统计能量谱极值转折点数、峰度和偏度的极值；根据能量谱主方向数及其统计曲线上极值转折点数，判断海况的复杂度；由海浪谱的峰度和偏度的极值判别海表分布，结合海浪空间特征、海水深度，确定海况类型及成因，实现海表状态分析，具体过程如图4(j>2的整数)所示。

图4 海状态分析

在海表状态分析前，首先，对数据进行防溢处理，保证影像切块的尺寸是2的p(p>0的整数)次方倍；其次，快速傅里叶变换的开窗尺寸影响海浪谱分布的表达，需实验确定；同时，为了更好表述海浪谱能量的分布规律，对其进行均值滤波，排除高频波的干扰，使谱能规律更显著，统计峰度、偏度，绘制海浪能量谱曲线图时，用未滤波数据，以全面描述海浪谱特征。

2 结果与分析

2.1 实验结果

实验数据图1～图3影像的海浪能量谱(图5)、峰度和偏度(图6)，海浪能量谱特征和海浪谱的峰度、偏度值及特征是海表状态判断的依据。

图5 海浪能量谱及统计图

图6 海浪谱的峰度和偏度

2.2 实验分析

1)海况复杂度。海表状态是由波向不同、波高各异和类型不一的海浪构成的海波系统。当海表仅由单一波向的风浪或涌浪构成，是单波系统海况，由两个及以上波向不同的海浪构成，是双(多)波系统海况。因此，海浪波向数是度量海波复杂度的参数之一。海浪能量谱中蕴含的海浪的波向信息与其主方向一致，所以主方向数即波向数，是判断海波系统复杂度的依据。但影像中的薄云等噪声使海浪能量谱含伪主方向，误判海波数，因此，本文提出由海浪能量谱的主方向数结合其曲线极值转折点数，确定海表海波系统的复杂度。

能量谱图中能量集中方向为一个，且其统计曲线上极值转折点也只为一个，那么该海况是单波系统，如图5(a1)、图5(b1)、图5(d1)所示。但图5(b1)中在约30°方向出现能量少许集中，是因为图2(a)影像含7%薄云造成海浪能量谱中出现伪主方向。

能量谱图中主方向数多于一个，同时其统计曲线上极值转折点数也为多个，则该海况是双或多波系统，如图5(c1)和图5(e1)所示。图5(c1)中在约120°和135°两个方向能量聚集明显，约45°方向能量少许集中，而图5(c2)中，仅在影像灰度值为60和100附近出现极值转折点，所以该海况是双波系统，原因是数据中含7%的云量所致。图5(e1)中海浪谱能量集中在约90°、45°和135° 3个方向，图5(e2)中统计曲线呈近似正弦曲线，分别在影像灰度值为50、75和100 3处有极值转折点，因此是复杂海况。

实验结果显示，光学影像中少量的云量引起海浪能量谱中能量轻度集中，而能量谱曲线上未出现极值点，因此，如果能量谱中主方向数与极值转折点数相同，由主方向数判断海况复杂度；相反，则根据极值转折点数，按能量集中程度由大到小选择极值转折点数个主方向作为海浪波向，从而有效排除薄云的干扰，较准确实现海况复杂度的判断，为进一步确定海浪谱类型、海表状态分析提供理论基础。

2)海表状态。海表状态的分布由海浪谱的峰度和偏度的极值判断，结合海况复杂度，海浪空间特征，分析海表状态所蕴含的海浪类型，兼顾海水深度数据，分析海浪类型的成因，比较与表1中相应区域内有效波高和波龄适宜的海况是否一致。

图6(a1)和图6(a2)中海浪谱的峰度和偏度极值都小于0.5，显示海表状态呈弱非线性分布；图5(a1)中该区域是单波系统海况，结合图1中海浪峰(谷)条带较短、分布无规律，且海水深于1 000 m，地形对海表状态无影响，因此，该海表状态是单波风浪，与统计表1中该区域海浪的有效波高1.71 m、波龄32是年轻海浪所对应的海况一致[19]。

图6(b1)、图6(b2)和图6(c1)、图6(c2)是同景不同区域影像海浪谱的峰度和偏度。图6(b1)和图6(c1)的峰度极值分别是13和15，图6(b2)和图6(c2)的偏度极值分别是2.69和3，这表明两区域的海表非线性特征突出。由图5(b1)知图6(b1)对应的海况是伴有弱交叉谱的单系统海浪，由图5(c1)知图6(b2)的海况是交叉谱海况。虽两区域影像由于含7%的薄云造成海浪谱的峰度极值变大，但图6(b2)中交叉谱海况的波能相互累积也引起峰度极值增大。此外，300～600 m和50～200 m的水深对海表状态的影响造成深水区有效波高小于浅水区，与表1中两区域对应的有效波高分别是0.68 m和0.70 m相符。其次，图2中海浪条带的条纹长度长于图1中海浪条带长度，且分布规律较明显，符合波龄85成熟海浪的特征，因此该区域近岸是涌浪(图2(a))，随离岸距离增大海表状态转为风浪和涌浪的混合浪(图2(b))。

图6(d1)和图6(d2)中峰度和偏度的极值分别接近4和2，显示该区域海表非线性特征明显，是因为海水深度为20～30 m(图3(a))，海底地形造成海表出现了折射和反射现象。光学影像中，海浪是长条带，分布比较规律，是典型的长涌浪单波海况(图5(d1))，与表1中显示该区域海浪的有效波高0.58 m、波龄180充分发展海浪的海况一致[19]。

图6(e1)和(e2)中峰度和偏度极值分别是-0.5和0.8，表明海表非线性分布，负峰度值显示该区域海表风力小，波浪波高不大，但波陡较大，与有效波高0.45 m相符(表1)。光学影像中(图3(b))，海浪条带长度很短，分布杂乱无规律，因此是由多个波向不同的风浪构成的典型近岸复杂风浪海况，与波龄45(表1)发展中海浪的海况一致。

实验结果表明，海浪谱的峰度和偏度受海浪类型、复杂度和成因影响。单波海况中，风浪的峰度和偏度极值都小于涌浪的峰度和偏度的极值，且对峰度的影响小于对偏度的影响；多波海况中，涌浪和风浪混合海况的峰度和偏度极值大于若干波向不同的风浪组成的混合海况的峰度和偏度极值，且对峰度的影响大于对偏度的影响。

3)影响因素分析。

(1)开窗尺寸。当快速傅里叶变换的窗口尺寸过小时，海浪谱分布规律难以正确表达；当窗口尺寸过大时，难以分析海浪谱能量聚集规律，影响海况复杂度的判断，因此，合适的窗口尺寸，才能完整表达海浪谱分布特征和能量集聚规律。本文从同景影像中，分别裁剪32像素×32像素、64像素×64像素、128像素×128像素、256像素×256像素、512像素×512像素的切块，提取它们的海浪谱的能量谱、峰度和偏度，绘制曲线图。实验结果表明，随开窗尺寸的增大，海浪能量谱、峰度和偏度曲线形状由渐变、相似、突变到周期性，结果显示128像素×128像素的窗口尺寸最优。

(2)海水深度。海表状态在沿岸区，严重地受海底地形、岛屿等影响。选择与图1同景的海水深度为10～300 m、编号1～4的4块影像(图7)，探索海水深度对海表状态的影响。

图7 海水深度影响探索实例

提取编号1～4影像的海浪能量谱、海浪谱的峰度和偏度，绘制它们的统计图曲线，统计海浪谱类型、峰度和偏度的极值(表2)。

表2 同景深度不同影像的海表状态参数

表2显示随海水深度的增加，海表状态的非线性呈现先强后弱，对海浪谱的峰度影响较大，偏度影响较小，但成因各异。

编号1影像的海况是单波风浪，海表的非线性分布是由于海底地形产生的折射造成的；编号2受海底地形影响产生浅化和反射现象，使海表的非线性特征显著，峰度和偏度的极值分别为4.25和1.40；编号3和编号4影像的海表非线性特征强是交叉谱能量叠加引起的，但也受海底地形的影响，编号3的海水深度浅于编号4，因此，对海表状态的影响前者大于后者。

3 结束语

基于光学成像原理，根据海表亮度场与海表状态的数学函数关系，以海浪在空间域和频率域的特征为基础，研究海浪能量谱及海浪谱的峰度和偏度，结合海水深度，实现海表变化状态的分析，为实时获取海表状态数据提供一种新思路。

首先，由海浪能量谱的主方向结合其统计图曲线上极值转折点数，确定海波系统的复杂度，解决数据中含噪声(薄云等)形成的伪主方向引起的海波系统误判；然后，根据海浪谱的峰度和偏度的极值，结合它们的曲线形态，分析海表分布；最后，由海波系统复杂度、海表分布、海浪表象的空间特征及实验区的海水深度，确定海浪类型及成因，实现海表状态分析。分析结果与表1中相应区域统计的有效波高、波龄等参数相符的海况一致，节约获取海表状态基础数据的成本，加深光学卫星在海浪监测中的应用。

但是，两个及以上波向构成的双(多)波系统海况中，所蕴含的每一单波系统对海表状态的贡献能力及关系是构建海表变化状态波参数模型的基础，能降低模型的构建难度，有待深入研究。