李漢東, 陳澤宗,2*, 趙 晨, 陳 曦, 陳小莉

(1.武漢大學(xué)電子信息學(xué)院， 武漢 430072； 2.武漢大學(xué)地球空間信息技術(shù)協(xié)同創(chuàng)新中心， 武漢 430079)

船載微波多普勒雷達(dá)從散射回波中獲取海面動(dòng)態(tài)信息，實(shí)現(xiàn)了走航觀測海浪的新手段[1]。目前對(duì)海面散射回波的研究大多建立在靜態(tài)雷達(dá)平臺(tái)上，而船載雷達(dá)平臺(tái)下天線接收的回波信號(hào)中不僅包含海態(tài)信息，也包含了平臺(tái)運(yùn)動(dòng)信息，改變了海面電磁散射的計(jì)算模型，給海面回波的研究帶來了新的挑戰(zhàn)。

目前海面電磁散射的計(jì)算方法主要分為數(shù)值法和近似法[2]。矩量法等數(shù)值法的計(jì)算精度高，但存在計(jì)算量大的缺點(diǎn)；小斜率近似法[3]是一種廣泛應(yīng)用的高精度近似方法，通過保留表面斜率級(jí)數(shù)的各階展開項(xiàng)獲取散射振幅的近似值，適用于復(fù)合尺度的海面電磁散射問題，具有計(jì)算簡便、適用性廣且精度較高等優(yōu)點(diǎn)。

船體在航行中容易受到海浪和海流等外部環(huán)境影響而產(chǎn)生橫搖、縱搖等六自由度運(yùn)動(dòng)，導(dǎo)致回波多普勒譜頻移和譜寬發(fā)生改變，進(jìn)而影響到海態(tài)參數(shù)的反演。近年來，中外學(xué)者提出了一系列實(shí)時(shí)或非實(shí)時(shí)的船舶運(yùn)動(dòng)預(yù)測算法。Fossen等[4]提出了動(dòng)態(tài)船舶運(yùn)動(dòng)預(yù)測的條帶理論公式；Ueng等[5]提出了實(shí)時(shí)預(yù)測的船舶運(yùn)動(dòng)仿真有效計(jì)算模型；Sandaruwan等[6]提出了高性能船舶姿態(tài)預(yù)測算法。

現(xiàn)基于船載微波多普勒雷達(dá)系統(tǒng)，根據(jù)船載雷達(dá)探測海面時(shí)的工作特點(diǎn)，研究走航狀態(tài)下微波電磁波與海表面的作用機(jī)理。通過二維時(shí)變海面建模與船舶運(yùn)動(dòng)建模，采用小斜率近似法計(jì)算海面電磁散射幅度，仿真得到船載雷達(dá)平臺(tái)下的海面散射回波，并分析不同運(yùn)動(dòng)方式(勻速運(yùn)動(dòng)、變速運(yùn)動(dòng)、搖晃運(yùn)動(dòng))對(duì)多普勒譜特征的影響。根據(jù)空間海浪反演算法計(jì)算得到海浪譜并與理論值對(duì)比分析，驗(yàn)證該海面回波模擬方法的正確性。

1 時(shí)變海面建模

有向海浪譜S(ω,θ)描述了海面組成諧波分量在頻域和方向上的分布，可以表示為無向海浪譜S(ω)與方向分布函數(shù)G(ω,θ)的乘積。采用Pierson-Moscowitz譜和Longuet-Higgins方向函數(shù)的表達(dá)式為

(1)

(2)

S(ω,θ)=S(ω)G(ω,θ)

(3)

式中：S(ω)為無向海浪普；α為無因次常量，α=0.008 1；β=0.74；g為重力加速度；U19.5表示海面上方19.5 m處的平均風(fēng)速；ω為海浪角頻率；G(ω,θ)為方向分布函數(shù)；S(ω,θ)為有向海浪譜；s為角散系數(shù)；θ為方向角；θ0為主效浪方向；Γ(*)為伽瑪函數(shù)。

根據(jù)線性波浪理論，時(shí)變海面的波面位移能夠表示為1組不同振幅、頻率和相位的余弦波的疊加[7]。采用線性疊加法計(jì)算得到波面位移h(x,y,t)為

(4)

式(4)中：(x,y)表示海面離散點(diǎn)坐標(biāo)；t為時(shí)間；M、N分別為頻率和方向維度的離散點(diǎn)數(shù)；ωi和ki分別為角頻率與波數(shù)；θj為方向角；εij為均勻分布在(0, 2π)間的隨機(jī)相位。

2 船舶運(yùn)動(dòng)建模

2.1 船舶運(yùn)動(dòng)模型

