高永麗, 王際朝, 孫國(guó)棟, 張 坤, 姜向陽(yáng), 王 寧

DCM模擬中基于OPSA方法的參數(shù)敏感性分析

高永麗1, 王際朝1, 孫國(guó)棟2, 張 坤3, 姜向陽(yáng)4, 王 寧4

(1. 中國(guó)石油大學(xué)(華東), 山東 青島 266580; 2. 中國(guó)科學(xué)院大氣物理研究所, 北京 100029; 3. 中國(guó)科學(xué)院海洋研究所, 山東 青島 266071; 4. 山東省海洋資源與環(huán)境研究院, 山東 煙臺(tái) 264006)

深層葉綠素最大值(deep chlorophyll maximum, DCM)現(xiàn)象是海洋與湖泊中普遍存在的生態(tài)現(xiàn)象。對(duì)其進(jìn)行數(shù)值模擬時(shí), 參數(shù)不確定性是導(dǎo)致模擬結(jié)果出現(xiàn)誤差的重要原因?；谝粋€(gè)經(jīng)典海洋生態(tài)模式(nutrients-phytoplankton model, NP), 本文通過(guò)最優(yōu)參數(shù)敏感性分析(optimization parameter sensitivity analysis, OPSA)方法探討了模式參數(shù)不確定性對(duì)DCM模擬的影響。研究表明, 背景場(chǎng)渾濁度、垂向湍流擴(kuò)散系數(shù)、浮游植物營(yíng)養(yǎng)鹽含量和硝酸鹽再循環(huán)系數(shù)為模式中的敏感參數(shù), 它們的擾動(dòng)將導(dǎo)致DCM模擬發(fā)生顯著改變。進(jìn)一步, 設(shè)計(jì)觀測(cè)系統(tǒng)模擬試驗(yàn)評(píng)估了消除敏感參數(shù)誤差DCM模擬的改進(jìn)程度。結(jié)果顯示, 去除4個(gè)敏感參數(shù)誤差DCM模擬平均改進(jìn)了56.83%, 約是去除不敏感參數(shù)誤差平均改進(jìn)程度(4.51%)的13倍。而且, 去除敏感參數(shù)誤差模擬改進(jìn)的穩(wěn)定性更好, 變異系數(shù)僅為9.44%, 去除不敏感參數(shù)誤差模擬改進(jìn)的變異系數(shù)達(dá)到了14.76%, 穩(wěn)定性較差。據(jù)此, 可優(yōu)先發(fā)展與敏感參數(shù)直接相關(guān)的動(dòng)力過(guò)程參數(shù)化方案, 或在有限的觀測(cè)資源下優(yōu)先對(duì)敏感參數(shù)展開目標(biāo)觀測(cè), 進(jìn)而為提高DCM模擬與預(yù)測(cè)提供科學(xué)指導(dǎo)。

DCM; 參數(shù)敏感性; 最優(yōu)擾動(dòng); OPSA方法; 目標(biāo)觀測(cè)

在海洋、湖泊等大部分水體中, 葉綠素垂向分布的最大值并不出現(xiàn)在水體表面, 而是出現(xiàn)在水面下一定深度的地方[1], 這就是DCM(deep chlorophyll maximum)現(xiàn)象。早在1965年, Yentsch[2]就發(fā)現(xiàn)在印度洋相對(duì)穩(wěn)定的水體結(jié)構(gòu)中, 物理與生物的耦合作用使得浮游植物在真光層底部出現(xiàn)最大值。隨后, 更多研究證實(shí)了該現(xiàn)象在其他大洋、湖泊, 甚至混合作用較弱的河口區(qū)也普遍存在[3]。DCM現(xiàn)象體現(xiàn)了浮游植物對(duì)光和營(yíng)養(yǎng)鹽的適應(yīng)性, 研究DCM現(xiàn)象對(duì)探究海洋生態(tài)系統(tǒng)的結(jié)構(gòu)與功能具有十分重要的意義, 因此一直是國(guó)內(nèi)外海洋生態(tài)動(dòng)力學(xué)領(lǐng)域關(guān)注的熱點(diǎn)之一[4-5]。

隨著高性能計(jì)算的快速發(fā)展, 數(shù)值模式逐漸成為海洋環(huán)流及海洋生態(tài)動(dòng)力學(xué)研究的重要工具[6-7]。與觀測(cè)相比, 數(shù)值模擬結(jié)果通常存在偏差。其中, 模式參數(shù)誤差是導(dǎo)致模擬結(jié)果不確定性的主要原因之一。海洋生態(tài)模式中含有大量與描述物理、生物化學(xué)過(guò)程相關(guān)的參數(shù)。這些參數(shù)的預(yù)設(shè)值只有很少一部分是由實(shí)際觀測(cè)獲得, 其余大部分是通過(guò)經(jīng)驗(yàn)、半經(jīng)驗(yàn)或者統(tǒng)計(jì)方法得到的估計(jì)值, 因此具有較大的不確定性[8]。通過(guò)參數(shù)敏感性分析探究模式中參數(shù)不確定性影響到底如何, 識(shí)別相對(duì)敏感的參數(shù), 并且集中有限的人力物力資源優(yōu)先對(duì)敏感參數(shù)展開觀測(cè), 對(duì)提高數(shù)值模擬技巧是非常重要的。

