李 淵，李云梅，呂 恒，吳傳慶，王珊珊，王永波

(1:江蘇省地理信息資源開發(fā)與利用協(xié)同創(chuàng)新中心，南京 210023)(2:環(huán)保部衛(wèi)星環(huán)境應(yīng)用中心，北京 100029)

太湖葉綠素a同化系統(tǒng)敏感性分析*

李 淵1，李云梅1**，呂 恒1，吳傳慶2，王珊珊1，王永波1

(1:江蘇省地理信息資源開發(fā)與利用協(xié)同創(chuàng)新中心，南京 210023)(2:環(huán)保部衛(wèi)星環(huán)境應(yīng)用中心，北京 100029)

太湖葉綠素a同化系統(tǒng)對于不同參數(shù)的敏感性將直接影響到該系統(tǒng)能否精確的估算太湖葉綠素a的濃度分布.利用2009年4月21日環(huán)境一號衛(wèi)星(HJ-1B CCD2)影像數(shù)據(jù)反演太湖葉綠素a濃度場信息.以此作為背景場信息，結(jié)合基于集合均方根濾波的太湖葉綠素a同化系統(tǒng)，分析和評價(jià)了樣本數(shù)目、同化時(shí)長、背景場誤差、觀測誤差和模型誤差對于同化系統(tǒng)性能的影響.結(jié)果表明：從計(jì)算成本、系統(tǒng)運(yùn)行時(shí)間和同化效果等方面分析，當(dāng)集合樣本數(shù)目達(dá)到30～40左右時(shí)同化系統(tǒng)取得了較好的結(jié)果；同化系統(tǒng)對于背景場誤差的估計(jì)變化不是很敏感，即初始場的估計(jì)是否準(zhǔn)確對于同化系統(tǒng)的性能影響不是很大；同化系統(tǒng)對于模型誤差和觀測誤差的變化較為敏感，不同的測試點(diǎn)位由于水體動力學(xué)性質(zhì)不一，其敏感性的表現(xiàn)形式有所差異；利用數(shù)據(jù)同化方法可以有效地估算太湖葉綠素a濃度.

太湖；集合均方根濾波；數(shù)據(jù)同化；葉綠素a；敏感性分析

太湖是我國第三大淡水湖，是長江中下游地區(qū)最典型的淺水湖泊，全湖水面面積2338 km2，湖泊長度69 km，平均寬度34 km，平均水深1.89m，平底水淺是太湖湖盆的一個(gè)顯著特點(diǎn)[1].近年來太湖水體富營養(yǎng)化嚴(yán)重，藍(lán)藻水華頻繁暴發(fā)，嚴(yán)重危害了區(qū)域的可持續(xù)發(fā)展.葉綠素a濃度作為水體富營養(yǎng)化和水質(zhì)評價(jià)的重要指標(biāo)，一直是水環(huán)境監(jiān)測的主要參數(shù).目前，對其進(jìn)行監(jiān)測的方法主要有地面采樣分析和遙感監(jiān)測.常規(guī)采樣監(jiān)測主要存在取樣點(diǎn)少和覆蓋面積小的缺點(diǎn)，難以反映太湖區(qū)域全面、動態(tài)的信息.利用遙感技術(shù)對葉綠素a濃度進(jìn)行監(jiān)測，有著速度快、范圍廣、相對成本低的優(yōu)勢.但是受衛(wèi)星的時(shí)間分辨率和天氣狀況的影響，遙感反演數(shù)據(jù)難以連續(xù)、動態(tài)地反映水體組分參數(shù).

數(shù)據(jù)同化方法提供了將多源數(shù)據(jù)融合，結(jié)合模型和觀測信息，進(jìn)行數(shù)據(jù)連續(xù)模擬、預(yù)測的手段.數(shù)據(jù)同化技術(shù)已經(jīng)廣泛應(yīng)用于氣象、土壤、水文、海洋、生態(tài)等多個(gè)領(lǐng)域. Pathmathevan等[2-3]進(jìn)行了土壤水分同化試驗(yàn)；Moradkhani等[4]基于Kalman濾波數(shù)據(jù)同化方法對水文模型的狀態(tài)、參數(shù)同時(shí)進(jìn)行估計(jì)；Gregg[5]利用三維全球海洋數(shù)值模型同化SeaWiFS海洋葉綠素?cái)?shù)據(jù)，有效地提高了模型結(jié)果精度；Gu等[6]利用數(shù)據(jù)同化方法重構(gòu)了高質(zhì)量的MODIS NDVI數(shù)據(jù).

數(shù)據(jù)同化方法主要分為兩類：變分法和順序同化法[7].近年來，以卡爾曼濾波為代表的順序數(shù)據(jù)同化方法在很多領(lǐng)域都得到了廣泛應(yīng)用.而基于集合論和統(tǒng)計(jì)估計(jì)理論的集合卡爾曼濾波具有程序設(shè)計(jì)相對簡單、不需要伴隨或切線性算子、可以應(yīng)用于非線性系統(tǒng)等優(yōu)點(diǎn)，克服與彌補(bǔ)了傳統(tǒng)卡爾曼濾波中的缺陷.集合卡爾曼濾波的方法首先是由Evensen在1994年引入的[8]，由于其相對四維變分而言計(jì)算相對簡單，可操作性強(qiáng)、計(jì)算成本相對較低，近年來引起了人們的廣泛關(guān)注.該方法在氣象、海洋方面都有廣泛的應(yīng)用[9].近年來，為了避免采樣誤差，確定性方法隨之產(chǎn)生，這些方法主要包括集合均方根濾波、集合調(diào)整濾波、集合變換卡爾曼濾波等[10-15].

