誤差系數(shù)對(duì)變分同化的質(zhì)量影響研究*

蔣金利 諸雪征 顧 進(jìn) 韓朝帥，2 馬 巖，3 李廣峰

（1.陸軍防化學(xué)院 北京 102205）（2.軍事科學(xué)院防化研究院 北京 102205）（3.武警四川總隊(duì)機(jī)動(dòng)第二支隊(duì) 康定 626000）

1 引言

數(shù)據(jù)同化算法是一種最初來(lái)源于數(shù)值天氣預(yù)報(bào)，為其提供初始場(chǎng)的數(shù)據(jù)處理技術(shù)，現(xiàn)已廣泛應(yīng)用于大氣海洋領(lǐng)域。在積分描寫動(dòng)力系統(tǒng)演變過(guò)程的數(shù)學(xué)模式的同時(shí)，不斷吸收觀測(cè)資料，給出系統(tǒng)狀況的整體估計(jì)，且觀測(cè)資料僅僅通過(guò)少量的觀測(cè)站便可獲取，大大降低經(jīng)濟(jì)成本，并有效提高了數(shù)值天氣預(yù)報(bào)的準(zhǔn)確度與可靠性。

國(guó)外早在1922年，Richadson就首次將同化的思想用于數(shù)值天氣預(yù)報(bào)的主觀分析中［1］。當(dāng)時(shí)他將觀測(cè)數(shù)據(jù)以手工的方式插值到網(wǎng)格點(diǎn)，以此作為預(yù)報(bào)的初始場(chǎng)，因?yàn)榉绞竭^(guò)于簡(jiǎn)單，未做其他誤差處理，導(dǎo)致最后預(yù)報(bào)失敗。我國(guó)在該領(lǐng)域最早的研究是1958年顧震潮學(xué)者提出數(shù)值天氣預(yù)報(bào)不是簡(jiǎn)單的初值問(wèn)題［2］。20世紀(jì)60年代我國(guó)著名氣象學(xué)家丑紀(jì)范創(chuàng)造性地提出將變分法和泛函分析引入數(shù)值天氣預(yù)報(bào)領(lǐng)域，該理論比國(guó)外同行的研究早了近十年。

目前主流的數(shù)據(jù)同化算法主要分為兩類（按照算法與模型之間的關(guān)聯(lián)機(jī)制分類）：

1）以卡爾曼濾波、集合卡爾曼濾波等算法為代表的順序數(shù)據(jù)同化算法。

2）以三維變分、四維變分等算法為代表的連續(xù)數(shù)據(jù)同化算法。

其中，自從20世紀(jì)90年代以來(lái)，隨著計(jì)算機(jī)性能的大幅度提升，變分同化算法在實(shí)際應(yīng)用中得到快速發(fā)展，尤其是在數(shù)值天氣預(yù)報(bào)領(lǐng)域，已經(jīng)被國(guó)際上多個(gè)業(yè)務(wù)中心作為業(yè)務(wù)算法進(jìn)行使用［3］。中國(guó)國(guó)家氣象局［4］，美國(guó)國(guó)家大氣研究中心［5～6］，英國(guó)國(guó)家天氣預(yù)報(bào)中心［7］以及歐洲中尺度天氣預(yù)報(bào)中心［8］都將變分同化算法作為核心業(yè)務(wù)算法實(shí)現(xiàn)了業(yè)務(wù)化。在數(shù)值天氣預(yù)報(bào)領(lǐng)域的成功實(shí)踐證明了變分同化算法的可行性和正確性。目前同化算法除了氣象領(lǐng)域之外，在海洋［9～11］，大氣污染物預(yù)測(cè)［12～15］，地質(zhì)土壤［16～18］等領(lǐng)域都得到了較大的發(fā)展。

2 變分同化基本原理

假設(shè)n維向量xt表示某一物質(zhì)擴(kuò)散的真實(shí)狀態(tài)，n維向量xb表示背景場(chǎng)狀態(tài)，m維向量y表示觀測(cè)向量，H表示觀測(cè)算子，是觀測(cè)值與真實(shí)狀態(tài)之間的紐帶函數(shù)。同時(shí)，我們用εb表示背景場(chǎng)的誤差，εo表示觀測(cè)場(chǎng)的誤差，那么我們便可以得到：

假設(shè)誤差準(zhǔn)確，且估計(jì)無(wú)偏，那么εb、εo的數(shù)學(xué)期望為0，則背景誤差協(xié)方差矩陣為

觀測(cè)誤差協(xié)方差矩陣為

在數(shù)據(jù)同化中，我們一般認(rèn)為背景與觀測(cè)是兩個(gè)相互獨(dú)立的事件。假設(shè)Pb（x）表示背景場(chǎng)的概率分布，Po（x）表示觀測(cè)值的概率分布，且兩者相互獨(dú)立，那么這兩個(gè)概率分布的乘積就是它們的聯(lián)合概率密度函數(shù)，也就是我們所求分析場(chǎng)的概率分布，即：

