陳冠宇，艾未華，程玉鑫，戈書(shū)睿，袁凌峰

(1. 國(guó)防科技大學(xué)氣象海洋學(xué)院，江蘇 南京 211101； 2. 中國(guó)華陰兵器試驗(yàn)中心，陜西 華陰 714200； 3. 海軍海洋水文氣象中心，北京 100071)

基于星載SAR數(shù)據(jù)和模式資料的海面風(fēng)場(chǎng)變分融合方法研究

陳冠宇1，艾未華1，程玉鑫2，戈書(shū)睿1，袁凌峰3

(1. 國(guó)防科技大學(xué)氣象海洋學(xué)院，江蘇 南京 211101； 2. 中國(guó)華陰兵器試驗(yàn)中心，陜西 華陰 714200； 3. 海軍海洋水文氣象中心，北京 100071)

為實(shí)現(xiàn)合成孔徑雷達(dá)數(shù)據(jù)與數(shù)值預(yù)報(bào)模式資料融合，提高海面風(fēng)場(chǎng)精度和業(yè)務(wù)化運(yùn)用水平，提出了一種基于星載SAR數(shù)據(jù)與模式資料的變分融合方法。其研究思路是采用二維連續(xù)小波變換提取SAR圖像中高精度風(fēng)條紋風(fēng)向，結(jié)合地球物理模型函數(shù)求解海面風(fēng)場(chǎng)的經(jīng)向分量和緯向分量，然后采用Kriging插值方法將數(shù)值預(yù)報(bào)模式風(fēng)速插值到SAR海面風(fēng)場(chǎng)覆蓋區(qū)域，得到SAR風(fēng)速觀測(cè)算子，由此構(gòu)建SAR風(fēng)場(chǎng)與模式風(fēng)場(chǎng)融合的代價(jià)函數(shù)，并采用變分方法求解分析風(fēng)場(chǎng)，最終得到融合后的海面風(fēng)場(chǎng)結(jié)果。仿真分析結(jié)果表明，變分融合后的海面風(fēng)速和風(fēng)向結(jié)果更接近于理想值，尤其在SAR海面風(fēng)場(chǎng)覆蓋區(qū)域更為明顯。選取ENVISAT/ASAR資料和與其時(shí)空匹配的歐洲中期天氣預(yù)報(bào)中心模式風(fēng)場(chǎng)資料開(kāi)展實(shí)例驗(yàn)證，結(jié)果表明融合后的海面風(fēng)場(chǎng)結(jié)果比模式風(fēng)場(chǎng)更加接近于浮標(biāo)觀測(cè)結(jié)果。

合成孔徑雷達(dá)； 小波變換； 變分融合方法； 海面風(fēng)場(chǎng)

引言

海面風(fēng)場(chǎng)是海洋上層運(yùn)動(dòng)的主要?jiǎng)恿?lái)源，它與海洋中絕大多數(shù)的物理過(guò)程密切相關(guān)，是研究海洋動(dòng)力環(huán)境的重要參數(shù)，在海洋災(zāi)害監(jiān)測(cè)、海洋資源的開(kāi)發(fā)利用等方面亦具有重要意義[1]，海面風(fēng)場(chǎng)的探測(cè)已成為當(dāng)前海洋學(xué)研究中的一個(gè)熱點(diǎn)。利用衛(wèi)星遙感手段進(jìn)行海洋觀測(cè)的技術(shù)已經(jīng)趨于成熟，尤其是海面風(fēng)場(chǎng)的觀測(cè)[2]。星載合成孔徑雷達(dá)(synthetic aperture radar，SAR)具有全天時(shí)、全天候、高分辨率(數(shù)十米至數(shù)米)等特點(diǎn)，尤其適用于海岸帶和島嶼區(qū)的觀測(cè)，可以彌補(bǔ)散射計(jì)、微波輻射計(jì)等載荷在海面風(fēng)場(chǎng)探測(cè)方面的不足[3-4]。相對(duì)于其他探測(cè)海面風(fēng)場(chǎng)的衛(wèi)星載荷，SAR具備反演獲取高精度海面風(fēng)場(chǎng)的特點(diǎn)，擁有廣闊的應(yīng)用前景[5]。與微波散射計(jì)不同，星載合成孔徑雷達(dá)僅能單視角方向?qū)Ｓ^測(cè)，因此海面風(fēng)場(chǎng)反演需先獲得風(fēng)向再反演風(fēng)速，其風(fēng)向獲取依賴(lài)于圖像風(fēng)條紋或外部風(fēng)向，其中外部風(fēng)向的時(shí)空分辨率不易與SAR風(fēng)場(chǎng)匹配，而基于風(fēng)條紋的風(fēng)向反演方法精度高，但統(tǒng)計(jì)表明僅有44%左右的SAR圖像包含風(fēng)條紋信息[6-7]；全極化SAR能夠獨(dú)立反演海面風(fēng)場(chǎng)，但對(duì)交叉極化弱信號(hào)的絕對(duì)輻射定標(biāo)要求較高，在具體工程實(shí)現(xiàn)方面存在較大困難[8]，限制了SAR海面風(fēng)場(chǎng)反演的業(yè)務(wù)化應(yīng)用。

