機載波譜儀海浪譜反演方法及其驗證

李秀仲, 何宜軍, 孟俊敏, 張振華

(1. 南京信息工程大學, 江蘇 南京 210044; 2. 國家海洋局第一海洋研究所, 山東 青島 266061; 3. 北京遙測技術研究所, 北京 100076)

機載波譜儀海浪譜反演方法及其驗證

李秀仲1, 何宜軍1, 孟俊敏2, 張振華3

(1. 南京信息工程大學, 江蘇 南京 210044; 2. 國家海洋局第一海洋研究所, 山東 青島 266061; 3. 北京遙測技術研究所, 北京 100076)

為了驗證波譜儀反演二維海浪譜的功能, 根據海浪波譜儀的信號形成機制, 總結了機載波譜儀反演海浪的流程。利用機載波譜儀回波數據, 通過自相關和互相關兩種功率譜估計方法, 反演了二維海浪譜。最后通過與浮標測量的二維海浪譜進行對比, 驗證了該機載波譜儀探測二維海浪譜的有效性。結果表明, 無論采用自相關函數還是互相關函數進行功率譜估計, 得到的主波波長和有效波高與實際二維海浪譜基本一致。互相關函數法得到的交叉譜能去除180°模糊現象, 其在計算有效波高時相對于自相關函數會稍微偏小。在計算斜率方差時可以采用5°~12°入射角范圍的后向散射系數進行公式擬合, 因此定標與否并不影響最后的二維海浪譜結果, 未來星載波譜儀只有靠多波束聯合才能實現。

波譜儀; 功率譜估計; 二維海浪譜; 180°模糊

海況信息的獲取有著重要的業務化需求(如海洋氣象預報、導航、近海和沿海活動等)和科研需求(海浪動力分析、大氣和海洋邊界層的相互作用、區域和整體的海洋和大氣系統耦合, 電磁信號與海面相互作用等)。過去幾十年, 雖然海浪預報已經進步很大[1], 但是在高海況條件、極值事件(如臺風、風暴潮等)方面還有很大的改進余地; 在提供準確參數預報方面不僅可以預報有效波高, 還可以改進峰值波長和波向預報。目前應用的海浪數值模式在預測二維海浪譜方面面臨很大缺點[1], 比如海浪激發源和耗散項的參數化以及采用的初始場存在誤差較大等。此外, 由于數值模型離散化的限制, 在海浪能量的方向分布上, 開闊海域時會產生較大誤差, 這樣就會導致波浪在靠近海岸時的演化結果中產生更大的預報誤差。

目前唯一星載測量海浪的雷達系統是合成孔徑雷達(Synthetic Aperture Radar, SAR)。自從Seasat衛星上搭載了SAR后[2], SAR就一直比較昂貴并且在測量海浪方向譜時并不簡單易用。為了去除SAR這種限制, 人們發展了其它類型雷達: Walsh等[3]1985開發了Ka波段的海面等值線雷達(surface contour radar, SCR), Jackson等[4-5]1985年開發了Ku波段雷達海浪波譜儀(Radar Ocean Wave Spectrometer, ROWS), Hauser等[6]1992年開發了C波段雷達測波儀(Radar pour l’Etude du Spectre des Surfaces par AnalyseCirculaire, RESSAC), Hauser等[7]2003年開發了C波段極化雷達觀測系統(Système de Télédétection pour l’Observation par Radar de la Mer, STORM), 以及Plant等[8]2005年開發了X波段相干真實孔徑雷達(coherent real aperture radar, CORAR)。這些機載測量雷達比方向浮標提供了更多的海浪方向性信息[9]。

中法海洋衛星(China–France Oceanography Satellite, CFOSAT)是一個由中國和法國空間機構正在準備的新衛星。它是為了同步海浪譜和海面風場觀測而設計的。該衛星上將搭載一個Ku波段的散射計測量風速, Ku波段的波譜儀(Surface Waves Investigation and Monitoring, SWIM)來測量海面浪場。CFOSAT是一個極軌衛星, 飛行高度為520 km。SWIM雷達能向海面發射6個不同入射角的波束: 0°、2°、4°、6°、8°和10°, 每個波束天線孔徑都是2°,天線旋轉速度是5.6 rad/min。SWIM能提供如下的地球物理參數: 0°~10°入射角的標準化后向散射系數(normalized radar cross-section, NRCS); 星下點的有效波高和風速; 6°、8°、10°波束能提供海浪方向譜。