將數(shù)據(jù)同化模型引入太湖水環(huán)境監(jiān)測中，可以利用遙感反演的水質(zhì)參數(shù)信息，結(jié)合太湖水體動力學(xué)數(shù)值模型，模擬水質(zhì)參數(shù)的擴(kuò)散和運(yùn)輸，從而連續(xù)、動態(tài)地反映太湖水質(zhì)情況.但是，模型的適用性、運(yùn)行效率等都有待于進(jìn)一步檢驗(yàn).本研究以葉綠素a濃度為監(jiān)測指標(biāo)，基于集合均方根濾波和太湖葉綠素?cái)U(kuò)散模型，利用2009年4月21日環(huán)境一號衛(wèi)星的葉綠素a濃度場反演信息，研究和分析樣本數(shù)目、同化時(shí)長、背景場誤差、觀測誤差和模型誤差對同化系統(tǒng)性能的影響.同時(shí)為同化系統(tǒng)在將來的實(shí)際應(yīng)用尋求最佳的參數(shù)組合.

1 數(shù)據(jù)與方法

1.1 研究區(qū)與樣點(diǎn)分布

圖1 太湖樣點(diǎn)分布Fig.1 The distribution of sample sites and test sites in Lake Taihu

以太湖為實(shí)驗(yàn)區(qū)，2009年4月21日在太湖布設(shè)了26個(gè)采樣點(diǎn)(圖1)，對葉綠素a濃度及水面光譜進(jìn)行了野外實(shí)地觀測.此外，本課題組在太湖湖心區(qū)布設(shè)有浮標(biāo)和平臺水質(zhì)監(jiān)測站點(diǎn)，因此在太湖湖心區(qū)隨機(jī)選取了16個(gè)同化系統(tǒng)敏感性分析的測試點(diǎn)(圖1).以這些測試點(diǎn)為例，分析和評價(jià)樣本數(shù)目、同化時(shí)長、背景場誤差、觀測誤差和模型誤差對同化系統(tǒng)性能的影響.

1.2 地面實(shí)驗(yàn)數(shù)據(jù)與方法

野外實(shí)測數(shù)據(jù)包括水面光譜的采集和26個(gè)樣點(diǎn)的葉綠素a濃度測量.

水面的反射光譜采用美國ASD公司生產(chǎn)的ASD Handheld Spectroradiometer便攜式光譜輻射計(jì)測量，其波段范圍為350～1050nm，光譜分辨率為2nm.為了減少水體鏡面反射和船體陰影的影響，更好地提取出反映水體信息的離水輻亮度Lw和遙感反射率Rrs，采用唐軍武等[16]提出的關(guān)于內(nèi)陸Ⅱ類水體水面以上光譜測量的方法.測量時(shí)，天空晴朗無云，湖面平靜，待船停穩(wěn)后，在甲板開闊處(距水面1m左右)分別測量標(biāo)準(zhǔn)灰板、遮擋標(biāo)準(zhǔn)灰板、水體和天空光的光譜信息，以上4個(gè)參數(shù)各測量10條光譜曲線，在測量水面反射光譜的同時(shí)記錄各測點(diǎn)的GPS坐標(biāo)和當(dāng)時(shí)的風(fēng)速、風(fēng)向以及時(shí)間.遙感反射率通常利用式(1)進(jìn)行計(jì)算[17]：

(1)

在光譜采集的同時(shí)采集表層水樣，低溫冷藏帶回實(shí)驗(yàn)室測量葉綠素濃度.葉綠素濃度的測量采用常規(guī)的化學(xué)分析方法，用0.45μm的GF/F濾膜過濾，90%的熱乙醇提取，然后利用分光光度計(jì)檢測，葉綠素濃度獲取的詳細(xì)步驟見文獻(xiàn)[18].

1.3 衛(wèi)星數(shù)據(jù)及預(yù)處理

環(huán)境一號衛(wèi)星是由我國于2008年9月發(fā)射升空的，包括兩顆光學(xué)衛(wèi)星(HJ-1A衛(wèi)星和HJ-1B衛(wèi)星)和一顆雷達(dá)衛(wèi)星(HJ-1C衛(wèi)星)，擁有光學(xué)、紅外、超光譜多種探測手段，具有大范圍、全天候、全天時(shí)、動態(tài)的環(huán)境和災(zāi)害監(jiān)測能力[19].伴隨著環(huán)境一號衛(wèi)星的發(fā)射升空和遙感技術(shù)的不斷發(fā)展，眾多學(xué)者利用環(huán)境一號衛(wèi)星在水體環(huán)境遙感監(jiān)測方面取得了大量的研究成果[20-22].本文利用2009年4月21日環(huán)境一號衛(wèi)星(HJ-1B CCD2)影像數(shù)據(jù)，結(jié)合野外實(shí)測葉綠素a濃度反演太湖葉綠素a濃度場信息.