許多研究證明多源數(shù)據(jù)的融合有利于提高氣象、海洋等產(chǎn)品的探測(cè)精度，探索將SAR探測(cè)資料反演的海面風(fēng)場(chǎng)數(shù)據(jù)與其他風(fēng)場(chǎng)資料進(jìn)行融合，有望提高風(fēng)速和風(fēng)向精度[9]。目前，已有學(xué)者對(duì)微波散射計(jì)、微波輻射計(jì)和微波高度計(jì)海面風(fēng)場(chǎng)資料進(jìn)行了融合研究，但針對(duì)SAR數(shù)據(jù)與其他風(fēng)場(chǎng)資料融合的相關(guān)研究較少。Anctil et al.[10]對(duì)高度計(jì)風(fēng)速與大西洋船測(cè)風(fēng)場(chǎng)進(jìn)行了資料融合，Chin et al.[11]采用結(jié)合多分辨率分析的插值方法將ERS-1衛(wèi)星散射計(jì)風(fēng)場(chǎng)數(shù)據(jù)與歐洲中期天氣預(yù)報(bào)中心模式風(fēng)場(chǎng)數(shù)據(jù)進(jìn)行融合，Atlas et al.[12]通過(guò)變分方法融合輻射計(jì)資料SSM/I海面風(fēng)速、常規(guī)觀測(cè)和歐洲中期天氣預(yù)報(bào)中心分析值，生成了變分分析風(fēng)場(chǎng)，美國(guó)國(guó)家大氣研究中心(NCAR)數(shù)據(jù)支持部門(mén)(DSS)開(kāi)發(fā)了QuikSCAT散射計(jì)觀測(cè)值和NCEP再分析資料的融合產(chǎn)品Q(chēng)SCAT-NCEP混合風(fēng)場(chǎng)資料[13]，但針對(duì)SAR數(shù)據(jù)與其他風(fēng)場(chǎng)資料融合的相關(guān)研究較少。

本文開(kāi)展基于星載SAR數(shù)據(jù)與模式資料的海面風(fēng)場(chǎng)變分融合方法研究，利用小波分析的時(shí)-頻局部化特性，反演出風(fēng)條紋區(qū)域高精度的風(fēng)向和風(fēng)速，并通過(guò)變分方法與數(shù)值預(yù)報(bào)模式風(fēng)場(chǎng)進(jìn)行融合，獲得分析風(fēng)場(chǎng)。分別仿真研究模式風(fēng)場(chǎng)存在誤差和SAR后向散射截面存在誤差兩種情況下的結(jié)果，對(duì)比融合得到的分析風(fēng)場(chǎng)與理想風(fēng)場(chǎng)的誤差，結(jié)果表明兩種情況下均是變分融合的結(jié)果更接近于理想風(fēng)場(chǎng)，尤其在風(fēng)條紋區(qū)域更為明顯；選取ENVISAT/ASAR數(shù)據(jù)與歐洲中期天氣預(yù)報(bào)中心模式風(fēng)場(chǎng)資料進(jìn)行海面風(fēng)場(chǎng)的變分融合試驗(yàn)，結(jié)果表明融合后的海面風(fēng)場(chǎng)結(jié)果較模式風(fēng)場(chǎng)更加接近浮標(biāo)觀測(cè)結(jié)果。仿真分析與試驗(yàn)結(jié)果均表明，文章提出的基于星載SAR數(shù)據(jù)與模式資料的海面風(fēng)場(chǎng)變分融合方法能夠有效提高海面模式風(fēng)場(chǎng)的精度。

1 變分融合方法

1.1 基于小波的SAR海面風(fēng)場(chǎng)反演方法

小波分析因其出色的時(shí)-頻局部化特性被譽(yù)為信號(hào)處理中的“數(shù)學(xué)顯微鏡”[14]。通過(guò)理論分析和實(shí)際的反演試驗(yàn)，Mexican-hat小波在SAR海面風(fēng)場(chǎng)反演方面取得了較好的效果[15-16]。二維Mexican-hat小波具有較好的時(shí)域和頻域局部化能力及信號(hào)能量集中特性，對(duì)于信噪比較高、風(fēng)條紋清晰的SAR圖像，該方法能夠獲取較高的海面風(fēng)向反演精度[17]。在空間-頻率域，二維Mexican-hat小波變換可表示為：

(1)

式中:k表示二維空間-頻率域的變量；·表示向量?jī)?nèi)積。

Mexican-hat小波法反演風(fēng)向的過(guò)程是首先對(duì)SAR圖像進(jìn)行二維連續(xù)Mexican-hat小波變換，獲取不同尺度下的小波能量譜圖像，提取風(fēng)條紋信息。然后對(duì)能量譜圖像進(jìn)行二維FFT變換，計(jì)算SAR圖像中風(fēng)條紋的波數(shù)譜。最后將二維波數(shù)譜峰值的連線做垂線，對(duì)其進(jìn)行風(fēng)向去模糊后得到海面風(fēng)向[18]。針對(duì)Mexican-hat小波法反演風(fēng)向存在的風(fēng)向180°模糊問(wèn)題，采用數(shù)值預(yù)報(bào)模式風(fēng)向?qū)ζ溥M(jìn)行去除。文獻(xiàn)[19]的研究表明，在二維連續(xù)小波分解尺度為8時(shí)，海面風(fēng)向反演獲取的效果較好。

Hersbach et al.[20]經(jīng)大量試驗(yàn)建立了SAR后向散射截面與海面風(fēng)場(chǎng)之間關(guān)系的CMOD地球物理模型：

σ0(V,θ,φ)=
B0(V,θ)(1+B1(V,θ)cosφ+B2(V,θ)cos2φ)α

(2)

其中B0、B1、B2是風(fēng)速V和入射角θ的函數(shù)，φ為風(fēng)向，α為常數(shù)。將反演出的風(fēng)向作為地球物理模型函數(shù)的相對(duì)風(fēng)向輸入，同時(shí)輸入定標(biāo)后的后向散射系數(shù)和入射角，通過(guò)迭代計(jì)算得到相應(yīng)的海面風(fēng)速，最終獲得SAR反演海面風(fēng)場(chǎng)信息。