然而, 圍繞模式參數(shù)的敏感性分析是一個(gè)復(fù)雜的問(wèn)題。早期研究中, 人們普遍使用的是數(shù)據(jù)同化方法, 例如非線性優(yōu)化技術(shù)[9], 伴隨方法[10]和弱約束參數(shù)估計(jì)方法[11]等。此類方法通過(guò)調(diào)整參數(shù)值使得模式輸出與觀測(cè)數(shù)據(jù)的偏差達(dá)到最小, 并以此確定最優(yōu)參數(shù)值。在此基礎(chǔ)上, 根據(jù)參數(shù)誤差與模式輸出變化之間的關(guān)系, 分析參數(shù)的敏感性。這些方法需要提供狀態(tài)變量和相關(guān)參數(shù)的初猜值, 且未能考慮參數(shù)間的相互作用, 因此具有一定的局限性。Chu等[12]發(fā)展了一種基于方差的非線性方法, 通過(guò)定義的敏感性指標(biāo)刻畫了浮游植物生物量對(duì)12個(gè)模式參數(shù)的相對(duì)敏感性, 并利用二階敏感性指標(biāo)分析了參數(shù)兩兩組合的非線性相互作用對(duì)模擬結(jié)果的影響。近年來(lái), Mu等[13]基于對(duì)海洋模式中初始誤差的研究, 又提出了條件非線性最優(yōu)參數(shù)擾動(dòng)方法。該方法考慮了參數(shù)擾動(dòng)間的非線性相互作用, 可識(shí)別那些在目標(biāo)時(shí)刻對(duì)所關(guān)注的自然事件的模擬或預(yù)報(bào)產(chǎn)生最大影響的重要參數(shù)或參數(shù)組合[14-15]。

事實(shí)上, 在一些數(shù)值模式(特別是生態(tài)模式)中, 模式參數(shù)可能多達(dá)10~50個(gè)。在此情形下, 使用以上方法識(shí)別敏感參數(shù), 并考察參數(shù)擾動(dòng)間的非線性相互作用時(shí), 將耗費(fèi)相當(dāng)大的計(jì)算量。為此, Wang等[16]提出了優(yōu)化參數(shù)敏感性分析方法(OPSA), 該方法旨在辨別出“累計(jì)”影響都非常小的參數(shù)擾動(dòng)(不敏感參數(shù)), 進(jìn)而識(shí)別出相對(duì)敏感的參數(shù)。該方法在識(shí)別敏感參數(shù)時(shí)能大大節(jié)約計(jì)算資源, 目前已成功應(yīng)用于黑潮延伸體模態(tài)轉(zhuǎn)變過(guò)程等研究[17]。

本文基于OPSA方法和經(jīng)典海洋生態(tài)模式(nut-rients-phytoplankton model), 識(shí)別出四個(gè)相對(duì)敏感的參數(shù), 并利用敏感性試驗(yàn)驗(yàn)證了去除敏感參數(shù)誤差確實(shí)能有效提高模擬技巧, 最后討論了基于敏感參數(shù)優(yōu)先實(shí)施目標(biāo)觀測(cè)的可靠性。

1 方法與模式介紹

1.1 OPSA方法

假設(shè)以()為狀態(tài)變量的數(shù)值模式如下:

()=M()(0) (1)

則去除擾動(dòng)中參數(shù)p的擾動(dòng)p(1≤≤)后,()的改變量如下:

假設(shè)擾動(dòng)=(1,2, …,p)的約束范圍為ε, 構(gòu)建如下最優(yōu)化問(wèn)題:

利用最優(yōu)化算法, 如遺傳算法(differential evolution)、譜梯度算法(spectral gradient algorithm)等求解上述最優(yōu)化問(wèn)題, 即得導(dǎo)致模式誤差出現(xiàn)最大值J的一類參數(shù)擾動(dòng)(p,p,…, p)。

接下來(lái), 考察多個(gè)參數(shù)同時(shí)發(fā)生擾動(dòng)引起的模式誤差。令:

則參數(shù)擾動(dòng)p與p(1≤≤)導(dǎo)致的“聯(lián)合誤差”為:

類似地, 可定義三個(gè)或更多參數(shù)同時(shí)發(fā)生擾動(dòng)導(dǎo)致的模式誤差。又因?yàn)?

, (8)

故:

由公式(9)可知, 考慮參數(shù)擾動(dòng)間的非線性相互作用時(shí), 任意兩個(gè)參數(shù)擾動(dòng)導(dǎo)致的最大模式誤差一定不超過(guò)兩個(gè)單參數(shù)擾動(dòng)引起的最大誤差之和。因此, 在用OPSA方法進(jìn)行敏感性分析時(shí), 只需求解針對(duì)單參數(shù)擾動(dòng)構(gòu)建的非線性優(yōu)化系統(tǒng), 多參數(shù)擾動(dòng)導(dǎo)致的模式誤差可用相應(yīng)單參數(shù)擾動(dòng)導(dǎo)致的誤差累積來(lái)度量, 這就省去了優(yōu)化多參數(shù)擾動(dòng)的計(jì)算量, 大大節(jié)約了計(jì)算資源。

1.2 NP模式

NP模式是描述海洋生態(tài)系統(tǒng)動(dòng)力過(guò)程的一維垂向數(shù)值模式, 不考慮物質(zhì)的水平運(yùn)輸, 常用來(lái)研究海洋和湖泊中浮游植物的垂向分布趨勢(shì)[18-19], 特別是DCM現(xiàn)象的模擬。Huisman等[20]在研究DCM振蕩與混沌現(xiàn)象時(shí)指出, 雖然NP模式只是一個(gè)理論模式, 但它呈現(xiàn)了許多與真實(shí)海洋環(huán)境中DCM現(xiàn)象一致的特征。例如, 模式模擬的DCM深度在100 m左右, 且DCM所在位置以上的營(yíng)養(yǎng)鹽濃度接近于0、以下的營(yíng)養(yǎng)鹽濃度是線性增加的, 這些表現(xiàn)與Klausmeier和Litchman對(duì)觀測(cè)資料的分析得出的結(jié)論是一致的[21]。因此, 雖然海洋生態(tài)系統(tǒng)的數(shù)值模式越來(lái)越呈現(xiàn)出復(fù)雜化的趨勢(shì), NP模式仍然經(jīng)常被用來(lái)研究葉綠素的垂向分布及其對(duì)物理過(guò)程的響應(yīng)機(jī)制[22]。

