基于功率譜獲取海浪波長波向算法的改進

陶 法, 馬舒慶

（成都信息工程學院電子工程學院,四川成都 610225）

1 引言

隨著海洋遙感技術的不斷發展,基于Seasat上的合成孔徑雷達（SAR）獲取海浪遙感圖像和無人機航拍獲取海浪圖像,對海洋和海洋氣象發揮了越來越重要的作用。準確有效地獲取海浪波長和波向,不僅對海浪預報、海洋工程,而且對反演海面風場同樣具有重要意義[1]。文中主要針對利用無人機航拍海浪圖像獲取海浪波長和波向算法進行了研究,無人機航拍海浪圖如圖1所示。

2 海浪模型

根據Louguet-Higgins的理論,將平穩海況下的海浪視為平穩的具有各態歷經性的隨機過程,波動可看作是無限多個振幅不等、頻率不等、初相位不等的簡單余弦波的疊加而成。忽略波與波之間的相互作用,可得到三維隨機海浪波高方程[2]:

其中,m,n分別為頻率分割數和方向分割數,km為波數,Φn為相對于x軸的傳播方向,x,y分別表示平行和垂直于主導波傳播方向的笛卡爾坐標,εmn為隨機相位角,Amn為不同頻率的簡單正弦波的振幅,其表達式為:

其中,S（ωm,Φn）為海浪方向譜,且 S（ω,Φ）=S（ω）G（Φ）,Δωm 、ΔΦn分別為ωm 、Φn的增量,S（ω）為頻率譜函數,即不同頻率間隔內的組成波提供的能量,代表海浪能量相對于組成波頻率的分布,又稱為能量譜,G（Φ）為方向函數。在方向函數G（Φ）一定的條件下海浪的振幅 Amn與頻率譜函數S（ω）成正比關系,根據Neumann于1952年得到的一般表達式為[3]:

圖1 截取達里諾爾湖的無人機航拍海浪圖

其中,指數 p取 4～6,q為2～4,A和B中包含風要素（風速、風時、風距）或波高要素（波高、周期）。

基于Neumann理論Pierson-Moscowitz 1964年提出[3]:

其中 α=8.1×10-3,β=0.74,α g2=0.78。從理論上講,S（ω）分布在 ω=0～ ∞,實際上其顯著部分集中在一定頻帶內,其P-M譜圖如圖2所示。

3 基于重心坐標的改進算法

目前從海浪圖像中提取海浪波長的方法主要有:二維傅里葉變換法（2D-FFT）[4]、一維傅里葉變換法（1D-FFT）[4]、基于聚類和hough變換法[5]、以及基于轉動慣量[6]等方法。

3.1 2D-FFT法

方法對海浪圖像進行二維傅氏變換,尋找代表海浪主頻率的最大點的位置,從而求出海浪傳播方向和波長。

其中L,θ代表海浪的波長和方向。xm,ym是功率譜圖中代表海浪主頻率最大值點的坐標,N代表海浪圖像的采樣數,Δx代表采樣間隔。

2D-FFT法比較直觀,但主要缺點有:由于噪聲的存在和主頻不明顯,極易產生誤差,而用平滑法去除噪聲又會引起主頻位置的偏移;在計算海浪波長和方向時只用到最大值點,大部分求出來的頻率分量都沒有用到。

圖2 P-M譜圖

3.2 基于轉動慣量法

假設海浪功率譜圖是一個有許多離散粒子組成的薄平面,那么波浪的方向就是使這些粒子圍繞某一軸的轉動慣量最小的軸的方向,設該軸為y=kx,通過點（x0,y0）到直線的距離公式可得到:

那么其轉動慣量為:

令轉動慣量最小即dJ/dk=0,可求出:

其中k1×k2=-1。k1,k2分別代表使轉動慣量最大和最小的值,并舍去最大的值。再沿著這條軸的方向在某一領域尋找最大值,該最大值點的坐標就是海浪主波峰的位置,從而可以計算出波浪的波長和波向。

方法在求取方向值時較方便,但求取最大值點時同樣受到噪聲點和主頻不明顯的影響。

3.3 對2D-FFT方法的改進

從圖2中可以看出在中心頻率點處S（ω）最大,即由頻率為中心頻率的組成波提供的能量最大,可認為該頻率為海浪的主頻率。但在實際中海浪不是單一頻率組成而且海浪的方向也不是單一方向的,由公式（1）可知每一個頻率組成波對海浪的波長和波向都起一定的作用,其作用的大小由其系數值決定。所以僅以主頻率計算海浪波長和波向存在誤差。