在CFOSAT的準備階段, 北京遙測技術研究所開發了機載波譜儀, 并在山東半島以南的黃海區域進行了校飛試驗, 以檢驗SWIM測量海浪方向譜的能力。本文針對機載波譜儀70圈的回波數據, 通過總結國內外反演算法, 反演了二維海浪方向譜, 并與浮標實測的二維海浪譜對比, 檢驗了波譜儀測量海浪二維譜的有效性, 為星載波譜儀的反演算法的開發奠定了基礎。

圖1 處理的波譜儀校飛試驗數據所在軌跡Fig. 1 Tracks of the processed data obtained using an airborne sea-wave spectrometer

表1 機載波譜儀主要參數Tab. 1 Main parameters of the airborne sea-wave spectrometer

1 機載波譜儀校飛介紹

機載波譜儀的校飛試驗位置選在黃海北部, 圖1中三角形和六角形分別為2014年6月12日10: 21和2014年6月21日11: 07浮標所在位置(分別是122.84°E、36.62°N和122.61°E、36.62°N), 浮標每半小時給出一個海浪譜結果。校飛時無風速計實測風速, 但根據浮標所測海浪波高較小可知, 兩個架次的風速在2~5 m/s變化。三角形附近的粗曲線為飛機2014年6月12日11: 6: 20~11: 11: 20共301 s的飛行軌跡, 也是本文處理的30圈回波數據的軌跡, 始末位置相差15.62 km, 起始位置與浮標相差7.17 km;細直線為飛機2014年6月21日10: 48: 58~10: 55: 37共400 s的飛行軌跡, 是本文處理的40圈的回波數據軌跡, 始末位置相差28.54 km, 與六角形所代表的浮標最近距離約為11.58 km。

對于機載波譜儀, 由于飛行高度不可能達到衛星的高度, 因此采用一個大波束, 以增加足印面積。該機載波譜儀的主要參數在表1中給出。

圖2顯示了此次機載波譜儀校飛的幾何示意圖,根據3 dB波束寬度, 每個脈沖在海面上形成的等天線增益值線可以簡化成一個橢圓形狀。根據傾斜調制的原理, 雷達接收圖中橢圓內圓環里的海面反射能量, 當波浪局地傾斜導致的海面法線對著雷達時海面回波強, 反之當波浪局地傾斜導致的海面法線偏離雷達時, 雷達接收的回波弱。因此雷達信號變化的幅值與海浪局地傾角有關。根據Jackson等[4-5]、Hauser等[10], 海面斜率與后向散射系數的關系為:

圖2 機載波譜儀觀測幾何Fig. 2 Geometry of observation of the airborne sea-wave spectrometer

式中, p為海面斜率概率密度函數, a為定標參數,x為波面高度, q為入射角。因此獲取海面斜率的雷達調制信號, 求取調制信號功率譜, 可得某一方向的海浪斜率譜通過雷達天線360°水平面旋轉, 即可得到二維方向譜, 即公式(3)中的

式中，Ly為波束方位向3 dB寬度,為調制譜, k為波數, j為方位角。

圖3模擬了機載波譜儀校飛1 min形成的足跡。波譜儀在海面形成的足跡與以下幾個參數有關: 飛機飛行高度、飛機飛行速度、天線旋轉速度以及脈沖重復頻率。為了顯示清晰, 該圖只顯示了5 Hz的脈沖重復頻率結果, 飛行高度為3 000 m、飛機飛行速度為60 m/s、天線轉速是6 rad/min。考慮到準鏡面散射理論適用范圍以及其水平分辨率因素, 圖中展示了入射角5°~12°時的波束足印, 水平徑向距離為375 m, 短半軸長則為214 m。