圖2 反演與實(shí)測葉綠素a濃度對比Fig.2 Comparison of retrieved and measured chlorophyll-a concentrations

圖像預(yù)處理部分包括4步，第1步：輻射定標(biāo)；第2步：幾何校正；第3步：大氣校正；第4步：水體提取.其中大氣校正參考Freitas在2009年提出的校正方法.首先計(jì)算大氣頂層的遙感反射率(RTOA)，然后在地面實(shí)測點(diǎn)位中尋找在時(shí)間和空間上的準(zhǔn)同步點(diǎn)位.利用該點(diǎn)位的實(shí)測遙感反射率(RBOA)，計(jì)算同步的RTOA與RBOA之間的差值，認(rèn)為該差值是大氣貢獻(xiàn)率.假設(shè)在研究區(qū)太湖上空大氣狀況均一，從而最終將該點(diǎn)的大氣貢獻(xiàn)率應(yīng)用到整個(gè)太湖.具體校正方法參考文獻(xiàn)[23].

1.4 初始背景場的生成

本研究將使用HJ-1衛(wèi)星的葉綠素a濃度反演結(jié)果作為基于風(fēng)生流的葉綠素a擴(kuò)散模型(見1.5節(jié))的初始濃度場及太湖葉綠素a同化系統(tǒng)(見2.2節(jié))的初始背景場.為此，首先利用經(jīng)過數(shù)據(jù)預(yù)處理的環(huán)境一號衛(wèi)星(HJ-1B CCD2)影像數(shù)據(jù)，結(jié)合地面實(shí)測葉綠素a濃度，構(gòu)建葉綠素a濃度反演模型.在滿足精度需求的條件下，遵循模型易于構(gòu)建和易于業(yè)務(wù)化運(yùn)行的原則，通過模型的構(gòu)建與遴選，最終確定了單波段反演模型.構(gòu)建的單波段反演模型為：

y=-42.67lnx-139.54

(2)

式中，x為環(huán)境一號衛(wèi)星(HJ-1B)第4波段(760～900nm)的遙感反射率，y為葉綠素a濃度(μg/L).其中，野外實(shí)測26個(gè)樣點(diǎn)的葉綠素a濃度數(shù)據(jù)，隨機(jī)選取18個(gè)樣點(diǎn)數(shù)據(jù)用來構(gòu)建反演模型，剩余8個(gè)點(diǎn)用來進(jìn)行模型精度驗(yàn)證.建模R2為0.938，驗(yàn)證數(shù)據(jù)與反演結(jié)果的散點(diǎn)圖可以看出(圖2)，驗(yàn)證數(shù)據(jù)較好地分布在1 ∶1線兩側(cè)，而且RMSE為4.1985μg/L，平均相對誤差為29.3%，精度滿足實(shí)驗(yàn)需求.基于單波段反演模型，結(jié)合HJ-1B衛(wèi)星影像數(shù)據(jù)，反演了2009年4月21日太湖葉綠素a濃度，結(jié)果如圖3所示.

圖3 太湖葉綠素a濃度分布Fig.3 The distribution of chlorophyll-a concentration in Lake Taihu

1.5 基于風(fēng)生流的葉綠素a擴(kuò)散模型

假設(shè)湖水為均勻不可壓的流體，垂直方向服從靜水壓力分布，采用笛卡爾左手直角坐標(biāo)系，x軸和y軸位于湖水的平均水平面上，x軸向東為正，y軸向北為正，z軸向上為正.流體動力學(xué)方程為[24-25]：

(3)

(4)

(5)

p=ρwg(η+z)

(6)

將連續(xù)方程垂向積分可得到自由面方程：

(7)

式中，u、v、w分別為x、y、z軸上的流速分量；η為垂直方向上湖面相對于平均水面的高度；ρw是水體密度；Az和Ah分別是垂直和水平的渦粘系數(shù)；f為柯氏力，為7.23×10-5；g為重力加速度；p為水的壓強(qiáng).

葉綠素a濃度的控制方程為：

(8)

式中，c為葉綠素a濃度(μg/L).該模型的數(shù)值解法參考文獻(xiàn)[24-25]，可以有效地估算和預(yù)測葉綠素a濃度.

圖4 葉綠素同化系統(tǒng)Fig.4 Framework of chlorophyll-a assimilation

2 葉綠素?cái)?shù)據(jù)同化模型構(gòu)建

2.1 數(shù)據(jù)同化方案

本研究中的太湖葉綠素a數(shù)據(jù)同化系統(tǒng)主要由數(shù)據(jù)同化算法(基于集合均方根濾波的數(shù)據(jù)同化方法)、葉綠素?cái)U(kuò)散模擬模型、數(shù)據(jù)(參數(shù)、觀測數(shù)據(jù)、輸出數(shù)據(jù))及誤差分析4部分組成(圖4).

2.2 集合均方根濾波

