基于變分多參數(shù)正則化方法融合衛(wèi)星散射計資料

張凱峰,鄧婉月,王挺,王慧鵬,項杰,4*,宋清濤,劉春霞

基于變分多參數(shù)正則化方法融合衛(wèi)星散射計資料

張凱峰1,鄧婉月2,王挺1,王慧鵬3,項杰1,4*,宋清濤5,劉春霞6

(1.解放軍理工大學(xué) 氣象海洋學(xué)院,江蘇 南京211101;2.新疆維吾爾自治區(qū)環(huán)境監(jiān)測總站,新疆 烏魯木齊830011;3.75839部隊,廣東 廣州510510;4.南京大學(xué) 中尺度災(zāi)害性天氣教育部重點實驗室,江蘇 南京210023;5.國家衛(wèi)星海洋應(yīng)用中心,北京100081;6.中國氣象局廣州熱帶海洋氣象研究所,廣東 廣州510080)

基于常規(guī)三維變分同化(3DVAR)思想和反問題中的正則化技術(shù),提出了適用于風(fēng)場融合的帶正則化約束項的3DVAR方法,在南海海域開展數(shù)據(jù)融合試驗,同時采用模型函數(shù)方法確定合理的正則化參數(shù),針對一次臺風(fēng)個例進行了Quik SCAT散射計海面風(fēng)場數(shù)據(jù)和華南中尺度模式海面風(fēng)場數(shù)據(jù)的融合試驗,結(jié)果表明采用帶正則化約束的3DVAR融合方法,明顯消除了常規(guī)3DVAR方法融合風(fēng)場時帶來的虛假信息,融合后分析風(fēng)場以及渦度場和散度場分布均勻,結(jié)構(gòu)清晰,氣旋中心顯著,且分析場中觀測起主導(dǎo)作用;采用信號自由度(DFS)方法對融合方法進行定量評估,發(fā)現(xiàn)相對常規(guī)3DVAR方法,帶正則化約束的3DVAR融合系統(tǒng)中觀測數(shù)據(jù)提供的DFS較多,同時提高了觀測場對分析場的影響;基于獨立觀測資料對融合結(jié)果進行檢驗發(fā)現(xiàn)相對華南中尺度模式和常規(guī)3DVAR方法的統(tǒng)計結(jié)果,帶正則化約束的3DVAR方法得到的風(fēng)場具有最小的均方根誤差和最大的相關(guān)系數(shù)。

3DVAR;正則化;Quik SCAT;信號自由度;模型函數(shù)

1 引言

在數(shù)值天氣預(yù)報(NWP)中,初始場是最重要的前提條件,業(yè)務(wù)NWP中心通常利用資料同化方法來獲得一個最優(yōu)的初始場,常用的方法有三維變分同化(3DVAR)、四維變分同化(4DVAR)等。觀測資料的多少是數(shù)值預(yù)報關(guān)注的重點。近年來,衛(wèi)星遙感技術(shù)的發(fā)展為全球的觀測提供了一種有效的技術(shù)手段,高時空分辨率的觀測數(shù)據(jù)豐富了原有的觀測網(wǎng)絡(luò),特別是資料稀疏的海洋區(qū)域,不過,由于極軌衛(wèi)星觀測時間一般與區(qū)域模式初始場時間不一致,因此如何合理有效地同化衛(wèi)星觀測資料對于提高數(shù)值模式初始場的精度以及數(shù)值預(yù)報的準確率十分重要,尤其是臺風(fēng)的分析預(yù)報。已有研究表明同化衛(wèi)星遙感資料能夠改進臺風(fēng)初始場精度,提高臺風(fēng)強度以及路徑的預(yù)報準確率[1-4]。

海面風(fēng)場是區(qū)域和全球海洋環(huán)流的主要動力,調(diào)節(jié)海氣之間水汽、熱量以及物質(zhì)的交換,影響著海洋天氣形勢的發(fā)展演變過程,因此獲取高時空分辨率的海面風(fēng)場具有重要的研究和應(yīng)用價值。目前獲取海面風(fēng)場產(chǎn)品主要依賴衛(wèi)星遙感,主要傳感器有微波輻射計、散射計以及高度計,其中只有散射計可以同時提供風(fēng)速、風(fēng)向,觀測海面風(fēng)場的散射計有SASS、ESCAT、NSCAT、QuikSCAT 以及近年來陸續(xù)發(fā)射的ASCAT、HY-2A、OSCAT、Rapid SCAT 等,為全球海洋環(huán)境的檢測提供了豐富的數(shù)據(jù),其中QuikSCAT得到了廣泛的應(yīng)用,不僅有原始的軌道刈幅數(shù)據(jù),還有經(jīng)過處理的規(guī)則網(wǎng)格化數(shù)據(jù),方便了數(shù)據(jù)的利用。衛(wèi)星觀測數(shù)據(jù)的日益增加使得風(fēng)場數(shù)據(jù)的融合變得越來越重要,許多學(xué)者對于海面風(fēng)場融合方法也做了相關(guān)的研究,主要包括統(tǒng)計方法和變分方法,基于統(tǒng)計插值方法開展風(fēng)場融合的研究較多,Tang和Liu[5]采用連續(xù)訂正方法對ERS-1散射計風(fēng)場和數(shù)值天氣預(yù)報模式風(fēng)場進行融合,得到了分辨率為1°,12 h一次的風(fēng)場產(chǎn)品;Royle等[6]基于分層Bayes方法融合了NSCAT散射計和NWP風(fēng)場,Perrie等[7]采用最優(yōu)插值法(OI)融合了ERS-2、NSCAT、ERS-2高度計和數(shù)值模式風(fēng)場,得到了西北大西洋天氣尺度高分辨率的網(wǎng)格化海面風(fēng)場數(shù)據(jù)。Atlas[8]在1993年就利用變分法融合了SSM/I風(fēng)速資料和ECMWF風(fēng)場數(shù)據(jù),為之后CCMP(Cross-Calibrated Multi-platfor m)數(shù)據(jù)集的發(fā)布提供了基礎(chǔ),Chao等[9]在加州中部海岸開展了Quik SCAT與區(qū)域數(shù)值模式風(fēng)場數(shù)據(jù)的融合,通過與獨立觀測資料進行對比得出融合風(fēng)場具有最小的誤差和最大的相關(guān)性。目前國外已經(jīng)發(fā)布了相關(guān)的海面風(fēng)場融合產(chǎn)品,主要包括CCMP數(shù)據(jù)集和QSCAT-NCEP混合風(fēng)場資料,為研究和應(yīng)用提供了長時間序列的海面風(fēng)場數(shù)據(jù)。國內(nèi)對于風(fēng)場融合也做了一些相關(guān)的研究,為了獲取較高時空分辨率的風(fēng)場,凌征等[10]采用Cressman方法將Quik SCAT風(fēng)場與氣象站風(fēng)場數(shù)據(jù)融合,為研究我國近海環(huán)流的特征提供了參考,蔣興偉和宋清濤[11]利用低通濾波方法融合了3種散射計資料,包括NSCAT、Quik SCAT和Seawinds,通過濾波的方式對散射計刈幅掃描數(shù)據(jù)進行處理,有效消除了衛(wèi)星數(shù)據(jù)反演過程中帶來的奇異值,并且保留了豐富的中小尺度信息,獲得了逐月的風(fēng)場數(shù)據(jù),空間分辨率為0.25°×0.25°。在我國HY-2A衛(wèi)星發(fā)射后,國內(nèi)學(xué)者對其搭載的散射計風(fēng)場產(chǎn)品也做了相關(guān)的融合研究[12-13]。而在南海海域,由于其地形復(fù)雜,島嶼較多,近海風(fēng)場具有十分明顯的季風(fēng)特征和地形特點,國外開發(fā)的全球海面風(fēng)場數(shù)據(jù)產(chǎn)品未必適合我國臨近海域的研究和應(yīng)用需要,因此,充分利用多源衛(wèi)星遙感風(fēng)場數(shù)據(jù),開展我國近海海域海面風(fēng)場數(shù)據(jù)融合研究十分必要[14]。

