基于資源三號光學影像的海洋表面狀態分析

袁策，田根，熊祖強

(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中相應區域統計的有效波高、波齡等參數相符的海況一致，節約獲取海表狀態基礎數據的成本，加深光學衛星在海浪監測中的應用。

但是，兩個及以上波向構成的雙(多)波系統海況中，所蘊含的每一單波系統對海表狀態的貢獻能力及關系是構建海表變化狀態波參數模型的基礎，能降低模型的構建難度，有待深入研究。