不同矩陣分解方法對(duì)海洋數(shù)據(jù)同化的影響*

管志斌，肖俊敏，季統(tǒng)凱，洪學(xué)海+，譚光明，馬 巖

1.中國(guó)礦業(yè)大學(xué)（北京）機(jī)電與信息工程學(xué)院，北京 100083

2.中國(guó)科學(xué)院 計(jì)算技術(shù)研究所 高性能計(jì)算機(jī)研究中心，北京 100190

1 引言

當(dāng)今時(shí)代，對(duì)海洋生態(tài)系統(tǒng)和海洋動(dòng)力學(xué)的研究早已進(jìn)入衛(wèi)星時(shí)代，科研人員可以通過(guò)數(shù)據(jù)傳輸衛(wèi)星，獲取分布于各個(gè)海洋不同區(qū)域內(nèi)，大量傳感器所收集的海洋中的多尺度信息。但是，隨著大量的海洋觀測(cè)數(shù)據(jù)的獲取，不可避免地出現(xiàn)了一些亟待解決的難題，即如何將這些海洋觀測(cè)數(shù)據(jù)或者數(shù)據(jù)結(jié)果，與相應(yīng)的數(shù)據(jù)模型進(jìn)行結(jié)合，從而使得其更接近于現(xiàn)實(shí)世界中實(shí)際的海洋狀態(tài)。

海洋數(shù)據(jù)同化是一種能夠?qū)⒑Ｑ笥^測(cè)數(shù)據(jù)與相應(yīng)的海洋數(shù)值模型進(jìn)行融合的有效方法。海洋數(shù)據(jù)同化方法主要是將收集到的海洋觀測(cè)數(shù)據(jù)與理論數(shù)值模型相結(jié)合，取兩者的優(yōu)點(diǎn)，從而得到更接近于實(shí)際狀態(tài)的結(jié)果[1-2]。高效的數(shù)據(jù)同化方法，不僅可以從大量的海洋觀測(cè)數(shù)據(jù)中快速地提取出有價(jià)值的信息，而且可以通過(guò)提高數(shù)據(jù)分析的精度和預(yù)測(cè)技巧，使人類加深對(duì)海洋生態(tài)系統(tǒng)和海洋動(dòng)力學(xué)等相關(guān)領(lǐng)域信息的理解和認(rèn)識(shí)，從而更好地處理海洋與人類社會(huì)關(guān)系。

海水溫度和鹽度是諸多海洋水文信息中兩個(gè)最基本、最重要的信息尺度。溫度是度量海水冷熱程度的基本物理量；鹽度是對(duì)海水含鹽量的一種衡量尺度，是海洋水體的最重要理化特性之一。隨著Argo（array for real-time geostrophic ocean-ography）浮標(biāo)的全球海洋大范圍覆蓋以及遙感觀測(cè)衛(wèi)星的連續(xù)性觀測(cè)，使得溫鹽等觀測(cè)數(shù)據(jù)被大量獲取，如何有效地利用這些海洋觀測(cè)數(shù)據(jù)，并對(duì)不同海洋尺度信息進(jìn)行融合，對(duì)海洋信息研究具有極其重要的意義[3-4]。

通過(guò)海洋數(shù)據(jù)同化，將觀測(cè)數(shù)據(jù)與數(shù)值模型進(jìn)行融合，并從大量海洋數(shù)據(jù)中提取有價(jià)值的信息。在實(shí)際中，預(yù)報(bào)誤差的形成，主要是由于模型欠缺以及預(yù)見(jiàn)性誤差，而數(shù)據(jù)同化可以確保系統(tǒng)的預(yù)報(bào)具有穩(wěn)健的正確性[5-6]。海洋數(shù)據(jù)同化方法，不僅可以結(jié)合原始的海洋觀測(cè)數(shù)據(jù)與數(shù)值模型的優(yōu)點(diǎn)，使預(yù)測(cè)結(jié)果更接近于客觀的實(shí)際海洋狀態(tài)，而且能夠使人類更深入地研究并理解海洋狀態(tài)[7-9]。