3.3.1 算法實現

將波浪圖像進行二維傅立葉變換得到頻譜圖并通過頻移使其低頻部分位于圖像中心,（為了計算方便可以保留頻譜中大于某一閾值的值,其他值設為0）。由于頻譜是中心對稱的所以只需選取頻譜圖的一半,如圖3所示。

設像素點為（m×n）的頻譜圖中每一點坐標代表圖像的頻率分量可設為 f（x,y）,其值代表圖像中該頻率分量所占的比重可設為 A（x,y）,則由所有頻率組成波提供的海浪功率譜為:

如果用某一特征頻率點為 f（x0,y0）來代表所有頻率點的話,則由該點計算出海浪功率譜為:

利用合力矩定理可以計算出該點坐標:

圖3 提取前1024個最大值點的圖像頻譜圖

即將圖像的功率譜看成密度不均勻的薄片平面,對其求重心,所得重心的坐標即是所要求的特征頻率點。

3.3.2 仿真結果

選取閾值為32點到1024點得出到中心點的距離和方向值,結果如圖4所示。

從圖中可以看到,隨著取樣點的增多計算出的方向角漸漸變小,距中心點距離漸漸變大,但變化平穩。其原因在于,隨著采樣點的增多圖像的高頻分量會相應地增多,其方向角變化趨勢與海浪紋理細節的變化趨勢相一致;距中心點距離變大和P-M譜中隨著海浪譜的增大其主頻率點向低頻方向漂移相一致。

其中閾值為32點所求的重心坐標為:（x0=5,y0=13）,方向角θg=71°,其空域仿真結果如圖5、6所示。

圖4 基于重心法的距離值Dg和方向值θg

圖5 基于重心點的空域海浪圖

圖6 截取閾值為32點的空域海浪圖

對2D-FFT方法的改進算法的優點為:所求的f（x0,y0）點是所有頻率點的重心擬合。不會造成大部分頻率點用不上的問題;該重心坐標是唯一的,不會產生主頻率點不明顯等問題;能夠根據選取的閾值合理地反映實際海浪的波長和波向。

4 分析與比較

圖7 基于重心法、最大值點的距離值 Δ D和方向值Δ θ

圖8 基于重心法和轉動慣量法的方向值

在對圖1加入不同比重的鹽椒噪聲條件下分別基于重心法、轉動慣量法、最大值法獲取到中心點的距離（Dg、Dm）、方向值（θg、θI、θm）然后作比較,結果如圖 9、10 所示 。

從圖中可以看到:與通過最大值點求取的相應值作比較,結果是在一定范圍內方向增加6°,距離增加1個像素;與基于轉動慣量法所測得的方向值數據基本一致;基于最大值法求取的距離值Dm、方向值θm受加入鹽椒噪聲影響明顯;加入鹽椒噪聲后基于重心法和轉動慣量法獲取的距離值Dg、方向值（θg、θI）基本保持不變。

圖9 加入鹽椒噪聲下重心法與最大值法的距中心點距離值

圖10 加入鹽椒噪聲下重心法、轉動慣量法、最大值法的方向值

5 結束語

利用重心坐標法對基于圖像功率譜法獲取最大值點計算海浪波長、波向算法加以改進,并通過大量的試驗與最大值點和轉動慣量法求取的數據對比,驗證了該方法計算方便、抗噪聲能力強、物理概念清晰,能夠更合理地表示海浪波長和波向,并有效地解決了基于二維傅里葉變換帶來的主頻不明顯和大部分求出來的頻率分量都沒有用到等缺點。

[1] 陳艷玲,黃城,丁曉利.ERS22 SAR反演海洋風矢量的研究[J].地球物理學報,2007,50（6）:1688-1694.

[2] 俞韋修.隨機波浪及其工程應用[M].大連:大連理工大學出版社,2003.

[3] 文圣常,余宙文.海浪理論與計算原理[M].北京:科學出版社,1984.

[4] 孫京生,劉征凱.利用SeasatSAR遙感圖象分析海浪的一種新方法[J].電子科學學刊,1988,10（6）.

[5] 劉政凱,畢亞雷,黃潤恒.HOUGH變換在海浪遙感圖像處理中的應用[J].環境遙感,1989,6（1）.

[6] Tang Huiming,Xu Shengrong.Detection of Direction and Wavelength of Ocean Wave by Power Spectrum of Ocean Wave Image[M].Pattern Recognition,1992:164-166.

[7] Rafael C,Gonzalez.Digital Image Processing[M].Prentice Hall,2002.