1.2 變分融合

本文采用變分方法對(duì)SAR風(fēng)場(chǎng)與模式風(fēng)場(chǎng)進(jìn)行融合。首先，將SAR海面風(fēng)場(chǎng)矢量(觀測(cè)風(fēng)場(chǎng))與模式風(fēng)場(chǎng)矢量(背景風(fēng)場(chǎng))分解為經(jīng)向分量U和緯向分量V(即水平與垂直分量)，分別利用Kriging插值法將分解的背景風(fēng)場(chǎng)插值到SAR風(fēng)場(chǎng)所在區(qū)域上得到SAR風(fēng)速觀測(cè)算子H[21]，然后利用變分融合方法求解分析風(fēng)場(chǎng)。

X=(V1,V2,…，VN)T

則：

觀測(cè)算子H可表示為：

這里μ為極小化處理時(shí)的拉格朗日乘子，

c(xi,xj)為背景風(fēng)場(chǎng)中網(wǎng)格點(diǎn)xi與xj之間的協(xié)方差函數(shù)，xk為觀測(cè)風(fēng)場(chǎng)區(qū)域上的插值點(diǎn)。本文采用指數(shù)型插值函數(shù)：

其中，L為插值點(diǎn)與觀測(cè)點(diǎn)之間的距離，C0、C1、a均為常數(shù)。

觀測(cè)風(fēng)場(chǎng)與背景風(fēng)場(chǎng)采用二維變分(2D-Var)的方法進(jìn)行融合計(jì)算，其在數(shù)學(xué)上可以描述為使如下定義的目標(biāo)函數(shù)最?。?/p>

J=Jq+Jm

(3)

其中，觀測(cè)風(fēng)場(chǎng)的目標(biāo)函數(shù)Jq和背景風(fēng)場(chǎng)的目標(biāo)函數(shù)Jm定義如下：

(4)

(5)

其中，U和V分別是融合得到的經(jīng)向風(fēng)速和緯向風(fēng)速(即分析風(fēng)速)，q和m代表觀測(cè)風(fēng)場(chǎng)和背景風(fēng)場(chǎng)，Q和M分別是觀測(cè)風(fēng)場(chǎng)和背景風(fēng)場(chǎng)的誤差協(xié)方差矩陣。將(4)和(5)式帶入(3)式可得：

(6)

(7)

利用最優(yōu)化算法求解(6)和(7)式便可得出融合后的分析風(fēng)速。

1.3 海面風(fēng)場(chǎng)變分融合流程

本文提出的基于星載SAR數(shù)據(jù)與模式資料的海面風(fēng)場(chǎng)變分融合方法，首先在SAR影像線性紋理特征明顯的區(qū)域(風(fēng)條紋區(qū)域)采用二維連續(xù)小波變換提取風(fēng)條紋信息，但是存在180°的風(fēng)向模糊問(wèn)題，對(duì)此可以利用模式資料的風(fēng)向?qū)ζ溥M(jìn)行風(fēng)向去模糊得到SAR海面風(fēng)向[22]。然后將反演出的風(fēng)向作為相對(duì)風(fēng)向輸入地球物理模型函數(shù)進(jìn)一步求解SAR海面風(fēng)速，并將風(fēng)速分解為經(jīng)向分量和緯向分量。利用Kriging插值法將分解的背景風(fēng)場(chǎng)插值到SAR風(fēng)場(chǎng)所在區(qū)域上得到SAR風(fēng)速觀測(cè)算子，構(gòu)建SAR風(fēng)速與背景風(fēng)速融合的代價(jià)泛函，最后通過(guò)變分方法求解分析風(fēng)速得出融合結(jié)果。具體反演流程如圖1所示。

圖1 海面風(fēng)場(chǎng)的變分融合流程圖Fig.1 Flow chart of variational fusion of sea surface wind field

2 仿真研究及誤差分析

為深入地分析融合效果，考查誤差對(duì)融合結(jié)果的影響，引入隨機(jī)擾動(dòng)，分別針對(duì)模式風(fēng)場(chǎng)存在誤差、SAR風(fēng)場(chǎng)無(wú)誤差和SAR后向散射截面存在誤差、模式風(fēng)場(chǎng)無(wú)誤差兩種情況，對(duì)變分融合結(jié)果進(jìn)行敏感性分析，將變分融合后的分析風(fēng)場(chǎng)與理想風(fēng)場(chǎng)進(jìn)行比對(duì)，驗(yàn)證變分方法融合SAR風(fēng)場(chǎng)和模式風(fēng)場(chǎng)的有效性。仿真模擬的SAR風(fēng)場(chǎng)區(qū)域及理想風(fēng)場(chǎng)探測(cè)點(diǎn)位置如圖2所示，圖中黑色箭頭代表理想風(fēng)場(chǎng)，白點(diǎn)代表SAR風(fēng)場(chǎng)坐標(biāo)位置，理想風(fēng)場(chǎng)經(jīng)向分量ui=8 m·s-1、緯向分量vi=6 m·s-1，分辨率為20 km。左側(cè)灰色區(qū)域?yàn)锳區(qū)，表示SAR風(fēng)場(chǎng)覆蓋區(qū)域。利用Kriging插值方法將理想風(fēng)場(chǎng)插值到A區(qū)上構(gòu)造理想SAR風(fēng)場(chǎng)，分辨率為30 km?；疑珔^(qū)域右側(cè)的區(qū)域?yàn)锽區(qū)(參考區(qū)域)，將參考區(qū)上的融合結(jié)果與A區(qū)的融合結(jié)果進(jìn)行比對(duì)，以考察SAR風(fēng)場(chǎng)覆蓋區(qū)域和非SAR風(fēng)場(chǎng)覆蓋區(qū)域的融合效果，這主要是考慮到并非SAR圖像上都存在風(fēng)條紋，在實(shí)際應(yīng)用中SAR圖像中僅部分區(qū)域包含風(fēng)條紋。整體風(fēng)場(chǎng)的經(jīng)向距離取值范圍為[-100 km，100 km]，緯向距離取值范圍為[-100 km，100 km]。