本文的結(jié)構(gòu)如下：第1章對(duì)海洋數(shù)據(jù)同化的基本概念進(jìn)行簡(jiǎn)單介紹；第2章中，在闡述海洋數(shù)據(jù)同化基本流程和數(shù)據(jù)來(lái)源之后，將分別介紹集合卡爾曼濾波同化方法，可有效節(jié)省計(jì)算開(kāi)銷的集合最優(yōu)插值，以及目前集合最優(yōu)插值中所采用的基于奇異值分解進(jìn)行矩陣求逆的方法；第3章將嘗試對(duì)同化過(guò)程中分析場(chǎng)的求解算法進(jìn)行改進(jìn)，主要集中在矩陣求逆運(yùn)算，對(duì)比LU分解求逆、Choleskey分解求逆、QR分解求逆以及奇異值分解（singular value decomposition，SVD）求逆；在第4章中，將給出同化程序基于四種不同矩陣分解的性能測(cè)試結(jié)果，并對(duì)數(shù)值結(jié)果進(jìn)行對(duì)比分析；第5章對(duì)全文進(jìn)行總結(jié)，并給出展望。

2 數(shù)據(jù)同化方法

2.1 海洋數(shù)據(jù)同化

在對(duì)海洋溫鹽數(shù)據(jù)（Argo溫鹽廓線數(shù)據(jù)）進(jìn)行同化的過(guò)程中，由于EnOI（ensemble optimal interpolation）所具有的不可忽略的計(jì)算便捷性等優(yōu)點(diǎn)，EnOI常常被作為實(shí)際應(yīng)用中的同化方法[10]。

具體的流程如圖1所示。

Fig.1 Flow chart of ocean data assimilation圖1 海洋數(shù)據(jù)同化流程圖

（1）需要獲取Argo溫鹽廓線數(shù)據(jù)，并進(jìn)行適當(dāng)?shù)馁|(zhì)量控制，使用已有的背景場(chǎng)數(shù)據(jù)，執(zhí)行同化命令，先對(duì)Argo數(shù)據(jù)進(jìn)行同化，包括Argo海流數(shù)據(jù)的同化、Argo溫度數(shù)據(jù)的同化、Argo鹽度數(shù)據(jù)同化，完成Argo數(shù)據(jù)的同化之后，得到更新的背景場(chǎng)[11]。

（2）在獲取海面高度計(jì)異常數(shù)據(jù)（sea level anomaly，SLA）之后，同樣先進(jìn)行質(zhì)量控制，再使用同化Argo得到的背景場(chǎng)，對(duì)高度計(jì)數(shù)據(jù)進(jìn)行同化，得到更新的背景場(chǎng)。

（3）在完成海表面溫度數(shù)據(jù)（sea surface temperature，SST）的質(zhì)量控制之后，使用（2）中更新過(guò)的背景場(chǎng)，來(lái)同化海表面溫度數(shù)據(jù)，完成海洋數(shù)據(jù)同化任務(wù)[12]。

2.2 數(shù)據(jù)來(lái)源

Argo是全球范圍內(nèi)的3 800個(gè)左右的自由漂移剖面浮標(biāo)，被放置在不同海洋的多個(gè)區(qū)域內(nèi)，用于測(cè)量海洋上層2 000 m左右范圍內(nèi)的溫度和鹽度的信息[13]。使用這些浮標(biāo)，可以使得其連續(xù)觀測(cè)的海洋數(shù)據(jù)，如海洋溫度、鹽度、流速等，在收集后數(shù)小時(shí)內(nèi)即可公開(kāi)并使用。

如圖2所示，Argo浮標(biāo)會(huì)自動(dòng)潛入距海平面2 km深處的等密度層上，并隨深層海流保持中性漂移，在經(jīng)過(guò)預(yù)定時(shí)間之后，會(huì)自動(dòng)上浮，并在上浮的過(guò)程中，利用Argo浮標(biāo)自身攜帶的各種溫度、鹽度等傳感器進(jìn)行連續(xù)的剖面測(cè)量。Argo浮標(biāo)給出的是離散的剖面觀測(cè)數(shù)據(jù)，包括海洋特定位置上的溫度、鹽度、壓力傳感器的聯(lián)合測(cè)量值[14-15]。當(dāng)浮標(biāo)到達(dá)海平面后，會(huì)自動(dòng)搜索用于數(shù)據(jù)傳輸?shù)男l(wèi)星，在實(shí)現(xiàn)自身定位后，會(huì)使用數(shù)據(jù)傳輸衛(wèi)星系統(tǒng)，將收集到的連續(xù)時(shí)間內(nèi)的測(cè)量數(shù)據(jù)，傳送到衛(wèi)星地面接收站[16]。當(dāng)浮標(biāo)收集的測(cè)量數(shù)據(jù)全部傳輸完畢之后，會(huì)再次自動(dòng)下沉，執(zhí)行下一個(gè)循環(huán)過(guò)程，進(jìn)行數(shù)據(jù)的收集操作。