在NP模式描述的非線性系統(tǒng)中, 浮游植物生長(zhǎng)受到光的強(qiáng)度與營(yíng)養(yǎng)鹽濃度的限制, 浮游植物()和營(yíng)養(yǎng)鹽()的動(dòng)力學(xué)過(guò)程滿足如下的一維反應(yīng)擴(kuò)散方程:

其中,代表水柱的深度,代表時(shí)間,代表垂向湍流擴(kuò)散系數(shù)。和分別代表浮游植物死亡率和下沉速率,和分別代表浮游植物細(xì)胞的營(yíng)養(yǎng)鹽含量和死掉的浮游植物參與營(yíng)養(yǎng)鹽再循環(huán)的系數(shù),H和H分別代表營(yíng)養(yǎng)鹽和光的半飽和常數(shù)。諸參數(shù)值見表1。

表1 NP模式中的參數(shù)設(shè)置

另外, 來(lái)自太陽(yáng)輻射的入射光在水面下的傳播滿足Lambert-Beer’s方程[4]:

其中,in代表入射光在水面的光強(qiáng)。

假設(shè)浮游植物和營(yíng)養(yǎng)鹽在模式的上下邊界都沒(méi)有通量交換, 且營(yíng)養(yǎng)鹽在水柱的底部保持不變:

此處,N表示營(yíng)養(yǎng)鹽在水柱底部Z處的值, 本文中=10 mmol/m3。

2 參數(shù)敏感性分析

2.1 試驗(yàn)設(shè)置

為了構(gòu)建求解最優(yōu)參數(shù)擾動(dòng)的非線性優(yōu)化系統(tǒng), 首先進(jìn)行以下試驗(yàn)設(shè)置。

2.1.1 選取參考態(tài)

設(shè)定水柱的深度是300 m, 空間步長(zhǎng)是0.05 m, 時(shí)間步長(zhǎng)是0.1 h, 模式中所有參數(shù)值見表1。采用向前差分格式將NP模式的控制方程組[公式(10)]離散化, 并且積分20 a。

在積分時(shí)間超過(guò)8 a時(shí), 模式達(dá)到平衡態(tài)。此時(shí), 隨著表層營(yíng)養(yǎng)鹽的消耗, 在水面下一定深度處,

圖1 一維NP模式中營(yíng)養(yǎng)鹽和浮游植物的時(shí)空演變及系統(tǒng)達(dá)到平衡態(tài)后的垂直分布

來(lái)自水柱底部的營(yíng)養(yǎng)鹽與來(lái)自水面的入射光達(dá)到平衡(圖1 a)。浮游植物在此大量繁殖, NP模式呈現(xiàn)出穩(wěn)定的DCM特征(圖1 b)。圖1 c表明, 在平衡態(tài)下, 營(yíng)養(yǎng)鹽從下往上逐漸擴(kuò)散, 從DCM所在位置到水柱底部是線性增加的, 表層營(yíng)養(yǎng)鹽幾乎接近于0。同時(shí), 浮游植物生物量的最大值出現(xiàn)在水面下84 m, 浮游植物的生存區(qū)間主要集中在水面下50~ 150 m內(nèi)(圖1 d), 這與許多研究中呈現(xiàn)的DCM模擬結(jié)果是一致的[14, 23]?？梢? DCM現(xiàn)象是浮游植物對(duì)光和營(yíng)養(yǎng)鹽適應(yīng)性競(jìng)爭(zhēng)的結(jié)果, 研究DCM現(xiàn)象對(duì)理解海洋生態(tài)系統(tǒng)動(dòng)力過(guò)程有著極為重要的意義。為考察參數(shù)擾動(dòng)引起的模擬誤差大小, 在下面的敏感性分析中, 將積分終止時(shí)刻的平衡態(tài)作為擾動(dòng)的參考態(tài)。

2.1.2 擾動(dòng)約束及目標(biāo)函數(shù)

另外, 為了度量該擾動(dòng)引起的DCM變化, 將水柱上各點(diǎn)浮游植物變化量的2范數(shù)定義為目標(biāo)函數(shù)[參考公式(4)], 具體表達(dá)式如下:

其中,為水柱的深度。

2.1.3 閾值

在用OPSA方法進(jìn)行參數(shù)敏感性分析時(shí), 求解每個(gè)參數(shù)對(duì)應(yīng)的優(yōu)化問(wèn)題[公式(5)]后, 可對(duì)所得的最大模擬誤差進(jìn)行排序:1,ε≤2,ε≤…≤10,ε, 其中J,ε(1≤≤)為相應(yīng)參數(shù)擾動(dòng)在約束=0.1下所引起的最大目標(biāo)函數(shù)值。接下來(lái), 為了確定哪些參數(shù)較為敏感, 需要確定閾值。假設(shè)1,+2,+…J,≤,1,+2,+…J,+J,>(+1≤10), 則認(rèn)為1,、2,…J,對(duì)應(yīng)的參數(shù)擾動(dòng)引起的模擬誤差的最大值仍然較小, 相應(yīng)的參數(shù)被視為不敏感參數(shù), 而J+1,、J+2,…10,對(duì)應(yīng)的參數(shù)則為相對(duì)敏感參數(shù)。