圖2 模擬的背景風(fēng)速探測(cè)點(diǎn)位置和SAR風(fēng)場(chǎng)反演區(qū)域Fig.2 Location of detection point of simulated background wind speed and retrieval area of SAR wind field

2.1 背景模式風(fēng)場(chǎng)存在誤差時(shí)風(fēng)場(chǎng)融合結(jié)果的敏感性分析

當(dāng)理想風(fēng)場(chǎng)的風(fēng)速為10 m·s-1左右時(shí)，背景風(fēng)場(chǎng)的風(fēng)速均方根誤差一般不超過(guò)1 m·s-1。由于選取的理想風(fēng)場(chǎng)的風(fēng)速均為中等風(fēng)速(10 m·s-1)，所以對(duì)該理想風(fēng)場(chǎng)添加ε=1 m·s-1的隨機(jī)誤差(記ε為最大擾動(dòng)振幅)，進(jìn)而得到擾動(dòng)風(fēng)場(chǎng)。將擾動(dòng)風(fēng)場(chǎng)作為背景場(chǎng)，然后引入模擬的SAR風(fēng)場(chǎng)，經(jīng)過(guò)變分融合得到分析風(fēng)場(chǎng)，具體研究流程如圖3所示。

圖3 背景風(fēng)場(chǎng)存在誤差時(shí)的研究流程圖Fig.3 Research flow chart of background wind field with error

利用Kriging插值計(jì)算得到的理想SAR風(fēng)速為10 m·s-1，背景風(fēng)速最小為9.15 m·s-1，最大為10.94 m·s-1。通過(guò)比對(duì)融合后的分析風(fēng)場(chǎng)與模式風(fēng)場(chǎng)的差值(圖4)可知，SAR風(fēng)場(chǎng)覆蓋的A區(qū)域融合后分析風(fēng)場(chǎng)的變化較為明顯，而融合后B區(qū)域的分析風(fēng)場(chǎng)相較于背景模式風(fēng)場(chǎng)變化不大。經(jīng)計(jì)算得出的A區(qū)域的分析風(fēng)速在9.63～10.57 m·s-1之間振蕩，擾動(dòng)明顯得到了抑制，即采用變分方法融合SAR風(fēng)場(chǎng)與模式風(fēng)場(chǎng)是有效的。

為考察本文方法在不同擾動(dòng)強(qiáng)度、不同區(qū)域內(nèi)的融合效果，采用文獻(xiàn)[23]的方法，針對(duì)融合結(jié)果對(duì)整體以及A區(qū)和B區(qū)的風(fēng)場(chǎng)進(jìn)行敏感性試驗(yàn)分析。取ε=0.5、1.0、2.0 m·s-1，分別將背景以及分析風(fēng)場(chǎng)和理想風(fēng)場(chǎng)進(jìn)行對(duì)比，計(jì)算得到各自對(duì)應(yīng)的均方根誤差，結(jié)果如表1-3所示，其中ebu為背景風(fēng)場(chǎng)在經(jīng)向分量上的均方根誤差，ebv為背景風(fēng)場(chǎng)在緯向分量上的均方根誤差， eau為分析風(fēng)場(chǎng)在經(jīng)向分量上的均方根誤差，eav為分析風(fēng)場(chǎng)在緯向分量上的均方根誤差。綜合分析表1-3可知，隨著ε的減小，分析風(fēng)場(chǎng)的均方根誤差隨之減小。當(dāng)ε=0.5、1.0、2.0 m·s-1時(shí)，分析風(fēng)場(chǎng)在經(jīng)向分量和緯向分量方向上的均方根誤差均比背景風(fēng)場(chǎng)在經(jīng)向分量和緯向分量方向上的均方根誤差小，說(shuō)明基于星載SAR數(shù)據(jù)與模式資料的變分融合方法對(duì)模式數(shù)據(jù)有積極效應(yīng)。由于B區(qū)背景風(fēng)場(chǎng)距離SAR風(fēng)場(chǎng)區(qū)域較遠(yuǎn)，不能使SAR風(fēng)速對(duì)背景風(fēng)場(chǎng)產(chǎn)生有效影響，因此，對(duì)于均方根誤差，背景和分析風(fēng)場(chǎng)的結(jié)果幾乎相同，結(jié)果如表3所示；A區(qū)中SAR風(fēng)場(chǎng)覆蓋的背景風(fēng)場(chǎng)受到較大的影響，從均方根誤差來(lái)看，分析風(fēng)場(chǎng)均明顯小于背景風(fēng)場(chǎng)(表2)，融合效果明顯，這些與圖4所得結(jié)果一致。統(tǒng)計(jì)表明：在A區(qū)，利用SAR風(fēng)場(chǎng)與背景風(fēng)場(chǎng)進(jìn)行變分融合后，背景風(fēng)場(chǎng)的均方根誤差能夠降低17%左右。整體風(fēng)場(chǎng)在經(jīng)向和緯向分量上的均方根誤差介于A區(qū)與B區(qū)風(fēng)場(chǎng)在經(jīng)向和緯向分量上的均方根誤差之間。由此可知，采用本文的方法進(jìn)行風(fēng)場(chǎng)融合是有效的。