Fig.2 Flow chart ofArgo data collection圖2 Argo數(shù)據(jù)收集流程圖

2.3 集合卡爾曼濾波

集合卡爾曼濾波（ensemble Kalman filter，EnKF），是一種利用蒙特卡羅的短期集合預(yù)報(bào)方法，來(lái)估計(jì)預(yù)報(bào)誤差協(xié)方差的四維集合同化方法；EnKF是由Evensen于1994年真正地應(yīng)用，并首次推廣到數(shù)據(jù)同化領(lǐng)域；EnKF將集合預(yù)報(bào)（ensemble forecast）的思想與卡爾曼濾波（Kalman filter，KF）相結(jié)合，總結(jié)兩者的優(yōu)點(diǎn)，通過(guò)對(duì)一組有限的預(yù)測(cè)集合進(jìn)行估計(jì)，得到背景場(chǎng)誤差協(xié)方差，在此基礎(chǔ)上，再進(jìn)行最優(yōu)分析，計(jì)算出具有復(fù)雜結(jié)構(gòu)的基于流的，隨空間、時(shí)間變化的預(yù)報(bào)誤差統(tǒng)計(jì)量[17-18]。

由于背景誤差協(xié)方差是決定數(shù)據(jù)同化質(zhì)量的關(guān)鍵因素之一，因此背景誤差協(xié)方差矩陣的建立與構(gòu)造就顯得尤為重要[19]。

集合卡爾曼濾波（EnKF）繼承了卡爾曼濾波（KF）的優(yōu)點(diǎn)，是一種基于迭代更新的計(jì)算方法[20]。

集合卡爾曼濾波中“分析場(chǎng)”的表達(dá)式如下：

其中，Af是預(yù)測(cè)矩陣，稱為背景場(chǎng)矩陣，A′是擾動(dòng)預(yù)測(cè)矩陣，D′=D-HAf，D是擾動(dòng)測(cè)量矩陣，?e是擾動(dòng)測(cè)量誤差矩陣；而H是測(cè)量算子，表示真實(shí)模型狀態(tài)值與測(cè)量值之間的映射關(guān)系。

令：

由式（1）可知，在Aa的計(jì)算過(guò)程中最耗費(fèi)時(shí)間的部分是對(duì)M求逆。為了得到M-1，傳統(tǒng)的方法是直接對(duì)M這個(gè)m×m矩陣，進(jìn)行特征分解，使得：

而由此所產(chǎn)生的計(jì)算開(kāi)銷則為m2，特別的，當(dāng)m特別大時(shí)，這種開(kāi)銷將無(wú)法承受。因此，可以使用其他更有效的分解算法，先將矩陣M分解成子矩陣的乘積，且這些矩陣的求逆過(guò)程具有較小的計(jì)算開(kāi)銷。

由于集合擾動(dòng)和測(cè)量誤差之間是不相關(guān)的（與預(yù)測(cè)誤差也是不相關(guān)的），可以得到：

此時(shí)對(duì)矩陣M的求逆過(guò)程，可以轉(zhuǎn)換為先對(duì)Q進(jìn)行求逆計(jì)算，從而降低計(jì)算開(kāi)銷。

通過(guò)EnKF可以得到同化過(guò)程的最優(yōu)解，但是EnKF使用的是大量的預(yù)報(bào)樣本和測(cè)量數(shù)據(jù)，這在數(shù)值計(jì)算上是較為低效的；因此，本文所采用的是，在數(shù)值計(jì)算上極其高效的，相比較于EnKF，可以得到次優(yōu)解的集合最優(yōu)插值方法。