卡爾曼濾波是一種典型的順序數(shù)據(jù)同化方法，1960年，Kalman[26]針對隨機(jī)過程狀態(tài)估計(jì)提出了卡爾曼濾波的思想.Evensen等在1994年提出了集合卡爾曼濾波算法[8]，克服了傳統(tǒng)卡爾曼濾波的不足，并且具有計(jì)算成本低、運(yùn)算效率高等優(yōu)點(diǎn)[27].但是在標(biāo)準(zhǔn)的集合卡爾曼濾波中，觀測值是被當(dāng)做隨機(jī)變量處理的，所以會在觀測值上加擾動.近年來，為了避免采樣誤差，確定性方法隨之產(chǎn)生. Whitaker等[13]提出確定性觀測法的集合均方根濾波就是其中之一.該方法不需要對觀測值進(jìn)行擾動，減少了誤差的引入.該算法利用蒙特卡羅方法的思想，用符合高斯分布的一組隨機(jī)變量(設(shè)數(shù)目為N)去代表隨機(jī)動態(tài)預(yù)報(bào)中的概率密度函數(shù)，通過向前積分，計(jì)算下一個(gè)時(shí)刻狀態(tài)總體的概率密度，并得到該時(shí)刻的統(tǒng)計(jì)特性(如均值與協(xié)方差).

具體算法如下：

② 計(jì)算卡爾曼增益矩陣Kn：

(9)

(10)

(11)

(12)

在求解增益矩陣的過程中，采用了通過樣本來計(jì)算增益矩陣的方法，這樣可以達(dá)到降低計(jì)算機(jī)存儲和提高運(yùn)算效率的目的.在這里觀測算子H可以不是線性的.

(13)

(14)

(15)

(16)

④ 進(jìn)行預(yù)報(bào)，也就是進(jìn)行狀態(tài)的時(shí)間更新：

(17)

式中，ηn-1(k)～N(0, Qn-1), k=1, 2, …, K， M代表葉綠素a擴(kuò)散模擬模型.

⑤ 回到步驟②，如此循環(huán)往復(fù)，直至沒有新的觀測資料加入或已滿足同化時(shí)長需求.

2.3 誤差統(tǒng)計(jì)

對于同化后的結(jié)果，采用兩種誤差分析方法：(1) 均方根誤差(rootmeansquareerror, RMSE)；(2) 相對誤差(relativeerror, RE).

(18)

(19)

式中,N為整個(gè)觀測時(shí)間，obsi為i時(shí)刻的“真值”，Xi為i時(shí)刻的葉綠素a擴(kuò)散模擬值或同化結(jié)果.

3 敏感性分析

從集合均方根濾波同化算法的設(shè)計(jì)過程可知，同化系統(tǒng)的性能會受到背景場誤差、觀測誤差、樣本數(shù)目和模型誤差等因素的影響[28].為了研究這些參數(shù)對于同化結(jié)果的影響，我們將通過調(diào)整參數(shù)的變化，來觀察這些參數(shù)對于同化系統(tǒng)性能的影響.目前，太湖葉綠素a同化系統(tǒng)基于Matlab軟件實(shí)現(xiàn)，水體動力學(xué)模型是基于Fortran語言封裝，在系統(tǒng)運(yùn)行時(shí)，由Matlab調(diào)用水體動力學(xué)模型.

本次實(shí)驗(yàn)共分為5部分，分別是樣本數(shù)目、同化時(shí)長、背景誤差、觀測誤差和模型誤差的敏感性分析.敏感性分析實(shí)驗(yàn)利用2009年4月21日環(huán)境一號衛(wèi)星的葉綠素a濃度反演結(jié)果作為葉綠素a擴(kuò)散模型的初值，結(jié)合葉綠素a擴(kuò)散模型，在模擬時(shí)長內(nèi)每小時(shí)輸出一次模擬結(jié)果，并將該結(jié)果作為理論意義上的“真值”.因此，當(dāng)我們把敏感性分析實(shí)驗(yàn)的同化結(jié)果與上述結(jié)果進(jìn)行對比分析后，即可判斷不同參數(shù)及其組合對太湖葉綠素a同化系統(tǒng)性能的影響.為了避免隨機(jī)性誤差因素的影響，關(guān)于每個(gè)變量的同化實(shí)驗(yàn)都重復(fù)了10次，以消除隨機(jī)性因素的影響.最終的誤差分析采用平均均方根誤差(ARMSE)，即以10次的均方根誤差的均值為評價(jià)標(biāo)準(zhǔn).

3.1 模型誤差敏感性分析