圖4 融合后的分析風(fēng)場(chǎng)與模式風(fēng)場(chǎng)的誤差(a.經(jīng)向分量U，b.緯向分量V)Fig.4 Error of analysis wind field and model wind field after fusion(a.meridional component U,b.zonal component V)

表1不同振幅擾動(dòng)下整體區(qū)域背景風(fēng)場(chǎng)和分析風(fēng)場(chǎng)相對(duì)理想風(fēng)場(chǎng)的均方根誤差

Table 1 Root mean square error of overall regional background wind field and analysis wind field compared to ideal wind field under different amplitude disturbances m·s-1

表2不同振幅擾動(dòng)下A區(qū)背景風(fēng)場(chǎng)和分析風(fēng)場(chǎng)相對(duì)理想風(fēng)場(chǎng)的均方根誤差

Table 2 Root mean square error of background wind field in area A and analysis wind field compared to ideal wind field under different amplitude disturbances m·s-1

表3不同振幅擾動(dòng)下B區(qū)背景風(fēng)場(chǎng)和分析風(fēng)場(chǎng)相對(duì)理想風(fēng)場(chǎng)的均方根誤差

Table 3 Root mean square error of background wind field in area B and analysis wind field compared to ideal wind field under different amplitude disturbances m·s-1

2.2 SAR后向散射截面存在誤差時(shí)風(fēng)場(chǎng)融合結(jié)果的敏感性分析

下面針對(duì)SAR后向散射截面存在誤差的情況，對(duì)變分融合結(jié)果進(jìn)行敏感性分析，具體研究流程如圖5所示。對(duì)理想的SAR后向散射截面σi添加隨機(jī)擾動(dòng)得到加擾動(dòng)的SAR后向散射截面σ0。Freeman初步給出了SAR反演海面風(fēng)場(chǎng)、海浪等海洋環(huán)境要素時(shí)的輻射定標(biāo)精度，認(rèn)為SAR探測(cè)風(fēng)速的誤差小于20%時(shí)，其絕對(duì)定標(biāo)精度、相對(duì)定標(biāo)精度均為1 dB[24]。當(dāng)風(fēng)速為10 m·s-1時(shí)，假定入射角為30°，相對(duì)風(fēng)向?yàn)?0°，代入地球物理模型求得σi=-11.8 dB，對(duì)σi添加最大為1 dB的擾動(dòng)，反演的SAR風(fēng)速與背景風(fēng)速之差最大值可達(dá)2.3 m·s-1，統(tǒng)計(jì)得到SAR風(fēng)速的均方根誤差為1.06 m·s-1。假設(shè)對(duì)背景風(fēng)場(chǎng)不作擾動(dòng)，即(ub,vb)=(ui,vi)，在這種情況下利用變分方法進(jìn)行風(fēng)場(chǎng)融合后，通過(guò)比對(duì)分析風(fēng)場(chǎng)與理想風(fēng)場(chǎng)的誤差(圖6)可知，分析風(fēng)場(chǎng)由于計(jì)入SAR后向散射截面隨機(jī)擾動(dòng)而有小幅度擾動(dòng)，但不明顯，其原因是在變分方法中分析風(fēng)場(chǎng)中的信息只有15%左右來(lái)自觀測(cè)場(chǎng)的貢獻(xiàn)，其余85%左右來(lái)自背景場(chǎng)的貢獻(xiàn)[25]。SAR后向散射截面受擾動(dòng)后，針對(duì)整體風(fēng)場(chǎng)以及A區(qū)風(fēng)場(chǎng)和B區(qū)風(fēng)場(chǎng)，計(jì)算各個(gè)風(fēng)場(chǎng)經(jīng)向分量和緯向分量相對(duì)于理想風(fēng)場(chǎng)經(jīng)向分量和緯向分量的均方根誤差，所得結(jié)果如表4所示。由表4可知，當(dāng)SAR后向散射截面存在誤差時(shí)，覆蓋區(qū)域的均方根誤差稍大一些，遠(yuǎn)離覆蓋區(qū)域的均方根誤差稍小一些，整體的風(fēng)速誤差僅在厘米每秒級(jí)。

圖5 SAR后向散射截面存在誤差時(shí)的研究流程圖Fig.5 Research flow chart of SAR backward scattering cross section with error

圖6 融合后的分析風(fēng)場(chǎng)與理想風(fēng)場(chǎng)的誤差(a.經(jīng)向分量U，b.緯向分量V)Fig.6 Error of analysis wind field and ideal wind field after fusion(a.meridional component U,b.zonal component V)

表4σi受擾動(dòng)下背景風(fēng)場(chǎng)和分析風(fēng)場(chǎng)相對(duì)理想風(fēng)場(chǎng)的均方根誤差

Table 4 Root mean square error of background wind field and analysis wind field compared to ideal wind field under σi disturbed m·s-1

3 實(shí)例研究及驗(yàn)證

3.1 試驗(yàn)數(shù)據(jù)

本文采用的SAR數(shù)據(jù)是兩軌ENVISAT/ASAR 的寬幅成像模式垂直極化數(shù)據(jù)；采用的背景風(fēng)場(chǎng)數(shù)據(jù)是與SAR數(shù)據(jù)時(shí)空匹配的歐洲中期天氣預(yù)報(bào)中心再分析資料；浮標(biāo)數(shù)據(jù)為NDBC(National Data Buoy Center)數(shù)據(jù)，其觀測(cè)高度為海表面5 m，通過(guò)風(fēng)場(chǎng)剖面能量法則將浮標(biāo)風(fēng)速轉(zhuǎn)化為海表面10 m等效風(fēng)速。