根據(jù)最大似然估計(jì)原理，極大似然估計(jì)值就是使Pa(x)取得最大值時(shí)的變量值，也就是所求的分析值。在實(shí)際應(yīng)用中，如果不知道實(shí)際的氣象條件，擴(kuò)散狀態(tài)，一般假設(shè)為各向同性，其概率則要服從高斯分布，概率密度函數(shù)也應(yīng)為高斯函數(shù)，則：

兩者的概率密度相乘，便可以得到分析值的概率密度函數(shù)。由于我們所要求的是使聯(lián)合概率密度函數(shù)取極大值時(shí)的函數(shù)值，并不是為了去極大值。所以，一般來(lái)說(shuō)，為了簡(jiǎn)單計(jì)算，我們一把將聯(lián)合概率密度函數(shù)取對(duì)數(shù)值。同時(shí)，消去負(fù)號(hào)，使求極大值問(wèn)題轉(zhuǎn)變?yōu)榍髽O小值問(wèn)題，也是簡(jiǎn)便計(jì)算的常用方法之一。以此為基礎(chǔ)，構(gòu)造出的目標(biāo)函數(shù)為

上式的含義是表示分析場(chǎng)與背景場(chǎng)的距離以及分析場(chǎng)與觀測(cè)場(chǎng)的距離。該函數(shù)的值越小，說(shuō)明距離越小，說(shuō)明分析值與觀測(cè)值之間的“距離”越小。

式中：x表示分析值，xb為背景場(chǎng)，B為背景誤差協(xié)方差，R為觀測(cè)誤差協(xié)方差，y為觀測(cè)值，H為觀測(cè)算子。

3 誤差協(xié)方差矩陣

誤差是整個(gè)同化過(guò)程中不可避免的。數(shù)學(xué)角度看，可以把同化過(guò)程理解為最小化誤差的過(guò)程。誤差在同化開始前主要體現(xiàn)在背景誤差協(xié)方差矩陣B上。由于實(shí)際的操作過(guò)程中，不管是所設(shè)定的背景場(chǎng)數(shù)據(jù)，亦或是由各種型號(hào)的偵察裝備所得到的觀測(cè)場(chǎng)數(shù)據(jù)，甚至是模型本身都不是絕對(duì)準(zhǔn)確的，都是有誤差的。但是由于很難確定“真值”xtrue到底是多少，就無(wú)法得知背景誤差和觀測(cè)誤差到底是多少。同時(shí)當(dāng)誤差一直存在于同化的整個(gè)過(guò)程中時(shí)，整個(gè)分析過(guò)程都是帶有誤差的，這個(gè)分析誤差是上述三種誤差在同化過(guò)程中相互作用的結(jié)果。

一般情況，誤差協(xié)方差矩陣假設(shè)觀測(cè)誤差與背景誤差不相關(guān)，即：

式中：εb為背景誤差；εo為觀測(cè)誤差，T為矩陣轉(zhuǎn)置。

在實(shí)際分析中，這也是完全合理的假設(shè)，背景場(chǎng)參量與觀測(cè)場(chǎng)參量相互獨(dú)立。

3.1 背景誤差協(xié)方差矩陣

背景誤差協(xié)方差矩陣表示預(yù)報(bào)誤差的概率分布函數(shù)，是一個(gè)高斯型的函數(shù)分布。同化系統(tǒng)所產(chǎn)生的分析場(chǎng)參量大概只有15%來(lái)自觀測(cè)資料，而剩余的85%均來(lái)自背景場(chǎng)［19］。在物理意義上，背景場(chǎng)是同化前該地域的狀態(tài)場(chǎng)，它將同化循環(huán)中的觀測(cè)數(shù)據(jù)不斷融合進(jìn)最新的分析場(chǎng)。

背景誤差協(xié)方差矩陣在同化中具有很重要的作用：首先它具有信息傳遞的作用，尤其是在觀測(cè)資料短缺的地區(qū)，分析增量的形狀幾乎完全由協(xié)方差決定，矩陣的相關(guān)性決定了從觀測(cè)點(diǎn)到周圍空間的信息傳播；其次具有信息平滑的作用，在觀測(cè)資料稠密地區(qū)，它的相關(guān)性也決定了信息的平滑程度；它還能體現(xiàn)各個(gè)參量之間的平衡關(guān)系，比如高度位勢(shì)和風(fēng)場(chǎng)之間的平衡。

與卡爾曼濾波數(shù)據(jù)同化方法［20～21］不同的是，變分法理論上并不能很好地處理背景誤差協(xié)方差矩陣B，因此背景項(xiàng)的設(shè)立就顯得尤為重要［22～23］。