2.4 集合最優(yōu)插值

集合最優(yōu)插值（EnOI），是在集合卡爾曼濾波（EnKF）思想的基礎(chǔ)之上提出的一種次優(yōu)的集合同化方法，記為 EnOI[21]。

與EnKF的不同之處在于，EnOI只需要傳入一個(gè)樣本數(shù)據(jù)，即EnOI只需要利用一組有限的、靜態(tài)的模式狀態(tài)的樣本集合，就可以完成對(duì)背景誤差協(xié)方差的估計(jì)與構(gòu)建，并實(shí)現(xiàn)對(duì)背景誤差的流依賴性的刻畫(huà)，而EnKF的使用過(guò)程需要多個(gè)預(yù)測(cè)樣本，因此，在計(jì)算資源的開(kāi)銷上，EnOI要遠(yuǎn)遠(yuǎn)小于EnKF，更加適用于大氣海洋數(shù)據(jù)的研究[22-23]。EnOI分析的計(jì)算方程與式（1）相似，可以通過(guò)如下方程進(jìn)行計(jì)算：

其中：

G=αHA′A′THT+?e?eT，′=(d-Hθf(wàn))，θf(wàn)是包含多個(gè)變量的背景場(chǎng)，θa是相應(yīng)的分析場(chǎng)，參數(shù)α∈(0,1]是額外引入的權(quán)重因子，作用在集合中的測(cè)量值上；很顯然，式（8）每次只需要計(jì)算一個(gè)模型狀態(tài)。

EnOI中，所使用的模型狀態(tài)樣本數(shù)據(jù)，來(lái)自于多年積分的模式結(jié)果，在客觀世界中，由于不同時(shí)間段的氣候差異較大，這會(huì)影響到分析結(jié)果，因此引入?yún)?shù)α，來(lái)降低氣候差異所帶來(lái)的負(fù)面影響[24-25]。

類似于式（5），則：

則集合最優(yōu)插值中的求解方程轉(zhuǎn)化為：

EnOI是一種非常受關(guān)注的方法，它可以有效地節(jié)省計(jì)算時(shí)間。當(dāng)不隨時(shí)間變化的靜態(tài)樣本集合被創(chuàng)建之后，在分析過(guò)程中[26]，由于僅對(duì)一個(gè)模型狀態(tài)進(jìn)行更新，因此同化過(guò)程中最后的計(jì)算開(kāi)銷會(huì)從O(N2)降低到O(nN)，其中N表示測(cè)量值向量的數(shù)量，n表示測(cè)量值向量中存儲(chǔ)的測(cè)量值的數(shù)量。

在EnOI中求解矩陣的逆所用的傳統(tǒng)方法是奇異值分解，在下一節(jié)，將介紹奇異值分解方法[27]在EnOI中的具體實(shí)現(xiàn)。

2.5 奇異值分解

定理1（奇異值分解（SVD））設(shè)A是秩為r的m×n實(shí)值矩陣，則存在m階的正交矩陣Um×m=(u1,u2,…,um)和n階的正交矩陣Vn×n=(v1,v2,…,vn)，及對(duì)角矩陣Σ，使得：

此時(shí)，根據(jù)奇異值分解（SVD）定理，可先對(duì)式（10）進(jìn)行SVD分解，得到：

其中，U和VT均為正定矩陣，Σ是對(duì)角矩陣，其對(duì)角線上存儲(chǔ)的是矩陣的奇異值，并且按遞減的順序在對(duì)角矩陣Σ的對(duì)角線上進(jìn)行排序。進(jìn)一步，根據(jù)式（9）可知：

其中，Λ-1=Σ-1(Σ-1)T。因此，式（8）可以轉(zhuǎn)化為SVD分解的形式，如下：

在實(shí)際應(yīng)用EnOI進(jìn)行海洋數(shù)據(jù)同化的過(guò)程中，通過(guò)實(shí)驗(yàn)了解到，EnOI算法中的使用SVD分解進(jìn)行矩陣求逆的過(guò)程相對(duì)耗時(shí)較大；因此，本文的主要目的是通過(guò)對(duì)矩陣求逆的過(guò)程進(jìn)行優(yōu)化，從而減少矩陣求逆過(guò)程的運(yùn)行時(shí)間；所采用的方法是，使用其他的矩陣分解求逆方法進(jìn)行替代，從而對(duì)EnOI的計(jì)算過(guò)程進(jìn)行優(yōu)化。在下一章中，將介紹其他矩陣分解求逆方法。