3.2 結(jié)果分析

首先在SAR影像紋理特征明顯區(qū)域采用二維小波法反演海面風(fēng)向，利用CMOD_IFR2模型求解海面風(fēng)速并分解為經(jīng)向分量和緯向分量。以圖8中第一軌SAR圖像浮白框區(qū)域SAR子圖像(圖7a)為例，說(shuō)明基于Mexican-hat小波變換的海面風(fēng)場(chǎng)反演過(guò)程，先對(duì)SAR圖像進(jìn)行二維連續(xù)小波變換獲得小波能量譜，再對(duì)其進(jìn)行二維FFT變換計(jì)算風(fēng)條紋的波數(shù)譜，波數(shù)譜峰值連線(實(shí)線)的垂線(虛線)就是所求風(fēng)向(圖7b)。但圖7b中所得風(fēng)向存在180°模糊，將模式風(fēng)向作為比對(duì)去模糊風(fēng)向。根據(jù)模式風(fēng)向判斷，2011年5月24日08時(shí)左右探測(cè)區(qū)域?yàn)闁|南風(fēng)，去除西北向的模糊風(fēng)向，可得風(fēng)向的反演結(jié)果為101.1°，結(jié)合地球物理模型函數(shù)求解出海面風(fēng)速為12.2 m·s-1。

圖7 小波法反演的SAR子圖像(a)和風(fēng)向(b)Fig.7 SAR sub image(a) and wind direction(b) inversed by the wavelet method

圖8為第一軌SAR數(shù)據(jù)反演風(fēng)場(chǎng)和2011年5月24日06時(shí)歐洲中期天氣預(yù)報(bào)中心再分析風(fēng)場(chǎng)資料的探測(cè)點(diǎn)位置及浮標(biāo)位置圖。圖中箭頭方向代表風(fēng)向，箭頭長(zhǎng)度代表風(fēng)速。SAR資料探測(cè)時(shí)間為2011年5月24日07：52。白色箭頭代表SAR風(fēng)場(chǎng)，經(jīng)緯網(wǎng)格點(diǎn)分辨率為0.3°×0.3°；黑色箭頭代表模式風(fēng)場(chǎng)，經(jīng)緯網(wǎng)格點(diǎn)分辨率為0.25°×0.25°；白點(diǎn)代表浮標(biāo)位置，浮標(biāo)號(hào)為46078，觀測(cè)時(shí)間為07:50。計(jì)算風(fēng)條紋區(qū)域SAR反演風(fēng)向和風(fēng)速并分解為經(jīng)向分量和緯向分量，以此作為觀測(cè)風(fēng)場(chǎng)，并以模式風(fēng)場(chǎng)作為背景場(chǎng)，經(jīng)過(guò)變分融合得到分析風(fēng)場(chǎng)，如圖9所示，圖中箭頭方向表示風(fēng)向，箭頭長(zhǎng)度表示風(fēng)速，圖中左上角為風(fēng)速標(biāo)尺，標(biāo)尺長(zhǎng)度代表風(fēng)速為10 m·s-1。

圖10中SAR數(shù)據(jù)探測(cè)時(shí)間為2011年8月23日20：55，與之時(shí)空匹配的2011年8月23日18時(shí)歐洲中期天氣預(yù)報(bào)中心再分析風(fēng)場(chǎng)資料的探測(cè)點(diǎn)位置及浮標(biāo)位置如圖所示，浮標(biāo)號(hào)為46075，觀測(cè)時(shí)間為20:50。圖11為經(jīng)過(guò)變分融合得到的分析風(fēng)場(chǎng)。

圖8 SAR反演風(fēng)場(chǎng)和數(shù)值預(yù)報(bào)模式風(fēng)場(chǎng)Fig.8 Wind field retrieval results of SAR and wind field of numerical forecasting model

圖10 SAR反演風(fēng)場(chǎng)和數(shù)值預(yù)報(bào)模式風(fēng)場(chǎng)Fig.10 Wind field retrieval results of SAR and wind field of numerical forecasting model

圖9 變分融合后的分析風(fēng)場(chǎng)Fig.9 Analysis wind field after variational fusion

圖11 變分融合后的分析風(fēng)場(chǎng)Fig.11 Analysis wind field after variational fusion

由圖11可知，變分融合后的分析風(fēng)場(chǎng)與模式風(fēng)場(chǎng)的趨勢(shì)基本保持一致，將浮標(biāo)點(diǎn)附近的SAR反演風(fēng)速、分析風(fēng)速和模式風(fēng)速同浮標(biāo)探測(cè)值進(jìn)行比對(duì)分析，結(jié)果如表5所示。

表5浮標(biāo)點(diǎn)附近的SAR反演風(fēng)速、分析風(fēng)速和模式風(fēng)速

Table 5 SAR retrieval wind speed， analysis wind speed and model wind speed near the buoy m·s-1