由于風(fēng)場融合問題本質(zhì)上屬于反問題[15-17],而反問題通常是不適定的,因此,開展風(fēng)場數(shù)據(jù)融合的研究既具有理論意義,也具有應(yīng)用價值,求解不適定問題的一個有效途徑是正則化方法[16]。已有研究表明,在大氣海洋領(lǐng)域,將正則化思想引入到變分同化系統(tǒng)里可以很好的解決不適定問題,通過利用L曲線或者模型函數(shù)方法選取最優(yōu)的單正則化參數(shù),可以提高計算效率,得到了較好的結(jié)果[18-20]。本文在常規(guī)3DVAR的基礎(chǔ)上,利用反問題中的正則化技術(shù)開展對南海海面風(fēng)場數(shù)據(jù)融合的研究,將Quik SCAT散射計海面風(fēng)場數(shù)據(jù)與華南中尺度模式(Guang Zhou Mesoscale Model,GZMM)輸出的海面風(fēng)場數(shù)據(jù)進行融合,并且采用模型函數(shù)方法確定合理的雙正則化參數(shù),以得到最佳的融合風(fēng)場,然后利用信號自由度方法計算帶正則化約束的3DVAR融合系統(tǒng)中的信息容量,探討分析場對觀測場的敏感性,最后通過西沙自動氣象站獨立觀測資料對融合結(jié)果進行驗證,評估融合效果。

2 研究區(qū)域以及數(shù)據(jù)

2.1 研究區(qū)域

南海海域及其周邊存在較多的島嶼,由于散射計風(fēng)場反演數(shù)據(jù)的精度受陸地的干擾較大,為了減少這種影響,選擇遠離陸地的區(qū)域,范圍為10.12°～17.92°N,110.12°～117.92°E(圖1中的紅色方框區(qū)域)。在該區(qū)域內(nèi)的永興島上,中國科學(xué)院南海海洋研究所于2008年建成了西沙海洋觀測站,其中包括無人自動氣象觀測站(AWS,圖1中的五角星位置),AWS的海面風(fēng)場數(shù)據(jù)可用來進行精度驗證。

2.2 數(shù)據(jù)介紹

2.2.1 獨立觀測數(shù)據(jù)

獨立觀測數(shù)據(jù)來自西沙永興島上的無人自動氣象站(AWS)提供的距離海面10 m處的風(fēng)場,時間范圍為2008年8-12月,時間分辨率為2 min,包括風(fēng)速、風(fēng)向,其中風(fēng)向采用16方位角表示,方位序號0～15表示N-NNW,2009年6月之后以度數(shù)來表示,所有數(shù)據(jù)均經(jīng)過質(zhì)量控制,剔除明顯錯誤。

2.2.2 Quik SCAT散射計資料

圖1 研究區(qū)域(紅色方框)和AWS位置(五角星)Fig.1 Research area(red box)and AWS location(five-pointed star)

QuikSCAT為太陽同步軌道海洋科學(xué)衛(wèi)星,于1999年6月19日發(fā)射,攜帶了測量海面風(fēng)場的散射計,運行周期為101 min,刈幅寬度為1 800 k m,工作頻率為13.4 GHz,測量距離海面10 m高度的風(fēng)場,空間分辨率為25 k m,每天可觀測的范圍覆蓋了全球90%的海洋。美國國家航空航天局噴氣推進動力實驗室(NASA/JPL)發(fā)布了兩種Quik SCAT L2B風(fēng)場產(chǎn)品,一種是標準的風(fēng)場產(chǎn)品,另一種是改進了模糊剔除技術(shù)得到的風(fēng)場產(chǎn)品[21]。Ebuchi等[22]利用太平洋以及大西洋浮標資料對Quik SCAT風(fēng)場進行驗證,得出散射計風(fēng)場與浮標一致性高,風(fēng)速、風(fēng)向均方根誤差分別為1.00 m/s和22.4°,精度較高。本文利用改進的模糊剔除技術(shù)的Quik SCAT風(fēng)場產(chǎn)品進行南海海面風(fēng)場融合研究。

2.2.3 模式風(fēng)場數(shù)據(jù)