3 矩陣分解求逆

在海洋數(shù)據(jù)同化應(yīng)用中，傳統(tǒng)的用于對(duì)式（10）進(jìn)行矩陣求逆運(yùn)算的方法是SVD分解；在同化程序中，對(duì)矩陣進(jìn)行分析，發(fā)現(xiàn)矩陣是一個(gè)稠密的對(duì)稱正定矩陣，矩陣維度是150×150；而矩陣求逆過(guò)程，需要迭代執(zhí)行394×750次；因此，可以嘗試采用其他的矩陣分解方法進(jìn)行替代，來(lái)實(shí)現(xiàn)海洋溫鹽數(shù)據(jù)同化過(guò)程中的矩陣求逆。

本章將分別介紹矩陣的LU分解[28]、Choleskey分解[29]以及QR分解[30]；然后，將這些分解方法分別應(yīng)用到同化過(guò)程的矩陣求逆過(guò)程中。

3.1 矩陣LU分解

定理2（LU分解）給定一個(gè)可逆矩陣A，則存在一個(gè)下三角矩陣L和一個(gè)上三角矩陣U，使得：

LU分解是將矩陣A通過(guò)初等變換L，變成一個(gè)上三角矩陣U，本質(zhì)上就是高斯消元的另一種表達(dá)形式。LU分解的計(jì)算時(shí)間復(fù)雜度為O(2n3/3)。

在此，將矩陣進(jìn)行LU分解，得到：

于是，可以得到式（8）的LU分解形式，如下：

3.2 矩陣Choleskey分解

定理3（Choleskey分解）給定一個(gè)n×n維的對(duì)稱正定矩陣A，則存在一個(gè)主對(duì)角線元素嚴(yán)格正定的下三角矩陣L，滿足：

其中，L?是矩陣L的共軛轉(zhuǎn)置矩陣。Choleskey分解的計(jì)算時(shí)間復(fù)雜度為O(n3/3)。

在這里，將矩陣進(jìn)行Choleskey分解，如下：

由此，模型分析值計(jì)算式（8）的Choleskey分解形式，如下：

3.3 矩陣QR分解

定理4（QR分解）給定一個(gè)n×n維的可逆矩陣A，如果矩陣A的列向量線性無(wú)關(guān)，則A可以分解為：

其中，Q是n×n正交矩陣，R是非奇異的上三角矩陣。

將矩陣的QR分解應(yīng)用到矩陣的求逆過(guò)程中，令：

則分析場(chǎng)的求解式（8）變?yōu)槿缦滦问剑?/p>

本章所介紹的矩陣LU分解、Choleskey分解、QR分解的主要思想是將待求逆矩陣分解為容易進(jìn)行求逆的三角矩陣或正定矩陣的乘積；通過(guò)對(duì)分解后的矩陣進(jìn)行求逆，再進(jìn)行相應(yīng)的矩陣乘積變換之后，即可得到最初的待求逆矩陣的逆；從時(shí)間復(fù)雜度上來(lái)看，三角矩陣或正定矩陣的求逆過(guò)程的時(shí)間復(fù)雜度要遠(yuǎn)小于最初稠密矩陣求逆過(guò)程的時(shí)間復(fù)雜度；在下一章中，將給出具體的實(shí)驗(yàn)結(jié)果，以及對(duì)實(shí)驗(yàn)結(jié)果的分析和總結(jié)。

4 實(shí)驗(yàn)結(jié)果與分析