當(dāng)船舶在海面上航行時(shí)，建立圖1所示的船體坐標(biāo)系G-XsYsZs。由于受到海洋環(huán)境的影響，船體會(huì)沿三坐標(biāo)軸分別產(chǎn)生平行運(yùn)動(dòng)與旋轉(zhuǎn)運(yùn)動(dòng)，統(tǒng)稱為六自由度運(yùn)動(dòng)。

圖1 船體六自由度運(yùn)動(dòng)Fig.1 Six degrees of freedom motion of hull

當(dāng)船舶沿固定航向在海面上航行時(shí)，海浪作用于船體使姿態(tài)角改變。由于橫蕩、縱蕩和垂蕩等平行運(yùn)動(dòng)不影響搖晃速度，在航向固定時(shí)艏搖對(duì)船體影響較小，因此重點(diǎn)研究橫搖和縱搖運(yùn)動(dòng)。根據(jù)船舶水動(dòng)力理論，該類運(yùn)動(dòng)可看作是受遭遇海浪激勵(lì)與流體阻力影響的阻尼運(yùn)動(dòng)，在平衡位置上下不斷振動(dòng)，因此會(huì)在天線上疊加時(shí)變的搖晃速度。根據(jù)阿基米德定律，相對(duì)于靜止海面的波面凸起會(huì)對(duì)船底接觸面施加等效于相同體積海水質(zhì)量的轉(zhuǎn)動(dòng)作用力F，引起船體姿態(tài)角改變(圖2)。

圖2 海浪引起船體搖晃運(yùn)動(dòng)Fig.2 Shipshaking caused by waves

當(dāng)船舶沿固定航向行駛時(shí)，將船體垂直投影到二維海平面上，并劃分為1 m×1 m的網(wǎng)格，如圖3所示。將船體形狀簡化為長(L)、寬(W)、高(D)的立方體，則投影邊界框內(nèi)近似包含L×W個(gè)采樣點(diǎn)。根據(jù)時(shí)變海面位移計(jì)算每個(gè)采樣點(diǎn)當(dāng)前時(shí)刻的波高(hi，j)與轉(zhuǎn)動(dòng)作用力，減去海水流體阻力并求和得到瞬時(shí)合力，根據(jù)流體動(dòng)力學(xué)公式計(jì)算得到船體搖晃運(yùn)動(dòng)加速度，進(jìn)而得到搖晃運(yùn)動(dòng)速度與瞬時(shí)姿態(tài)角。

圖3 船體搖晃計(jì)算示意圖Fig.3 Schematic diagram for calculation of ship shaking

由圖3可知，搖晃矢量rm為轉(zhuǎn)動(dòng)徑向矢量在平行于船艏(船舷)方向的分量，搖晃力矩Mm為搖晃矢量與轉(zhuǎn)動(dòng)作用力的矢量積，表達(dá)式為

(5)

式(5)中：ρsea為海水密度；Aunit為網(wǎng)格面積;n為方向矢量。搖晃阻力(R)與作用力(F)表達(dá)式分別為

R=bIωm

(6)

(7)

式中：b為阻尼系數(shù)；ωm為搖晃角速度；K為取決于船舶類型的搖晃系數(shù)；轉(zhuǎn)動(dòng)慣量(I)表達(dá)式為

(8)

根據(jù)流體動(dòng)力學(xué)公式計(jì)算得

am=F/I,ωm=ωm+amΔt,θm=ωmΔt

(9)

式(9)中：αm為搖晃角加速度；θm為姿態(tài)角；Δt為時(shí)間步長。

2.2 船舶運(yùn)動(dòng)仿真結(jié)果分析

以2018年9月南海海試實(shí)驗(yàn)為例，船載微波多普勒雷達(dá)被安裝在“實(shí)驗(yàn)1號(hào)”科考船的船頭[8-9]。設(shè)定船舶運(yùn)動(dòng)仿真參數(shù)，海面10 m風(fēng)速分別為8、10 m/s，計(jì)算得到3 min內(nèi)5.14 m/s航速下的搖晃運(yùn)動(dòng)速度，如圖4所示。

由圖4可知，海況越高時(shí)，搖晃運(yùn)動(dòng)幅度越劇烈，角速度為1～2(°)/s，周期為8～12 s。

圖5給出了2018年9月11日一段3 min姿態(tài)傳感器實(shí)測數(shù)據(jù)，期間海面10 m風(fēng)速為8 m/s。由于傳感器安裝位置存在傾斜角，導(dǎo)致橫縱搖角存在固定偏移量。傳感器測得的橫縱搖角幅度變化范圍同樣在1～2(°)/s以內(nèi)，周期為8～12 s，與仿真結(jié)果基本一致，驗(yàn)證了該船舶運(yùn)動(dòng)仿真算法的有效性。