同化實(shí)驗(yàn)參數(shù)設(shè)置為：葉綠素a擴(kuò)散模型誤差：最小值為1%，步長1%，最大值為10%；同化時(shí)長12h；觀測誤差1%；背景場誤差30%；樣本數(shù)目為25.模型誤差主要是由于模型本身無法精確描述和刻畫水體動力學(xué)的變化而引起的系統(tǒng)誤差，同時(shí)由于邊界條件、模型參數(shù)等的不確定也會引起模型誤差.圖5反映了模型誤差對于同化系統(tǒng)的影響.不同點(diǎn)位對于模型誤差敏感性是不同的，出現(xiàn)最小的ARMSE的模型誤差的值各不相同(圖5).在集合卡爾曼濾波的時(shí)間更新方程中，如果忽略模型誤差一項(xiàng)會導(dǎo)致低估預(yù)報(bào)誤差協(xié)方差矩陣[29]，從而引起預(yù)報(bào)的偏差.但本身對于模型誤差的估計(jì)也是很難的一件事情.從敏感性分析的結(jié)果中可以看出，大部分的點(diǎn)位當(dāng)模型誤差為4%～6%時(shí)，ARMSE最小.這也說明，對于葉綠素a擴(kuò)散模型而言，當(dāng)模型誤差在4%～6%之間的時(shí)候，同化效果較好.當(dāng)模型誤差較小時(shí)，會導(dǎo)致對預(yù)報(bào)誤差協(xié)方差矩陣的低估；當(dāng)模型誤差較大時(shí)，又會導(dǎo)致預(yù)報(bào)精度下降.

圖5 模型誤差敏感性分析結(jié)果Fig.5 Sensitivity analysis results of the model error

圖6 樣本數(shù)目敏感性分析結(jié)果Fig.6 Sensitivity analysis results of the number of the ensemble

3.2 樣本數(shù)目敏感性分析

同化實(shí)驗(yàn)參數(shù)設(shè)置為：樣本數(shù)目：最小值為10，步長5，最大值為50；背景場誤差30%；觀測誤差1%；葉綠素a擴(kuò)散模型誤差為5%；同化時(shí)長為12h.樣本數(shù)目的多少決定了通過這些樣本是否能夠更加準(zhǔn)確地反映狀態(tài)變量的空間分布，而樣本數(shù)目過多又會增加系統(tǒng)運(yùn)行的時(shí)間，所以分析樣本數(shù)目對同化系統(tǒng)的敏感性就變得十分重要.考慮到時(shí)間的計(jì)算成本，目前實(shí)驗(yàn)中最大樣本數(shù)目設(shè)計(jì)為50.樣本的敏感性分析實(shí)驗(yàn)結(jié)果可以看出，隨著集合樣本數(shù)目的增加，ARMSE總體上呈現(xiàn)逐漸減小的趨勢(圖6).當(dāng)樣本數(shù)目在10～25附近時(shí)，ARMSE會出現(xiàn)不同程度的輕微波動，當(dāng)集合數(shù)目大于30后，ARMSE總體減小的趨勢就很明顯了.從系統(tǒng)運(yùn)行時(shí)間和同化效果來看，當(dāng)集合樣本數(shù)目在30～40左右時(shí)，系統(tǒng)的運(yùn)行效率和同化效果都取得了較好的效果.關(guān)于集合樣本數(shù)目敏感性分析這一結(jié)果也符合了集合均方根濾波通過蒙特卡洛方法解決非線性問題的原理與方法.即大量樣本可以更加準(zhǔn)確地反映狀態(tài)變量的空間分布，更加準(zhǔn)確地反映和刻畫狀態(tài)變量在真實(shí)空間中的分布.但大量樣本也會導(dǎo)致同化系統(tǒng)運(yùn)行時(shí)間過長，增加計(jì)算成本.

3.3 同化時(shí)長敏感性分析

同化實(shí)驗(yàn)參數(shù)設(shè)置為：同化時(shí)長：最小值為2h，步長2h，最大值為12h；背景場誤差30%；觀測誤差1%；葉綠素a擴(kuò)散模型誤差為5%；樣本數(shù)目為25.同化時(shí)長決定了同化觀測數(shù)據(jù)的周期長短與同化結(jié)果之間的關(guān)系，較長的同化時(shí)長需要更多的觀測數(shù)據(jù)，同時(shí)系統(tǒng)的運(yùn)行時(shí)間也會更長，但卻能得到較好的結(jié)果.同化時(shí)長的確定往往與觀測數(shù)據(jù)的采集頻率、數(shù)據(jù)處理系統(tǒng)的配置等因素有關(guān).同化時(shí)長的敏感性分析結(jié)果可以看出(圖7)，隨著同化時(shí)長的增加，ARMSE逐漸減小.其中，10#點(diǎn)位出現(xiàn)了一個(gè)上升的波動，其余點(diǎn)位都呈現(xiàn)整體下降的趨勢. 10#點(diǎn)位在2～6h平均均方根誤差上升，說明在該點(diǎn)位該同化時(shí)長的情況下，同化濾波出現(xiàn)了濾波發(fā)散的不穩(wěn)定情況.說明在某些時(shí)段，同化時(shí)長因素對不同湖區(qū)的影響具有差異性.但整體趨勢還是隨著同化時(shí)長的增加，ARMSE逐漸減小.

圖7 同化時(shí)長敏感性分析結(jié)果Fig.7 Sensitivity analysis results of the assimilation time

3.4 背景場誤差敏感性分析