采用控制變量法，測(cè)試5種不同時(shí)間段的數(shù)據(jù)分別在不同進(jìn)程數(shù)量下，采用SVD分解求逆、LU分解求逆、Choskey分解求逆、QR分解求逆之后，數(shù)據(jù)同化程序的運(yùn)行時(shí)間的對(duì)比結(jié)果，以及本文所采用的三種矩陣分解方法，使同化程序性能提升倍數(shù)的比較。其中所使用的數(shù)據(jù)，來(lái)源于不同時(shí)期的海洋Argo溫鹽廓線曲線數(shù)據(jù)、海表面溫度數(shù)據(jù)（SST）、高度計(jì)數(shù)據(jù)（SLA）；不同的時(shí)期包含有：2012_003，2012_005，2012_007，2012_033，2012_035，即2012年的第3、5、7、33、35天的數(shù)據(jù)。使用的進(jìn)程數(shù)（NP）分別是2、4、6、8、10、12、14、16。

表1～表5（a）給出的是，針對(duì)5種不同的數(shù)據(jù)（2012_003，2012_005，2012_007，2012_033，2012_035），分別基于矩陣SVD分解、LU分解、Choleskey分解、QR分解來(lái)實(shí)現(xiàn)矩陣求逆之后，在不同進(jìn)程數(shù)下（NP=2,4,6,8,10,12,14,16），數(shù)據(jù)同化程序的運(yùn)行時(shí)間（elapsed time）對(duì)比，斜體數(shù)值表示耗時(shí)最長(zhǎng)且性能最差，加粗?jǐn)?shù)值表示耗時(shí)最少且性能最優(yōu)。表1～表5（b）展示的是，相比SVD分解而言，分別采用LU分解、Choleskey分解、QR分解之后，在不同進(jìn)程數(shù)下（NP=2,4,6,8,10,12,14,16），同化程序的性能提升倍數(shù)。表1～表5清晰地展示了LU分解、Choleskey分解、QR分解相比于SVD分解，能夠使得集合最優(yōu)插值的計(jì)算效率提升至少兩倍以上。

從表1～表5中可以明顯地看出，傳統(tǒng)的SVD分解進(jìn)行矩陣求逆，程序運(yùn)行所需的時(shí)間要遠(yuǎn)多于LU分解、Choleskey分解和QR分解，即總體耗時(shí)：SVD＞QR＞LU＞Choleskey；在進(jìn)程數(shù)為2時(shí)，SVD分解求解所需要的時(shí)間最多，是其他3種分解方法進(jìn)行求解所需時(shí)間的3倍左右；隨著進(jìn)程數(shù)目的增加，4種分解方法進(jìn)行求解所需的時(shí)間都逐漸降低，但采用LU分解、Choleskey分解、QR分解之后，程序的運(yùn)行時(shí)間依然明顯小于傳統(tǒng)方法所采用的SVD分解。

Table 1(a) Comparison of runtime on assimilation program(AP)(Data1)表1(a) 同化程序的運(yùn)行時(shí)間對(duì)比（數(shù)據(jù)1） s

Table 1(b) Multiple of performance enhancement on assimilation program(AP)(Data1)表1(b) 同化程序的性能提升倍數(shù)（數(shù)據(jù)1）

Table 2(a) Comparison of runtime on assimilation program(AP)(Data2)表2(a)同化程序的運(yùn)行時(shí)間對(duì)比（數(shù)據(jù)2） s

Table 3(a) Comparison of runtime on assimilation program(AP)(Data3)表3(a) 同化程序的運(yùn)行時(shí)間對(duì)比（數(shù)據(jù)3） s

Table 4(a) Comparison of runtime on assimilation program(AP)(Data4)表4(a) 同化程序的運(yùn)行時(shí)間對(duì)比（數(shù)據(jù)4） s

此外，從表1～表5（b）中可以看到，在并行運(yùn)行下，分別采用LU分解、Choleskey分解、QR分解之后，同化程序在性能上，相比較于SVD分解提升了至少兩倍左右；而LU分解和Choleskey分解的提升效果最為明顯。

Table 2(b) Multiple of performance enhancement on assimilation program(AP)(Data2)表2(b) 同化程序的性能提升倍數(shù)（數(shù)據(jù)2）

Table 3(b) Multiple of performance enhancement on assimilation program(AP)(Data3)表3(b) 同化程序的性能提升倍數(shù)（數(shù)據(jù)3）

Table 4(b) Multiple of performance enhancement on assimilation program(AP)(Data4)表4(b) 同化程序的性能提升倍數(shù)（數(shù)據(jù)4）