由表5可以看出，SAR反演風(fēng)速與浮標(biāo)風(fēng)速最為接近，模式風(fēng)速與之相比誤差較大。這一方面是由于基于小波的SAR反演風(fēng)場(chǎng)具有精度高、誤差小的特點(diǎn)，另一方面SAR圖像探測(cè)時(shí)間與浮標(biāo)探測(cè)時(shí)間幾乎相同，而模式風(fēng)場(chǎng)的探測(cè)時(shí)間與浮標(biāo)觀測(cè)時(shí)間相差較多，其中浮標(biāo)2處的模式風(fēng)速與浮標(biāo)觀測(cè)值相比誤差很大，這可能與模式風(fēng)速探測(cè)時(shí)間同浮標(biāo)觀測(cè)時(shí)間相差近3 h有關(guān)。通過(guò)變分融合得到的分析風(fēng)速較模式風(fēng)速更加靠近浮標(biāo)風(fēng)速，其中浮標(biāo)1處的分析風(fēng)速與浮標(biāo)觀測(cè)值十分接近，浮標(biāo)2處的分析風(fēng)速相較于模式風(fēng)速也得到改善，證明本文提出的風(fēng)場(chǎng)融合方法有效。

4 小結(jié)

本文提出的基于星載SAR數(shù)據(jù)與模式資料的海面風(fēng)場(chǎng)變分融合方法，主要是利用SAR反演海面風(fēng)場(chǎng)精度高、區(qū)域廣的特點(diǎn)，通過(guò)小波變換法提取較高精度的風(fēng)條紋風(fēng)向，結(jié)合地球物理模型函數(shù)反演風(fēng)速并分解為經(jīng)向分量和緯向分量，以數(shù)值預(yù)報(bào)模式資料作為背景風(fēng)場(chǎng)，構(gòu)建變分融合的代價(jià)泛函，最后利用變分方法求解分析風(fēng)場(chǎng)，實(shí)現(xiàn)合成孔徑雷達(dá)數(shù)據(jù)反演海面風(fēng)場(chǎng)與數(shù)值預(yù)報(bào)模式風(fēng)場(chǎng)資料的融合，以期提高區(qū)域性海面模式風(fēng)場(chǎng)的整體精度。仿真研究表明，在背景模式風(fēng)場(chǎng)存在誤差、SAR反演風(fēng)場(chǎng)無(wú)誤差的情況下，SAR風(fēng)場(chǎng)對(duì)背景風(fēng)場(chǎng)的調(diào)整有積極作用，能夠使背景風(fēng)場(chǎng)的均方根誤差降低17%左右，遠(yuǎn)離SAR風(fēng)場(chǎng)覆蓋區(qū)域效果次之。在SAR后向散射截面存在擾動(dòng)、背景模式風(fēng)場(chǎng)無(wú)誤差的情況下，融合后的分析風(fēng)場(chǎng)受擾動(dòng)的幅度較小，說(shuō)明在一定的輻射定標(biāo)精度內(nèi)SAR后向散射截面誤差對(duì)風(fēng)場(chǎng)融合結(jié)果的影響相對(duì)較小，利用本文方法融合的風(fēng)場(chǎng)具有較強(qiáng)的抗噪性。但該方法需要SAR圖像中存在一定區(qū)域的風(fēng)條紋信息，以便提取較高精度的風(fēng)條紋風(fēng)向及風(fēng)速信息與背景風(fēng)場(chǎng)進(jìn)行融合，進(jìn)而提升海面風(fēng)場(chǎng)融合的整體精度。

[1] Shimada T, Sawada M, Sha W, et al. Low-level easterly winds blowing through the Tsugaru Strait, Japan. Part I: Case study and statistical characteristics based on observations[J]. Mon Wea Rev, 2010, 138(10):3806-3821.

[2] Sheng Z. The estimation of lower refractivity uncertainty from radar sea clutter using the Bayesian-MCMC method[J]. Chin Phys B, 2013,22(2):580-585.

[3] Liang Z. New GMF+RAIN model based on rain rate and application in typhoon wind retrieval[J]. Acta Phys Sinica, 2010, 59(10):7478-7490.

[4] Martin S. An introduction to ocean remote sensing[M]. Cambridge: Cambridge University Press, 2014: 401-435.

[5] Portabella M, Stoffelen A, Johannessen J A. Toward an optimal inversion method for synthetic aperture radar wind retrieval[J]. J Geophys Res, 2002, 107(C8):1-13.

[6] Apel J R. An improved model of the ocean surface wave vector spectrum and its effects on radar backscatter[J]. J Geophys Res, 1994, 99(C8):16269-16291.

[7] Attema E P W, Duchossois G, Kohlhammer G. ERS-1/2 SAR land applications: Overview and main results[C]// Geoscience and remote sensing symposium proceedings, 1998. IGARSS’98. 1998 IEEE International. 1998:1796-1798.

[8] Zhang B, Perrie W, Vachon P W, et al. Ocean vector winds retrieval from C-Band fully polarimetric SAR measurements[J]. IEEE Trans Geosci Remote Sens, 2012, 50(11):4252-4261.

[9] Grody N C, Weng F. Microwave emission and scattering from deserts: Theory compared with satellite measurements[J]. IEEE Trans Geosci Remote Sens, 2008, 46(2):361-375.

[10] Anctil F, Donelan M A, Forristall G Z, et al. Deep-water field evaluation of the NDBC-SWADE 3-m discus directional buoy[J]. J Atmos Oceanic Technol, 1993, 10(1):97.

[11] Chin T M, Milliff R F, Large W G. Basin-scale, high-wavenumber sea surface wind fields from a multiresolution analysis of scatterometer data[J]. J Atmos Oceanic Technol, 1998, 15(3):741.

[12] Atlas R, Hoffman R N, Ardizzone J, et al. A new cross-calibrated, multi-satellite ocean surface wind product[C]// Geoscience and Remote Sensing Symposium, 2008. IGARSS 2008. IEEE International. 2008:I-106-I-109.