本文采用的模式數(shù)據(jù)為華南中尺度模式輸出的海面風(fēng)場數(shù)據(jù)。GZMM由中國氣象局廣州熱帶海洋氣象研究所開發(fā),預(yù)報模式基于GRAPES_TMM(GRAPES Tropical Mesoscale Model)的非靜力平衡模式,同時加入了適應(yīng)本地的物理過程、地形處理及適合南海的邊界層技術(shù)方案、熱帶對流技術(shù)方案、水物質(zhì)技術(shù)方案。水平方向采用等距的經(jīng)-緯格點和Arakawa-C格式,模式的預(yù)報范圍:10.0°～35.0°N,95.0°～130.0°E。本文所用數(shù)據(jù)范圍為2008年9-12月,模式水平分辨率為0.12°,每日兩次(00、12 UTC)預(yù)報,分別發(fā)布84 h時效的預(yù)報,海面10 m風(fēng)場是其業(yè)務(wù)產(chǎn)品之一,GZMM并沒有同化Quik SCAT風(fēng)場數(shù)據(jù),保證了數(shù)據(jù)之間的獨立性。

2.2.4 數(shù)據(jù)質(zhì)量控制

在數(shù)據(jù)融合之前,必須對Quik SCAT數(shù)據(jù)進行質(zhì)量控制,以防引進不良觀測信息影響融合結(jié)果。首先,根據(jù)wvc_qual_flag標識,剔除覆蓋陸地以及冰面的數(shù)據(jù);Quik SCAT風(fēng)速設(shè)計范圍為3～30 m/s,剔除該范圍之外的數(shù)據(jù)。第二,對于降雨標識數(shù)據(jù),在分析臺風(fēng)時,由于臺風(fēng)周圍降雨數(shù)據(jù)較多,剔除過多觀測數(shù)據(jù)難以辨別臺風(fēng)中心的位置,并且會影響融合后臺風(fēng)位置以及強度的變化,在融合時保留這些數(shù)據(jù)[5],非極端天氣條件考慮降雨數(shù)據(jù)的剔除。第三,由于極軌衛(wèi)星接收資料的連續(xù)性,為了充分利用觀測資料,融合時間窗口選擇3 h。

3 研究方法

3.1 方法介紹

3.1.1 常規(guī)3DVAR方法

3DVAR的目的是融合觀測數(shù)據(jù)與背景場?；谝韵录僭O(shè):(1)觀測及背景誤差服從均值為0的正態(tài)分布;(2)觀測誤差和背景誤差不相關(guān),則3DVAR問題歸結(jié)為如下常規(guī)的目標函數(shù)的極小化(稱為常規(guī)的3DVAR):

其中,x和xb分別為分析場和背景場;y是觀測場;H是觀測算子;R和B分別為觀測誤差協(xié)方差矩陣和背景誤差協(xié)方差矩陣; ( )T和 ( )-1分別表示矩陣的轉(zhuǎn)置和逆。

3.1.2 帶正則化約束的3DVAR方法

同化問題本身屬于反問題,而B的估計不可避免地含有誤差,因此可以引入正則化參數(shù)α來修正B;同時,對于風(fēng)場融合來說,引入的觀測風(fēng)場信息可能在融合后的風(fēng)場中產(chǎn)生虛假的空間梯度,因此,引入先驗的平滑項作為正則化項進行約束[9,23],本文選取風(fēng)場的渦度梯度和散度梯度作為正則化項加入到目標函數(shù)中,表達式為:

式中,d A=d x d y,β為正則化參數(shù),則常規(guī)的目標函數(shù)(1)變?yōu)?

故對于本文要研究的海面風(fēng)場融合來說,

式中,U和V分別為風(fēng)矢量的緯向和經(jīng)向分量;下標q和m分別代表Quik SCAT和模式;α、β為正則化參數(shù),ζ和D分別為海面風(fēng)場的渦度和散度。為了區(qū)別,J～的極小化問題稱為帶正則化約束的3DVAR。在本文中,H取為線性插值算子(雙線性插值)。對Jr進行離散處理可得:

式中,W1=L1TL1,W2=L2TL2,L1和L2為二階偏導(dǎo)離散后的系數(shù)矩陣(詳細推導(dǎo)見附錄)。本文采用有限內(nèi)存的擬牛頓法(L-BFGS)實現(xiàn)目標函數(shù)的最小化。通過計算可得目標函數(shù)(4)的梯度:?J～=,其中,

3.1.3 誤差協(xié)方差矩陣的統(tǒng)計

觀測誤差協(xié)方差矩陣及背景誤差協(xié)方差矩陣的估計是資料同化的一個重要環(huán)節(jié)。本文中的B矩陣采用Parish和Derber[24]提出的NMC方法進行計算:

對于區(qū)域模式來說,T1和T2一般為12 h和24 h[25]。

觀測誤差假定是不相關(guān)的,則R為一對角矩陣,觀測誤差方差為對角線元素,則R=σ2I,根據(jù)Quik-SCAT與AWS對比結(jié)果(圖2)可知,U、V的誤差方差分別為:

圖2 Quik SCAT與AWS風(fēng)速U、V分量散點圖Fig.2 The U and V component scatter plot bet ween QuikSCAT and AWS

3.1.4 正則化參數(shù)的選取

選取合適的正則化參數(shù)是求解反問題時采用正則化方法的關(guān)鍵,常用的選取方法主要有Morozov偏差原理[26]或吸收Morozov偏差原理、L曲線準則[27]、廣義交叉檢驗準則[28]等,其中Morozov偏差原理或吸收Morozov偏差原理是一種十分有效的參數(shù)選取方法,然而利用數(shù)值方法直接求解關(guān)于正則化參數(shù)的非線性方程時往往計算量很大,而且只具有局部收斂性,結(jié)果并不令人滿意。為了有效求解這樣的非線性方程,Kunisch和Zou[29]提出了模型函數(shù)方法,一些學(xué)者將該方法進行了進一步的研究和發(fā)展[30-31],并且將該方法應(yīng)用到多個正則化參數(shù)的選擇。下面針對本文要研究的海面風(fēng)場融合給出一般性的描述。從目標函數(shù)(3)可以得到:

則吸收Morozov偏差方程為:

式中,γ＞1、κ＞1為吸收系數(shù);c≥1為常數(shù);σ為誤差水平;式(9)等價于:

上述方程中Fα,β()沒有顯示的表達式,且關(guān)于α,β是非線性的,因此,采用線性模型函數(shù)mkα,β()=Tk+Ckα+Dkβ近似表示Fα,β(),Tk、Ck、Dk為待定參數(shù),代入式(10),則

確定雙正則化參數(shù)的流程圖如圖3所示。

圖3 模型函數(shù)計算雙參數(shù)流程Fig.3 Two parameter calculation flow based on model function

3.1.5 信號自由度的計算

信號自由度(DFS)表示一個反演系統(tǒng)(如3DVAR系統(tǒng))能從觀測數(shù)據(jù)中獲得的獨立狀態(tài)變量的數(shù)量,同時,它也可以表示同化系統(tǒng)中分析場對觀測數(shù)據(jù)的敏感性,可以用于評估同化或者反演系統(tǒng)[32]。對于常規(guī)3DVAR方法來說,最優(yōu)分析場xa可以表示為:

式中,xb是背景場,y是觀測場,H是觀測算子,K是最優(yōu)卡爾曼增益矩陣,

H是H在xb附近線性化的切線性算子H=H'(xb)。

同化系統(tǒng)式(12)的DFS定義為:

式中,tr{ · }表示 { · }的跡。

而對于帶正則化約束的3DVAR模型來說,最優(yōu)分析xa表示為:

當(dāng)H為線性算子即H=H時,帶正則化約束的3DVAR模型DFS為:

由于本文研究構(gòu)造的融合系統(tǒng)中,矩陣維數(shù)較小,可以直接進行DFS的計算。

3.2 方案設(shè)計流程

方案實施流程圖如圖4所示。首先運用觀測場和背景場分別統(tǒng)計觀測誤差協(xié)方差矩陣和背景誤差協(xié)方差矩陣,并確定觀測算子;然后采用擬牛頓法分別求解得到常規(guī)3DVAR方法和帶正則化約束的3DVAR方法的最優(yōu)分析場,借助模型函數(shù)方法確定合理的雙正則化參數(shù),計算信號自由度,對融合系統(tǒng)進行評估;最后利用自動氣象站獨立觀測數(shù)據(jù)檢驗融合效果。為了方便,文中將常規(guī)3DVAR方法和帶正則化約束的3DVAR分別簡稱為常規(guī)方法和正則化方法。

圖4 方案設(shè)計流程圖Fig.4 The flow chart for program design

4 數(shù)值試驗及結(jié)果

4.1 臺風(fēng)個例分析

以2008年第16號臺風(fēng)“米克拉”為例進行融合試驗,分析正則化方法的有效性。9月28日,在南海中部海面生成的熱帶低壓,于29日2時正式編報并于8時加強為熱帶風(fēng)暴,中心風(fēng)力8級,風(fēng)速達18 m/s。以2008年9月28日12時12 h預(yù)報場為背景場,9月28日Quik SCAT L2B升軌資料為觀測場,融合結(jié)果為29日00時分析場,參照中國氣象局熱帶氣旋資料中心最佳路徑數(shù)據(jù)集(www.typhoon.gov.cn),29日00時熱帶風(fēng)暴中心位置位于16.3°N,111.2°E。圖5為觀測場和背景場的環(huán)流分布形勢,圖5a為Quik SCAT軌道刈幅觀測數(shù)據(jù),可以看出,圖中有一個明顯的氣旋環(huán)流,但氣旋中心不太明顯,中心位置大約在16°N,111°E附近,氣旋周圍風(fēng)速較大,最大風(fēng)速在氣旋的西南側(cè),達到28 m/s,位于背景場(圖5b)的西北側(cè);圖5b為背景場的環(huán)流分布形勢,圖中氣旋性環(huán)流明顯且中心顯著,中心位于15.7°N,112.8°E附近,最大風(fēng)速達到16 m/s,位于氣旋中心的東南側(cè)。

圖5 2008年9月28日QuikSCAT觀測場(a)與背景場(b)環(huán)流形勢Fig.5 Quik SCAT observation field(a)and background field(b)circulation distribution in September 28,2008

為了顯示正則化約束的作用,做如下試驗:β=0和β≠0分別研究單參數(shù)以及雙參數(shù)對融合結(jié)果的影響。令β=0分析調(diào)節(jié)α對融合結(jié)果的影響,設(shè)計了3組試驗,即α=1,α=0.1,α=0.01分析融合結(jié)果,如圖6所示。

定義均方根誤差(RMSE):

式中,y為Quik SCAT觀測場;H為觀測算子;當(dāng)x=xb時,即為觀測場與背景場之間風(fēng)速的RMSE;x=xa時為觀測場與分析場之間風(fēng)速的RMSE;N為觀測個數(shù)。通過計算融合衛(wèi)星散射計前觀測場與背景場(O-B)以及融合后觀測場與分析場(OA)之間的差異,分析觀測信息對分析場的影響。

表1為觀測與背景和觀測與分析之間的均方根偏差RMSE,以及相應(yīng)的Jo與Jb的數(shù)值,后者反映了分析場與觀測場以及背景場偏差的大小。當(dāng)Jo＞Jb時,分析場主要由背景場確定(即分析場接近背景場),當(dāng)Jo＜Jb時,分析場主要由觀測場確定。

當(dāng)α=1時,即為常規(guī)方法,所得到的分析場如圖6a所示。從圖6a可以看出,融合Quik SCAT散射計觀測資料之后,分析場的熱帶氣旋中心位置及強度有了明顯的改變(相對背景場來說),氣旋環(huán)流中心接近于觀測,并且最大風(fēng)速變大,但是,分析場(融合后的風(fēng)場)分布不均勻,風(fēng)場渦度場和散度場分布散亂,不規(guī)則,存在許多異常的小范圍渦度低值區(qū)和散度高值區(qū),并沒有體現(xiàn)出氣旋中心,可見融合Quik SCAT風(fēng)場的同時,引入了虛假的觀測信息,增加了分析場的噪聲;另一方面,雖然O-A的RMSE為2.49 m/s小于O-B的RMSE值,為3.15 m/s,但是,Jo＞Jb(表1),這表明分析場接近于背景場,背景場起主要作用。隨著α減小,O-A的RMSE也逐漸減小,Jo與Jb的值變化顯著:Jo的值從大變小,Jb的值從小變大,并且,熱帶氣旋周圍風(fēng)速大值區(qū)域逐漸變大,渦度高值區(qū)和散度低值區(qū)變化顯著;當(dāng)α=0.01時,O-A的RMSE減小為1.98 m/s,Jo的數(shù)值為3 386.535明顯小于Jb的值(18 623.01),這表明分析場更接近于觀測場,觀測場起主要作用,但從渦度場和散度場的分布可以看出,通過調(diào)整正則化參數(shù)α,分析場的渦度場及散度場散亂現(xiàn)象沒有改變,虛假噪聲并沒有得 到有效的抑制,難以確定渦度、散度中心。