由此可見, 閾值的選取對(duì)于用OPSA方法識(shí)別敏感參數(shù)是非常關(guān)鍵的。在實(shí)際應(yīng)用中, 往往是根據(jù)所研究的具體問(wèn)題來(lái)設(shè)置的?？紤]到不同水體中的海洋生態(tài)系統(tǒng), 浮游植物濃度有著較大差異[24]。且即使是同一水體中, 不同年份的同一季節(jié), 因氣候和水文環(huán)境差異導(dǎo)致的浮游植物生物量也有很大不同[25]。當(dāng)水體富營(yíng)養(yǎng)化時(shí), 還會(huì)出現(xiàn)浮游植物過(guò)量繁殖的水華現(xiàn)象, 這時(shí)浮游植物生物量可能達(dá)到平時(shí)的幾倍[26]。因此, 為了考察參數(shù)不確定性對(duì)DCM模擬的影響程度, 本文中假設(shè)水柱各點(diǎn)浮游植物生物量的變化最大不超過(guò)參考態(tài)的一倍, 則計(jì)算可得相應(yīng)的目標(biāo)函數(shù)值為8.441 2×108, 以此作為來(lái)判斷累計(jì)參數(shù)擾動(dòng)引起的模擬誤差是否達(dá)到閾值。

2.2 最優(yōu)化試驗(yàn)

本節(jié)中的最優(yōu)化試驗(yàn)是如下設(shè)計(jì)的: 首先分別計(jì)算對(duì)所有參數(shù)全部疊加擾動(dòng)以及在此基礎(chǔ)上逐一去掉其中某個(gè)參數(shù)擾動(dòng)所對(duì)應(yīng)的模擬誤差, 然后將二者的差作為目標(biāo)函數(shù), 對(duì)所有參數(shù)擾動(dòng)進(jìn)行優(yōu)化(優(yōu)化維數(shù)為10維), 即得到相應(yīng)的最優(yōu)擾動(dòng)及最大模擬誤差。這里, 為了逐個(gè)求解每個(gè)參數(shù)擾動(dòng)導(dǎo)致的最大模擬誤差, 共進(jìn)行了10次優(yōu)化試驗(yàn)。

其中, 尋找使目標(biāo)函數(shù)達(dá)到極大值的最優(yōu)參數(shù)擾動(dòng), 需要求解非線性最優(yōu)化問(wèn)題[公式(5)]。為此, 選擇譜投影梯度算法[27](spectral projected gradient, SPG)。該方法除需要提供約束條件和目標(biāo)函數(shù)外, 還需要目標(biāo)函數(shù)關(guān)于參數(shù)擾動(dòng)的梯度信息(根據(jù)梯度定義計(jì)算), 具體計(jì)算過(guò)程可參考文獻(xiàn)[28]。另外, 為使模式重新達(dá)到平衡態(tài), 優(yōu)化時(shí)間為3年。且考慮到NP模式在迭代500次左右已逐漸達(dá)到平衡態(tài), 因此迭代次數(shù)設(shè)置為1 000次。

為考察單參數(shù)的敏感性, 求解最優(yōu)化試驗(yàn)所得的最優(yōu)參數(shù)擾動(dòng)見表2。其中, 參數(shù)3(浮游植物死亡率)的不確定性影響最大時(shí), 各參數(shù)擾動(dòng)全部位于約束邊界。其余參數(shù)對(duì)應(yīng)的最優(yōu)擾動(dòng)中, 各參數(shù)擾動(dòng)并沒(méi)有全部達(dá)到邊界值。由此可見, 多參數(shù)同時(shí)發(fā)生擾動(dòng)引起的動(dòng)力過(guò)程變化之間確實(shí)存在非線性相互作用。在進(jìn)行參數(shù)敏感性分析時(shí), 必須考慮這種非線性相互作用帶來(lái)的影響, 才能對(duì)非線性系統(tǒng)進(jìn)行合理的量化分析。

表2 考察單參數(shù)敏感性所得的最優(yōu)參數(shù)擾動(dòng)

表3是求解上述優(yōu)化試驗(yàn)得到的目標(biāo)函數(shù)值?？梢? 參數(shù)1(硝酸鹽半飽和常數(shù))的不確定性影響最小, 參數(shù)9(背景場(chǎng)渾濁度)的不確定性影響最大。各目標(biāo)函數(shù)值的變化范圍從0.192到5.335, 最小值與最大值相差近30倍。因此, 在數(shù)值模擬中參數(shù)擾動(dòng)的不確定性影響相差很大, 對(duì)它們進(jìn)行敏感性分析, 從而確定最優(yōu)參數(shù)化方案和觀測(cè)方案是非常有必要的。接下來(lái), 將表3中的目標(biāo)函數(shù)值按從小到大的順序進(jìn)行排序(圖2 a), 依次為:1、4、10、2、3、8、5、6、7和9, 然后逐個(gè)進(jìn)行疊加并與閾值比較。

表3 考察單參數(shù)敏感性所得的目標(biāo)函數(shù)值

圖2 各參數(shù)的目標(biāo)函數(shù)值及累加的(按從小到大的順序)目標(biāo)函數(shù)值排序

注: 其中代表設(shè)定的閾值; (b)中橫軸刻度“…P”表示從前一參數(shù)對(duì)應(yīng)的目標(biāo)函數(shù)值累加到P對(duì)應(yīng)的目標(biāo)函數(shù)值

