• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

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

    2023-07-29 07:23:02高永麗王際朝孫國棟姜向陽
    海洋科學(xué) 2023年5期
    關(guān)鍵詞:營養(yǎng)鹽擾動敏感性

    高永麗, 王際朝, 孫國棟, 張 坤, 姜向陽, 王 寧

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

    高永麗1, 王際朝1, 孫國棟2, 張 坤3, 姜向陽4, 王 寧4

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

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

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

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

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

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

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

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

    1 方法與模式介紹

    1.1 OPSA方法

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

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

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

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

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

    接下來, 考察多個參數(shù)同時發(fā)生擾動引起的模式誤差。令:

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

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

    , (8)

    故:

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

    1.2 NP模式

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

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

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

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

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

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

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

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

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

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

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

    2.1.1 選取參考態(tài)

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

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

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

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

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

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

    其中,為水柱的深度。

    2.1.3 閾值

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    其中,和分別為數(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), 且波動較小, 變異系數(shù)也在合理的范圍之內(nèi)。去除不敏感參數(shù)誤差, DCM模擬的改進(jìn)則很有限, 變異系數(shù)已接近15%。因此, 先識別敏感參數(shù), 再優(yōu)先改進(jìn)敏感參數(shù)誤差, 對于提高DCM模擬技巧是可靠的。

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

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

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

    3 結(jié)論與討論

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

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

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

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

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

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

    作為中國三大河口三角洲之一的黃河三角洲, 在環(huán)渤海地區(qū)發(fā)展中具有重要的戰(zhàn)略地位。該區(qū)域生態(tài)系統(tǒng)獨(dú)具特色, 處于大氣、河流、海洋與陸地的交接帶, 是世界上典型的河口濕地生態(tài)系統(tǒng)。基于OPSA方法, 識別敏感參數(shù), 優(yōu)先發(fā)展與敏感參數(shù)直接相關(guān)的動力過程參數(shù)化方案, 可有效提高黃河三角洲生態(tài)系統(tǒng)的模擬與預(yù)報技巧。在此基礎(chǔ)上開展目標(biāo)觀測研究, 科學(xué)指導(dǎo)黃河三角洲區(qū)域生態(tài)環(huán)境監(jiān)測網(wǎng)絡(luò)的建設(shè), 對于推動區(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ù)的不確定性對深層葉綠素最大值現(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ù)值模式在海洋觀測設(shè)計中的重要作用: 回顧與展望[J]. 中國科學(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

    國家自然科學(xué)基金項目(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ù)報性的研究, E-mail: kzhang@qdio.ac.con

    (本文編輯: 楊 悅)

    猜你喜歡
    營養(yǎng)鹽擾動敏感性
    Bernoulli泛函上典則酉對合的擾動
    (h)性質(zhì)及其擾動
    釔對Mg-Zn-Y-Zr合金熱裂敏感性影響
    涼水河子河營養(yǎng)鹽滯留能力評估
    小噪聲擾動的二維擴(kuò)散的極大似然估計
    AH70DB鋼焊接熱影響區(qū)組織及其冷裂敏感性
    焊接(2016年1期)2016-02-27 12:55:37
    瓊東海域冬季、夏季營養(yǎng)鹽結(jié)構(gòu)特征及其對浮游植物生長的影響
    2012年冬季南海西北部營養(yǎng)鹽分布及結(jié)構(gòu)特征
    用于光伏MPPT中的模糊控制占空比擾動法
    如何培養(yǎng)和提高新聞敏感性
    新聞傳播(2015年8期)2015-07-18 11:08:24
    久久久久久久国产电影| 午夜福利视频精品| www.熟女人妻精品国产 | videos熟女内射| 91精品国产国语对白视频| av又黄又爽大尺度在线免费看| av福利片在线| 中文字幕人妻熟女乱码| 日韩精品免费视频一区二区三区 | 亚洲精品456在线播放app| 黄色毛片三级朝国网站| 国产精品女同一区二区软件| av一本久久久久| 欧美精品国产亚洲| 午夜av观看不卡| 亚洲经典国产精华液单| 国产精品国产av在线观看| 亚洲少妇的诱惑av| 丝瓜视频免费看黄片| 大片免费播放器 马上看| 精品午夜福利在线看| 久久99蜜桃精品久久| 中文字幕人妻丝袜制服| kizo精华| av在线老鸭窝| 免费看av在线观看网站| a级毛片黄视频| 91精品国产国语对白视频| 亚洲情色 制服丝袜| av片东京热男人的天堂| 国产永久视频网站| 欧美激情极品国产一区二区三区 | 蜜臀久久99精品久久宅男| 免费大片黄手机在线观看| 一级爰片在线观看| 啦啦啦啦在线视频资源| videossex国产| 欧美另类一区| 少妇被粗大猛烈的视频| 男女边摸边吃奶| 亚洲一码二码三码区别大吗| 大陆偷拍与自拍| 女人久久www免费人成看片| 欧美激情极品国产一区二区三区 | 亚洲精品国产色婷婷电影| 免费久久久久久久精品成人欧美视频 | 赤兔流量卡办理| 久久久久久久亚洲中文字幕| 97在线人人人人妻| 少妇精品久久久久久久| 亚洲av在线观看美女高潮| 国产精品国产三级国产专区5o| 熟女人妻精品中文字幕| 亚洲欧美一区二区三区国产| 欧美 亚洲 国产 日韩一| 久久人人爽av亚洲精品天堂| 日韩制服骚丝袜av| 中文字幕人妻熟女乱码| 少妇人妻 视频| 日本色播在线视频| 久久99一区二区三区| 精品视频人人做人人爽| 日韩制服丝袜自拍偷拍| 高清毛片免费看| 久久99热6这里只有精品| 久久精品久久久久久久性| 天堂俺去俺来也www色官网| 水蜜桃什么品种好| 秋霞在线观看毛片| 九色亚洲精品在线播放| 大香蕉97超碰在线| 国产片特级美女逼逼视频| 熟女人妻精品中文字幕| 中文乱码字字幕精品一区二区三区| 水蜜桃什么品种好| 成年人午夜在线观看视频| 91精品国产国语对白视频| 在线观看美女被高潮喷水网站| 欧美日韩视频高清一区二区三区二| 观看美女的网站| 2018国产大陆天天弄谢| 日韩中文字幕视频在线看片| 亚洲美女视频黄频| 26uuu在线亚洲综合色| 大香蕉久久成人网| 男女啪啪激烈高潮av片| 一级毛片黄色毛片免费观看视频| 看十八女毛片水多多多| 捣出白浆h1v1| 老司机影院毛片| 亚洲人成77777在线视频| 欧美激情 高清一区二区三区| 日本欧美视频一区| 国产无遮挡羞羞视频在线观看| 好男人视频免费观看在线| 天美传媒精品一区二区| 一级片'在线观看视频| 精品国产一区二区久久| 大香蕉久久成人网| 欧美日韩综合久久久久久| 午夜福利,免费看| 人妻 亚洲 视频| 国产又爽黄色视频| 老熟女久久久| 婷婷色av中文字幕| 日韩一区二区三区影片| 色婷婷久久久亚洲欧美| videossex国产| 中文字幕av电影在线播放| 在线观看国产h片| 啦啦啦在线观看免费高清www| 日韩av免费高清视频| 一级毛片 在线播放| 一级a做视频免费观看| 欧美日韩视频高清一区二区三区二| 久久青草综合色| 国精品久久久久久国模美| 久久久久人妻精品一区果冻| 2018国产大陆天天弄谢| 日本猛色少妇xxxxx猛交久久| 两个人看的免费小视频| 欧美日韩国产mv在线观看视频| av国产久精品久网站免费入址| 欧美日韩国产mv在线观看视频| 2021少妇久久久久久久久久久| 秋霞在线观看毛片| 老熟女久久久| 亚洲欧美色中文字幕在线| 国产在线一区二区三区精| 建设人人有责人人尽责人人享有的| 三级国产精品片| 午夜精品国产一区二区电影| 午夜福利视频在线观看免费| 午夜免费观看性视频| 中国美白少妇内射xxxbb| 午夜免费男女啪啪视频观看| av视频免费观看在线观看| 赤兔流量卡办理| 18+在线观看网站| 国产一区二区三区综合在线观看 | 国产一区二区在线观看av| 日韩制服骚丝袜av| 2018国产大陆天天弄谢| 久久ye,这里只有精品| 久久综合国产亚洲精品| videossex国产| 免费日韩欧美在线观看| 国产 精品1| 欧美丝袜亚洲另类| 极品少妇高潮喷水抽搐| 亚洲国产看品久久| 成人午夜精彩视频在线观看| 人妻 亚洲 视频| 日韩一区二区视频免费看| 国产又色又爽无遮挡免| 女性生殖器流出的白浆| 亚洲精品aⅴ在线观看| 国产精品.久久久| 精品人妻熟女毛片av久久网站| 捣出白浆h1v1| 国产免费现黄频在线看| 中文字幕人妻丝袜制服| 亚洲精品aⅴ在线观看| 蜜臀久久99精品久久宅男| 在线观看免费日韩欧美大片| 只有这里有精品99| 免费播放大片免费观看视频在线观看| 五月伊人婷婷丁香| 久久久亚洲精品成人影院| av网站免费在线观看视频| 在线天堂最新版资源| 大话2 男鬼变身卡| 母亲3免费完整高清在线观看 | 欧美3d第一页| 丝袜在线中文字幕| 女的被弄到高潮叫床怎么办| 久久久久久久国产电影| 精品人妻一区二区三区麻豆| 免费看不卡的av| av网站免费在线观看视频| 日日撸夜夜添| 九色成人免费人妻av| 在线 av 中文字幕| 中文字幕亚洲精品专区| 国产av一区二区精品久久| 大陆偷拍与自拍| 久久精品国产自在天天线| 看非洲黑人一级黄片| 亚洲伊人久久精品综合| 亚洲精品乱久久久久久| 中文精品一卡2卡3卡4更新| √禁漫天堂资源中文www| 97在线视频观看| 青春草亚洲视频在线观看| 日本av手机在线免费观看| 国产精品久久久久久av不卡| 一区二区三区乱码不卡18| 婷婷色综合www| 水蜜桃什么品种好| 免费黄网站久久成人精品| 亚洲国产日韩一区二区| 亚洲,欧美精品.| 亚洲成色77777| 日本-黄色视频高清免费观看| 国产欧美日韩综合在线一区二区| 日本午夜av视频| 亚洲综合色惰| 人人妻人人添人人爽欧美一区卜| 2021少妇久久久久久久久久久| 女的被弄到高潮叫床怎么办| 9热在线视频观看99| 亚洲五月色婷婷综合| 中文字幕精品免费在线观看视频 | 不卡视频在线观看欧美| 午夜福利网站1000一区二区三区| 99久久中文字幕三级久久日本| 91久久精品国产一区二区三区| 久久国产亚洲av麻豆专区| 春色校园在线视频观看| 国产精品熟女久久久久浪| 女的被弄到高潮叫床怎么办| 国产不卡av网站在线观看| 久久久久久久久久人人人人人人| 亚洲一区二区三区欧美精品| 亚洲av成人精品一二三区| 777米奇影视久久| 亚洲精品久久午夜乱码| a 毛片基地| 久久精品国产a三级三级三级| 赤兔流量卡办理| 精品福利永久在线观看| 黑人高潮一二区| 国产免费一区二区三区四区乱码| 麻豆乱淫一区二区| 国产一区有黄有色的免费视频| 日韩,欧美,国产一区二区三区| 母亲3免费完整高清在线观看 | 欧美激情国产日韩精品一区| 国产永久视频网站| 777米奇影视久久| 99久久综合免费| 成人二区视频| 日韩电影二区| 国产精品.久久久| 最近最新中文字幕大全免费视频 | 大香蕉久久成人网| 亚洲国产av影院在线观看| 国产男人的电影天堂91| 亚洲精品av麻豆狂野| 中文乱码字字幕精品一区二区三区| 晚上一个人看的免费电影| 久久狼人影院| 久久精品人人爽人人爽视色| 亚洲精品乱码久久久久久按摩| 少妇人妻久久综合中文| 亚洲精品av麻豆狂野| 亚洲人与动物交配视频| 乱人伦中国视频| 国产不卡av网站在线观看| 夫妻性生交免费视频一级片| 1024视频免费在线观看| 菩萨蛮人人尽说江南好唐韦庄| 国产激情久久老熟女| 国产成人a∨麻豆精品| 人妻少妇偷人精品九色| 免费黄网站久久成人精品| 一区二区日韩欧美中文字幕 | 蜜臀久久99精品久久宅男| 国产成人91sexporn| 成人漫画全彩无遮挡| 最后的刺客免费高清国语| 亚洲精品第二区| 乱码一卡2卡4卡精品| 午夜福利,免费看| 精品熟女少妇av免费看| 亚洲天堂av无毛| 男女高潮啪啪啪动态图| 国产精品久久久久成人av| 99久久人妻综合| 久久久久精品人妻al黑| 男女边摸边吃奶| 免费av中文字幕在线| 丝袜美足系列| 国产精品国产三级国产专区5o| 亚洲三级黄色毛片| 中文字幕人妻丝袜制服| 亚洲av.av天堂| 又黄又粗又硬又大视频| 国产一区有黄有色的免费视频| 久久精品人人爽人人爽视色| 亚洲天堂av无毛| 色视频在线一区二区三区| 在线观看www视频免费| 日韩在线高清观看一区二区三区| 亚洲一区二区三区欧美精品| 一级片'在线观看视频| 国产精品久久久久成人av| 欧美亚洲日本最大视频资源| 久久这里有精品视频免费| 日本黄大片高清| 人人妻人人澡人人看| 美女xxoo啪啪120秒动态图| 水蜜桃什么品种好| 天天影视国产精品| 久久精品人人爽人人爽视色| 国产一区二区三区av在线| 伦理电影免费视频| 亚洲五月色婷婷综合| 两个人免费观看高清视频| 久久久亚洲精品成人影院| 丰满乱子伦码专区| 国产成人aa在线观看| 午夜福利乱码中文字幕| 大话2 男鬼变身卡| 天堂中文最新版在线下载| 亚洲国产av新网站| 美女内射精品一级片tv| 熟妇人妻不卡中文字幕| 九草在线视频观看| 亚洲美女视频黄频| 午夜老司机福利剧场| tube8黄色片| 成人毛片60女人毛片免费| 日韩大片免费观看网站| 在线观看免费高清a一片| 日韩一区二区三区影片| 国产成人午夜福利电影在线观看| freevideosex欧美| 国产亚洲最大av| 亚洲精品中文字幕在线视频| 久久亚洲国产成人精品v| 日韩av在线免费看完整版不卡| 国产亚洲欧美精品永久| 亚洲av免费高清在线观看| 国产一区有黄有色的免费视频| 亚洲成人av在线免费| 亚洲伊人久久精品综合| 嫩草影院入口| 99久久人妻综合| 成人亚洲精品一区在线观看| 久久精品久久久久久噜噜老黄| 久久午夜福利片| 97在线人人人人妻| 永久免费av网站大全| av又黄又爽大尺度在线免费看| 美国免费a级毛片| 亚洲伊人色综图| 色94色欧美一区二区| 日本色播在线视频| 黄片播放在线免费| 久久久精品区二区三区| 免费看不卡的av| 亚洲av中文av极速乱| 色哟哟·www| 久久久久精品性色| 巨乳人妻的诱惑在线观看| 精品熟女少妇av免费看| av免费在线看不卡| 2021少妇久久久久久久久久久| 老熟女久久久| 免费女性裸体啪啪无遮挡网站| 少妇的逼好多水| 97精品久久久久久久久久精品| 亚洲欧美成人精品一区二区| 国产精品国产三级专区第一集| 欧美少妇被猛烈插入视频| 97精品久久久久久久久久精品| 日韩成人伦理影院| 最黄视频免费看| 国产在视频线精品| 日韩成人伦理影院| 国产一区二区三区av在线| 亚洲高清免费不卡视频| 成人综合一区亚洲| 这个男人来自地球电影免费观看 | 新久久久久国产一级毛片| 高清黄色对白视频在线免费看| 国国产精品蜜臀av免费| 免费看av在线观看网站| 熟女电影av网| 成人影院久久| 十八禁高潮呻吟视频| 欧美xxⅹ黑人| 久久国产亚洲av麻豆专区| 亚洲第一av免费看| 69精品国产乱码久久久| 秋霞在线观看毛片| 欧美最新免费一区二区三区| 国产一区亚洲一区在线观看| 国产欧美日韩一区二区三区在线| 午夜福利网站1000一区二区三区| 亚洲精品久久午夜乱码| h视频一区二区三区| 亚洲天堂av无毛| av不卡在线播放| 18+在线观看网站| 99热国产这里只有精品6| 精品少妇内射三级| 91成人精品电影| 男女边摸边吃奶| 精品卡一卡二卡四卡免费| 久久久精品免费免费高清| 一区在线观看完整版| 日韩欧美精品免费久久| 一区二区日韩欧美中文字幕 | 在线观看三级黄色| 侵犯人妻中文字幕一二三四区| 亚洲成av片中文字幕在线观看 | av有码第一页| 十分钟在线观看高清视频www| 日韩 亚洲 欧美在线| 99re6热这里在线精品视频| 菩萨蛮人人尽说江南好唐韦庄| 亚洲四区av| 亚洲国产精品999| 成人国产麻豆网| 9热在线视频观看99| a级片在线免费高清观看视频| 香蕉国产在线看| 午夜91福利影院| 欧美变态另类bdsm刘玥| 亚洲av成人精品一二三区| 日韩 亚洲 欧美在线| videosex国产| 亚洲精品色激情综合| 天堂俺去俺来也www色官网| www.av在线官网国产| 亚洲综合色惰| 最近的中文字幕免费完整| 免费av不卡在线播放| 多毛熟女@视频| 黑人巨大精品欧美一区二区蜜桃 | 日日爽夜夜爽网站| 精品人妻一区二区三区麻豆| 中文天堂在线官网| 久久久久久久精品精品| 国产高清三级在线| 一区二区三区精品91| 久久精品国产综合久久久 | 18禁裸乳无遮挡动漫免费视频| 最近的中文字幕免费完整| 大香蕉久久网| 免费不卡的大黄色大毛片视频在线观看| 欧美成人午夜精品| 香蕉丝袜av| 色视频在线一区二区三区| 高清在线视频一区二区三区| 纵有疾风起免费观看全集完整版| 国产日韩一区二区三区精品不卡| 在线看a的网站| 中国三级夫妇交换| 激情视频va一区二区三区| 女性生殖器流出的白浆| 97在线视频观看| 男女午夜视频在线观看 | 女性被躁到高潮视频| 大香蕉久久网| 久久ye,这里只有精品| 国产深夜福利视频在线观看| 高清黄色对白视频在线免费看| 欧美日韩视频高清一区二区三区二| 欧美精品高潮呻吟av久久| av播播在线观看一区| 婷婷色综合www| 日本黄色日本黄色录像| 午夜久久久在线观看| 亚洲av在线观看美女高潮| freevideosex欧美| av国产久精品久网站免费入址| 亚洲欧美成人精品一区二区| 女人精品久久久久毛片| 久久久精品区二区三区| 男女免费视频国产| 嫩草影院入口| 免费看不卡的av| 18禁裸乳无遮挡动漫免费视频| 免费黄色在线免费观看| 日韩熟女老妇一区二区性免费视频| 大片电影免费在线观看免费| 久久久久久人妻| 一区二区三区精品91| 久久国内精品自在自线图片| 欧美人与善性xxx| 一边亲一边摸免费视频| 晚上一个人看的免费电影| 22中文网久久字幕| 9191精品国产免费久久| 人妻一区二区av| 精品午夜福利在线看| 午夜精品国产一区二区电影| 亚洲av中文av极速乱| 久久国内精品自在自线图片| 国产成人精品福利久久| 欧美精品av麻豆av| 国产色婷婷99| 久久久久久伊人网av| 中文字幕制服av| 一级,二级,三级黄色视频| 亚洲精品国产av成人精品| 精品少妇黑人巨大在线播放| 国产精品三级大全| 在线观看一区二区三区激情| 国产一区二区三区av在线| 亚洲欧美日韩另类电影网站| 大话2 男鬼变身卡| av不卡在线播放| 亚洲人成77777在线视频| 久久韩国三级中文字幕| 18+在线观看网站| 黑人高潮一二区| 精品99又大又爽又粗少妇毛片| 91国产中文字幕| 色婷婷久久久亚洲欧美| 成年人免费黄色播放视频| 青青草视频在线视频观看| 国产又色又爽无遮挡免| 大码成人一级视频| 十八禁高潮呻吟视频| 国产精品一国产av| 亚洲人与动物交配视频| 熟女电影av网| 91精品伊人久久大香线蕉| 免费看光身美女| 亚洲欧洲日产国产| 国产高清三级在线| 少妇人妻久久综合中文| 免费黄网站久久成人精品| 大香蕉久久成人网| av女优亚洲男人天堂| 精品国产露脸久久av麻豆| 国产精品久久久久久av不卡| 99国产精品免费福利视频| 成人影院久久| 人妻 亚洲 视频| 久久久国产一区二区| 又黄又爽又刺激的免费视频.| 亚洲国产av影院在线观看| 成人二区视频| 超色免费av| 亚洲五月色婷婷综合| 亚洲av日韩在线播放| 欧美日韩国产mv在线观看视频| 久久久久久久久久久免费av| 少妇猛男粗大的猛烈进出视频| 亚洲成av片中文字幕在线观看 | 精品少妇黑人巨大在线播放| 伦理电影免费视频| 男女无遮挡免费网站观看| 最近最新中文字幕免费大全7| 香蕉精品网在线| 国产黄频视频在线观看| 最近的中文字幕免费完整| 亚洲国产精品一区二区三区在线| 国产1区2区3区精品| 尾随美女入室| 国产片内射在线| 国产免费又黄又爽又色| 色吧在线观看| 精品亚洲乱码少妇综合久久| 两性夫妻黄色片 | 久久精品熟女亚洲av麻豆精品| 菩萨蛮人人尽说江南好唐韦庄| 激情视频va一区二区三区| 视频在线观看一区二区三区| 色婷婷久久久亚洲欧美| 精品国产露脸久久av麻豆| 亚洲高清免费不卡视频| 亚洲中文av在线| 视频中文字幕在线观看| 成人午夜精彩视频在线观看| 黑人高潮一二区| av在线app专区| 亚洲国产精品专区欧美| 日韩制服丝袜自拍偷拍| 国产熟女午夜一区二区三区| 激情视频va一区二区三区| 日韩av在线免费看完整版不卡| 国产精品久久久久久精品电影小说| 一级毛片黄色毛片免费观看视频| 91精品伊人久久大香线蕉| av有码第一页| 国产av国产精品国产| 精品少妇黑人巨大在线播放| 色视频在线一区二区三区| 高清视频免费观看一区二区| 成人18禁高潮啪啪吃奶动态图| 午夜免费男女啪啪视频观看| 9色porny在线观看| 国产探花极品一区二区| 亚洲精品av麻豆狂野| 女性生殖器流出的白浆| 亚洲高清免费不卡视频| 日本免费在线观看一区| 久久午夜福利片| 人人妻人人澡人人爽人人夜夜| 久久韩国三级中文字幕| 看十八女毛片水多多多| 精品卡一卡二卡四卡免费| 国产精品一区二区在线观看99| 精品卡一卡二卡四卡免费| 香蕉丝袜av| 五月天丁香电影| 久久久久久久久久久久大奶| 亚洲人与动物交配视频| 超色免费av| 欧美亚洲日本最大视频资源| 欧美日韩精品成人综合77777| 久久久久久久国产电影| www.色视频.com| 午夜影院在线不卡| 多毛熟女@视频| 尾随美女入室| av在线播放精品| 国产成人精品久久久久久| 欧美国产精品va在线观看不卡|