基于變分多參數正則化方法融合衛星散射計資料

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

基于變分多參數正則化方法融合衛星散射計資料

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

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

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

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

1 引言

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

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

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

2 研究區域以及數據

2.1 研究區域

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

2.2 數據介紹

2.2.1 獨立觀測數據

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

2.2.2 Quik SCAT散射計資料

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

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

2.2.3 模式風場數據

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

2.2.4 數據質量控制

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

3 研究方法

3.1 方法介紹

3.1.1 常規3DVAR方法

3DVAR的目的是融合觀測數據與背景場。基于以下假設:(1)觀測及背景誤差服從均值為0的正態分布;(2)觀測誤差和背景誤差不相關,則3DVAR問題歸結為如下常規的目標函數的極小化(稱為常規的3DVAR):

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

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

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

式中,d A=d x d y,β為正則化參數,則常規的目標函數(1)變為:

故對于本文要研究的海面風場融合來說,

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

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

3.1.3 誤差協方差矩陣的統計

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

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

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

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

3.1.4 正則化參數的選取

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

則吸收Morozov偏差方程為:

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

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

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

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

3.1.5 信號自由度的計算

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

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

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

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

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

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

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

由于本文研究構造的融合系統中,矩陣維數較小,可以直接進行DFS的計算。

3.2 方案設計流程

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

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

4 數值試驗及結果

4.1 臺風個例分析

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

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

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

定義均方根誤差(RMSE):

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

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

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

圖6 α=1(a,d,g),α=0.1(b,e,h)以及α=0.01(c,f,i)時分析風場(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的數值變化Tal.1 The RMSE of O-B and O-A as well as the value of J o and J b whenβ=0

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

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

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

4.2 敏感性分析

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

式中,p為觀測數量。這里針對本文的情形計算相應的OI,計算結果見表2。

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

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

5 融合效果評估

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

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

基于獨立的西沙自動站觀測數據,對融合效果進行評估,通過以下幾個統計量進行檢驗,統計結果見表3。定義A為對比數據,B為參考數據,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):

相關系數(R):

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

表3 統計對比結果Tab.3 The results of statistical comparison

6 結論和討論

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

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

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

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

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

[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] 凌征,王桂華,陳大可,等.中國近海風場融合[C]//中國“數字海洋”論壇論文集.天津:國家海洋信息中心,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] 蔣興偉,宋清濤.海洋衛星微波遙感技術發展現狀與展望[J].科技導報,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] 齊亞琳,林明森.數據融合技術在海洋二號衛星數據中的應用[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微波散射計風場與NCEP再分析風場的全球海洋表面混合風場[J].海洋預報,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] 項杰,王慧鵬,王春明,等.南海海面風場變分融合的初步研究[J].熱帶氣象學報,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] 肖庭延,于慎根,王彥飛.反問題的數值解法[M].北京:科學出版社,2003.Xiao Tingyan,Yu Shengen,Wang Yanfei.Numerical Solution of Inverse Problem[M].Beijing:Science Press,2003.

[16] 劉繼軍.不適定問題的正則化方法及應用[M].北京:科學出版社,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] 姜祝輝,黃思訓,杜華棟,等.利用變分結合正則化方法對高度計風速資料調整海面風場的研究[J].物理學報,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] 趙延來,黃思訓,杜華棟,等.正則化方法同化多普勒天氣雷達資料及對降雨預報的影響[J].物理學報,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] 曹小群,黃思訓,張衛民,等.區域三維變分同化中背景誤差協方差的模擬[J].氣象科學,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.

附錄

關于正則化項Jr的處理:

海面風場的渦度ζ和散度D,表達式如下:

式中,u,v分別為單個網格點的緯向和經向風速,對Jr進行處理可得:

式(A2)可以簡化為:

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

式中,Δs為單個網格面積,W1=L1TL1,W2=L2TL2,L1和L2為二階偏導離散后的系數矩陣,分別為:

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

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

張凱峰,鄧婉月,王挺,等.基于變分多參數正則化方法融合衛星散射計資料[J].海洋學報,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。

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

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