圖2 b表明, 參數(shù)1、4、10、2、3和8所對(duì)應(yīng)的目標(biāo)函數(shù)值累加后仍未達(dá)到給定的閾值, 即所允許的最大模式誤差。由公式(5)可知, 這6個(gè)參數(shù)同時(shí)發(fā)生擾動(dòng)導(dǎo)致的模式誤差不超過(guò)它們各自擾動(dòng)導(dǎo)致的模式誤差之和, 因此將這6個(gè)參數(shù)被視作相對(duì)不敏感參數(shù)。其他參數(shù)5(垂直湍流擴(kuò)散系數(shù))、6(浮游植物營(yíng)養(yǎng)鹽含量)、7(硝酸鹽再循環(huán)系數(shù))和9(背景場(chǎng)渾濁度)則為敏感參數(shù)。

事實(shí)上, 在海洋生態(tài)系統(tǒng)中, 垂直湍流擴(kuò)散(5)除了對(duì)浮游植物的接觸率、攝食率和生長(zhǎng)率有直接影響外, 還在很大程度上影響著營(yíng)養(yǎng)鹽的輸送[7], 是對(duì)生態(tài)系統(tǒng)起著關(guān)鍵作用的物理過(guò)程。浮游植物營(yíng)養(yǎng)鹽含量(6)與硝酸鹽再循環(huán)系數(shù)(7)則體現(xiàn)了系統(tǒng)中浮游植物死亡后營(yíng)養(yǎng)鹽的再次利用。而背景場(chǎng)渾濁度(9)影響的是光的入射率, 反映了光對(duì)該生態(tài)系統(tǒng)的控制作用。可見, 敏感參數(shù)既包含物理參數(shù), 又有描述生態(tài)動(dòng)力過(guò)程的參數(shù), 直接或間接影響著營(yíng)養(yǎng)鹽和光這兩大海洋生態(tài)系統(tǒng)的主要因素, 在很大程度上影響著模擬結(jié)果的準(zhǔn)確性。

2.3 觀測(cè)系統(tǒng)模擬試驗(yàn)(observing system simulation experiment, OSSE)

在上一小節(jié)中, 通過(guò)優(yōu)化得到相應(yīng)的最優(yōu)擾動(dòng)及其所導(dǎo)致的最大模擬誤差, 并由此識(shí)別出敏感參數(shù)為5、6、7和9。那么去除敏感參數(shù)誤差, DCM模擬是否一定會(huì)得到改進(jìn), 具體改進(jìn)程度如何呢？為此, 本小節(jié)將通過(guò)觀測(cè)系統(tǒng)模擬試驗(yàn)(OSSE)來(lái)回答這一問(wèn)題。

如圖3所示, OSSE試驗(yàn)包含真值試驗(yàn), 控制試驗(yàn)和同化試驗(yàn)。具體地, 真值試驗(yàn)即參數(shù)均未發(fā)生擾動(dòng)時(shí)的DCM模擬(參考態(tài)); 對(duì)真值試驗(yàn)中的10個(gè)參數(shù)疊加100組隨機(jī)誤差并計(jì)算由此導(dǎo)致的模擬偏差即為控制試驗(yàn); 最后計(jì)算去除其中部分參數(shù)誤差(敏感/不敏感)后的模擬偏差作為同化試驗(yàn)。在此基礎(chǔ)上, 用下式度量DCM模擬的改進(jìn)程度:

圖3 OSSE試驗(yàn)設(shè)計(jì)圖

其中,1為控制試驗(yàn)導(dǎo)致的模擬誤差,2為同化試驗(yàn)導(dǎo)致的模擬誤差。此處, 共進(jìn)行了兩組同化試驗(yàn), 包括去除敏感參數(shù)誤差和去除不敏感參數(shù)誤差(圖4)。這里,值為正表示去掉部分參數(shù)誤差后, DCM模擬得到了改進(jìn), 而值為負(fù), 則可能是由參數(shù)擾動(dòng)間的非線性相互作用導(dǎo)致的。若去除的部分參數(shù)誤差跟剩余的參數(shù)誤差之間是相互“抑制”的關(guān)系, 去除這部分參數(shù)誤差后, 導(dǎo)致的模擬結(jié)果就可能更差。

圖4 觀測(cè)系統(tǒng)模擬試驗(yàn)結(jié)果

Fig. 4 Results of OSSE experiments: removing the perturbations sensitive/insensitive parameters

據(jù)圖4 a, 當(dāng)去除敏感參數(shù)誤差時(shí), 100組隨機(jī)試驗(yàn)中有95組的DCM模擬都得到了改善, 且改善均值為56.83%; 當(dāng)去除不敏感參數(shù)誤差時(shí), 100組隨機(jī)試驗(yàn)中只有64組的DCM模擬得到了改善, 改善均值僅為4.51%(圖4 b)。由此可見, 識(shí)別敏感參數(shù)后, 有針對(duì)性地去除敏感參數(shù)誤差, 確實(shí)能有效提高DCM模擬水平。

僅考慮模擬改進(jìn)程度的大小, 不能保證去除敏感參數(shù)隨機(jī)誤差后DCM模擬一定會(huì)得到改善, 還需考察去除敏感參數(shù)誤差后模擬改進(jìn)效果的穩(wěn)定性。為此, 使用變異系數(shù)(coefficient of variation)進(jìn)一步考察OSSE試驗(yàn)所得兩組數(shù)據(jù)的相對(duì)離散程度, 公式如下:

其中,和分別為數(shù)據(jù)的標(biāo)準(zhǔn)差和均值。

由圖5, 去除敏感參數(shù)誤差所得DCM模擬的改進(jìn)程度, 其變異系數(shù)為9.44%; 而去除不敏感參數(shù)誤差所得DCM模擬的改進(jìn)程度, 其變異系數(shù)為14.76%。由此可見, 去除模式中敏感參數(shù)的誤差, 大多數(shù)情況下DCM模擬都得到了改進(jìn), 且波動(dòng)較小, 變異系數(shù)也在合理的范圍之內(nèi)。去除不敏感參數(shù)誤差, DCM模擬的改進(jìn)則很有限, 變異系數(shù)已接近15%。因此, 先識(shí)別敏感參數(shù), 再優(yōu)先改進(jìn)敏感參數(shù)誤差, 對(duì)于提高DCM模擬技巧是可靠的。

圖5 改進(jìn)敏感/不敏感參數(shù)誤差所得深層葉綠素最大值模擬的變異系數(shù)

進(jìn)一步, 將去除敏感/不敏感參數(shù)擾動(dòng)的所得的浮游植物生物量垂向分布取均值, 考察DCM位置與強(qiáng)度的變化。由圖6a可見, 與疊加所有參數(shù)擾動(dòng)相比, 去除敏感參數(shù)擾動(dòng)后DCM位置與強(qiáng)度都發(fā)生了很大改變, 且已與參考態(tài)非常接近。而去除不敏感參數(shù)擾動(dòng)后DCM并沒(méi)有較大變化(圖6b中紅線與藍(lán)線幾乎重合), 仍然遠(yuǎn)離參考態(tài)。可見, 識(shí)別具有較大不確定性的敏感參數(shù)對(duì)于改進(jìn)DCM模擬是非常重要的, 去除敏感參數(shù)誤差DCM模擬的改進(jìn)非常明顯。這在含有大量參數(shù)尚未經(jīng)由觀測(cè)明確的海洋生態(tài)數(shù)值模式中, 顯得尤為重要。

圖6 分別改善敏感參數(shù)誤差/不敏感參數(shù)誤差后, 浮游植物的垂向分布

3 結(jié)論與討論

本文基于OPSA方法探討了DCM數(shù)值模擬中的參數(shù)敏感性:

1) 求解所建立的非線性優(yōu)化系統(tǒng), 發(fā)現(xiàn)最優(yōu)擾動(dòng)中單個(gè)參數(shù)的擾動(dòng)并不總是位于邊界上, 不同參數(shù)擾動(dòng)間確實(shí)存在非線性相互作用。

2) 針對(duì)DCM模擬, 根據(jù)OPSA方法識(shí)別的不敏感參數(shù)為:1、4、10、2、3和8, 敏感參數(shù)為:5(垂直湍流擴(kuò)散系數(shù))、6(浮游植物營(yíng)養(yǎng)鹽含量)、7(硝酸鹽再循環(huán)系數(shù))和9(背景場(chǎng)渾濁度)。

3) OSSE試驗(yàn)表明, 對(duì)模式參數(shù)疊加100組隨機(jī)擾動(dòng), 去除敏感參數(shù)擾動(dòng)后有95組的DCM模擬都得到了改進(jìn), 平均改進(jìn)程度為56.83%; 而去除不敏感參數(shù)擾動(dòng), 僅有64組的DCM模擬得到了改進(jìn), 平均改進(jìn)程度僅為4.51%。因此識(shí)別敏感參數(shù)對(duì)改進(jìn)DCM模擬有著重要意義。

次表層的浮游植物貢獻(xiàn)了水體中絕大部分初級(jí)生產(chǎn)力, 因此人們對(duì)浮游植物的垂向分布非常關(guān)注。在用數(shù)值模式對(duì)DCM現(xiàn)象進(jìn)行數(shù)值模擬時(shí), 參數(shù)擾動(dòng)對(duì)模擬結(jié)果具有很大影響?；贠PSA方法識(shí)別的敏感參數(shù)為5、6、7和9, 這與使用條件非線性最優(yōu)參數(shù)擾動(dòng)方法(CNOP-P)識(shí)別的敏感參數(shù)組合是一致的[14]。三個(gè)直接影響營(yíng)養(yǎng)鹽分布的參數(shù)(垂向湍流擴(kuò)散、浮游植物營(yíng)養(yǎng)鹽含量、硝酸鹽再循環(huán)系數(shù))與光的限制參數(shù)(背景場(chǎng)渾濁度)是影響浮游植物垂向分布的最主要因素。本研究使用OPSA方法共求解了10次非線性優(yōu)化系統(tǒng)(10維), 而使用條件非線性最優(yōu)擾動(dòng)方法則進(jìn)行了210次最優(yōu)化試驗(yàn)(4維)才識(shí)別出四個(gè)敏感參數(shù), 前者大大節(jié)約了計(jì)算量。事實(shí)上, 隨著數(shù)值模式復(fù)雜程度越來(lái)越高, 模式中參數(shù)個(gè)數(shù)也越來(lái)越多, 使用OPSA方法能以較小的計(jì)算量識(shí)別敏感參數(shù)的優(yōu)勢(shì)將更加凸顯。