變分公式中的B是這樣定義的：

式中xtrue代表真值，εb為背景誤差，很多時(shí)候B矩陣也表示為

不難看出，B矩陣是正定對(duì)稱矩陣，而且特征值均為非負(fù)值。

假設(shè)背景誤差為向量(e1，e2，e3)，其背景誤差協(xié)方差矩陣可寫成：

3.2 觀測(cè)誤差協(xié)方差矩陣

觀測(cè)誤差協(xié)方差矩陣包含關(guān)于觀測(cè)誤差的統(tǒng)計(jì)信息。對(duì)觀測(cè)系統(tǒng)來(lái)說(shuō)，不同位置上的觀測(cè)誤差在統(tǒng)計(jì)上是相互獨(dú)立的。觀測(cè)誤差一般包括觀測(cè)算子、儀器和代表性誤差。

觀測(cè)算子誤差主要是因?yàn)閷?duì)兩種不同種類的數(shù)據(jù)之間聯(lián)系的理解不足導(dǎo)致，主要出現(xiàn)在非直接觀測(cè)中。

代表性誤差一般來(lái)源于兩方面：一是模型不完美，得不到完全準(zhǔn)確的測(cè)量；二是有限的模型分辨率，這類誤差的減小可以利用高密度觀測(cè)，提高模式分辨率。

儀器誤差是出廠時(shí)帶有的，一般會(huì)帶有數(shù)據(jù)參數(shù)，考慮設(shè)備老化等實(shí)際情況可以對(duì)參數(shù)進(jìn)行適當(dāng)調(diào)整。除此之外，也鮮有好的解決辦法。

一般而言，對(duì)R的處理較簡(jiǎn)化，在觀測(cè)之間相互獨(dú)立情況下，將其簡(jiǎn)化為對(duì)角矩陣，方便計(jì)算。

4 變分同化誤差系數(shù)影響分析

在變分同化模塊中，同化質(zhì)量受到多種不同因素的影響，其本身各類參數(shù)的影響最大。在上一節(jié)的基礎(chǔ)上，通過(guò)數(shù)據(jù)模擬方法［24］分析背景場(chǎng)參量和觀測(cè)場(chǎng)參量對(duì)同化質(zhì)量的影響。

簡(jiǎn)單起見，在同一高度層進(jìn)行仿真試驗(yàn)，只考慮同一種化學(xué)物質(zhì)一個(gè)平面層內(nèi)的連續(xù)擴(kuò)散，不考慮不同觀測(cè)設(shè)備所獲取的數(shù)據(jù)格式的不統(tǒng)一，簡(jiǎn)化觀測(cè)算子為單位矩陣或單位向量。

設(shè)定統(tǒng)一的試驗(yàn)先驗(yàn)條件：界定危害區(qū)域?yàn)?0km×10km，風(fēng)向?yàn)?°，風(fēng)速固定為3m/s，大氣穩(wěn)定度為D級(jí)，擴(kuò)散參數(shù)以GB/T為準(zhǔn)，迭代次數(shù)1000次。假定同化窗口0時(shí)刻為擴(kuò)散發(fā)生4小時(shí)之后，以模型首次預(yù)報(bào)值加上隨機(jī)誤差生成背景場(chǎng)，以模型后續(xù)預(yù)報(bào)值作為真值，將預(yù)報(bào)值加上20%隨機(jī)誤差生成觀測(cè)值。

4.1 背景誤差系數(shù)對(duì)同化的影響分析

保證以上條件不變，設(shè)定觀測(cè)均方差為1，選取算法為四維變分?jǐn)?shù)據(jù)同化算法，簡(jiǎn)單分析背景場(chǎng)對(duì)同化結(jié)果的影響。

首次試驗(yàn)，設(shè)定初始背景協(xié)方差矩陣為二維單位矩陣，同化結(jié)果走勢(shì)如圖1所示。

圖1 首次試驗(yàn)同化結(jié)果圖

圖中可以看出，在20次迭代的時(shí)候，分析值明顯比觀測(cè)值和背景值更加接近真值，起到同化效果。隨著迭代次數(shù)的不斷增多，各條曲線間的差異便不是那么容易可以看出來(lái)了。在200次迭代的同化圖中，不難看出，基本已經(jīng)達(dá)到了同化效果，并且具有收斂性。這與圖2的誤差以及方差走勢(shì)圖也是相一致的。

圖2 首次同化誤差變化圖

從圖中可以明顯看出在200次迭代后，誤差及方差基本收斂。

第二次試驗(yàn)，調(diào)整初始背景場(chǎng)誤差系數(shù)，設(shè)置為首次的0.5倍大小，再次進(jìn)行數(shù)值模擬試驗(yàn)。結(jié)果如圖3所示。

圖3 第二次背景參數(shù)影響分析試驗(yàn)

