基于特征正交分解與離散經驗插值的淺水波模式降階解法

陶蘇林， 李雨鴻， 周 寅， 劉 垚

(1.南京信息工程大學應用氣象學院,南京 210044；2.南京大橋機器有限公司技術開發(fā)中心,南京 211101；3.遼寧省氣象科學研究所,沈陽 110166；4.94857部隊61分隊,蕪湖 241007；5.寧夏氣象防災減災重點實驗室,銀川 750002)

描述實際大氣問題建立的動力學模型通常是高維或無限維的非線性系統(tǒng),模型有效性受限于離散格式較多的自由度[1]。尋求既可保證模擬精度又可降低計算量的模型降階技術是區(qū)別于并行計算的有效思路。構造降階模型基本思路是以較精確的數學描述刻畫原數值模型的主要動力學特征,并且這種數學描述的階數和計算代價遠低于原數值模型[2]。

特征正交分解方法(proper orthogonal decomposition method, POD)提供已知大量數據在最小二乘意義下的一組正交基,以實現給定數據的最優(yōu)低維逼近,其降階效能在諸多研究中得到證明[3—6]。POD降階模型可較好模擬流場運動細節(jié)、刻畫較強的非線性特征,但對非線性項的計算代價依然依賴于原始系統(tǒng)的維度m。離散經驗插值方法(discrete empirical interpolation method, DEIM)則進一步給出了非線性項降階處理思路,即在精選的n個(n<

1 淺水波模式

引入Parrett和Cullen測試使用的一維淺水波模式[9]。該模式形式簡單,但其運行依然可以達到極為復雜的程度,在測試氣象問題方面十分有效,且有典型性。該模式刻畫均質不可壓縮無黏旋轉淺流體的基本運動,忽略層結,描述大氣與海洋運動的重要特征:①流體是旋轉的,科氏加速度是重要的；②運動的形態(tài)比是小量。Lorenc[10]基于該模式討論了有關非線性系統(tǒng)的資料同化問題。此外,Griffith等[11]也利用該模式構建討論變分同化方案在估算系統(tǒng)誤差分量與時間相關的模式誤差分量研究中的有效性。

該模式包含科氏項和剛性底部地形邊界以描述淺層大氣運動,其控制方程為[9]

{?u?t+u?u?x+?φ?x-fv+g?H?x=0

?v?t+u?v?x+fu=0

?φ?t+u?φ?x+φ?u?x=0

式(1)中,x∈Ω=(0,2πL)且t∈(0,T)；u=u(x,t)且v=v(x,t) 分別表示緯向和經向風速；f為科氏參數；φ=gη(x,t)為重力場,其中g為重力加速度,且η(x,t)>0為底部以上大氣流體的平均深度。剛性底邊界由H=H(x)定義。

式(1)雖然形式簡單,但為非線性,并且通過該框架可探究忽略層結情況下旋轉作用和底部地形對大氣運動的影響。淺水波模式[式(1)]的離散格式參見文獻[12],該格式通過引入擴散項來抑制二階有限差分格式在地形躍變(如山峰)附近產生的虛假震蕩。Parrett和Cullen[9]通過與Glimm法得到的參考解比較,證明使用的離散格式可以提供足夠精確的解。Griffith[13]也根據Parrett和Cullen提出的包含地形條件的案例測試了該離散格式,并得到了與Parrett和Cullen研究一致的結果。

若定義tk時刻的風場和位勢場：

{uk=(uk1,uk2,…,ukj)T

vk=(vk1,vk2,…,vkj)T

φk=(φk1,φk2,…,φkj)T

并且：

{mk=(uk○φk)

nk=(vk○φk)

于是淺水波模式[式(1)]離散格式可寫成:

(4)

式(4)中,f1(uk-1,mk-1)、f2(φk-1)和f3(vk-1,mk-1)為非線性項,具體形式為

(5)

其他矩陣形式如下:

(6)

(7)

(8)

(9)

(10)

(11)

(12)

取科氏參數f=7.292×10-5s-1,擴散系數K=5.0×102m2·s-1,重力加速度g=9.8 m·s-2,空間尺度坐標L=9.55×105m,時間尺度坐標T=1.2×105s。利用這些參數設置使得淺水波模式[式(4)]可以模擬30°N的地球大氣運動[11]。另外,模式剛性底邊界定義為

(13)

式(13)中,a=10Δx。

令空間格點數J=100,積分時間步數N=100；在[0,2πL]×[0,T]范圍內,淺水波模式[式(4)]的初始值條件定義為

(14)

式(14)中，大氣流體平均深度為ηm=1 m,并且j=1,2,…,J。因此,數值試驗中大氣流體運動是由中心位置的地形開始發(fā)展的。此外,為保證離散格式[式(4)]的穩(wěn)定性,時間步長Δt應遵循CFL(courant-friedrichs-lewy)條件。

2 POD/DEIM降階解法

給定Hilbert空間L2內一個自由度為m的非線性系統(tǒng):

(15)

式(15)中，空間x∈Ω,時間t∈[0,T],解y(x,t)=[y(x1,t),y(x2,t),…,y(xm,t)]T∈Rm×1,系數矩陣M,K∈Rm×m以及關于y(x,t)的非線性函數F={F[y(x1,t)],F[y(x2,t)],…,F[y(xm,t)]}T(l ?R時F:l→R)。定義該模型的緊湊形式為

(16)

則構建該模型的POD/DEIM降階模型步驟如下。

(1)根據一定離散方法構建原始數學物理模型的數值格式,將原始模型的連續(xù)形式轉變?yōu)橛善⒎址匠虡嫵傻碾x散形式,即確定模式。

(2)積分確定模式獲得約定時段的模型數值解,按照約定取樣方式選取瞬像(snapshot)[14]構成瞬像集合Y。對瞬像集合Y進行奇異值分解,獲得降序排列的特征值及其對應的特征向量。這些特征向量構成正交基函數,進一步用于構建POD基函數。

(3)選取前r個正交基函數構成一個子空間,重構模式變量yPOD。

(4)將重構解代入第1步的確定模式,并根據Galerkin投影變換進一步將確定模式投影至由POD基函數張成的降維空間,獲得原始模式的降階形式

(17)

(5)進一步利用DEIM方法處理POD降階模型[式(17)]中的非線性項。

(6)求解常微分方程,得到降階系數。返回第3步,利用r個正交基函數重構模式變量yPOD。

2.1 降階基向量選取

2.1.1 獲得瞬像樣本

進而基于數據集①～⑤構建新的樣本數據集:{mk=0,m1,…,mN}, {nk=0,n1,…,nN}, {φk=0,φ1,…,φN}。其中:k=0,1,6,11,…,N-4時刻的狀態(tài)變量取自數據集①相應時刻的狀態(tài)變量；k=2,7,12,…,N-3時刻的狀態(tài)變量取自數據集②；k=3,8,13,…,N-2時刻的狀態(tài)變量取自數據集③；k=4,9,14,…,N-1時刻的狀態(tài)變量取自數據集④；k=5,10,15,…,N時刻的狀態(tài)變量取自數據集⑤。

(18)

(19)

2.1.2 特征值分解

對矩陣A、B和C分別進行奇異值分解,可以獲得按降序排列的特征值：

(20)

并定義：

(21)

2.1.3 確定POD基向量

選擇r滿足：r=argmin{I(r):I(r)≥γ0}且γ0∈[0,1]則可獲得與這r個特征值對應的一組正交基向量,并構成降階空間:

(22)

2.2 淺水波模式的POD降階模型

利用式(22)提供的POD基函數近似狀態(tài)變量mk、nk和φk:

(23)

式(23)中,αm,k、αn,k和αφ,k∈Rr×1為待求解系數。因此,只要獲得所有時刻的αm,k、αn,k和αφ,k,就可以得到所有時間層的狀態(tài)變量。

將式(23)代入淺水波模式[式(4)]，可得：

(24)

利用(Φm)T、(Φn)T和(Φφ)T分別左乘式(24),則可得到淺水波模式的POD降階模型:

(25)

該降階模型的初始條件為

(26)

周期邊界條件通過矩陣Mi(i=1,2,…,5)設定,具體見文獻[12]。因此,利用式(25)求解淺水波模式的POD降階模型,獲得系數αm,k、αn,k和αφ,k,從而重構狀態(tài)變量。

2.3 淺水波模式的POD/DEIM降階模式

(27)

式(27)中，Pf1、Pf2和Pf3根據DEIM算法確定。并且可以得到：

(28)

因此,淺水波模式的POD/DEIM降階形式為

{αm,k=(Φm)TM1MPOD1∈Rr×J(+Φmαm,k-2)+

(Φm)TM2MPOD2∈Rr×J(+Φnαn,k-1)+

(Φm)TM3MPOD3∈Rr×J(+Φφαφ,k-1)+

αn,k=(Φn)TM2POD2∈Rr×J(+Φnαn,k-2)+

αφ,k=(Φφ)TM2POD2∈Rr×J(+Φφαφ,k-2)+

(Φφ)TM5MPOD5∈Rr×J(+Φmαm,k-1)-(Φφ)T

式(26)為該降階模型[式(29)]的初始條件,周期邊界條件通過矩陣Mi(i=1,2,…,5)設定。并且,由式(29)可知,POD/DEIM降階模型非線性項的計算僅利用s(s≤J)個空間格點。因此,POD/DEIM降階模型的計算代價較POD降階模型有了進一步減小。

2.4 降階模型精度驗證

(30)

3 算例與分析

令FULL-MODEL表示全階模型、POD-ROM表示POD降階模型、POD/DEIM-ROM表示POD/DEIM降階模型。算例結果基于CPU型號為Intel酷睿i3—2100(主頻3.10 GHz)、內存2 GB以及安裝Windows XP Service Pack 3的計算機獲得。

3.1 降階模型計算效率

設計8種數值試驗討論不同空間格點數時降階模型計算效率。設定降階模型需捕獲全階模型狀態(tài)變量99.8%以上的能量,即對于狀態(tài)變量m、n和φ均需滿足I(r)≥99.8%；設定[0,2πL]區(qū)間內格點數分別為150、200、250、300、350、400、450和500,積分步數為格點數2倍,瞬像樣本數為積分步數1/2。

隨著格點數增加,FULL-MODEL、POD-ROM 及POD/DEIM-ROM這三種模型CPU耗時均呈增加趨勢,且FULL-MODEL的CPU耗時增幅最大(表1)?？臻g格點數相同時,三種模型耗費CPU時間排序為:FULL-MODEL > POD-ROM > POD/DEIM-ROM耗時。由此可見,經過非線性項處理的POD/DEIM-ROM模型計算效率最高。同時,隨著格點數及積分步數的增加,POD-ROM以及POD/DEIM-ROM所需POD基函數及非線性項基函數均呈增加趨勢。并且,為了捕獲全階模型超過99.8%的能量,狀態(tài)變量φ較m和n需選擇更多基函數,非線性項f1比f2和f3也需使用更多基函數。

利用DEIM方法處理非線性項后,降階模型可捕獲的能量有所下降。POD/DEIM-ROM數值解(狀態(tài)變量u、v和φ)與FULL-MODEL數值解的均方根誤差均大于POD-ROM與FULL-MODEL數值解的均方根誤差(表2)。但隨格點數增加,兩種降階模型數值解誤差未明顯增加。POD-ROM模擬的各個狀態(tài)變量(u、v和φ)與FULL-MODEL數值解的均方根誤差變化在空間格點數增至300后趨于平緩?？臻g格點數小于300時,狀態(tài)變量u和φ的均方根誤差先減后增,狀態(tài)變量v的均方根誤差則呈逐漸增加趨勢。POD/DEIM-ROM模擬的狀態(tài)變量u和φ與FULL-MODEL的均方根誤差在格點數少于400時逐漸增加,之后有所減少；而狀態(tài)變量v與FULL-MODEL的均方根誤差在格點數為150的數值試驗中達最小,其他趨勢與u和φ基本一致。

POD-ROM數值解(u、v和φ)與FULL-MODEL的相關系數總體達0.995以上,其中u和v相關系數達0.999(表3)。POD/DEIM-ROM數值解相關系數均低于POD-ROM與FULL-MODEL數值解相關系數,表明非線性項的進一步近似處理將丟失更多的模式能量。因此,使用DEIM方法進一步處理模式〗非線性項后,需增加非線性項基函數維數以提高POD/DEIM-ROM捕獲全階模型的能量。

表1 格點數不同時POD-ROM和POD/DEIM-ROM所需CPU時間和POD基函數

表2 格點數不同時POD-ROM和POD/DEIM-ROM與FULL-MODEL的均方根誤差

表3 格點數不同時POD-ROM、POD/DEIM-ROM與FULL-MODEL的相關系數

此外,從POD/DEIM-ROM模擬的各個狀態(tài)變量與FULL-MODEL的相關系數來看,水平風場v分量數值解效果最好,其次是u分量,位勢場φ模擬效果最差。隨格點數增加,POD-ROM和POD/DEIM-ROM數值解(u、v)與FULL-MODEL數值解相關系數無明顯變化,但POD/DEIM-ROM數值解φ相關系數明顯下降,且當格點數增至500時,相關系數僅為0.508。

3.2 瞬像選擇與影響

瞬像選擇是降階模型應用的重要步驟。瞬像樣本維數是降階模型的重要參數,特別是使用DEIM方法處理非線性項后,還需考慮非線性項的瞬像選擇,相當于增加了降階模型參數。討論降階模型對瞬像樣本選擇的敏感性,對于指導降階模型的應用有很重要的意義。

根據狀態(tài)變量和非線性項瞬像維數的不同組合,設計了36種敏感試驗。設定狀態(tài)變量m、n和φ均滿足I(r)≥99.8%；設定[0,2πL]區(qū)間內格點數為100,積分步數為格點數的3倍；狀態(tài)變量m、n和φ以及非線性項f1、f2和f3的瞬像樣本數分別取10、20、50、100、200和300,并設定狀態(tài)變量m、n和φ的瞬像維數相等,非線性項f1、f2和f3的瞬像樣本數相等。

圖1 狀態(tài)變量和非線性項瞬像維數不同時POD-ROM和POD/DEIM-ROM與FULL-MODEL數值解的均方根誤差Fig.1 The RMSE of the POD-ROM and POD/DEIM-ROM in different snapshot dimension of the state variable and nonlinear term with FULL-MODEL

圖2 狀態(tài)變量和非線性項瞬像維數不同時POD-ROM和POD/DEIM-ROM與FULL-MODEL數值解的相關系數Fig.2 The correlation coefficient of the POD-ROM and POD/DEIM-ROM in different snapshot dimension of the state variable and nonlinear term with FULL-MODEL

僅從CPU耗時來看,狀態(tài)變量m、n和φ與非線性項f1、f2和f3的瞬像維數不同組合對POD-ROM和POD/DEIM-ROM計算效率沒有明顯影響。但從降階模型模擬的狀態(tài)變量與全階模型狀態(tài)變量的均方根誤差以及相關系數來看,不同組合的瞬像維數對降階模型數值解影響很明顯(圖1、圖2)。POD/DEIM-ROM數值解與FULL-MODEL狀態(tài)變量的均方根誤差較大,其對應的相關系數也較低,POD-ROM結果卻與之相反。當狀態(tài)變量瞬像維數為200時,POD-ROM和POD/DEIM-ROM數值解與FULL-MODEL數值解均方根誤差均較大。并且,POD-ROM數值解均方根誤差較其他類型的瞬像維數組合明顯增大,對應的相關系數也明顯減小,POD/DEIM-ROM對瞬像選擇也有類似敏感性。因此,本試驗當區(qū)間[0,2πL]內空間格點數為100、積分步數為300時,POD-ROM和POD/DEIM-ROM狀態(tài)變量瞬像維數應避免設置為200。

3.3 DEIM插值點選擇與影響

利用POD/DEIM降階理論處理非線性項的關鍵是確定合理的DEIM插值點。為了討論不同數量DEIM插值點對POD/DEIM-ROM數值解的計算精度影響,對數值試驗作如下設置:①設定[0,2πL]區(qū)間內格點數為100,積分步數為300；②設定狀態(tài)變量m、n和φ以及非線性項f1、f2和f3的瞬像樣本數為100；③根據狀態(tài)變量(m、n和φ)瞬像矩陣POD譜(圖3)確定選擇狀態(tài)變量前10個POD基函數,即令r=10；④設定非線性項f1、f2和f3的DEIM插值點分別為90、80、70、60、50、40、30、20和10,以測試POD/DEIM-ROM模擬效率。

隨著DEIM插值點維數減少,POD/DEIM-ROM降階模型耗費的CPU時間呈明顯減少趨勢,但該降階模型數值解u、v和φ與FULL-MODEL數值解均方根誤差和相關系數無明顯變化(表4)。因此,在數值試驗條件下,淺水波模式的POD/DEIM-ROM數值解受DEIM插值點維數變化影響不明顯,但計算效率得到了很大提高。此外,從各個狀態(tài)變量模擬結果來看,選擇相同數量DEIM插值點時,POD/DEIM-ROM對位勢場φ的模擬結果最差,均方根誤差最大,相關系數最小,水平風場u和v分量模擬結果較好。

圖3 狀態(tài)變量(m、n和φ)與非線性項(f1、f2和f3)瞬像矩陣POD譜Fig.3 The POD spectrum corresponding to the solution snapshots for the state variables (m, n and φ) and nonlinear terms (f1, f2 and f3)

表4 DEIM插值點維數不同時POD/DEIM-ROM所需CPU時間及其與FULL-MODEL的均方根誤差和相關系數

全階模型FULL-MODEL預報的狀態(tài)變量u、v和φ如圖4。狀態(tài)變量瞬像維數為10,非線性項DEIM插值點維數為10時的POD/DEIM-ROM數值解u、v和φ則如圖5。

狀態(tài)變量瞬像維數為10,非線性項DEIM插值點維數為10時的POD/DEIM-ROM數值解u、v和φ與FULL-DODEL數值解均方根誤差如圖6所示。

POD/DEIM-ROM在前向積分中期和后期模擬結果最大值逐漸偏離全階模型結果。并且,隨著POD/DEIM-ROM積分時刻增加,狀態(tài)變量u、v和φ均方根誤差均呈明顯上升趨勢。而非線性項DEIM插值點維數在10～90間變化時,POD/DEIM-ROM模擬結果未發(fā)生明顯變化,表明不同DEIM插值點維數對POD/DEIM-ROM數值解結果影響相對較小。

圖4 全階模型預報的狀態(tài)變量Fig.4 The forecasted state variable by full model

圖5 降階模型預報的狀態(tài)變量Fig.5 The forecasted state variable by reduced-order model

圖6 降階模型與全階模型預報的狀態(tài)變量均方根誤差Fig.6 RMSE of the state variable between the reduced-order model and full model

4 討論

隨著對大氣運動過程的深入理解,現有數值天氣預報模式不斷得到更新完善。比如,在原先的中尺度模式中嵌套小尺度模式,以提高模擬分辨率(特別是空間分辨率)。但這些復雜的數值模式需要消耗大量計算資源。應用并行計算技術是解決這個問題的重要手段,如朱小謙[15]研究了高分辨率譜模式的可擴展并行算法,提高了模式運行效率。但如果可以直接對原始數值模式進行降階處理,則可簡化很多額外算法研究。

采取直接對原始數值模式進行降階處理的思路,對旋轉大氣中有限區(qū)域淺水波模式進行降階處理。通過POD方法獲得淺水波模式的POD-ROM,并利用DEIM方法進一步對非線性項進行處理,獲得淺水波模式的POD/DEIM-ROM,試驗結果表明POD/DEIM降階算法是有效的。取狀態(tài)變量瞬像矩陣前10個POD基函數,POD-ROM就可以捕獲全階模式99.8%以上的能量,并且CPU耗時明顯減少,該優(yōu)勢在空間格點變得稠密時尤其明顯。這些結果與Cao等[16]、Du等[17]和Dimitriu等[18]等研究結論相同。而通過DEIM方法對淺水波模式的3個非線性項進行近似處理后,POD/DEIM-ROM的CPU耗時進一步減少。

但是,降階模型數值解與全階模式數值解仍然存在一定偏差。特別是經過DEIM處理后的降階模型因為丟失了全階模式的部分非線性(但可能是非常重要的)信息,其數值解與全階模式數值解均方根誤差明顯大于POD-ROM。因此,為了實現模型降階,同時又保持較高的計算精度,需考慮對降階模型進行訂正。目前,有相關研究取得了初步成果。如,Du等[17]通過在瞬像矩陣創(chuàng)建過程應用SobolevH1范數來訂正POD降階模型；Wang等[8]通過利用正則化技術對二維Burgers方程降階模型進行訂正,提高了POD/DEIM降階模型精度。本文還未涉及降階模型訂正,將在未來研究中展開討論。

此外,還需注意POD/DEIM方法雖然進一步提高了降階模型的計算效率,同時也引入了更多的可變參數。這些參數的設置對于降階模型的模擬精度影響不容忽視。3.2與3.3節(jié)試驗表明淺水波模式的POD-ROM和POD/DEIM-ROM對瞬像和DEIM插值點維數等可變參數的敏感性。因此,如何選擇可變參數的最優(yōu)取值將是POD/DEIM降階模型應用中的重要研究問題。

5 結論

利用特征正交分解方法(POD)和離散經驗插值方法(DEIM)建立了旋轉大氣中有限區(qū)域淺水波模式的降階模型,并研究了降階模型(ROM)的計算效率以及瞬像和DEIM插值點的選擇與影響。研究得出如下結論。

(1) POD/DEIM-ROM計算效率明顯高于POD-ROM和FULL-MODEL,并且可以捕獲FULL-MODEL超過99.8%的能量,特別當空間格點數量明顯增加時,POD/DEIM-ROM計算效率最高。

(2) 狀態(tài)變量和非線性項瞬像樣本數變化對降階模型CPU耗時影響不大,但對數值解質量影響比較明顯。POD/DEIM-ROM數值解質量受DEIM插值點維數變化影響不明顯,但隨DEIM插值點數量減小,POD/DEIM-ROM的CPU耗時明顯減少。

實際應用中,數值模式的POD/DEIM-ROM對瞬像和DEIM插值點維數有較強的敏感性。這些可變參數的最優(yōu)取值問題將是未來研究難點。

致謝:感謝南京信息工程大學數學與統(tǒng)計學院王曰朋教授和應用氣象學院申雙和教授對本研究的支持和指導。