圖3 機載波譜儀探測足跡Fig. 3 Footprint detected by the sea-wave spectrometer

2 機載海浪譜反演方法

本文根據波譜儀的工作原理, 在國內外波譜儀研究基礎上, 采用以下流程進行了二維海浪譜的反演。

圖4中, 回波數據是雷達接收的不同時刻電磁波能量, 天線增益是由實驗室實測后的二維天線增益經過飛機姿態角校正后得到的增益[11], r代表斜距, m( r)是沿斜距的調制信號, 由本文公式(4)中的組成, x是水平距離,是沿水平方向的調制信號。此處海浪斜率譜為由調制譜根據公式(3)轉換得到斜率譜為標準化過程,j為方位角。

2.1 計算0σ

波譜儀獲取回波能量之后, 為了獲取海浪的無量綱的標準化后向散射系數0s, 采用雷達方程為

式中,rP和tP是接收和發射的能量, R是斜距,tG和rG是天線發射和接收的增益, l是單位為米的載頻電磁波長。雷達接收的截面(Radar Cross-Section, RCS)s用公式(5)表達

圖4 機載波譜儀海浪譜反演方法Fig. 4 Retrieving method of sea-wave spectrum using an airborne spectrometer

S是電磁波波束照射的面積, 這個面積能夠用一個矩形的寬度dx和長度dy來表示, 并且是方位向的波束寬度, 所以。將此公式代入公式(4), 并令可得到公式(6)

2.2 計算斜率方差和調制系數

基于小入射角下NRCS, 有一個簡便的方法能用來計算斜率方差(mean square slope, MSS)。Barrick[12]和Valenzuela[13]導出了有限粗糙面的電磁散射公式, 總結出小入射角范圍是0°~15°時準鏡面散射是海面散射的主要部分。準鏡面散射公式是此處天底點的菲尼爾反射系數,是二維海面斜率概率密度分布函數, zx表示距離向海面斜率, zy表示方位向海面斜率。正比于滿足條件的微小面元的數量, 所以也可以寫成。如果是高斯分布, 那么等式(8)變成

此處us和cs分別是海面斜率在逆風向和側風向的標準偏差。MSS是方位角j處的斜率方差。如果海面斜率概率分布函數(Probability Density Function, PDF)是高斯的, 那么我們可以得到[14-15]

所以, 當0s通過式(7)獲取后, 利用方位角j處的回波對上式擬合, 可以獲取MSS。此處選擇5°入射角作為擬合上限是考慮了地面分辨率隨著入射角減小而變差這個因素。Freilich等[16]假設入射角0°~18°時降雨雷達散射能量在海面風速0~20 m/s范圍內主要是由準鏡面散射機制產生的, 其中較大入射角時布拉格散射機制變得相對重要了。此處選擇12°入射角作為上限并假設此時準鏡面散射占主要部分是合適的。需要注意的是, 我們獲取的是系數的負倒數, 所以0s是否進行了定標并沒有關系。

根據海面斜率分布近似高斯函數形式的結論,則通過公式(2)調制系數計算公式變為

由于海浪譜是由信號調制譜得來的, 因此計算海浪譜之前, 需要通過調制信號計算調制譜。采用兩種方法估計調制譜, 第一種是自相關函數功率譜估計方法

該方法是采用同一信號進行自相關得到的功率譜, 在進行譜估計之前, 需要對調制信號進行加窗處理, 以消除由于信號截斷造成的頻譜泄露問題。第二種是交叉譜估計方法

j￠為同圈內與j有一定時間間隔的方位角, 交叉譜計算時, 主要針對相鄰的兩個回波, 以達到去除噪聲的目的。

這兩種方法求取的海浪譜都存在180°模糊的問題, 本文借鑒SAR消除180°模糊的方法[17], 根據公式(13)中求取的互相關函數, 則觀測到的海浪傳播的徑向軌道速度為

2.4 海浪參數的確定

有效波高(significant wave height, SWH)HSW代表了海浪能量的大小, 其計算公式為

主波波數kp為F(k,j)最大值所對應的波數, 主波波長則為主波波向為方向譜取最大值處對應的方向j。在求取浮標所測海浪譜主波波長和主波波向時, 對于浮標提供的頻率譜, 需要轉換成波數譜。根據