同化實(shí)驗(yàn)參數(shù)設(shè)置為：背景場誤差：最小值為20%，步長2%，最大值為38%；同化時(shí)長12h；觀測誤差1%；葉綠素a擴(kuò)散模型誤差為5%；樣本數(shù)目為25.背景場誤差是在HJ-1衛(wèi)星反演結(jié)果基礎(chǔ)上加入了滿足高斯正態(tài)分布的誤差擾動.背景場誤差反映了對狀態(tài)變量空間的初始誤差估計(jì)，決定了對狀態(tài)空間的初始估計(jì).背景場誤差敏感性分析結(jié)果可以看出(圖8)，當(dāng)背景場誤差被低估的時(shí)候，即背景場誤差為20%～28%時(shí)，ARMSE基本都處于一個(gè)平穩(wěn)的波動狀態(tài)，沒有明顯的變化趨勢.當(dāng)背景場誤差被高估的時(shí)候，即背景場誤差為32%～38%時(shí)，不同的點(diǎn)位表現(xiàn)得有所不同.以3#、6#、11#、16#點(diǎn)位為代表，當(dāng)背景場誤差為32%～38%時(shí)，這些點(diǎn)位的ARMSE都呈上升趨勢.以4#、5#、8#、9#點(diǎn)位為代表，當(dāng)背景場誤差為32%～38%時(shí)，這些點(diǎn)位的平均均方根誤差仍表現(xiàn)為一平穩(wěn)波動狀態(tài).但就整體而言，這些點(diǎn)位在背景場誤差估計(jì)存在偏差的情況下，ARMSE變化不大.隨著背景場誤差的增大，ARMSE呈現(xiàn)一個(gè)有規(guī)律的波動狀態(tài).而這也符合順序同化的特點(diǎn)，即使初始場的狀態(tài)變量的估計(jì)存在較大誤差，同化系統(tǒng)經(jīng)過一段時(shí)間也會達(dá)到穩(wěn)定.

圖8 背景場誤差敏感性分析結(jié)果Fig.8 Sensitivity analysis results of the background error

3.5 觀測誤差敏感性分析

同化實(shí)驗(yàn)參數(shù)設(shè)置為：觀測誤差：最小值為1%，步長1%，最大值為10%；同化時(shí)長12h；背景場誤差30%；葉綠素a擴(kuò)散模型誤差為5%；樣本數(shù)目為25.觀測誤差是在模擬“真值”結(jié)果基礎(chǔ)上加入了滿足高斯正態(tài)分布的誤差擾動.觀測誤差主要包括儀器采集數(shù)據(jù)時(shí)的儀器誤差和代表性誤差.觀測誤差敏感性分析結(jié)果可以看出(圖9)，隨著觀測誤差的增加，ARMSE主要分為兩種情況.其一，以1#、2#、4#、5#、9#、16#點(diǎn)位為代表，隨著觀測誤差的增加，ARMSE也在增加，整體呈上升的趨勢.其二，剩余的點(diǎn)位隨著觀測誤差的增加，ARMSE圍繞某一均值處于一個(gè)波動狀態(tài).說明位于太湖不同區(qū)域的點(diǎn)位，由于水體動力學(xué)性質(zhì)的不同，其同化效果隨觀測誤差的變化而異.

圖9 觀測誤差敏感性分析結(jié)果Fig.9 Sensitivity analysis results of the observation error

4 太湖葉綠素a同化實(shí)驗(yàn)?zāi)M

模擬實(shí)驗(yàn)基于2009年4月21日環(huán)境一號衛(wèi)星的反演葉綠素?cái)?shù)據(jù)作為葉綠素?cái)U(kuò)散模型的“真實(shí)”初值.將該初值在葉綠素?cái)U(kuò)散模型下進(jìn)行12h模擬，每1小時(shí)輸出一次測試點(diǎn)位的模擬結(jié)果.同時(shí)在“真實(shí)”初值場上進(jìn)行30%的擾動，作為同化系統(tǒng)的初始背景場信息.利用之前分析的敏感性參數(shù)，采用觀測誤差為1%，模型誤差為5%，背景場誤差為30%，集合樣本數(shù)為50，同化時(shí)長12h參數(shù)設(shè)置同化系統(tǒng)，分析和評價(jià)太湖葉綠素a同化系統(tǒng)的同化效果.同化12h后的結(jié)果與“真實(shí)”初值在葉綠素a擴(kuò)散模型作用下12h的預(yù)測結(jié)果作相對誤差分析可以看出，全湖范圍內(nèi)大部分的區(qū)域都將誤差控制在30%以下(圖10).由于在數(shù)值計(jì)算時(shí)，將全湖離散為2813個(gè)離散的網(wǎng)格點(diǎn)，通過統(tǒng)計(jì)得到全湖有56.5%網(wǎng)格點(diǎn)的誤差控制在30%以下.即通過同化16個(gè)網(wǎng)格點(diǎn)數(shù)據(jù)，使得全湖56.5%的區(qū)域誤差得到了改善.從同化的16個(gè)點(diǎn)位數(shù)據(jù)來看，相對誤差分別為16.3%、4.7%、16.7%、13.7%、3.5%、18.4%、19.9%、15.1%、13.8%、25.9%、11.4%、17.2%、11.9%、20.6%、16.7%、18.1%.從單點(diǎn)的同化效果來看，精度平均提升了49.2%.但是，我們也發(fā)現(xiàn)在太湖西部部分區(qū)域同化效果不理想，一方面可能是由于目前同化系統(tǒng)還不夠完善，在系統(tǒng)運(yùn)行效率和計(jì)算成本方面有待進(jìn)一步改善；另一方面，由于模擬實(shí)驗(yàn)中數(shù)據(jù)同化時(shí)長在1～12h之內(nèi)，因此，未考慮藻類生長模型及水質(zhì)模型.但是，從長時(shí)間尺度來看，葉綠素a濃度變化還會受到水溫、營養(yǎng)鹽含量、光照等因素的影響.因此，為了提高數(shù)據(jù)同化精度，在將來進(jìn)一步的研究中還需要在基于風(fēng)生流的太湖水體動力學(xué)模型基礎(chǔ)上耦合藻類生長模型，從而更加全面準(zhǔn)確地模擬和估算太湖葉綠素a濃度.就目前實(shí)際情況而言，由于太湖東部水質(zhì)狀況較好，所以大部分的浮標(biāo)站點(diǎn)位于太湖中部和北部，其數(shù)目不少于20個(gè).因此，利用這些浮標(biāo)站點(diǎn)的連續(xù)動態(tài)觀測數(shù)據(jù)，利用基于集合均方根的數(shù)據(jù)同化方法對于太湖葉綠素a濃度的估算是有效可行的.