圖6 α=1(a,d,g),α=0.1(b,e,h)以及α=0.01(c,f,i)時分析風(fēng)場(a～c)、渦度場(d～f)和散度場(g～i)分布Fig.6 The distribution of wind field(a-c),vorticity(d-f)and divergence(g-i)whenα=1(a,d,g),α=0.1(b,e,h)andα=0.01(c,f,i)

表1 β=0時O-B和O-A的RMSE以及J o與J b的數(shù)值變化Tal.1 The RMSE of O-B and O-A as well as the value of J o and J b whenβ=0

為了克服這一現(xiàn)象,引入先驗的平滑約束項來調(diào)整風(fēng)場,即β≠0,通過調(diào)節(jié)目標函數(shù)J～中的正則化參數(shù)α、β,以得到最優(yōu)的融合風(fēng)場。利用模型函數(shù)方法可得合理的雙正則化參數(shù)分別為α=0.002,β=3.30,相應(yīng)的融合后風(fēng)場如圖7a所示。與常規(guī)方法得到的分析風(fēng)場相比,正則化方法得到的整個環(huán)流場分布均勻,渦度場(圖7b)和散度場(圖7c)變化顯著,明顯消除了圖6d～6i中渦度場和散度場出現(xiàn)的異常低值和高值區(qū)域,分布均勻,結(jié)構(gòu)清晰,可以辨認熱帶氣旋渦旋中心的大概位置,位于16.4°N,111.5°E附近,平滑約束項對風(fēng)場融合引進的虛假噪聲進行了有效的抑制,O-A的RMSE為1.79 m/s,分析場主要依賴于觀測場,達到了融合的目的。

試驗分析結(jié)果表明,不考慮平滑約束項即β=0時,通過調(diào)節(jié)單參數(shù)α改變分析場中觀測場和背景場的主導(dǎo)作用:α越小,分析場越接近于觀測場,觀測場起主要作用,背景場得到了有效地抑制,但分析風(fēng)場噪點較多,渦度場和散度場出現(xiàn)了異常的低值和高值區(qū)域,結(jié)構(gòu)異常散亂,熱帶氣旋中心也不顯著;在風(fēng)場融合系統(tǒng)中加入先驗的平滑約束后,通過確定合理的雙正則化參數(shù)對分析風(fēng)場進行調(diào)整,可以得出,正則化方法進一步改善了融合效果,有效消除了風(fēng)場融合時出現(xiàn)的噪點現(xiàn)象,明顯改善了分析風(fēng)場以及渦度和散度場的分布情況,風(fēng)場的渦度場和散度場分布更加均勻,結(jié)構(gòu)清晰,熱帶氣旋中心十分顯著,中心位置與實況更加接近,融合后氣旋的結(jié)構(gòu)和強度均優(yōu)于散射計風(fēng)場、模式風(fēng)場以及運用常規(guī)方法融合得到的結(jié)果,且觀測場與分析場風(fēng)速的均方根誤差最小,正則化方法有效融合了觀測信息。

圖7 正則化方法的分析風(fēng)場(a)、渦度場(b)和散度場(c)Fig.7 Analysis of wind field(a),vorticity field(b)and divergence field(c)of regularization method

4.2 敏感性分析

下面利用信號自由度(DFS)方法分別計算常規(guī)方法(α=1,β=0)和正則化方法(α=0.002、β=3.30)中U和V分量的DFS,考察分析場對觀測數(shù)據(jù)的敏感性(或觀測數(shù)據(jù)對分析場的影響),并且對比分析兩種方法之間的差異。Lupu等[33]指出,對于任意選定的資料,可以用觀測影響(OI)表述觀測數(shù)據(jù)對分析場的平均影響,即:

式中,p為觀測數(shù)量。這里針對本文的情形計算相應(yīng)的OI,計算結(jié)果見表2。

表2 信號自由度和觀測影響Tab.2 Signal freedomand observation influence

從表2可以看出,常規(guī)方法融合系統(tǒng)中,U風(fēng)速 和V風(fēng)速提供的DFS分別為37.284 2和34.754 8,觀測場對分析場的影響OI分別為3.25%和3.03%,觀測對分析的影響較小,而在正則化方法融合系統(tǒng)中,觀測提供的DFS增大,分析場中包含的觀測信息量明顯增加,U分量和V分量提供的DFS分別變?yōu)?5.562 7和85.082 9,觀測數(shù)據(jù)對分析場的影響相對常規(guī)方法也有了明顯的改變,其中U分量的觀測影響(OI)從3.25%增加到7.46%,V分量的OI從3.03%增加到7.42%,增加了大約4%,可見經(jīng)過正則化約束提高了融合系統(tǒng)的信息容量,增強了分析場對觀測數(shù)據(jù)的敏感性。

5 融合效果評估

利用常規(guī)方法以及正則化方法分別進行長時間的融合試驗,選取數(shù)據(jù)融合的時間范圍為2008年9月26日至11月25日共兩個月,將Quik SCAT散射計二級風(fēng)場的升軌和降軌數(shù)據(jù)(觀測場)分別與GZMM模式風(fēng)場數(shù)據(jù)(背景場)進行融合,利用公式(16)分別統(tǒng)計每天升軌和降軌的觀測與背景以及觀測與分析的均方根誤差,評估引入觀測信息對模式初始場的影響以及不同融合系統(tǒng)所得結(jié)果之間的差異,并且基于獨立的觀測數(shù)據(jù),對融合結(jié)果進行驗證,檢驗融合效果。