w是角頻率, f為波浪頻率, h和g分別為水深和重力加速度, 水深h在30 m左右, 根據上式可以得到每個頻率對應的波數。

3 二維海浪譜及其參數的驗證

浮標測量的二維海浪譜如圖5、圖6所示, 其中HSW分別為0.47, 0.6 m, 主波波長是83, 39 m, 波浪傳播方向是南偏東20°, 12°, 若正北為0°, 逆時針為正, 則浪向為200°, 192°。

圖5 2014年6月12日10: 21波浪騎士所測二維海浪譜Fig. 5 Two-dimensional spectrum detected from a sea wave rider at 10: 21 on June 12, 2014

圖6 2014年6月21日11: 07波浪騎士所測二維海浪譜Fig. 6 Two-dimensional spectrum detected from a sea wave rider at 11: 07 on June 21, 2014

由于本次校飛時海況較低, 因此海浪信息并不夠規則, 此處將這兩塊數據所得結果各自進行了平均。利用前面兩種功率譜計算方法得到兩種海浪譜反演結果如下圖, 對于第一塊數據, 自相關函數方法得到的結果顯示, 主波波長93 m, 波浪傳播方向為30°, 與浮標結果相反, 有效波高為0.46 m。而交叉譜方法得到的這3個參數分別為93 m、210°、0.41 m, 海浪傳播方向基本與浮標所測結果一致。

第2塊數據中, 自相關函數方法得到的主波波長、波浪傳播方向、有效波高分別為128 m, 270°和 0.57 m, 交叉譜方法得到的這3個參數分別為37 m、170°、0.55 m。很顯然, 該結果中自相關函數法由于噪聲的影響無法測得正確的結果。

圖7 2014年6月12日自相關函數法所得海浪譜Fig. 7 Sea-wave spectrum obtained from the auto-correlation method on June 12, 2014

圖8 2014年6月12日互相關函數法所得海浪譜Fig. 8 Sea-wave spectrum obtained from the cross-correlation method on June 12, 2014

圖9 2014年6月21日自相關函數法所得海浪譜Fig. 9 Sea-wave spectrum obtained from the auto-correlation method on June 21, 2014

圖10 2014年6月21日互相關函數法所得海浪譜Fig. 10 Sea-wave spectrum obtained from the cross-correlation method on June 21, 2014

這兩種功率譜估計方法所測主波波長和有效波高與浮標一致性都較好, 誤差在可接受范圍之內。根據交叉譜方法中相位差的方法, 基本可以去除180°模糊。圖7顯示方位向30°時能量大于210°, 這應該是由于在210°方位向時異常回波較多, 因此會有較多的回波被刪除掉, 導致結果顯示波浪傳播在30°方位向。圖8中顯示有幾圈180°模糊并沒有去除掉, 因此會在30°方位向時產生一個較小的峰值。交叉譜方法得到的有效波高較小, 這應該與互相關函數功率譜估計的過程中用到了相鄰兩束回波進行互相關,相關性強度必定會比自相關函數減弱, 但這樣做的優點是可以去除兩束回波中不相關的信息, 對去噪有一定幫助。

若都采用互相關函數法, 6月12日結果得到的海浪傳播方向10°的誤差比21日22°的誤差小, 但主波波長10 m的誤差大于21日2 m的誤差。12日和21日校飛時海面有效波高分別是0.47 m和0.6 m, 所以海況均很小, 若單獨處理出某幾圈的結果, 誤差會較大, 因此本文兩塊數據分別進行了30圈數據的平均和40圈數據的平均。通過本次結果處理可知, 交叉譜方法得到的海浪譜在去除180°模糊方面有較大優勢, 并且在去噪方面有一定優勢, 因此交叉譜比自相關函數譜更有優勢。由于此次校飛試驗所在的海況環境較低, 對于更高海況和更大的飛行高度尚沒有試驗, 以及風浪和涌浪混合時海浪反演的準確性無法得到驗證。

4 結論