從第二次的試驗(yàn)結(jié)果可以看出，背景初始誤差的變化會(huì)導(dǎo)致同化精度提高，但是并不影響同化的誤差以及方差的收斂性。

第三次試驗(yàn)，調(diào)整背景場(chǎng)誤差系數(shù)，設(shè)置為首次的0.2倍大小，再次進(jìn)行數(shù)值模擬試驗(yàn)。結(jié)果如圖4所示。

圖4 第三次背景誤差系數(shù)影響分析試驗(yàn)

從第三次試驗(yàn)的同化結(jié)果可以看出，背景誤差系數(shù)的改變會(huì)影響初始同化的精度，但是不會(huì)影響同化過(guò)程的收斂性。

圖5是三次背景場(chǎng)參數(shù)試驗(yàn)同化效果對(duì)比圖。

由圖5可以看出，改變初始的背景場(chǎng)誤差系數(shù)，會(huì)對(duì)同化的效果和精度產(chǎn)生影響。當(dāng)背景誤差系數(shù)減小的時(shí)候，同化精度會(huì)有提高。

圖5 背景誤差系數(shù)分析試驗(yàn)對(duì)比圖

但是，從以上三次試驗(yàn)分析來(lái)看，背景場(chǎng)誤差系數(shù)對(duì)同化結(jié)果造成影響的同時(shí)，并不會(huì)影響整個(gè)同化窗口內(nèi)的收斂性，三次試驗(yàn)基本都在200次迭代之后開始收斂。

4.2 觀測(cè)誤差系數(shù)對(duì)同化的影響分析

保持初始設(shè)置條件不變，選取算法為三維變分同化算法，設(shè)置B為對(duì)角矩陣，對(duì)角線元素為誤差方差，設(shè)定其為固定值0.4。這次試驗(yàn)分析同化后的兩個(gè)小時(shí)內(nèi)的同化情況。

由于觀測(cè)誤差協(xié)方差矩陣R為對(duì)角陣，首次試驗(yàn)設(shè)定觀測(cè)誤差均方差為0.1。

試驗(yàn)結(jié)果如圖6所示。

圖6 首次觀測(cè)參數(shù)影響分析試驗(yàn)

從首次試驗(yàn)的結(jié)果來(lái)看，由于背景誤差系數(shù)設(shè)定為固定數(shù)值，所以整體試驗(yàn)的收斂性較好，基本20次迭代后已成收斂狀態(tài)。同時(shí)，由于觀測(cè)誤差協(xié)方差矩陣的性質(zhì)為對(duì)角矩陣，且對(duì)角線元素也已經(jīng)設(shè)置，所以觀測(cè)場(chǎng)方差的圖形是一條直線。

第二次試驗(yàn)，保持其他參數(shù)不變，上調(diào)觀測(cè)誤差系數(shù)為0.25，結(jié)果如圖7所示。

圖7 第二次觀測(cè)參數(shù)影響分析試驗(yàn)

就第二次的試驗(yàn)結(jié)果來(lái)看，誤差系數(shù)加大，同化的穩(wěn)定性加大，同化的波動(dòng)加大，但是收斂性并沒有改變。

第三次試驗(yàn)，保持其他參數(shù)不變的基礎(chǔ)上，上調(diào)觀測(cè)誤差參數(shù)為0.30，結(jié)果如圖8所示。

圖8 第三次觀測(cè)響參數(shù)分析試驗(yàn)

第三次同化的結(jié)果來(lái)看，同化窗口內(nèi)波動(dòng)加大，收斂性不變。

圖9是三次同化效果的對(duì)比圖。

圖9 觀測(cè)參數(shù)分析試驗(yàn)對(duì)比圖

對(duì)比來(lái)看，觀測(cè)參數(shù)的變化會(huì)影響同化的質(zhì)量，誤差越大，波動(dòng)越大，效果越差。同時(shí)，參數(shù)改變不影響收斂性。

5 結(jié)語(yǔ)

通過(guò)對(duì)兩種不同的誤差參數(shù)的試驗(yàn)結(jié)果來(lái)看，背景誤差系數(shù)和觀測(cè)誤差系數(shù)的改變會(huì)影響同化的質(zhì)量，誤差系數(shù)越大，同化窗口內(nèi)的波動(dòng)越大，效果越差，這也是符合常識(shí)認(rèn)知的；但是，通過(guò)對(duì)比發(fā)現(xiàn)，背景場(chǎng)參數(shù)變化對(duì)同化的影響程度要高于觀測(cè)參數(shù)的影響，這也是為什么在實(shí)際工作中，認(rèn)為背景場(chǎng)更加重要的原因；本研究的結(jié)果也進(jìn)一步表明，不管是背景參數(shù)還是觀測(cè)參數(shù)的改變，誤差系數(shù)的大小都不影響同化的整體收斂性。