圖10 全湖相對誤差分析Fig.10 Relative error analysis of Lake Taihu

5 結(jié)論

1) 集合樣本數(shù)目越多、同化時(shí)長越長都可以越有效地降低同化系統(tǒng)的誤差，但也會增加系統(tǒng)的計(jì)算成本和運(yùn)算時(shí)間.當(dāng)集合樣本數(shù)目在30～40左右時(shí)，系統(tǒng)的運(yùn)行效率和同化效果都取得了較好的效果.

2) 背景場誤差對于同化系統(tǒng)的影響不是很大，即使背景場誤差估計(jì)偏差較大，同化系統(tǒng)的平均均方根誤差變化也不是很大.

3) 觀測誤差對于同化系統(tǒng)的影響較大，不同的點(diǎn)位呈現(xiàn)出的趨勢有所不同.有的點(diǎn)位隨著觀測誤差增大，同化系統(tǒng)的誤差也會增大；有的點(diǎn)位則出現(xiàn)圍繞某一均值上下波動.

4) 模型誤差對于基于集合均方根濾波的數(shù)據(jù)同化方法較為重要，準(zhǔn)確估計(jì)模型誤差對于提高同化系統(tǒng)的性能至關(guān)重要.對于葉綠素a擴(kuò)散模型而言，當(dāng)模型誤差在4%～6%之間時(shí)，同化效果較好.

[1] 秦伯強(qiáng),胡維平,陳偉民.太湖水環(huán)境演化過程與機(jī)理.北京:科學(xué)出版社,2004:1-296.

[2] Pathmathevan M, Koilke T, Li Xetal. A simplified land data assimilation scheme and its application to soil moisture experiments in 2002(SMEX02).WaterResourcesResearch, 2003,39(12):1341. doi: 10.1029/2003WR002124.

[3] Pathmathevan M, Koilke T, Li X. A new satellite-based data assimilation algorithm to determine spatial and temporal variations of soil moisture and temperature profiles.JournaloftheMeteorologicalSocietyofJapan, 2003,81(5):1111-1135.

[4] Moradkhani H, Sorooshian S, Gupta HVetal. Dual state-parameter estimation ofhydrological models using ensemble Kalman filter.AdvancesinWaterResources, 2005,28(2):135-147.

[5] Gregg WW. Assimilation of seawifs ocean chlorophyll data into a three-dimensional global ocean model.JournalofMarineSystems, 2008,69(3/4):205-225.

[6] Gu J, Li X, Huang CLetal. A simplified data assimilation method for reconstructing time-series MODIS NDVI data.AdvancesinSpaceResearch, 2009,44(4):501-509.

[7] Anderson LA, Robinson AR, Lozano CJ. Physical and biological modeling in the Gulf Stream region: I. Data assimilation methodology.Deep-seaResearchPartI, 2000,47:1787-1827.

[8] Evensen G. Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics.JournalofGeophysicalResearch, 1994,99(C5):10143-10162.

[9] Evensen G. The ensemble Kalman filter: Theoretical formulation and practical implementation.OceanDynamics, 2003,53(4):343-367.

[10] Burgers G, Leeuwen PJV, Evensen G. Analysis scheme in the ensemble Kalman filter.MonthlyWeatherReview, 1998,126(6):1719-1724.

[11] Houtekamer PL, Mitchell HL. Data assimilation using an ensemble Kalman filter technique.MonthlyWeatherReview, 1998,126(3):796-811.

[12] Houtekamer PL, Mitchell HL. A sequential ensemble Kalman filter for atmospheric data assimilation.MonthlyWeatherReview, 2001,129(1):123-137.

[13] Whitaker JS, Hamill TM. Ensemble data assimilation without perturbed observations.MonthlyWeatherReview, 2002,130(7):1913-1924.

[14] Anderson JL. An ensemble adjustment Kalman filter for data assimilation.MonthlyWeatherReview, 2001,129(12):2884-2903.