本文分析了二維海浪譜的近實時探測在業務化和科研方面的需求, 介紹了未來中法星星載波譜儀探測海浪譜的幾何特性, 根據機載波譜儀校飛試驗和機載波譜儀的信號形成機制, 總結了波譜儀反演海浪的流程, 并利用該流程處理了機載校飛的回波數據, 通過自相關函數和互相關函數兩種功率譜估計方法, 反演了二維海浪譜。最后根據浮標得到的二維海浪譜對機載波譜儀探測二維海浪譜的有效性進行了驗證, 并對比了兩種海浪譜計算方法。主要結論如下:

1) 工作在Ku 波段的機載海洋波譜儀是一種真實孔徑雷達系統, 波譜儀波束基于長波對海面微尺度波的傾斜調制, 通過方位向360°掃描測量海浪譜,通過本文的驗證, 這種方法探測二維海浪譜以及波向、波長和有效波高是有效的。

2) 根據波譜儀信號形成機制, 采用互相關函數進行功率譜估計在獲取主波波長、波向和有效波高時, 能與浮標所測結果較為一致且精度較高。交叉譜方法相對于自相關函數方法能去除海浪傳播的180°模糊, 因此交叉譜方法更有優勢。

3) 互相關函數法得到的交叉譜, 在計算有效波高時會偏小, 這與互相關函數功率譜估計的方法有關。由于互相關法進行時, 兩個回波信號必須有差別,因此相關強度會變弱, 導致結果總能量偏小。海浪信息越規則, 這兩種方法結論就會越一致。

4) 在計算斜率方差時采用了5°~12°入射角范圍進行公式擬合, 根據本文結論可知, 低海況時定標與否對最后反演的二維海浪譜結果不會產生決定性影響, 但該方法對于更高海況適用性尚未可知。未來星載波譜儀只有靠多波束聯合才能實現。

致謝: 非常感謝北京遙測技術研究所提供的機載波譜儀數據和國家海洋局第一海洋研究所提供的浮標數據。

[1] Cavaleri L. Wave modeling: Where to go in the future[J]. Bulletin of the American Meteorological Society, 2006, 87(2): 207-214.

[2] Gonzalez F I, Beal R C, Brown W E, et al. Seasat synthetic aperture radar: Ocean wave detection capabilities[J]. Science, 1979, 204(4400): 1418-1421.

[3] Walsh E J, Hancock III D W, Hines D E, et al. Directional wave spectra measured with the surface contour radar[J]. Journal of Physical Oceanography, 1985, 15(5): 566-592.

[4] Jackson F C, Walton W T, Peng C Y. A comparison of in situ and airborne radar observations of ocean wave directionality[J]. Journal of Geophysical Research: Oceans (1978–2012), 1985, 90(C1): 1005-1018.

[5] Jackson F C, Walton W T, Baker P L. Aircraft and satellite measurement of ocean wave directional spectra using scanning-beam microwave radars[J]. Journal of Geophysical Research: Oceans, 1985, 90(C1): 987-1004.

[6] Hauser D, Caudal G C, Rijckenberg G J, et al. RESSAC: A new airborne FM/CW radar ocean wave spectrometer[J]. Geoscience and Remote Sensing, 1992, 30(5): 981-995.

[7] Danièle H, Podvin T, Dechambre M, et al. Polarimetric measurements over the sea-surface with the airborne STORM radar in the context of the geophysical validation of the ENVISAT ASAR[J/OL]. [2015-09-17]. http: //earth.esa.int/workshops/polinsar2003/participants/ hauser5/POLINSAR-VALPARESO.pdf, 2003.

[8] Plant W J, Keller W C, Hayes K. Simultaneous measurement of ocean winds and waves with an airborne coherent real aperture radar[J]. Journal of Atmospheric and Oceanic Technology, 2005, 22(7): 832-846.

[9] Pettersson H, Graber H C, Hauser D, et al. Directional wave measurements from three wave sensors during the FETCH experiment[J]. Journal of Geophysical Research:Oceans (1978–2012), 2003, 108(C3): 209.