圖8為長時間序列的升軌和降軌融合結(jié)果中風(fēng)速的O-A和O-B的風(fēng)速RMSE變化趨勢,其中黑色虛線表示觀測場與背景場的RMSE,藍色虛線和紅色實線分別為觀測場與常規(guī)方法融合分析場以及觀測場與正則化方法分析場的RMSE,在一個成功的同化系統(tǒng)中,O-A的RMSE應(yīng)小于O-B。從圖8可以看出,數(shù)據(jù)融合前觀測場與背景場差異較大,RMSE在2～4 m/s之間波動,波動幅度較大,融合Quik SCAT散射計資料之后,RMSE變化顯著,O-A的RMSE明顯小于O-B,且大部分數(shù)值在2 m/s以下,表明融合后觀測場接近于分析場,觀測起主導(dǎo)作用,而利用正則化方法的融合效果更好,O-A的RMSE均小于常規(guī)方法的統(tǒng)計結(jié)果,升軌(圖8a)和降軌(圖8b)的結(jié)果一致。

基于獨立的西沙自動站觀測數(shù)據(jù),對融合效果進行評估,通過以下幾個統(tǒng)計量進行檢驗,統(tǒng)計結(jié)果見表3。定義A為對比數(shù)據(jù),B為參考數(shù)據(jù),N為匹配樣本量。

平均值(Mean):

圖8 升軌(a)和降軌(b)的O-A和O-B的RMSE變化序列Fig.8 The RMSEseries of O-A(a)and O-B(b)for ascend(a)and descend(b)respectively

平均偏差(Bias):

均方根誤差(RMSE):

相關(guān)系數(shù)(R):

從統(tǒng)計分析結(jié)果可以看出(表3),觀測誤差,即Quik SCAT風(fēng)速RMSE為1.68 m/s,與設(shè)計精度保持一致,U、V分量RMSE也均在2 m/s以內(nèi),與自動站一致性好,相關(guān)系數(shù)在0.92以上;背景場風(fēng)速相對自動站風(fēng)場偏差比較大,風(fēng)速RMSE為2.26 m/s,U、V分量均方根誤差分別為2.55 m/s和2.34 m/s,U分量相關(guān)性較差;利用常規(guī)的3DVAR方法融合后,分析場風(fēng)速以及風(fēng)速分量相對背景場統(tǒng)計結(jié)果有所改善,風(fēng)速RMSE為2.07 m/s,U分量相關(guān)性提高比較顯著,相關(guān)性達到了0.85;而正則化方法融合結(jié)果風(fēng)速以及風(fēng)速分量精度變化顯著,風(fēng)速RMSE為1.92 m/s,其中V分量精度的變化較大,RMSE減小到1.96 m/s,相關(guān)性達到0.94,與獨立觀測數(shù)據(jù)一致性高?？梢娙诤螿uikSCAT散射計觀測數(shù)據(jù)后,風(fēng)場的精度相對GZMM模式精度有了明顯的改善,其中采用正則化方法得到的風(fēng)場具有最小的均方根誤差和最大的相關(guān)系數(shù)。

表3 統(tǒng)計對比結(jié)果Tab.3 The results of statistical comparison

6 結(jié)論和討論

本文主要利用正則化方法在南海海域開展了衛(wèi)星散射計(Quik SCAT)和區(qū)域中尺度模式(GZMM)風(fēng)場資料的數(shù)據(jù)融合試驗,并采用模型函數(shù)方法選取合理的雙正則化參數(shù),針對一次臺風(fēng)個例以及長時間序列融合試驗驗證了正則化方法的有效性,并對融合結(jié)果進行了檢驗和評估,具體如下:

(1)個例融合結(jié)果發(fā)現(xiàn)常規(guī)方法得到的分析風(fēng)場中噪點較多,渦度場和散度場分布異常散亂,難以辨認氣旋中心,而采用正則化方法后,明顯消除了常規(guī)方法融合風(fēng)場時帶來的虛假信息,風(fēng)場以及渦度場和散度場中的噪點現(xiàn)象基本消失,分布均勻,結(jié)構(gòu)清晰,熱帶氣旋中心顯著,相對于模式與實況更加接近,且Quik-SCAT觀測場與融合后分析場的均方根誤差明顯小于其與背景場的統(tǒng)計結(jié)果,達到了融合的目的,同時正則化方法融合系統(tǒng)中觀測數(shù)據(jù)提供的DFS較多,觀測對分析的影響較強,經(jīng)過正則化約束后,觀測對分析的影響大約增加了4%,增強了分析對觀測的敏感性。

(2)采用西沙永興島自動站觀測資料對長時間序列的融合結(jié)果進行檢驗,得出GZMM模式U和V分量均方根誤差分別為2.55 m/s和2.34 m/s,相關(guān)性分別為0.78和0.89,常規(guī)方法U和V分量均方根誤差分別為2.28 m/s和2.22 m/s,相關(guān)性分別為0.85和0.91,正則化方法融合結(jié)果中分析場U和V分量均方根誤差分別為2.10 m/s和1.96 m/s,相關(guān)性分別為0.86和0.94,優(yōu)于GZMM模式和常規(guī)方法融合的統(tǒng)計結(jié)果,因此融合Quik SCAT觀測資料后,對背景場的影響比較明顯,U分量和V分量的精度和相關(guān)性都得到了改進,正則化方法得到的分析場具有最小的均方根誤差和最大的相關(guān)系數(shù)。

通過臺風(fēng)個例和長時間序列的融合試驗檢驗結(jié)果發(fā)現(xiàn)正則化方法是合理有效的,不過在臺風(fēng)個例融合試驗中為了體現(xiàn)正則化方法的優(yōu)勢,并沒有考慮降雨對散射計數(shù)據(jù)的影響,因此在后續(xù)的工作中需要改進質(zhì)量控制方案,對數(shù)據(jù)進行合理的處理,并進一步驗證正則化方法的普遍性,為開展我國周邊海域的風(fēng)場融合研究提供參考。