[15] Tippett MK, Anderson JL, Bishop CHetal. Ensemble square root filters.MonthlyWeatherReview, 2003,131(7):1485-1490.

[16] 唐軍武.海洋光學(xué)特性模擬與遙感模型[學(xué)位論文].北京:中國科學(xué)院遙感應(yīng)用研究所,1999:107-110.

[17] 唐軍武，田國良，汪小勇等.水體光譜測量與分析Ⅰ水面以上測量法.遙感學(xué)報(bào),2004,8(1)：37-44.

[18] Le CF, Li YM, Zha Yetal. Specific absorption coefficient and the phytoplankton package effect in Lake Taihu, China.Hydrobiologia, 2009,619(1):27-37.

[19] 王 橋,吳傳慶,厲 青.環(huán)境一號衛(wèi)星及其在環(huán)境監(jiān)測中的應(yīng)用.遙感學(xué)報(bào),2010,14(1):104-121.

[20] 徐祎凡,李云梅,王 橋等.基于環(huán)境一號衛(wèi)星多光譜影像數(shù)據(jù)的三湖一庫富營養(yǎng)化狀態(tài)評價(jià).環(huán)境科學(xué)學(xué)報(bào),2011,31(1):81-92.

[21] 曠 達(dá),韓秀珍,劉 翔等.基于環(huán)境一號衛(wèi)星的太湖葉綠素a濃度提取.中國環(huán)境科學(xué),2010,30(9)：1269-1273.

[22] 金 焰,張 詠,牛志春等.環(huán)境一號CCD數(shù)據(jù)在太湖藍(lán)藻水華遙感監(jiān)測中的應(yīng)用.環(huán)境監(jiān)測管理與技術(shù)，2010,22(5):53-54.

[23] Freitas FH. Spectral merging of MODIS/MERIS Ocean Colour Data to improve monitoring of coastal water processes [Dissertation]. Enschede: International Institute for Geo-Information Science and Earth Observation, 2009.

[24] Zhang Z, Song ZY, Lv GN. A new implicit scheme for solving 3-D shallow water flows.JournalofHydrodynamics(SeriesB), 2009,21(6):790-798.

[25] Zhang Z, Song ZY. Three-dimensional numerical modeling for wind-driven circulation and pollutant transport in a large scale lake. International Conference of Bioinformatics and Biomedical Engineering. Chengdu: IEEE, 2010:1-7.

[26] Kalman RE. A new approach to linear filtering and prediction problems.JournalofFluidsEngineering, 1960,82(1):35-45.

[27] 劉成思.集合卡爾曼濾波資料同化方案的設(shè)計(jì)和研究[學(xué)位論文].北京:中國氣象科學(xué)研究院,2005:13-14.

[28] 黃春林,李 新.土壤水分同化系統(tǒng)的敏感性試驗(yàn)研究.水科學(xué)進(jìn)展,2006,17(4):457-465.

[29] 曾忠一.大氣科學(xué)中的反問題(上冊).臺北:國立編譯館,2006.

Sensitivity analysis on Lake Taihu data assimilation scheme of chlorophyll-a concentration

LI Yuan1, LI Yunmei1, Lü Heng1, WU Chuanqing2, WANG Shanshan1& WANG Yongbo1

(1:JiangsuCenterforCollaborativeInnovationinGeographicalInformationResourceDevelopmentandApplication,Nanjing210023,P.R.China)(2:SatelliteEnvironmentApplicationCenter,MinistryofEnvironmentalProtection,Beijing100029,P.R.China)

Sensibility of the Lake Taihu chlorophyll-a assimilation system to different parameters directly control the accuracy of estimate the chlorophyll-a concentration distribution when using this assimilation system. We used multispectral data of Environmental Satellite 1(HJ-1), obtained on April 21st, 2009, combined withinsitudata to retrieve the concentration of chlorophyll-a in Lake Taihu. We developed a Lake Taihu chlorophyll-a data assimilation system based on ensemble square root Klaman filter(EnSRF) technique. Take the retrieved chlorophyll-a concentration of Lake Taihu as the initial background value, then combined with the data assimilation system to analyze the influence of the ensemble size, the assimilation time, the background error, the observation error and the model error on the assimilation system. The results indicate: taking the computing cost, time cost of system and the performance of the assimilation system into consideration, the assimilation system performs well when the ensemble size are 30-40; the assimilation system is not very sensitive to the accuracy of estimation of the background; both the observation and the model errors are very important for the performance of the system; different test stations have different water dynamic properties, so they have different performance; the estimation of chlorophyll-a concentration can be improved by using the data assimilation method.

Lake Taihu; ensemble square root filter; data assimilation; chlorophyll-a; sensitivity analysis

*國家自然科學(xué)基金項(xiàng)目(41271343)、江蘇省高校自然科學(xué)研究重大項(xiàng)目(11KJA170003)和中國科學(xué)院數(shù)字地球重點(diǎn)實(shí)驗(yàn)室開放基金項(xiàng)目(2012LDE009)聯(lián)合資助. 2014-03-03收稿；2014-06-12收修改稿.李淵(1985～)，男，博士研究生；E-mail: liyuannjnu@163.com.

**通信作者；E-mail:yunmei2009@gmail.com.