[13] Grody N C, Weng F. Microwave emission and scattering from deserts: Theory compared with satellite measurements[J]. IEEE Trans Geosci Remote Sens, 2008, 46(2):361-375.

[14] 張日偉, 嚴(yán)衛(wèi), 艾未華,等. 基于Gabor小波變換的機(jī)載SAR海面風(fēng)向反演方法[J]. 微波學(xué)報(bào), 2011, 27(5):79-83.

[15] 孔毅, 趙現(xiàn)斌, 艾未華,等. 基于墨西哥帽小波變換的機(jī)載SAR海面風(fēng)場(chǎng)反演[J]. 解放軍理工大學(xué)學(xué)報(bào)(自然科學(xué)版), 2011, 12(3):301-306.

[16] Leite G C, Ushizima D M, Medeiros F N, et al. Wavelet analysis for wind fields estimation.[J]. Sensors, 2010, 10(6):5994-6016.

[17] Ai W H. Ocean surface wind direction retrieval from multi-polarization airborne SAR based on wavelet[J]. Acta Phys Sinica, 2012, 61(14):321-328.

[18] 程玉鑫, 艾未華, 孔毅,等. 基于影像紋理特征和外部風(fēng)向的星載SAR海面風(fēng)場(chǎng)反演研究[J]. 海洋科學(xué), 2015, 39(12):157-164.

[19] 艾未華, 孔毅, 趙現(xiàn)斌. 基于小波的多極化機(jī)載合成孔徑雷達(dá)海面風(fēng)向反演[J]. 物理學(xué)報(bào), 2012, 61(14):465-473.

[20] Hersbach H, Stoffelen A, Haan S D. An improved C-band scatterometer ocean geophysical model function: CMOD5[J]. J Geophys Res Oceans, 2007, 112(C3):225-237.

[21] Jiang Z H. A new approach to adjusting sea surface wind using altimeter wind data by variational regularization method[J]. Acta Phys Sinica, 2010, 59(12):8968-8977.

[22] 張毅, 蔣興偉, 林明森,等. 合成孔徑雷達(dá)海面風(fēng)場(chǎng)反演研究進(jìn)展[J]. 遙感技術(shù)與應(yīng)用, 2010, 25(3):423-429.

[23] 姜祝輝, 黃思訓(xùn), 杜華棟,等. 利用變分結(jié)合正則化方法對(duì)高度計(jì)風(fēng)速資料調(diào)整海面風(fēng)場(chǎng)的研究[J]. 物理學(xué)報(bào), 2010, 59(12):8968-8977.

[24] Freeman A. SAR calibration: An overview[J]. IEEE Trans Geosci Remote Sens, 1992, 30(6):1107-1121.

[25] Comstock K K, Wood R, Yuter S E, et al. Reflectivity and rain rate in and below drizzling stratocumulus[J]. Quart J Roy Meteor Soc, 2004, 130(603):2891-2918.

Researchonoceanwindvariationalfusionmethodbasedonspace-borneSARandmodeldata

CHEN Guanyu1, AI Weihua1, CHENG Yuxin2, GE Shurui1, YUAN Lingfeng3

(1.CollegeofMeteorologyandOceanography,NationalUniversityofDefenseTechnology,Nanjing211101,China; 2.HuayinWeaponTestCenterofChina,Huayin714200,China; 3.NavyHydrometeorologicalCenter,Beijing100071,China)

Based on the space-borne SAR(synthetic aperture radar) data and model data, a variational fusion method is presented to implement the fusion of SAR data and numerical weather prediction model data, and to improve the accuracy of sea surface wind field and operational level. First, two-dimensional continuous wavelet transform is applied to extract wind direction of high precision wind stripe in SAR image, and geophysical model function is used to calculate meridional and zonal components of sea surface wind field. Then the Kriging interpolation method is used to interpolate the wind speed of numerical weather prediction model to the coverage area of SAR sea surface wind field, resulting in that the SAR wind speed observation operator is obtained and the cost function of SAR wind field combined with model wind field is established. The variational method is adopted to solve the analysis wind field. Finally, the sea surface wind field results after fusion are obtained. Simulation results show that through the variational fusion, the sea surface wind speed and wind direction results are closer to the ideal values, especially in the coverage area of SAR wind field. The ENVISAT/ASAR sounding data and the spatial temporally matched model wind data of ECMWF(European Center for Mediumrange Weather Forecasts) are selected to carry out the case verification, and the results show that the sea surface wind field after fusion is closer to the buoy observed result than the model wind field, so it is confirmed that the variational fusion method is effective.

synthetic aperture radar; wavelet transform; variational fusion method; sea surface wind field

陳冠宇,艾未華,程玉鑫,等.基于星載SAR數(shù)據(jù)和模式資料的海面風(fēng)場(chǎng)變分融合方法研究[J].海洋氣象學(xué)報(bào),2017,37(4):65-74.

Chen Guanyu,Ai Weihua,Cheng Yuxin,et al.Research on ocean wind variational fusion method based on space-borne SAR and model data[J].Journal of Marine Meteorology,2017,37(4):65-74.

10.19513/j.cnki.issn2096-3599.2017.04.008.(in Chinese)

P407.2

A

2096-3599(2017)04-0065-10

10.19513/j.cnki.issn2096-3599.2017.04.008

2017-07-19；

2017-08-16

國(guó)家自然科學(xué)基金項(xiàng)目(41475019，41375029)

陳冠宇(1994—)，男，碩士研究生，主要從事微波海洋遙感研究，guanyumail@126.com。

艾未華(1979—)，男，副教授，主要從事氣象海洋雷達(dá)技術(shù)研究，awhzjax@126.com。