值得一提的是，采用矩陣QR分解的同化程序，在性能提升效果上遠(yuǎn)不如LU分解和Choleskey分解。這一現(xiàn)象的主要原因是，采用QR分解之后，得到的是一個(gè)正交矩陣Q和一個(gè)上三角矩陣的乘積，而采用LU分解和Choleskey分解之后得到的是下三角矩陣和上三角矩陣的乘積；由于稠密正交矩陣Q求逆過(guò)程的時(shí)間復(fù)雜度要大于下三角矩陣求逆的時(shí)間復(fù)雜度，因此，造成了QR分解的效果不如LU分解和Choleskey分解。

Table 5(a) Comparison of runtime on assimilation program(AP)(Data5)表5(a) 同化程序的運(yùn)行時(shí)間對(duì)比（數(shù)據(jù)5） s

從同化程序運(yùn)行時(shí)間上來(lái)看可以了解：QR分解＞LU分解＞Choleskey分解，即采用矩陣的Choleskey分解之后，同化程序的運(yùn)行時(shí)間會(huì)被縮短到最低，大幅度縮減了同化程序的運(yùn)行時(shí)間，提高了程序計(jì)算效率。

圖3～圖5分別給出了4、8、16個(gè)進(jìn)程下的同化程序性能對(duì)比。

在圖3～圖5（a）中，同化程序分別采用SVD分解、LU分解、Choleskey分解以及QR分解，并針對(duì)5種不同數(shù)據(jù)進(jìn)行測(cè)試（2012_003，2012_005，2012_007，2012_033，2012_035），得到的運(yùn)行時(shí)間對(duì)比；圖3～圖5（b）給出的則是分別采用LU分解、Choleskey分解以及QR分解之后，同化程序的運(yùn)行時(shí)間對(duì)比。

Table 5(b) Multiple of performance enhancement on assimilation program(AP)(Data5)表5(b) 同化程序的性能提升倍數(shù)（數(shù)據(jù)5）

從圖3～圖5可以看出，在性能對(duì)比上，分別采用LU分解、Choleskey分解或QR分解進(jìn)行矩陣求逆的同化程序，其性能遠(yuǎn)遠(yuǎn)優(yōu)于采用SVD分解進(jìn)行矩陣求逆的同化程序。

綜上所述，采用本文3種矩陣分解算法的同化程序，相比較于采用SVD分解算法的同化程序，其性能提升至少兩倍左右；并且采用Choleskey分解能夠最大程度地減少同化程序運(yùn)行時(shí)間，且其擴(kuò)展性與LU分解和QR分解同樣優(yōu)異。因此，Choleskey分解值得被應(yīng)用于EnOI的矩陣求逆過(guò)程中，它將大幅度地提升目前海洋數(shù)據(jù)同化程序的計(jì)算效率。

Fig.3 Performance comparison of assimilation programs(AP)(NP=4)圖3 同化程序的性能對(duì)比(NP=4)

Fig.4 Performance comparison of assimilation programs(AP)(NP=8)圖4 同化程序的性能對(duì)比(NP=8)

Fig.5 Performance comparison of assimilation programs(AP)(NP=16)圖5 同化程序的性能對(duì)比(NP=16)

5 結(jié)束語(yǔ)

使用集合最優(yōu)插值（EnOI）進(jìn)行海洋數(shù)據(jù)同化時(shí)，由于EnOI的矩陣求逆過(guò)程所采用的SVD分解非常耗時(shí)，因此本文著眼于對(duì)矩陣求逆過(guò)程進(jìn)行優(yōu)化，嘗試使用矩陣LU分解、Choleskey分解、QR分解來(lái)替代SVD分解。通過(guò)并行程序?qū)嶒?yàn)結(jié)果的對(duì)比和分析可以看出，本文所嘗試的LU分解、Choleskey分解、QR分解，相比較于SVD分解，可以使同化程序的性能提升至少兩倍左右，并且Choleskey分解明顯優(yōu)于其他的矩陣分解方法，可以最大程度地降低時(shí)間復(fù)雜度，縮短同化程序運(yùn)行時(shí)間。

因此，在未來(lái)的工作中，可以考慮將矩陣Choleskey分解應(yīng)用到實(shí)際的工程項(xiàng)目中，以達(dá)到提高性能的目的。