此外, OSSE試驗(yàn)證實(shí)去除敏感參數(shù)擾動(dòng)能更好地改進(jìn)數(shù)值模擬技巧。因此, 一方面可借助參數(shù)敏感性分析優(yōu)先發(fā)展與敏感參數(shù)直接相關(guān)的動(dòng)力過(guò)程參數(shù)化方案, 有效節(jié)省計(jì)算時(shí)間與機(jī)時(shí)消耗。另一方面, 考慮到開展大范圍的海洋觀測(cè)往往消耗巨大, 也可基于參數(shù)敏感性分析識(shí)別的敏感參數(shù), 利用有限的觀測(cè)資源優(yōu)先對(duì)敏感參數(shù)展開觀測(cè)。這正是目標(biāo)觀測(cè)的思想[29]。進(jìn)而據(jù)此對(duì)模式進(jìn)行校正, 可使其提供更好的模擬與預(yù)報(bào)。

作為中國(guó)三大河口三角洲之一的黃河三角洲, 在環(huán)渤海地區(qū)發(fā)展中具有重要的戰(zhàn)略地位。該區(qū)域生態(tài)系統(tǒng)獨(dú)具特色, 處于大氣、河流、海洋與陸地的交接帶, 是世界上典型的河口濕地生態(tài)系統(tǒng)?；贠PSA方法, 識(shí)別敏感參數(shù), 優(yōu)先發(fā)展與敏感參數(shù)直接相關(guān)的動(dòng)力過(guò)程參數(shù)化方案, 可有效提高黃河三角洲生態(tài)系統(tǒng)的模擬與預(yù)報(bào)技巧。在此基礎(chǔ)上開展目標(biāo)觀測(cè)研究, 科學(xué)指導(dǎo)黃河三角洲區(qū)域生態(tài)環(huán)境監(jiān)測(cè)網(wǎng)絡(luò)的建設(shè), 對(duì)于推動(dòng)區(qū)域海洋經(jīng)濟(jì)高質(zhì)量發(fā)展具有重要意義。

[1] 倪曉波, 黃大吉. 海洋次表層葉綠素最大值的分布和形成機(jī)制研究[J]. 海洋科學(xué), 2006, 30(5): 58-64.

NI Xiaobo, HUANG Daji. Subsurface chlorophyll maximum: its temporal-spatial distribution and formation mechanism in the ocean[J]. Marine Sciences, 2006, 30(5): 58-64.

[2] YENTSCH C S. Distribution of chlorophyll and phaeo-phytin in the open ocean[J]. Deep-Sea Research and Oceanographic Abstracts, 1965, 12(5): 653-666.

[3] ABBOTT M K. Mixing and the dynamics of the deep chlorophyll maximum in the Lake Tahoe[J]. Limnology and Oceanography, 1984, 29(4): 862-878.

[4] EVANS G T, PARSLOW J S. A model of annual plankton cycles[J]. Biological Oceanographic, 1985, 3(3): 327-347.

[5] LONGHURST A R. Toward an ecological geography of the sea[M]. Amsterdam, Netherlands: Elsevier, 1998.

[6] OMLIN M, BRUN R, REICHERT P, et al. Biogeochemical model of Lake Zürich: sensitivity, identifiability and uncertainty analysis[J]. Ecological Modelling, 2001, 141(1/3): 105-123.

[7] JPLLIFF J K, KINDLE J C, SHULMAN I G, et al. Summary diagrams for coupled Hydrodynamic-Eco-system model skill assessment[J]. Journal of Marine Systems, 2009, 76(1/2): 64-82.

[8] KEVIN C R, LUKE A W, JORDAN S R, et al. Improving the precision of lake ecosystem metabolism estimates by identifying predictors of model uncertainty[J]. Limnology and Oceanography Methods, 2014, 12(5): 303-312.

[9] FASHAM M J R, EVANS G T. The use of optimization techniques to model marine ecosystem dynamics at the JGOFS station at 47°N 20°W[J]. Philosophical Transactions of the Royal Society of London B Biological Sciences, 1995, 348(1324): 203-209.

[10] EVENSEN G, DEE D P, SCHR?TER J, et al. Parameter estimation in dynamical models[J]. CHASSIGNET E P, VERRON J. Ocean modeling and parameterizations, NATO Science Series, 1998, 516: 373-398.

[11] LOSA S N, KIVMAN G A, RYABCHENKO V A, et al. Weak constraint parameter estimation for a simple ocean ecosystem model: what can we learn about the model and data?[J]. Journal of Marine Systems, 2004, 45(1/2): 1-20.

[12] CHU P C, IVANOV L M, MARGOLINA T M, et al. On non-linear sensitivity of marine biological models to parameter variations[J]. Ecological Modelling, 2007, 206(3/4): 369-382.

[13] MU M, DUAN W S, WANG Q, et al. An extension of conditional nonlinear optimal perturbation approach and its applications[J]. Nonlinear Processes Geophysics: 2010, 17(2): 211-220.

[14] GAO Y L, MU M, ZHANG K, et al. Impacts of parameter uncertainties on deep chlorophyll maximum simulation revealed by the CNOP-P approach[J]. Journal of Oceanology and Limnology, 2020, 38(5): 1382- 1393.

[15] SUN G D, MU M. The analyses of the net primary production due to regional and seasonal temperature differences in eastern China using the LPJ model[J]. Ecological Modelling: 2014, 289: 66-76.

[16] WANG Q, TANG Y. An optimization strategy for identifying parameter sensitivity in atmospheric and oceanic models[J]. Monthly Weather Review, 2017, 145(8): 3293-3305.

[17] WANG Q, PIERINI S, TANG Y, et al. Parameter sensitivity analysis of the short-range prediction of Kuroshio extension transition processes using an optimization approach[J]. Theoretical and Applied Climatology, 2019, 138(3/4): 1481-1492.