[10] Hauser D, Soussi E, Thouvenot E, et al. SWIMSAT: A real-aperture radar to measure directional spectra of ocean waves from space-main characteristics and performance simulation[J]. Journal of Atmospheric and Oceanic Technology, 2001, 18(3): 421-437.

[11] Li Xiuzhong, He Yijun, Zhang Biao, et al. The construction of a three-dimensional antenna gain matrix and its impact on retrieving sea surface mean square slope based on aircraft wave spectrometer data[J]. Journal of Atmospheric and Oceanic Technology, 2016, 33(4): 847-856.

[12] Barrick D E. Rough surface scattering based on the specular point theory[J]. Antennas and Propagation, 1968, 16(4): 449-454.

[13] Valenzuela G R. Theories for the interaction of electromagnetic and oceanic waves-A review[J]. Boundary-Layer Meteorology, 1978, 13(1-4): 61-85.

[14] Jackson F C, Walton W T, Hines D E, et al. Sea surface mean square slope from Ku-band backscatter data[J]. Journal of Geophysical Research: Oceans (1978–2012), 1992, 97(C7): 11411-11427.

[15] Hesany V, Plant W J, Keller W C. The normalized radar cross section of the sea at 10 incidence[J]. Geoscience and Remote Sensing, 2000, 38(1): 64-72.

[16] Freilich M H, Vanhoff B A. The relationship between winds, surface roughness, radar backscatter at low incidence angles from TRMM Precipitation Radar measurements[J]. J Atmos Ocean Technol, 2003, 20(4): 549-562.

[17] Engen G, Johnsen H. SAR-ocean wave inversion using image cross spectra[J]. Geoscience and Remote Sensing, 1995, 33(4): 1047-1056.

Received:Dec. 24, 2015

Retrieval method of an ocean wave spectrum using an airborne spectrometer and performing the validation

LI Xiu-zhong1, HE Yi-jun1, MENG Jun-min2, ZHANG Zhen-hua3
(1. Nanjing University of Science Information and Technology, Nanjing 210044, China; 2. First Institute of Oceanography, State Oceanic Administration, Qingdao 266061, China; 3. Beijing Research Institute of Telemetry, Beijing 100076, China)

spectrometer; power spectrum estimation; two-dimensional ocean wave spectra; 180° ambiguity

The ocean wave retrieval method is designed on the basis of the signal formation principle of an ocean wave spectrometer. By using the spectra estimation methods via the auto-correlation and cross-correlation functions, the two-dimensional ocean wave spectrum is obtained. Finally, after comparing the spectrum received from the buoy and that retrieved from the spectrometer, the effectiveness of detecting a two-dimensional spectrum from an airborne spectrometer is evaluated. We observe that in the environment of flight, the results of the methods using the auto-correlation and cross-correlation functions for retrieving ocean wave spectrum are consistent with that obtained from the buoy. From the cross spectrum, the ambiguity of 180° is excluded, although the significant wave height is smaller than that from auto-correlation method. When the sea-slope variance is calculated, the radar backscattering coefficients of the incidence angles at 5°–12° are fitted. Therefore, calibration of the radar backscattering coefficients is not required. Moreover, the future spaceborne spectrometer will be able to attain calibration of the radar backscattering coefficients using multibeam joints.

P731.22; TP732.1

A

1000-3096(2016)12-0123-08

10.11759/hykx20151224004

(本文編輯: 李曉燕)

2015-12-24;

2016-03-19

國家重點研發計劃子課題(2016YFC1401005); 江蘇省自然科學基金項目(BK2011008); 國家自然科學基金(41476158); 江蘇省高等教育優勢學科項目

[Foundation: Research and Development of the National Key Program Corpus, No.2016YFC1401005; Natural Science Foundation of Jiangsu Province, No.BK2011008; National Science Foundation of China, No.41476158; The Preponderant Discipline Project of High Education in Jiangsu Province]

李秀仲(1985-), 山東濟寧人, 博士研究生, 主要從事海洋微波遙感研究, E-mail: qdlixiuzhong@163.com; 何宜軍, 通信作者,教授, 博士生導師, E-mail: yjhe@nuist.edu.cn