3 海面回波模擬算法

小斜率近似法綜合了基爾霍夫近似法與微擾法，在掠入射角下也具有較高的精度，適用于復(fù)合尺度的海面電磁散射問題。當(dāng)微波段電磁波照射在二維海面h(r,t)=h(x,y,t)上時(shí)，海面電磁散射幾何模型如圖6所示，其中Ki為入射波矢量，Ks為散射波矢量，θi、φi、θs、φs分別為入射角、入射方位角、散射角和散射方位角。

根據(jù)圖6中的幾何關(guān)系可知：

Ki=k0-q0z，Ks=k1+q1z

(10)

式(10)中：k0和q0分別是Ki投影在海表面與z軸上的分量；k1和q1分別是Ks投影在海表面與z軸上的分量。根據(jù)一階小斜率近似法計(jì)算得到散射振幅表達(dá)式為

圖4 縱搖、橫搖運(yùn)動(dòng)仿真結(jié)果Fig 4 Simulation results of pitch and roll

圖5 實(shí)測姿態(tài)角數(shù)據(jù)Fig.5 Measured attitude angle data

圖6 海面電磁散射幾何模型Fig.6 Geometrical model of EM scattering from sea surface

(11)

G(r,h)exp[-j(k1-k0)r+j(q0+q1)h(r)]

(12)