致謝:感謝國家科技基礎(chǔ)條件平臺——國家地球系統(tǒng)科學(xué)數(shù)據(jù)共享平臺——南海及其鄰近海區(qū)科學(xué)數(shù)據(jù)中心(http://ocean.geodata.cn)提供數(shù)據(jù)支撐,感謝中國氣象局廣州熱帶海洋氣象研究所提供的GZMM(Guang Zhou Mesoscale Model)模式資料,感謝美國NASA噴氣推進實驗室(JPL)提供的Quik SCAT風(fēng)場資料。

[1] Zou Xiaolei,Xiao Qingnong.Studies on theinitialization and si mulation of a mature hurricane using a variational bogus data assi milation scheme[J].Journal of the At mospheric Sciences,2000,57(6):836-860.

[2] Chen S,Vandenberghe F,Petty G W,et al.Application of SSM/I satellite data to a hurricane si mulation[J].Quarterly Journal of the Royal Meteorological Society,2004,130(598):801-825.

[3] Chen Shuhua.The i mpact of assi milating SSM/I and QuikSCAT satellite winds on Hurricane Isidore si mulations[J].Monthly Weather Review,2007,135(2):549-566.

[4] Singh R,Pal P K,Kishtawal C M,et al.Thei mpact of variational assi milation of SSM/I and Quik SCAT satellite observations on the numerical si mulation of Indian Ocean tropical cyclones[J].Weat her and Forecasting,2008,23(3):460-476.

[5] Tang Wenqing,Liu W T.Objective Interpolation of Scatterometer Winds[M].California:JPL Publication,1996.

[6] Royle J A,Berliner L M,Wikle C K,et al.A hierarchical spatial model for constructing wind fields fro m scatterometer data in the labrador sea[M]//Case Studies in Bayesian Statistics.New York:Springer,1999:367-382.

[7] Perrie W,Dunlap E,Vachon P W,et al.Marine wind analysis from remotely sensed measurements[J].Canadian Journal of Remote Sensing,2002,28(3):450-465.

[8] Atlas R.A multiyear global surface wind velocity dataset using SSM/I wind observations[J].Bulletin of the American Meteorological Society,1996,77(5):869-882.

[9] Chao Yi,Li Zhijin,Kindle J C,et al.A high-resolution surface vector wind product for coastal oceans:Blending satellite scatterometer measurements with regional mesoscale at mospheric model si mulations[J].Geophysical Research Letters,2003,30(1):1013.

[10] 凌征,王桂華,陳大可,等.中國近海風(fēng)場融合[C]//中國“數(shù)字海洋”論壇論文集.天津:國家海洋信息中心,2008:90-94.Ling Zheng,Wang Guihua,Chen Dake,et al.China's offshore wind field fusion[C]//China's“Digital Ocean”Foru m Proceedings.Tianjin:National Mariti me Infor mation Centres,2008:90-94.

[11] 蔣興偉,宋清濤.海洋衛(wèi)星微波遙感技術(shù)發(fā)展現(xiàn)狀與展望[J].科技導(dǎo)報,2010,28(3):105-111.Jiang Xingwei,Song Qingtao.Satellite microwave measurements of the global oceans and f uture missions[J].Science&Technology Review,2010,28(3):105-111.

[12] 齊亞琳,林明森.數(shù)據(jù)融合技術(shù)在海洋二號衛(wèi)星數(shù)據(jù)中的應(yīng)用[J].航天器工程,2012,21(3):117-123.Qi Yalin,Lin Mingsen.Application of the data fusion techniquein HY-2 satellite data[J].Spacecraft Engineering,2012,21(3):117-123.

[13] 劉宇昕,張毅,王兆徽,等.基于ASCAT微波散射計風(fēng)場與NCEP再分析風(fēng)場的全球海洋表面混合風(fēng)場[J].海洋預(yù)報,2014,31(3):10-18.Liu Yuxin,Zhang Yi,Wang Zhaohui,et al.Global sea surface blended winds based on the scatterometer winds and NCEPreanalyzed winds[J].Marine Forecasts,2014,31(3):10-18.

[14] 項杰,王慧鵬,王春明,等.南海海面風(fēng)場變分融合的初步研究[J].熱帶氣象學(xué)報,2015,31(2):153-160.Xiang Jie,Wang Huipeng,Wang Chun ming,et al.Preli minary study on variational blending of surface vector windsin South China Sea[J].Journal of Tropical Meteorology,2015,31(2):153-160.

[15] 肖庭延,于慎根,王彥飛.反問題的數(shù)值解法[M].北京:科學(xué)出版社,2003.Xiao Tingyan,Yu Shengen,Wang Yanfei.Numerical Solution of Inverse Problem[M].Beijing:Science Press,2003.

[16] 劉繼軍.不適定問題的正則化方法及應(yīng)用[M].北京:科學(xué)出版社,2005.Liu Jijun..The Regularization Method and Application of Unsuitable Problem[M].Beijing:Science Press,2005.

[17] Wang Wanqiu,Xie Pingping.A multiplatfor m-merged(MPM)SST analysis[J].Journal of Cli mate,2007,20(9):1662-1679.

[18] 姜祝輝,黃思訓(xùn),杜華棟,等.利用變分結(jié)合正則化方法對高度計風(fēng)速資料調(diào)整海面風(fēng)場的研究[J].物理學(xué)報,2010,59(12):8968-8977.Jiang Zhuhui,Huang Sixun,Du Huadong,et al.A new approach to adjusting sea surface wind using altimeter wind data by variational regularization method[J].Acta Physica Sinica,2010,59(12):8968-8977.

[19] 趙延來,黃思訓(xùn),杜華棟,等.正則化方法同化多普勒天氣雷達資料及對降雨預(yù)報的影響[J].物理學(xué)報,2011,60(7):875-886.Zhao Yanlai,Huang Sixun,Du Huadong,et al.Regularization method of assi milating Doppler radar data anditsinfluence on precipitation forecast[J].Acta Physica Sinica,2011,60(7):875-886.

[20] Wang Sicheng,Huang Sixun,Xiang Jie,et al.Three-di mensional ionospheric tomography reconstruction using the model function approach in Tikhonov regularization[J].Journal of Geophysical Research:Space Physics,2016,121(12):12104-12115.

[21] Jet Propulsion Laboratory.Quik SCAT science data product user's manual,Version2.22001[EB/OL].[2016-10-09].daac,jpl.nasa.gov/quikscat/qscat doc.ht ml

[22] Ebuchi N,Graber H C,Caruso MJ.Evaluation of wind vectors observed by QuikSCAT/Sea Winds using ocean buoy data[J].Jour nal of At mospheric and Oceanic Technology,2002,19(12):2049-2062.

[23] Hoff man R N,Leidner S M,Henderson J M,et al.At wo-di mensional variational analysis met hod for NSCAT ambiguity removal:Methodology,sensitivity,and tuning[J].Journal of At mospheric and Oceanic Technology,2003,20(5):585-605.

[24] Parrish D F,Derber J C.The National Meteorological Center's spectral statistical-inter polation analysis system[J].Mont hly Weather Review,1992,120(8):1747-1763.

[25] 曹小群,黃思訓(xùn),張衛(wèi)民,等.區(qū)域三維變分同化中背景誤差協(xié)方差的模擬[J].氣象科學(xué),2008,28(1):8-14.Cao Xiaoqun,Huang Sixun,Zhang Wei min,et al.Modeling background error covariance in regional 3D-VAR[J].Scientia Meteorologica Sinica,2008,28(1):8-14.

[26] Kunisch K.On a class of damped Morozov principles[J].Computing,1993,50(3):185-198.

[27] Hansen P C,O’Leary D P.The use of t he L-curvein the regularization of discreteill-posed problems[J].SIA MJournal on Scientific Computing,1993,14(6):1487-1503.

[28] Golub G H,Heath M,Wahba G.Generalized cross-validation as a method for choosing a good ridge parameter[J].Technometrics,1979,21(2):215-223.

[29] Kunisch K,Zou J.Iterative choices of regularization parameters in linear inverse problems[J].Inverse Problems,1998,14(5):1247.

[30] Lu S,Pereverzev S V.Multi-parameter regularization and its nu merical realization[J].Nu merische Mathematik,2011,118(1):1-31.

[31] Wang Z.Multi-parameter Tikhonov regularization and model f unction approach to the damped Morozov principle for choosing regularization parameters[J].Journal of Computational and Applied Mathematics,2012,236(7):1815-1832.

[32] Rodgers C D.Inverse Methods for At mospheric Sounding:Theory and Practice[M].Singapore:World Scientific Publishing Co Pte Ltd.,2000.

[33] Lupu C,Gaut hier P,Laroche S.Evaluation of t hei mpact of observations on analysesin 3D-and 4D-Var based on infor mation content[J].Monthly Weather Review,2011,139(3):726-737.

附錄

關(guān)于正則化項Jr的處理:

海面風(fēng)場的渦度ζ和散度D,表達式如下:

式中,u,v分別為單個網(wǎng)格點的緯向和經(jīng)向風(fēng)速,對Jr進行處理可得:

式(A2)可以簡化為:

將積分離散化,利用前差法計算二階偏導(dǎo),最終可得:

式中,Δs為單個網(wǎng)格面積,W1=L1TL1,W2=L2TL2,L1和L2為二階偏導(dǎo)離散后的系數(shù)矩陣,分別為:

根據(jù)本文設(shè)計模型,分析風(fēng)場的分辨率為0.12°×0.12°,因此Δx≈Δy,Δs≈ Δx( )2,將吸收進系數(shù)β中,可以得到:

Blending satellite scatterometer data based on variational with multi-parameter regularization method

Zhang Kaifeng1,Deng Wanyue2,Wang Ting1,Wang Huipeng3,Xiang Jie1,4,Song Qingtao5,Liu Chunxia6

(1.Institute of Meteorology and Oceanogr aphy,PLA University of Science and Technology,Nanjing 211101,China;2.Xinjiang Uygur Autono mous Region Environ mental Monitoring Station,Uru mqi 830011,China;3.Unit No.75839 of PLA,Guangzhou 510510,China;4.Key Laboratory of Mesoscale Severe Weather(Nanjing University),Ministry of Education,Nanjing 210023,China;5.National Satellite Ocean Application Ser vice,State Oceanic Ad ministration,Beijing 100081,China;6.Guangzhou Institute of Tropical and Marine Meteorology,China Meteorological Ad ministration,Guangzhou 510080,China)

A 3DVAR method with regularization constraints is proposed to blend sea surface wind datain the South China Sea based on the traditional 3DVAR and regularization technology of the inverse problem,and the model function method which is used to deter mine the reasonable regularization parameters and then the blended experiments of the satellite scatterometer(Quik SCAT)and Guang Zhou Mesoscale Model(GZMM)sea surface wind field data are carried out for a typhoon case.Results show that when we use the regularization method for experiments,the false infor mation caused by the traditional 3DVAR is eli minated obviously and the noiseis al most disappeared,at the same ti me,the wind field and vorticity field as well as divergence field are distributed evenly,and the structure is clear,more i mportantly,it is clear that the cyclone center is remarkable,and observation is dramatic in the analysis field.Besides,the degrees of freedomfor signal(DFS)method is used to evaluate blended systems quantitatively,it is found that the regularized constraint 3DVAR system has a higher DFSand observation influence related to traditional 3DVAR.The blended results are tested based on theindependent observation data,it indicates that the result of regularized constraint 3DVAR method has the s mallest root mean square error and maxi mumcorrelation coefficient,which is better than the statistical result of GZMMand the conventional 3DVAR method.

3DVAR;regularization;Quik SCAT;the degrees of freedomfor signal;model function

P715.6

A

0253-4193(2017)12-0122-14

張凱峰,鄧婉月,王挺,等.基于變分多參數(shù)正則化方法融合衛(wèi)星散射計資料[J].海洋學(xué)報,2017,39(12):122-135,

10.3969/j.issn.0253-4193.2017.12.012

Zhang Kaifeng,Deng Wanyue,Wang Ting,et al.Blending satellite scatterometer data based on variational with multi-parameter regularization method[J].Haiyang Xuebao,2017,39(12):122-135,doi:10.3969/j.issn.0253-4193.2017.12.012

2017-01-23;

2017-03-16。

國家自然科學(xué)基金(41275113);全球變化與海氣相互作用專項。

張凱峰(1994—),男,山西省運城市人,主要從事海面風(fēng)場融合研究。E-mail:15380426538@163.com*通信作者:項杰,男,安徽省巢湖市人,教授,主要從事海面風(fēng)場融合研究。E-mail:xjieah@aliyun.com