邱吉東,張 彪,陳忠彪,何宜軍
(南京信息工程大學(xué) 海洋科學(xué)學(xué)院,江蘇 南京 210044)
基于X波段海洋雷達的風(fēng)浪聯(lián)合反演方法研究
邱吉東,張 彪,陳忠彪,何宜軍
(南京信息工程大學(xué) 海洋科學(xué)學(xué)院,江蘇 南京 210044)
提出了一種新的利用X波段海洋雷達聯(lián)合反演海面風(fēng)速與海浪譜的方法,該方法不需要額外的信息輸入來反演海浪譜。通過利用風(fēng)速與雷達后向散射強度的經(jīng)驗關(guān)系獲得海表風(fēng)速,然后將反演的風(fēng)速輸入風(fēng)浪譜,通過求解該模擬風(fēng)浪譜與雷達觀測圖像譜的約束函數(shù)的最小值來確定海浪譜。利用實驗數(shù)據(jù)對反演方法進行了驗證,風(fēng)速、有效波高、主波周期以及主波波向反演的均方根誤差分別為1.9m/s,0.4m,1.2 s和9.6°,證明了該方法的可行性。
X波段海洋雷達;風(fēng)速;海浪譜;聯(lián)合反演
X波段海洋雷達被應(yīng)用于海表物理參數(shù)反演已經(jīng)有近半個世紀的歷史。從1965年Wright從 X波段雷達圖像直接進行海浪傳播方向和波長的判讀[1],到1983年,Ziemer和Rosenthal對雷達圖像應(yīng)用二維傅里葉變換獲得海浪波數(shù)譜[2],并利用波數(shù)譜估算了海浪參數(shù),再到Y(jié)oung提出對雷達圖像序列進行三維傅立葉變換[3],利用X波段雷達海浪反演方法一直在不斷進步。
本世紀初X波段雷達在風(fēng)速反演方面上取得了飛速的發(fā)展,Heiko Dankert等提出利用神經(jīng)網(wǎng)絡(luò)算法和光流場算法估算風(fēng)速和風(fēng)向[4-5];Lund等提出用最小二乘法擬合得到的風(fēng)速與諧波函數(shù)的多項式關(guān)系[6-7];Vicen-Bueno等建立了基于雷達圖像灰度值的物理模型函數(shù)來反演風(fēng)速[8];W Huang提出基于譜分析的風(fēng)速反演方法[9]。上述方法在風(fēng)速和海浪的反演上都是單獨反演,并且需要分別定標(biāo)。
本文在現(xiàn)有研究的基礎(chǔ)上提出了一種新的風(fēng)速反演經(jīng)驗方程及風(fēng)浪聯(lián)合反演方法,只要在風(fēng)速反演上進行定標(biāo),對海浪的反演不需要額外的定標(biāo)。第二節(jié)將介紹風(fēng)速和海浪譜的反演方法;第三節(jié)將給出反演方法的結(jié)果驗證;最后對全文進行了總結(jié)與討論。
1.1 風(fēng)速的反演方法
首先對雷達圖像序列I(θ,r,t)求平均
式中:T為雷達圖像序列的時間總長度;θ為方位向;r為距離向;t為時間。
對雷達平均圖像的每一方位向求和:
式中:M為單位方位向上圖像元的總數(shù)。E反映了雷達回波能量在方位向上的總體水平。
將所有方位向上的E求平均得到S:
式中:Θ表示雷達圖像的方位向總數(shù);S反映了雷達圖像序列的總體后向散射回波的水平[9]。
海面后向散射強度與風(fēng)速存在一定的非線性關(guān)系[8,10],經(jīng)驗的物理模型函數(shù)可以確定風(fēng)速:
式中:a,b和c為待定系數(shù);u10為海表10 m風(fēng)速。通過擬合雷達圖像序列總體散射值S與風(fēng)速的關(guān)系,可以標(biāo)定a,b和c來反演風(fēng)速。
本文提出了新的經(jīng)驗物理模型函數(shù)來反演風(fēng)速,具體形式如下:
式中:a,b,c和d都是待定系數(shù)。
1.2 海浪的反演方法
1.2.1 海浪譜的數(shù)值仿真 本文中海浪頻譜選用了JONSWAP譜[11-13],其形式如下:
式中:g為重力加速度;ωp為譜峰頻率;α為能量尺度參數(shù),本文中取α=0.003;γ為譜升因子;σ為峰形參數(shù);γ和σ分別由下式確定:
式中:cp為波動傳播峰值速度。由于本研究試驗資料在近岸地區(qū)獲取,根據(jù)線性波動理論,在水深較淺的區(qū)域,波動傳播速度由水深決定,cp=為重力加速度,h為水深。
實際海面是二維的,所以還應(yīng)該考慮海浪能量在不同方向上的分布,這就需要引入方向譜。方向譜可以寫出如下形式:
其中方向函數(shù)G(ω,θ)選用光易型方向函數(shù):
式中:Γ為伽馬函數(shù);s為角擴散系數(shù),與頻率與風(fēng)速有關(guān):
Goda建議的峰值角擴散系數(shù)smax如表1所示[14]:
表1 峰值角擴散系數(shù)smax取值表
1.2.2 雷達圖像譜的提取 本文采用了Young等提出的利用傅里葉變換獲取雷達圖像譜的方法[3]。首先假定選擇分析的波場具有平穩(wěn)性的和各項同性,我們可以對其整個波場應(yīng)用三維傅立葉變換。對雷達灰度序列圖像I(x,y,t)三維傅立葉變換得到三維波數(shù)頻率雷達圖像譜f(kx,ky,ω),即
雷達圖像譜f(kx,ky,ω)的能量主要包括以下3部分:
(1)波浪成分,這部分能量占圖像譜的主要部分;(2)由斑點噪聲引起的背景噪聲能量;
(3)高次諧波能量,這是由于雷達的非線性成像機制引起的。
之前的研究通過利用調(diào)制傳遞函數(shù)的方法提取海浪成分[15-17],而調(diào)制傳遞函數(shù)又受到距離向、方位向的影響,單一的調(diào)制傳遞函數(shù)在反演海浪譜時具有一定的限制性,本文通過海浪譜仿真擬合的方法不需要利用到調(diào)制傳遞函數(shù)。
1.2.3 海浪譜的反演 海浪譜反演的流程圖如下:
圖1 海浪譜反演流程圖
本文中的海浪譜仿真只與3個參數(shù)有關(guān):(1)10m風(fēng)速u10;(2)峰值波數(shù)kp;(3)波浪主波傳播方向θp。10m風(fēng)速可以通過反演得到,所以只需要確定峰值波數(shù)kp和波浪主波傳播方向θp這兩個變量就可以確定風(fēng)浪譜。假設(shè)風(fēng)浪在風(fēng)向上是充分成長的,這一假設(shè)適用于大部分風(fēng)浪占主導(dǎo)的情況,但在涌浪占主導(dǎo)的情況下會出現(xiàn)較大誤差。所以在方向譜峰值角擴散系數(shù)Smax的選擇上,本文使用了風(fēng)浪情況下的Smax=10。
首先通過雷達圖像序列反演得到風(fēng)速和雷達圖像譜,然后將歸一化的仿真海浪譜和歸一化雷達圖像譜帶入約束函數(shù)。由于需要通過約束函數(shù)確定最接近雷達觀測圖像譜的仿真海浪譜,所以約束函數(shù)中必須包括歸一化的仿真海浪譜、歸一化的雷達圖像譜以及歸一化的兩譜之差,所以確定約束函數(shù)形式[18]如下:
式中:FI為歸一化雷達圖像譜;Fs為歸一化仿真海浪譜。
設(shè)定峰值波數(shù)kp的步長為0.001 rad/m,主波波向θp的步長為1°,通過迭代,求得J的最小值,其對應(yīng)的kp和θp就是需要的結(jié)果。此時的仿真海浪譜即可認為是有效海浪譜。最后通過對海浪譜進行計算,可以得到有效波高Hs、主波波長λ和主波周期T等參數(shù),如下所示:
2.1 實驗數(shù)據(jù)
為了驗證反演算法及測試儀器的性能,我們于2015年1月在福建平潭島進行了觀測實驗。平潭島在臺灣海峽以西,島嶼南北長29 km,東西寬19 km,是全國第五大島。實驗的時間在冬季,主要原因是這段時間寒潮的影響比較頻繁,因而經(jīng)常有高海況出現(xiàn),風(fēng)力最高可以達到25 m/s,風(fēng)向常年以北風(fēng)和東風(fēng)為主。
本文實驗所用的導(dǎo)航X雷達是基于日本FURUNO公司生產(chǎn)的船載導(dǎo)航X波段雷達改裝而成。天線的水平波束寬度是1.2°,電磁波的中心頻率是9.41GHz,對應(yīng)的電磁波波長為3.19 cm,極化方式為VV極化。天線的轉(zhuǎn)速是24 rot/min,即旋轉(zhuǎn)周期為2.5 s,用本系統(tǒng)觀測海浪時一般每組數(shù)據(jù)(即每個雷達圖像序列)記錄32幅雷達圖像,所以共需要80 s。
圖2為2015年1月15日17時00分所采集的圖像序列中的一幅圖像,其中黑框范圍內(nèi)為反演程序利用的256×256大小的數(shù)據(jù)區(qū)域,星號代表的是浮標(biāo)所處位置。
圖2 2015年1月15日17時00分雷達圖像,黑框內(nèi)為計算采樣范圍,星號為浮標(biāo)位置
2.2 結(jié)果分析
2.2.1 風(fēng)速的反演結(jié)果 本文對2015年1月在福建平潭試驗中獲取的810組風(fēng)速與雷達相對應(yīng)的有效數(shù)據(jù)進行分組,其中539組用于訓(xùn)練擬合經(jīng)驗?zāi)P停?71組用于驗證風(fēng)速反演結(jié)果。
圖3為分別利用方程(4)和方程(5)對訓(xùn)練數(shù)據(jù)組進行擬合的結(jié)果,黑色叉號為10 m風(fēng)速與雷達后向散射系數(shù)S的對應(yīng)散點,紅色實線為現(xiàn)有經(jīng)驗?zāi)P蛿M合線(S=226.7log(u10+0.75)+499),藍色實線為本文提出的新經(jīng)驗公式擬合線(S=106tanh (0.3292u10-1.862)+668.8)。從圖中可以明顯看出,本文提出的新經(jīng)驗公式的結(jié)果更優(yōu),包括相關(guān)系數(shù)更高,均方根誤差更低。這是由于在低風(fēng)速下,風(fēng)速對海面的粗糙度影響是有限的,這也就造成了在低風(fēng)速時,海面后向散射強度S隨風(fēng)速增長緩慢;而在風(fēng)速達到或超過15m/s時,海面后向散射強度S又會趨于飽和。
圖4為利用不同經(jīng)驗?zāi)P蛯︼L(fēng)速的反演結(jié)果,單位為m/s。藍色圈為新經(jīng)驗?zāi)P头囱萁Y(jié)果,紅色星為現(xiàn)有經(jīng)驗?zāi)P头囱萁Y(jié)果。其中現(xiàn)有研究的經(jīng)驗?zāi)P蚚8-9]反演結(jié)果的均方根誤差 (Root Mean Square Error,RMSE)為2.9 m/s,平均相對誤差(Mean Relative Error,MRE)為0.1;而本文提出的新經(jīng)驗?zāi)P头囱萁Y(jié)果的均方根誤差為1.9m/s,平均相對誤差為-0.06。由此可見,新經(jīng)驗?zāi)P偷姆囱萁Y(jié)果相對于現(xiàn)有模型有了較大程度上的改進。
圖3 10m風(fēng)速與雷達后向散射系數(shù)S擬合圖像
圖4 雷達反演風(fēng)速與實測風(fēng)速對比圖
2.2.2 海浪的反演結(jié)果 圖5(a)和圖5(b)分別為2015年1月8日5時0分在福建平潭島試驗獲取的雷達圖像譜和與之對應(yīng)的反演獲取的海浪譜。浮標(biāo)測量的有效波高為2.57m,主波波向為75°,主波周期為6.7 s;風(fēng)速儀測得10min平均風(fēng)速為17.81 m/s,10min平均風(fēng)向為84.17°。反演海浪譜計算的有效波高為1.81m,主波波向為85.38°,主波周期為5.62 s。
圖5 2015年1月8日5時0分平潭島試驗獲取的歸一化雷達圖像譜和歸一化反演海浪譜
圖6(a)和圖6(b)分別是2015年1月15日17時0分在福建平潭島試驗獲取的雷達圖像譜和與之對應(yīng)的反演獲取的海浪譜。浮標(biāo)測量的有效波高為2.06 m,主波波向為59.8°,主波周期為7 s;風(fēng)速儀測量的10min平均風(fēng)速為13.07m/s,10min平均風(fēng)向為65.75°。反演海浪譜計算的有效波高為2.02 m,主波波向為61.3°,主波周期為5.65 s。
圖6 2015年1月15日17時0分平潭島試驗獲取的歸一化雷達圖像譜和歸一化反演海浪譜
圖7(a)和7(b)分別為2015年1月6日16時0分在福建平潭島試驗獲取的雷達圖像譜和與之對應(yīng)的反演獲取的海浪譜。浮標(biāo)測得海浪有效波高為2.22m,主波波向77.6°,主波周期為7 s。從圖7(a)中可以看出,在圖像譜中,除了能量主要集中第一象限外,在第三象限仍有明顯的能量波峰存在,圖像譜能量分布存在多個峰值,反映出海浪并不是單一的風(fēng)浪,而本文所使用的參數(shù)化海浪譜為風(fēng)浪譜,導(dǎo)致這一類型的圖像譜對于仿真海浪譜的擬合存在非常大的誤差。在反演結(jié)果中有效波高僅為0.97m,主波波向為90°,主波周期為2.92 s。
圖7 2015年1月6日16時0分平潭島試驗獲取的歸一化雷達圖像譜和歸一化反演海浪譜
圖8~圖10分別是2015年1月6日至15日在福建平潭島試驗的99組有效數(shù)據(jù)的有效波高、主波周期和主波周期的觀測數(shù)據(jù)與反演數(shù)據(jù)的對比圖。風(fēng)浪的篩選標(biāo)準(zhǔn)為浮標(biāo)測量波向與風(fēng)速計測量風(fēng)向偏差在20°以內(nèi),風(fēng)速大于5 m/s。其中有效波高的均方根誤差為0.4m,平均相對誤差為0.1;主波周期的均方根誤差為1.2 s,平均相對誤差為-0.1;主波波向的均方根誤差為9.6°,平均相對誤差為0.05。
圖8 2015年1月6日至15日平潭試驗有效波高對比圖
圖9 2015年1月6日至15日平潭試驗主波周期對比圖
本文提出了一種新的利用X波段海洋雷達聯(lián)合反演海面風(fēng)速與海浪譜的方法,該方法不需要額外的信息輸入來反演海浪譜。在風(fēng)速反演部分,提出的經(jīng)驗物理模型相較于現(xiàn)有模型,有效提升了風(fēng)速反演的精度。在海浪反演部分,通過對比實驗數(shù)據(jù)和雷達反演數(shù)據(jù),對反演方法進行了驗證,風(fēng)速反演的均方根誤差達到1.9m/s,有效波高的均方根誤差為0.4 m,主波周期的均方根誤差為1.2 s,主波波向的均方根誤差為9.6°,證明了該方法的可行性。與現(xiàn)有研究相比,該方法還具有不需要對有效波高反演進行定標(biāo)的優(yōu)點。
在未來的工作中,將繼續(xù)尋找有效方法提高風(fēng)速反演的精度,尤其是低海況情況下風(fēng)速的反演和涌浪占主體或者混合浪的海浪譜的反演。另外需要采用多地點長時間序列的觀測數(shù)據(jù)對方法進行進一步驗證,并提高風(fēng)速和海浪譜的聯(lián)合反演精度與反演效率。
圖10 2015年1月6日至15日平潭試驗主波波向?qū)Ρ葓D
[1]WrightFF.WaveObservationsby Ship-Board Radar[J].Ocean SciOcean Eng,1965,1:506-514.
[2]Ziemer F,RosenthalW,Carlson H.MeasurementsofDirectionalWave Spectraby Ship Radar[C]//IAPSOSymp,General Assembly,Int Assoc Phys SciOceans,Hamburg,Germany,1983.
[3]Young I R,Rosenthal W,Ziemer F.A Three-Dimensional Analysis of Marine Radar Images for the Determination of Ocean Wave Directionality and Surface Currents[J].JournalofGeophysicalResearch Oceans,1985,90(C1):1049-1059.
[4]DankertH,RosenthalW.Ocean Surface Determination from X-Hand Radar-Image Sequences[J].Journal of Geophysical Research Oceans,2004,109(C4):1-11.
[5]王劍,段華敏.X波段雷達圖像提取海洋表面風(fēng)場[J].海洋技術(shù)學(xué)報,2010,29(3):5-8.
[6]Lund B,Graber H C,Romeiser R.Wind Retrieval From Shipborne Nautical X-Band Radar Data[J].IEEETransactionson Geoscience &Remote Sensing,2012,50(10):3800-3811.
[7]王慧,盧志忠.基于波數(shù)能量譜的海面風(fēng)向反演算法[J].華中科技大學(xué)學(xué)報:自然科學(xué)版,2014(12):96-100.
[8]Vicenbueno R,Horstmann J,TerrilE,etal.Real-TimeOceanWind VectorRetrieval from Marine Radar Image SequencesAcquired at Grazing Angle[J].JournalofAtmospheric&Oceanic Technology,2013,30(30):127-139.
[9]HuangW,Wang Y.A Spectra-Analysis-Based Algorithm forWind Speed Estimation From X-Band NauticalRadar Images[J].IEEE Geoscience&Remote Sensing Letters,2016,13(5):701-705.
[10]Chen Z,He Y,Zhang B,etal.Determination of Nearshore Sea SurfaceWind Vector from Marine X-Band Radar Images[J].Ocean Engineering,2015,96:79-85.
[11]Hasselmann D E,Dunckel M,Ewing JA.DirectionalWave Spectra Observed During JONSWAP 1973[J].Journal of Physical Oceanography,1980,10(8):1264-1280.
[12]崔利民.X波段雷達海浪與海流遙感機理及信息提取方法研究[D].青島:中國科學(xué)院研究生院(海洋研究所),2010.
[13]吳艷琴,吳雄斌,程豐,等.基于X波段雷達的海洋動力學(xué)參數(shù)提取算法初步研究[J].遙感學(xué)報,2007,11(06):817-825.
[14]Goda Y.Random Seasand Design ofMaritime Structures[M].World Scientific,2000:443.
[15]PlantW J.The Modulation Transfer Function:Conceptand Applications[M]//Radar Scattering from Modulated Wind Waves.Springer Netherlands,1989:155-172.
[16]Chen Z,He Y,Zhang B,etal.ANew Algorithm to RetrieveWave Parameters From Marine X-Band Radar Image Sequences[J].IEEE Transactionson Geoscience&Remote Sensing,2014,52(7):4083-4091.
[17]王淑娟.X波段雷達圖像的海浪信息提取[D].青島:中國海洋大學(xué),2007.
[18]Zhang B,LiX,PerrieW,etal.Synergistic MeasurementsofOceanWinds and Waves from SAR[J].JournalofGeophysicalResearch Oceans,2015,120(9):6164-6184.
Research on theWind-Wave Synergistic Retrieval Method Based on X-Band Marine Radar Observations
QIU Ji-dong,ZHANG Biao,CHEN Zhong-biao,HE Yi-jun
School ofMarine Sciences,Nanjing University of Information Science and Technology,Nanjing 210044,Jiangsu Province,China
This paper presents a new method for synergistic retrieval of ocean surface wind speeds and wave spectrum based on X-band marine radar observations.The proposed algorithm needs not any external data source to obtain wave information.The surface wind speeds are first estimated based on an empirical relationship between the wind speed and radar scattering intensity.Subsequently,the ocean surface wave spectrum is retrieved byminimizing the constraint function including the simulated wave spectrum and observed radar image spectrum.The retrieved wind and wave parameters are compared with thosemeasured by anemometers and buoys. The rootmean square error(RMSE)ofwind speed,significantwave height,and dominantwave period and wave direction is 1.9m/s,0.4m,1.2 s and 9.6°,respectively.
X-band marine radar;wind speed;wave spectrum;synergistic retrieval
P714;P225.1
A
1003-2029(2017)01-0001-06
10.3969/j.issn.1003-2029.2017.01.001
2016-08-15
江蘇省青年基金資助項目(BK20150905)
邱吉東(1991-),男,碩士研究生,主要研究方向為海洋微波遙感。Email:jidong.qiu@nuist.edu.cn