式中：Pinc為入射波功率；G(r,h)表示消除截?cái)嗾`差的錐形入射波[2]。將海水視為理想導(dǎo)體，計(jì)算得到不同極化方式下的散射系數(shù)矩陣B(k0,k1)。其中BHH代表水平極化，BVV代表垂直極化。

(13)

(14)

假設(shè)船載雷達(dá)平臺(tái)以徑向速度(v)向斜距為Ri的目標(biāo)距離元靠近，T為海面采樣時(shí)間間隔，則第i次采樣時(shí)回波信號(hào)接收時(shí)延(τi)表達(dá)式為

(15)

式(15)中：c為電磁波傳播速度。f0為入射波頻率。則船載雷達(dá)天線接收回波的散射振幅Aship為

Aship(k0,k1,t)=S(k0,k1,t)exp(-j2πf0τ)

(16)

式(16)中：根據(jù)散射振幅可推導(dǎo)得到多普勒譜Sdop為

(17)

采用能量均衡矩估計(jì)法計(jì)算多普勒譜的頻移fd和譜寬Bw，其中f表示頻率。

(18)

(19)

4 回波多普勒譜仿真與分析

船載雷達(dá)平臺(tái)下不同時(shí)刻天線照射的海面區(qū)域不同，需要根據(jù)實(shí)際運(yùn)動(dòng)情況選取每一時(shí)刻對(duì)應(yīng)的時(shí)變海面區(qū)域，并計(jì)算其電磁散射幅度。

設(shè)定入射電磁波頻率為2.85 GHz，海面10 m風(fēng)速為10 m/s，雷達(dá)采樣時(shí)間間隔設(shè)定為0.003 9 s，共采樣128點(diǎn)，計(jì)算0.5 s內(nèi)時(shí)變海面的電磁散射幅度和多普勒譜。為了減小海面隨機(jī)性的影響，方便多普勒譜的觀察與分析，仿真多普勒譜取100個(gè)隨機(jī)海面樣本求取均值。

圖7(a)為在勻速運(yùn)動(dòng)情況下船速分別為0、2.57、5.14 m/s時(shí)的海面回波多普勒譜，圖7(b)為多普勒譜頻移和譜寬隨船速變化。由圖7可知，隨著船速的增大，多普勒譜頻移呈線性增長趨勢(shì)，譜峰值逐漸減小，譜寬逐漸變大，且船速越大，譜寬增長越快。當(dāng)船速越大時(shí)，采樣時(shí)間間隔內(nèi)移動(dòng)的距離越大，參與計(jì)算的相鄰海面回波的相關(guān)性越弱，導(dǎo)致多普勒譜峰值降低與譜寬展寬。

圖7 不同船速的海面回波仿真(勻速)Fig 7 Sea echo simulation with different ship speeds (constant speed)

在實(shí)際的走航觀測中，船舶的航行并不總是理想的勻速運(yùn)動(dòng)，船速可能會(huì)在相干積累時(shí)間內(nèi)發(fā)生變化。假設(shè)雷達(dá)平臺(tái)的徑向速度v(t)是如式(20)所示的時(shí)變量。

v(t)=v0+at

(20)

式(20)中：v0為初始速度；a為恒定加速度。圖8(a)給出了初始速度是2.57 m/s，加速度分別是0、1、2 m/s2的海面回波多普勒譜，圖8(b)為多普勒譜頻移和譜寬隨船加速度變化。

由圖8可知，隨著船加速度的增大，多普勒譜頻移呈線性增長趨勢(shì)。相比勻速運(yùn)動(dòng)的情況，多普勒譜譜型對(duì)船速的變化十分敏感。當(dāng)船速改變時(shí)，譜峰值迅速減小，譜寬向速度變化側(cè)嚴(yán)重展寬，且加速度越大，譜峰值越小，譜寬增長越快。在存在初始船速的情況下，正加速度的譜寬變化曲線斜率更大，表明相同加速度的勻加速運(yùn)動(dòng)相比勻減速運(yùn)動(dòng)對(duì)譜寬的影響更明顯。

圖9分別為未添加搖晃速度、添加了海面10 m風(fēng)速為8、10 m/s時(shí)搖晃速度的海面回波多普勒譜。

計(jì)算得到頻移和譜寬如表1所示。

由表1可知，當(dāng)雷達(dá)天線添加搖晃速度時(shí)，頻移會(huì)發(fā)生改變，多普勒譜峰值減小，譜寬展寬，且海況越大時(shí)，搖晃速度越大，頻移變化量越大，譜寬展寬程度越劇烈。由于不同時(shí)刻的船體搖晃速度不同，頻移和譜寬變化量也會(huì)隨時(shí)間改變。

圖8 不同船加速度的海面回波仿真(變速)Fig.8 Sea echo simulation with different ship accelerations (variable speed)

表1 多普勒譜參數(shù)

5 基于海浪反演的模擬算法驗(yàn)證

根據(jù)二維時(shí)變海面不同距離元的相位差，擴(kuò)展得到探測范圍300～900 m共80 個(gè)距離元的距離多普勒譜，提取多普勒譜中心頻率并轉(zhuǎn)換為多距離元徑向速度序列。根據(jù)空間采樣的海浪反演算法，推導(dǎo)得到海浪譜，并反演海浪參數(shù)。圖10給出了不同均勻船速下反演海浪譜與理論靶譜的對(duì)比。

采用矩估計(jì)法對(duì)反演海浪譜作誤差分析，如表2所示。

表2 海浪譜誤差分析(勻速)

當(dāng)均勻船速增大時(shí)，有效浪高反演結(jié)果誤差增大，海浪譜相關(guān)系數(shù)減小。由圖10可知，不同船速下的反演海浪譜均在理論值附近上下波動(dòng)，整體趨勢(shì)和理論靶譜高度重合，誤差在可接受的范圍內(nèi)，驗(yàn)證了該海面回波模擬算法的正確性。

圖10 不同船速下的反演海浪譜Fig.10 Inversion wave spectrum with different ship speeds

圖11給出了初始速度為2.57 m/s、不同船加速度下反演海浪譜與理論靶譜的對(duì)比。

圖11 不同船加速度下的反演海浪譜Fig.11 Inversion wave spectrum of different ship acceleration

同樣采用矩估計(jì)法對(duì)反演海浪譜作誤差分析，如表3所示。

表3 海浪譜誤差分析(變速)

根據(jù)相干多普勒雷達(dá)原理，當(dāng)船加速度增大時(shí)，船速在采樣時(shí)間內(nèi)發(fā)生改變，天線接收的海面回波信號(hào)的相關(guān)性被削弱，導(dǎo)致多普勒譜峰值驟降與譜寬嚴(yán)重展寬，提取的徑向速度序列起伏更劇烈，反演海浪譜能量增加。由表3可知，有效浪高反演誤差明顯增大，海浪譜相關(guān)系數(shù)減小。

6 結(jié)論

基于時(shí)變海面建立了船舶運(yùn)動(dòng)模型，并通過與實(shí)測數(shù)據(jù)對(duì)比驗(yàn)證了模型的有效性。采用一階小斜率近似法計(jì)算海面電磁散射幅度，仿真得到了船載雷達(dá)平臺(tái)下的海面散射回波，并分析了不同運(yùn)動(dòng)方式(勻速運(yùn)動(dòng)、變速運(yùn)動(dòng)、搖晃運(yùn)動(dòng))對(duì)多普勒譜特征的影響。結(jié)果表明，船舶運(yùn)動(dòng)會(huì)導(dǎo)致多普勒譜譜型改變，船速和船加速度越大，頻移和譜寬變化越明顯。根據(jù)空間海浪反演算法計(jì)算得到海浪譜并與理論值對(duì)比分析，驗(yàn)證了該海面回波模擬算法的正確性，研究結(jié)果對(duì)船載雷達(dá)平臺(tái)下海浪反演與算法改進(jìn)具有重要的指導(dǎo)意義。