[18] 高永麗. 垂向湍流擴(kuò)散系數(shù)的不確定性對(duì)深層葉綠素最大值現(xiàn)象模擬的影響[J]. 海洋科學(xué), 2019, 43(2): 34-40.

GAO Yongli. A study of uncertainty related to the coefficient of vertical turbulence diffusion in ocean ecosystem model[J]. Marine Sciences, 2019, 43(2): 34-40.

[19] HODGES B A, RUDNICK D L. Simple models of steady deep maxima in chlorophyll and biomass[J]. Deep Sea Research Part I: Oceanographic Research Papers, 2004, 51(8): 999-1015.

[20] HUISMAN J, PHAM T N, KARL D M, et al. Reduced mixing generates oscillations and chaos in the oceanic deep chlorophyll maximum[J]. Nature, 2006, 439(7074): 322-325.

[21] KLAUSMEIER C A, LITCHMAN E. Algal games: The vertical distribution of phytoplankton in poorly mixed water columns[J]. Journal of Oceanology and Limnology, 2001, 46(8): 1998-2007.

[22] LICCARDO A, FIERRO A, IUDICONE D, et al. Response of the deep chlorophyll maximum to fluctuations in vertical mixing intensity[J]. Progress in Oceanography, 2013, 109: 33-46.

[23] WANG Q, MU M. A new application of conditional nonlinear optimal perturbation approach to boundary condition uncertainty[J]. Journal of Geophysical Research: Oceans, 2015, 120(12): 7979-7996.

[24] BOYCE D G, MARLON R L, BORIS W, et al. Global phytoplankton decline over the past century[J]. Nature, 2010, 466(29): 591-466.

[25] MCQUATTERS-GOLLOP A, REID P C, EDWARDS M, et al. Is there a decline in marine phytoplankton?[J]. Nature, 2011, 472(7342): E6-E7.

[26] BOYD P W, WATSON A J, LAW C S, et al. A mesoscale phytoplankton bloom in the polar Southern Ocean stimulated by iron fertilization[J]. Nature, 2000, 407(6805): 695-702.

[27] BIRGIN E G, MARTíNEZ J M, RAYDAN M, et al. Nonmonotone spectral projected gradient methods on convex sets[J]. SIAM Journal of Optimization, 2000, 10(4): 1196-1211.

[28] LIU X, MU M, WANG Q, et al. The nonlinear optimal triggering perturbation of the Kuroshio large meander and its evolution in a regional ocean model[J]. Journal of Physical Oceanography, 2018, 48(8): 1771-1786.

[29] 張坤, 穆穆, 王強(qiáng). 數(shù)值模式在海洋觀測(cè)設(shè)計(jì)中的重要作用: 回顧與展望[J]. 中國(guó)科學(xué): 地球科學(xué), 2021, 51(5): 653-665.

ZHANG Kun, MU Mu, WANG Qiang. Increasingly important role of numerical modeling in oceanic observation design strategy: A review[J]. Science China Earth Sciences, 2021, 51(5): 653-665.

Parameter sensitivity analysis of a biological ocean model based on OPSA method

GAO Yong-Li1, WANG Ji-chao1, SUN Guo-dong2, ZHANG Kun3, JIANG Xiang-yang4, WANG Ning4

(1. China University of Petroleum (East China), Qingdao 266580, China; 2. Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100029, China; 3. Institute of Oceanology, Chinese Academy of Sciences, Qingdao 266071, China; 4. Shandong Marine Resources and Environment Research Institute, Yantai 264006, China)

The deep chlorophyll maximum (DCM) is a common ecological phenomenon in oceans and lakes. Numerical simulation has emerged as an important tool for studying this phenomenon, and parameter uncertainty is a primary source of uncertainty in the simulations. By optimization parameter sensitivity analysis, the sensitivities of 10 parameters in the nutrient–phytoplankton model related to DCM simulation were investigated. The results revealed the sensitive parameters in the DCM simulation to be background turbidity, vertical turbulent diffusion, nutrient content of phytoplankton, and recycling coefficient of nitrate; perturbations in these parameters lead to considerable changes in the DCM. In addition, the observing system simulation experiment was designed to evaluate the improvement in DCM simulation while eliminating sensitive parameter errors. The results revealed that the average improvement in DCM simulation resulting from the removal of the sensitive parameter errors is 56.83%, which is approximately 13 times that obtained from the removal of the insensitive parameter errors (4.51%). Moreover, the coefficient of variation was calculated to examine the stability of simulation improvement. The values obtained were 9.44% for the removal of sensitive parameter perturbations and 14.76% for the removal of insensitive parameter perturbations, indicating decreased stability. This study suggests that prioritizing the parameterization scheme and target observation related to sensitive parameters may provide valuable insights for the advancement of DCM simulations and predictions.

DCM; parameter sensitivity; optimal perturbation; OPSA; target observation

May 17, 2022

731.26

A

1000-3096(2023)5-0139-10

10.11759/hykx20220517003

2022-05-17;

2022-07-18

國(guó)家自然科學(xué)基金項(xiàng)目(92158202, 41576015)

[The National Natural Science Foundation of China, Nos. 92158202, 41576015]

高永麗(1981—), 女, 山東膠州人, 漢族, 講師, 博士, 主要從事非線性最優(yōu)化、參數(shù)敏感性分析方面的教學(xué)與研究, E-mail: gaoyongli@upc.edu.cn; 張坤(1988—),通信作者, 男, 山東濟(jì)寧人, 副研究員, 主要從事海洋環(huán)流及其可預(yù)報(bào)性的研究, E-mail: kzhang@qdio.ac.con

(本文編輯: 楊 悅)