亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        一種新的模式傾向誤差估計(jì)算法及其在ENSO模擬中的應(yīng)用*

        2022-09-21 02:37:10高艷秋唐佑民張繼才
        海洋與湖沼 2022年5期
        關(guān)鍵詞:模態(tài)特征

        何 群 高艷秋 唐佑民, 4 張繼才

        一種新的模式傾向誤差估計(jì)算法及其在ENSO模擬中的應(yīng)用*

        何 群1, 2高艷秋2, 3①唐佑民2, 3, 4張繼才1

        (1. 浙江大學(xué)海洋學(xué)院 浙江舟山 316021; 2. 自然資源部第二海洋研究所衛(wèi)星海洋環(huán)境動(dòng)力學(xué)國家重點(diǎn)實(shí)驗(yàn)室 浙江杭州 310012; 3. 南方海洋科學(xué)與工程廣東省實(shí)驗(yàn)室(珠海) 廣東珠海 519082; 4. 河海大學(xué)海洋學(xué)院 江蘇南京 210098)

        氣候模式是我們理解、模擬和預(yù)報(bào)氣候演變的重要工具。然而即使是目前最先進(jìn)的耦合模式, 其模擬和預(yù)報(bào)與大氣/海洋的真實(shí)狀態(tài)相比, 仍存在較大偏差, 這是由于在模式的傾向方程中不可避免地存在系統(tǒng)性的誤差(傾向誤差)。因此, 減小模式傾向誤差對改進(jìn)模式的模擬和預(yù)報(bào)效果具有重要意義。該研究首先發(fā)展了一種新的計(jì)算模式傾向誤差的估計(jì)算法——基于局地集合變換卡爾曼濾波器(local ensemble transform kalman filter, LETKF)同化技術(shù)的傾向誤差估計(jì)算法。在此基礎(chǔ)上, 將新發(fā)展的算法應(yīng)用到Zebiak-Cane (ZC)模式, 通過同化海表面溫度異常(sea surface temperature anomaly, SSTA)數(shù)據(jù), 估計(jì)隨時(shí)空變化的傾向誤差, 并使用計(jì)算得到的傾向誤差訂正模式, 進(jìn)行積分模擬。結(jié)果表明: (1) 傾向誤差和ZC模式的模擬偏差具有高度相關(guān)性; (2) 訂正后的模式改善了對厄爾尼諾-南方濤動(dòng)(El Ni?o-Southern Oscillation, ENSO)的一些重要特征的模擬。這說明新發(fā)展的模式傾向誤差估計(jì)算法十分有效且在ENSO模擬中具有較好的應(yīng)用價(jià)值, 此外, 這種新的模式傾向誤差估計(jì)算法, 計(jì)算高效簡便, 可便捷地應(yīng)用于各模式中, 利于推廣。

        模式傾向誤差; 參數(shù)估計(jì); 局地集合變換卡爾曼濾波器; Zebiak-Cane模式

        厄爾尼諾-南方濤動(dòng)(El Ni?o-Southern Oscillation, ENSO)是短期氣候年際變化的最強(qiáng)信號(hào)之一(連濤等, 2017), 主要表現(xiàn)為熱帶太平洋海表面溫度每隔2~7 a的大范圍變暖或變冷, 其中暖事件稱為厄爾尼諾(El Ni?o), 冷事件稱為拉尼娜(La Ni?a)。ENSO影響著全球的溫度和降水, 它的發(fā)生往往會(huì)在全球多地造成嚴(yán)重的自然災(zāi)害, 如El Ni?o期間澳大利亞、印度尼西亞等地的炎熱干旱, 南美太平洋沿岸國家的暴雨洪澇, 以及我國的南澇北旱。而La Ni?a與El Ni?o作用相反, 當(dāng)La Ni?a事件發(fā)生時(shí), 南美洲地區(qū)旱情嚴(yán)重, 東南亞、澳大利亞等地區(qū)會(huì)出現(xiàn)暴雨, 甚至出現(xiàn)洪澇災(zāi)害(夏詠華, 2001; 王琳等, 2019)。因此, ENSO研究一直是大氣海洋領(lǐng)域的熱點(diǎn)和焦點(diǎn)。

        近幾十年里, ENSO的觀測、模擬和預(yù)報(bào)研究取得了長足進(jìn)展, 目前國際上已有20多個(gè)模式提供6~12個(gè)月的ENSO預(yù)報(bào)(https://iri.columbia.edu/our- expertise/climate/forecasts/enso/current/), 但它們普遍存在一定的預(yù)報(bào)誤差, 特別是對歷史上幾次El Ni?o事件的錯(cuò)誤預(yù)報(bào)嚴(yán)重制約了ENSO的實(shí)際預(yù)報(bào)技巧(穆穆等, 2017)。此外, 近些年來中部型El Ni?o事件頻發(fā)(Yeh, 2009), 給模式的模擬和預(yù)報(bào)帶來了新的挑戰(zhàn)。因此, 如何減小模式誤差, 提高模式對ENSO的模擬、預(yù)報(bào)能力, 是一項(xiàng)十分重要且前沿的工作。

        影響ENSO預(yù)報(bào)技巧的因素主要分為初始誤差和模式誤差兩種。與初始誤差相關(guān)的研究十分豐富且成熟: 一方面, 可通過數(shù)據(jù)同化方法優(yōu)化模式預(yù)報(bào)的初始條件(Chen, 1997, 1998), 另一方面, 可通過加強(qiáng)對重點(diǎn)區(qū)域的觀測, 獲取更有效的觀測資料, 優(yōu)化模式初始場(Mu, 2015; Hu, 2016)。除初始誤差外, 越來越多的研究開始關(guān)注模式誤差的影響(Qi, 2017), 且有研究表明, 對初始誤差的改進(jìn)僅能在提前幾個(gè)月的預(yù)報(bào)中體現(xiàn), 而對模式誤差的改進(jìn)能提高長期預(yù)報(bào)的技巧(Zheng, 2009)。

        模式誤差主要來源于物理內(nèi)核不匹配、物理方案近似和參數(shù)誤差三個(gè)方面(Wu, 2016)。鑒于改進(jìn)前兩個(gè)方面的難度, 許多學(xué)者通過改進(jìn)參數(shù)誤差來降低模式誤差。然而, 針對模式參數(shù)的優(yōu)化通常只能優(yōu)化特定的參數(shù)。因此, Roads (1987)提出在狀態(tài)趨勢方程中添加一個(gè)常數(shù)項(xiàng), 稱為恒定傾向誤差, 用于表示多種來源模式誤差的綜合效應(yīng), 并利用觀測信息估計(jì)恒定傾向誤差。依據(jù)這一思想, Moore等(1999)在隨機(jī)動(dòng)力系統(tǒng)的框架下, 把傾向誤差視為隨機(jī)誤差, 并利用隨機(jī)最優(yōu)理論來定量計(jì)算增長最快的隨機(jī)誤差, 即最優(yōu)強(qiáng)迫向量。這種最優(yōu)強(qiáng)迫向量能刻畫ENSO預(yù)報(bào)的最大不確定性, 被廣泛應(yīng)用在集合預(yù)報(bào)中。但Barkmeijer等(2003)認(rèn)為這種最優(yōu)強(qiáng)迫向量的計(jì)算量巨大, 不適用于現(xiàn)實(shí)中的高維非線性數(shù)值模式。為了解決這一限制, 他們提出了線性強(qiáng)迫奇異向量(forcing singular vector, FSV)方法, 通過求解線性最優(yōu)問題得到一個(gè)不隨時(shí)間變化的傾向誤差。在此基礎(chǔ)上, Duan等(2013)將FSV方法進(jìn)行擴(kuò)展, 發(fā)展出非線性強(qiáng)迫奇異向量(non-linear forcing singular vector, NFSV)方法, 基于變分思想計(jì)算得到隨時(shí)間和空間變化的最優(yōu)傾向誤差。他們(Duan, 2014)將NFSV方法應(yīng)用于著名的Zebiak-Cane (ZC)模式, 對1980~2004年間三次典型的東部(Eastern-Pacific, EP)型El Ni?o事件和三次典型的中部(Central-Pacific, CP)型El Ni?o事件進(jìn)行模擬, 結(jié)果表明, 使用NFSV方法校正的ZC模式, 不僅能夠重現(xiàn)EP型El Ni?o事件, 并且很好地再現(xiàn)了三次CP型El Ni?o事件?;诖? Tao等(2019)進(jìn)一步建立了NFSV型傾向誤差預(yù)報(bào)模式, 將預(yù)報(bào)模式與中間型海洋—大氣耦合模式(Intermediate Complex Model, ICM)耦合, 構(gòu)建了ENSO的NFSV-ICM預(yù)報(bào)模式, 并對1997~2016年間的ENSO事件進(jìn)行預(yù)報(bào)實(shí)驗(yàn)。與ICM模式相比, NFSV-ICM具有更高的預(yù)報(bào)技巧, 有效預(yù)報(bào)時(shí)效從6個(gè)月提高至12個(gè)月。

        如果視每個(gè)模式格點(diǎn)的傾向誤差是一個(gè)參數(shù), 那么傾向誤差估計(jì)實(shí)質(zhì)等價(jià)于資料同化中的參數(shù)估計(jì)。因此, 除了變分方法外, 目前常用于資料同化的集合卡爾曼濾波器及其變種在理論上都可以應(yīng)用于傾向誤差估計(jì)(吳新榮, 2013)。此外, 與四維變分方法相比, 卡爾曼濾波器無需發(fā)展伴隨模式, 計(jì)算較為簡單, 同時(shí)考慮了預(yù)報(bào)誤差的動(dòng)力演變, 可以顯式地提供預(yù)報(bào)的初始擾動(dòng)(沈浙奇等,2016)。基于此, 本研究擬發(fā)展一種基于集合卡爾曼濾波器的模式傾向誤差估計(jì)算法, 通過同化海表面溫度異常(sea surface temperature anomaly, SSTA)數(shù)據(jù)來計(jì)算模式傾向誤差, 并考察其對ENSO模擬的影響。

        本文分為三個(gè)部分。第一部分介紹了研究所用的模式、觀測資料和模式傾向誤差估計(jì)算法的具體內(nèi)容。第二部分為模式傾向誤差的估計(jì)及其在ENSO模擬中的應(yīng)用。最后一部分是總結(jié)和討論。

        1 模式與方法

        1.1 Zebiak-Cane模式

        ZC模式(Cane, 1985, 1986; Zebiak, 1987)是拉蒙特-多爾蒂地球觀測站(Lamont-Doherty Earth Observatory, LDEO)的業(yè)務(wù)化預(yù)報(bào)模式, 顯式地描寫了Bjerknes-Wyrtki理論的物理本質(zhì)(徐輝, 2006), 因成功預(yù)報(bào)出了1986~1987年的El Ni?o事件而著名, 被廣泛地應(yīng)用于ENSO研究(Chen, 1995, 2004, 2008; Gao, 2020)。

        ZC模式由海洋斜壓模式、海洋表層模式和大氣模式耦合而成, 其中海洋斜壓模式?jīng)Q定海洋對風(fēng)強(qiáng)迫的響應(yīng), 海洋表層模式用于計(jì)算SSTA, 而大氣模式用來模擬風(fēng)對SSTA的響應(yīng)(Perigaud, 1996; 岳彩軍等, 2004)。海洋模式中, SSTA的發(fā)展方程(徐輝, 2006)為

        我們在ZC模式的SSTA訂正方程中引入代表模式傾向誤差項(xiàng)的參數(shù)(,,):

        式中, Δ為模式積分步長, 為10 d。

        1.2 觀測資料

        觀測資料采用的是美國國家海洋和大氣管理局(National Oceanic and Atmospheric Administration, NOAA) ERSST v5數(shù)據(jù)集內(nèi)的月平均SSTA數(shù)據(jù), 數(shù)據(jù)分辨率為2°×2°, 區(qū)域范圍選取101.25°E~73.125°W, 29°S~29°N, 與ZC模式區(qū)域一致, 時(shí)間長度選取1950年1月至2020年12月, 共852個(gè)月。需要對選取的SSTA數(shù)據(jù)進(jìn)行預(yù)處理: 為了與ZC模式網(wǎng)格分辨率一致, 將數(shù)據(jù)插值到5.625°×2.000°的網(wǎng)格分辨率, 然后使用100個(gè)符合均值為0, 方差為0.1 °C(Huang, 2016)的高斯分布的隨機(jī)數(shù)對其進(jìn)行擾動(dòng)。

        1.3 LETKF框架下的傾向誤差估計(jì)算法

        當(dāng)集合卡爾曼濾波器(ensemble Kalman filter, EnKF)同化技術(shù)被用于多維海氣模式時(shí), 使用有限的集合樣本來估計(jì)背景誤差協(xié)方差矩陣會(huì)經(jīng)常導(dǎo)致偽相關(guān)問題(Tang, 2016), 而局地集合變換卡爾曼濾波器(local ensemble transform Kalman filter, LETKF)同化技術(shù)(Hunt, 2007)能克服遠(yuǎn)距離虛假相關(guān)這一缺點(diǎn), 在海洋和氣象領(lǐng)域有著廣泛的應(yīng)用。例如, Ruckstuhl等(2020)采用LETKF技術(shù)估計(jì)了COSMO-DE模式中的二維粗糙長度參數(shù), 改進(jìn)了模式對降水和云的預(yù)報(bào)。Gao等(2021)基于LETKF技術(shù)建立了ZC模式的同化系統(tǒng), 對模式的關(guān)鍵參數(shù)進(jìn)行估計(jì), 發(fā)現(xiàn)參數(shù)估計(jì)降低了模式誤差, 提高了模式對ENSO的預(yù)報(bào)效果。因此, 本文選取LETKF同化技術(shù)作為估計(jì)模式傾向誤差參數(shù)的方法。

        由于現(xiàn)實(shí)中沒有參數(shù)的直接觀測, 所以參數(shù)的估計(jì)是通過狀態(tài)向量附加技術(shù)(state vector augmentation technique, SVAT)實(shí)現(xiàn)的, 具體來說, 就是把參數(shù)視作狀態(tài)向量的一部分, 在進(jìn)行參數(shù)估計(jì)時(shí)使用狀態(tài)向量的觀測算子(吳新榮, 2013)。計(jì)算步驟如下:

        其中,表示將模式格點(diǎn)的變量投影到觀測空間上的線性觀測算子, 插值過程基于雙線性插值方法。

        (2) 計(jì)算卡爾曼增益矩陣

        (3) 計(jì)算分析場的集合平均:

        其中,obs代表狀態(tài)向量的觀測。

        (4) 為了更新每個(gè)集合成員, LETKF還需計(jì)算分析集合擾動(dòng)。通過一個(gè)變換矩陣, 可將背景場擾動(dòng)b轉(zhuǎn)換成分析場擾動(dòng)a, 得到分析集合的擾動(dòng)為

        基于以上計(jì)算步驟, 在給定待估參數(shù)初始場的前提下, 可計(jì)算得到參數(shù)在每一步的分析值, 取集合均值為最終分析結(jié)果。如此, 得到傾向誤差在不同模式格點(diǎn)和時(shí)間點(diǎn)的估計(jì)值。

        2 結(jié)果與分析

        2.1 模式傾向誤差的計(jì)算及分析

        2.1.1 傾向誤差估計(jì)初始場設(shè)計(jì) 使用LETKF同化技術(shù)對參數(shù)進(jìn)行估計(jì)時(shí)需要給出待估參數(shù)的初始值。由于參數(shù)代表的是模式的傾向誤差, 沒有實(shí)際觀測數(shù)據(jù), 因此, 我們將ZC模式只進(jìn)行狀態(tài)估計(jì)(only state estimate, OSE)時(shí)同化輸出的SSTA分析值與觀測值的偏差(error)即模擬誤差作為待估參數(shù)的初始場, 這種偏差一定程度上代表了模式誤差, 并使用一組符合均值為0, 方差為0.1 °C的高斯分布的隨機(jī)數(shù)對其進(jìn)行擾動(dòng), 生成100個(gè)數(shù)據(jù)集合, 作為參數(shù)的初始集合。同樣對初始時(shí)刻的SSTA進(jìn)行擾動(dòng), 生成狀態(tài)變量SSTA的初始集合。

        圖1a展示了SSTA觀測值和OSE實(shí)驗(yàn)分析值的Ni?o 3.4指數(shù)時(shí)間序列, 圖1b為觀測值減去實(shí)驗(yàn)分析值所得的差值。可以看出, 運(yùn)用LETKF同化技術(shù)得到的SSTA分析值和觀測值具有很好的相關(guān)性, 如表1所示, 相關(guān)系數(shù)高達(dá)0.97, 相比同化前0.86的相關(guān)系數(shù)有較大提高, 且平均絕對誤差由同化前的0.32降低到了0.15, 這說明了LETKF同化技術(shù)的有效性。但從圖1b可見, 分析值和觀測值之間仍存在一定的差異, 偏差多為負(fù)值, 表明模式常常高估于觀測, 尤其是在冷事件的模擬上, 差值最大達(dá)到了0.8 °C, 且冷事件模擬的偏差普遍高于暖事件的偏差。

        圖1 海表面溫度異常(sea surface temperature anomaly, SSTA)觀測值(紅)與OSE實(shí)驗(yàn)分析值(藍(lán))的Ni?o 3.4指數(shù)時(shí)間序列(a)以及兩者的差值(b)

        注: OSE: only state estimate, 僅進(jìn)行狀態(tài)估計(jì)實(shí)驗(yàn)

        表1 ZC模式同化實(shí)驗(yàn)效果對比

        Tab.1 Comparison in the effect of the assimilation experiment in ZC model

        2.1.2 傾向誤差計(jì)算 得到傾向誤差估計(jì)所需的初始場之后, 使用LETKF同化技術(shù)對傾向誤差進(jìn)行估計(jì)。通常, 基于卡爾曼濾波器的同化實(shí)驗(yàn)是從一個(gè)初始時(shí)刻出發(fā), 在模式的積分過程中將觀測數(shù)據(jù)同化進(jìn)模式, 對模擬結(jié)果進(jìn)行校正, 以輸出更為接近觀測狀態(tài)的分析值。同樣, 在進(jìn)行參數(shù)估計(jì)時(shí), 也是從某一初始時(shí)刻起, 通過同化觀測資料對參數(shù)進(jìn)行調(diào)整。例如, 以1950年1月為起始月, 基于該月的初始場, 同化SSTA數(shù)據(jù)對模式傾向誤差進(jìn)行估計(jì), 我們發(fā)現(xiàn), 隨著同化的進(jìn)行, 參數(shù)分析值在經(jīng)過一段時(shí)間(大約10 a)的調(diào)整后趨于穩(wěn)定。為了進(jìn)一步驗(yàn)證該結(jié)論, 基于1950年1月至2011年1月間共733個(gè)月的初始場, 計(jì)算對應(yīng)的模式傾向誤差。結(jié)果表明, 傾向誤差估計(jì)值均經(jīng)過約10 a的調(diào)整后趨于穩(wěn)定, 且不同起始月份得到的穩(wěn)定值不同。為了表述簡潔, 圖2僅展示了1950年1月至1951年1月間不同起始月份對應(yīng)的傾向誤差估計(jì)值。因此, 我們可以認(rèn)為, 模式傾向誤差與起始月份所給定的參數(shù)初始場之間存在一定的相關(guān)關(guān)系。從式(5)和式(6)也可以看出, 參數(shù)的分析值會(huì)受到初始值的影響。我們將穩(wěn)定后的估計(jì)值作為對應(yīng)于不同起始月的模式傾向誤差估計(jì)值, 并基于這些估計(jì)值展開分析。

        為了檢驗(yàn)?zāi)J絻A向誤差估計(jì)實(shí)驗(yàn)中同化的有效性, 將同化得到的傾向誤差估計(jì)值加入ZC模式, 重復(fù)OSE實(shí)驗(yàn), 發(fā)現(xiàn)SSTA分析值和觀測值的相關(guān)性較不含模式傾向誤差項(xiàng)的OSE實(shí)驗(yàn)略有提高, 由0.970提高到了0.972, 且平均絕對誤差進(jìn)一步由0.15降低至0.13, 體現(xiàn)了LETKF同化技術(shù)在估計(jì)模式傾向誤差上的有效作用。

        2.1.3 傾向誤差特征 由于Ni?o3.4區(qū)域是表征ENSO特征的一個(gè)關(guān)鍵區(qū)域, 所以圖3展示了Ni?o 3.4區(qū)域平均的傾向誤差估計(jì)值及其對應(yīng)的誤差(即觀測值和OSE實(shí)驗(yàn)分析值的差值)??梢钥吹? 傾向誤差估計(jì)值與誤差相關(guān)較好, 兩者的相關(guān)系數(shù)為0.71, 表明傾向誤差估計(jì)值可以在一定程度上代表模式誤差。

        圖2 分別以不同起始月份為起點(diǎn)得出的模式傾向誤差估計(jì)值在模式區(qū)域的平均

        注: 圖中所示起始月份為1950年1月至1951年1月

        圖3 誤差(黑)和模式傾向誤差估計(jì)值(藍(lán))在Ni?o 3.4區(qū)域的平均

        注:表示相關(guān)系數(shù)

        為了進(jìn)一步比較兩者在赤道地區(qū)的演變特征, 圖4展示了誤差(a)和傾向誤差估計(jì)值(b)在5°S~5°N區(qū)域平均的時(shí)間—經(jīng)度剖面圖。可以看到, 在時(shí)間和空間上, 兩場都有較好的對應(yīng)關(guān)系。當(dāng)誤差場中出現(xiàn)一個(gè)正(負(fù))差值時(shí), 相應(yīng)地, 在傾向誤差場也會(huì)出現(xiàn)正(負(fù))值, 兩場的相關(guān)系數(shù)達(dá)到了0.63。

        對誤差和模式傾向誤差估計(jì)值進(jìn)行經(jīng)驗(yàn)正交函數(shù)(empirical orthogonal function, EOF)分解, 從而可以了解兩場的時(shí)空演變特征。圖5為誤差和模式傾向誤差估計(jì)值EOF分解所得的4個(gè)主要特征模態(tài)的空間分布, 圖6則給出了對應(yīng)的時(shí)間系數(shù)。從圖5a可以看出, 誤差的第一特征模態(tài)呈東正西負(fù)的空間分布, 正值區(qū)主要出現(xiàn)在中部太平洋地區(qū), 并逐漸向東、向赤道外擴(kuò)散, 負(fù)值區(qū)在中、西部太平洋呈“馬蹄形”環(huán)繞正值區(qū)。傾向誤差估計(jì)值的第一特征模態(tài)(圖5e)也呈現(xiàn)了一個(gè)東正西負(fù)的分布, 正值主要集中在赤道附近的中、東部太平洋地區(qū), 負(fù)值在北太平洋中部最大, 向西、向南逐漸減小, 整體呈現(xiàn)出一個(gè)類似El Ni?o事件分布的空間型, 與誤差的第一特征模態(tài)具有相似的空間特征, 兩場的相關(guān)系數(shù)達(dá)到了0.72。同樣地, 從圖6a可見, 第一特征模態(tài)對應(yīng)時(shí)間系數(shù)的相關(guān)高達(dá)0.87, 說明在時(shí)空兩個(gè)維度, 兩場的演變都是較為一致的。

        進(jìn)一步比較第二、三、四特征模態(tài), 可以看到, 誤差的第二特征模態(tài)(圖5b)在西太平洋和北太平洋地區(qū)分別出現(xiàn)了較明顯的正值和負(fù)值, 這些顯著的特征在傾向誤差估計(jì)值的第二特征模態(tài)(圖5f)也有所體現(xiàn)。此外, 兩場的第三特征模態(tài)(圖5c和5g)空間分布高度相似, 自西向東依次呈現(xiàn)正、負(fù)、正、負(fù)的空間結(jié)構(gòu), 且兩場的相關(guān)系數(shù)達(dá)到了0.58, 圖6c對應(yīng)的時(shí)間序列相關(guān)為0.77。在第四特征模態(tài)(圖5d和5h)中, 誤差場的負(fù)值集中在近赤道中、東部太平洋地區(qū), 正值分布在西太平洋地區(qū)和赤道兩側(cè), 呈現(xiàn)了類似La Ni?a事件的空間型, 而在估計(jì)值場中, 近赤道西太平洋地區(qū)出現(xiàn)了負(fù)值, 其他地區(qū)的空間分布均與誤差場相似??偟膩碚f, 四個(gè)特征模態(tài)的空間場和對應(yīng)時(shí)間系數(shù)序列的相關(guān)均達(dá)到了0.5以上, 進(jìn)一步論證了模式傾向誤差估計(jì)值與誤差之間的高度相關(guān)關(guān)系, 這體現(xiàn)了LETKF同化技術(shù)用于估計(jì)模式傾向誤差的有效性。

        圖4 模擬誤差(a)和模式傾向誤差估計(jì)值(b)的時(shí)間-經(jīng)度分布圖

        圖5 模擬誤差(a, b, c, d)和傾向誤差估計(jì)值(e, f, g, h) EOF分解所得的4個(gè)主要特征模態(tài)空間分布

        注:EOF1表示EOF分解的第一特征模態(tài); EOF2表示EOF分解的第二特征模態(tài); EOF3表示EOF分解的第三特征模態(tài); EOF4表示EOF分解的第四特征模態(tài)

        圖6 模擬誤差和傾向誤差估計(jì)值EOF分解所得的4個(gè)主要特征模態(tài)對應(yīng)的時(shí)間系數(shù)

        注:表示相關(guān)系數(shù)

        2.2 考慮模式傾向誤差的ENSO模擬

        上文分析了模式傾向誤差估計(jì)值的時(shí)空特征, 揭示了其與模擬誤差之間的高度相關(guān)關(guān)系。接下來, 展開ENSO模擬實(shí)驗(yàn), 評(píng)估模式傾向誤差估計(jì)在改善模式模擬效果方面的作用。將2.1.3中模式傾向誤差EOF分解所得的4個(gè)主要空間模態(tài), 在每一個(gè)時(shí)間步乘以一個(gè)隨機(jī)的系數(shù)[~(-1,1)], 再進(jìn)行加和, 將合成的空間場加入模式傾向方程。為了與ZC模式進(jìn)行區(qū)分, 將帶有模式傾向誤差的ZC模式稱為F-ZC模式。自由積分F-ZC模式150 a, 積分結(jié)果用來考察模式傾向誤差估計(jì)對ENSO模擬的影響。

        2.2.1 赤道地區(qū)海溫距平場的時(shí)間–經(jīng)度分布特征 圖7為ZC模式和F-ZC模式自由積分的模擬結(jié)果沿赤道地區(qū)的時(shí)間-經(jīng)度分布圖??梢钥闯? ZC模式模擬的El Ni?o事件十分顯著, 絕大多數(shù)事件的峰值高達(dá)4 °C, 相對來說, 模擬的La Ni?a事件強(qiáng)度偏低, 出現(xiàn)不對稱的振蕩。F-ZC模式的模擬同樣準(zhǔn)確地模擬出了這種不對稱性, 且第2.1.1節(jié)中提到, ZC模式在冷事件的模擬上常高于觀測, 而F-ZC模式在一定程度上修正了這種高估, 其模擬的La Ni?a事件的強(qiáng)度和頻率均有所增加。此外, ZC模式模擬的暖事件在空間形態(tài)分布上均呈現(xiàn)出EP型El Ni?o事件這一單一形態(tài)特征, 而F-ZC模式不僅模擬出了EP型El Ni?o事件, 也模擬出了類似于CP型El Ni?o事件的空間型, 如第14 a、第77 a和第98 a的三次模擬, 盡管強(qiáng)度較低, 但模擬結(jié)果呈現(xiàn)出多樣性。

        2.2.2 Ni?o 3.4指數(shù)周期性的模擬 為了考察模式模擬ENSO周期的能力, 對Ni?o 3.4區(qū)SSTA時(shí)間序列進(jìn)行功率譜分析, 并將其與觀測進(jìn)行對比。圖8a、8b、8c分別是SSTA觀測值、ZC模式和F-ZC模式自由積分模擬結(jié)果的功率譜分析。整體上, 觀測和兩組模擬實(shí)驗(yàn)的周期均與通常所認(rèn)為的ENSO 2~7 a的主要周期一致。具體地, 觀測除存在3.63 a的顯著周期外, 也顯示了4.88 a的一個(gè)次要周期。而ZC模式和F-ZC模式的模擬結(jié)果都顯示出一個(gè)明顯的主周期, 分別為3.97和4.16 a, 在F-ZC模式中開展多次模擬實(shí)驗(yàn), 模擬結(jié)果的主周期均穩(wěn)定在4.16 a, 在觀測所涵蓋的范圍內(nèi)。且相比較之下, F-ZC模式模擬結(jié)果的主周期更加接近于觀測的次要周期。

        2.2.3 ENSO鎖相的模擬 除強(qiáng)度和周期外, ENSO的季節(jié)鎖相也是衡量ENSO模擬效果的一個(gè)重要指標(biāo)。盡管各個(gè)ENSO事件的振幅各不相同, 但ENSO相位非常相似, SSTA的峰值總是出現(xiàn)在冬季。為了更清楚地描述El Ni?o事件和La Ni?a事件的相位特征, 分別選取觀測、ZC模式模擬和F-ZC模式模擬中顯著的El Ni?o事件和La Ni?a事件進(jìn)行合成(中華人民共和國國家質(zhì)量監(jiān)督檢驗(yàn)檢疫總局等, 2017)。

        圖7 ZC模式(a, b, c)和F-ZC模式(d, e, f) 150 a自由積分模擬SSTA沿赤道緯圈的時(shí)間—經(jīng)度分布圖

        圖8 觀測(a)、ZC模式(b)和F-ZC模式(c)模擬結(jié)果在Ni?o3.4區(qū)的SSTA時(shí)間序列功率譜

        圖9a和9b分別繪制了暖事件和冷事件前后Ni?o3.4指數(shù)的時(shí)間演變, 其中, (0)年表示事件達(dá)到峰值的當(dāng)年。從觀測可以看到, 暖事件發(fā)生時(shí), Ni?o3.4區(qū)域逐漸變暖, 從當(dāng)年的春季開始進(jìn)入El Ni?o事件, 在冬季12月份達(dá)到峰值, 之后逐漸減弱, 直至次年夏季事件結(jié)束。ZC模式模擬輸出的El Ni?o事件峰值出現(xiàn)在當(dāng)年的秋季9月份, 相比較之下, F-ZC模式模擬的El Ni?o事件峰值出現(xiàn)的時(shí)間更加接近于觀測的峰值時(shí)間, 出現(xiàn)在冬季12月份。在事件強(qiáng)度方面, F-ZC模式模擬的事件強(qiáng)度整體上也是更接近于觀測的, 尤其是在達(dá)到峰值前的事件發(fā)展期, ZC模式在當(dāng)年年初就已經(jīng)出現(xiàn)了一個(gè)中等強(qiáng)度的El Ni?o事件, F-ZC模式則表現(xiàn)更好。在冷事件的模擬上, 觀測的峰值也出現(xiàn)在冬季12月份, 從當(dāng)年年初開始變冷后, 在夏季發(fā)生La Ni?a事件, 逐步發(fā)展變強(qiáng), 降至最低溫之后開始升溫, 直至次年夏季結(jié)束一次事件。但兩個(gè)模式的模擬實(shí)驗(yàn)均未表現(xiàn)出明顯的相位鎖定, 均在當(dāng)年的冬季和次年的夏季出現(xiàn)了兩個(gè)不同的峰值, 且模擬的冷事件強(qiáng)度和持續(xù)時(shí)間都偏低。從整體趨勢上看, F-ZC模式的模擬輸出更接近于觀測, 對次年夏季的第二次La Ni?a事件有一定程度的校正。重復(fù)開展多次模擬實(shí)驗(yàn), 在季節(jié)鎖相的模擬上, F-ZC模式表現(xiàn)均優(yōu)于ZC模式。

        2.2.4 海溫距平場時(shí)空演變特征 為了進(jìn)一步比較模式對主要空間模態(tài)的模擬效果, 分別對觀測的SSTA值、ZC模式和F-ZC模式模擬所得的SSTA分析值進(jìn)行EOF分解, 圖10給出了觀測及模式模擬的SSTA主要特征模態(tài)。由于自由積分模擬實(shí)驗(yàn)屬于理想實(shí)驗(yàn), 沒有對應(yīng)的真實(shí)時(shí)間, 因此, 這里僅比較主要空間模態(tài)的分布。

        圖9 合成的El Ni?o事件(a)和La Ni?a事件(b)在Ni?o 3.4區(qū)的SSTA演變

        注: (0)表示事件達(dá)到峰值的當(dāng)年; (+1)表示次年; (+2)表示第三年

        觀測的SSTA第一特征向量(圖10a)呈東正西負(fù)的分布, 正值區(qū)從美洲西岸逐漸向西延伸到日界線附近, 而負(fù)值區(qū)則從西太平洋地區(qū)呈“C字型”向東北及東南延伸。這一模態(tài)與El Ni?o事件發(fā)生時(shí)的熱帶太平洋SSTA分布型相似, 可稱為ENSO模態(tài), 它可以解釋SSTA總方差的72%。觀測的第二特征模態(tài)(圖10b)與第一特征模態(tài)恰好相反, 是一個(gè)西部為正值而東部為負(fù)值的模態(tài), 赤道西太平洋的正值最大并向東部延展, 赤道東太平洋的負(fù)值最大并向西部延展, 最終在中太地區(qū)相遇, 使從西太而來的正值擴(kuò)散至南、北太平洋, 呈“樹杈狀”分布。

        在模擬實(shí)驗(yàn)的SSTA場EOF分解中, ZC模式和F-ZC模式模擬的SSTA主要特征向量均與觀測較為接近??梢钥吹? 兩個(gè)模擬實(shí)驗(yàn)的第一特征模態(tài)(圖10c和10e)非常相似, 均為一個(gè)典型的中、東部型El Ni?o事件空間形態(tài), 正值區(qū)從東南太平洋逐漸向西發(fā)展。同樣, 和觀測類似, 第一特征模態(tài)的貢獻(xiàn)也是最大的, 分別達(dá)到了82%和78%, 這表明ZC模式能較好地模擬出ENSO事件最主要的形態(tài)特征, 同樣加入了模式傾向誤差參數(shù)的F-ZC模式仍保持住了這一優(yōu)勢。進(jìn)一步比較第二特征模態(tài)(圖10d和10f), 發(fā)現(xiàn)ZC模式和F-ZC模式的模擬都表現(xiàn)出了從東太地區(qū)向西延展的負(fù)值, 與觀測是一致的, 但兩者均未模擬出從西太起逐步向東傳播的正值。盡管如此, 與ZC模式相比, F-ZC模式弱化了在西太的負(fù)值, 強(qiáng)化了在中太的正值, 模擬出了更加接近觀測的第二特征模態(tài), 說明模式傾向誤差的加入對ENSO模擬效果產(chǎn)生了正面影響。重復(fù)開展多次模擬實(shí)驗(yàn), 實(shí)驗(yàn)結(jié)果穩(wěn)定一致。

        圖10 觀測(a, b)、ZC模式(c, d)和F-ZC模式(e, f)模擬的SSTA值EOF分解的前兩個(gè)主要特征模態(tài)空間分布

        注: EOF1表示EOF分解的第一特征模態(tài); EOF2表示EOF分解的第二特征模態(tài)

        3 結(jié)論

        本研究基于ZC模式, 延續(xù)Roads等(Roads, 1987; Moore, 1999; Barkmeijer,2003; Duan, 2013, 2014; Tao, 2019)將模式傾向誤差視作模式參數(shù)的思想, 利用LETKF同化技術(shù)對模式傾向誤差進(jìn)行估計(jì), 發(fā)展出一種新的模式傾向誤差的計(jì)算方法。在這基礎(chǔ)上, 基于1950~2020年的SSTA觀測數(shù)據(jù), 計(jì)算得到了對應(yīng)于1950~2011年的模式傾向誤差估計(jì)值。通過相關(guān)性分析和EOF分析發(fā)現(xiàn), 模式傾向誤差估計(jì)值和模擬誤差之間具有高度相關(guān)關(guān)系。進(jìn)一步地將傾向誤差估計(jì)項(xiàng)加入ZC模式, 重復(fù)開展了多次ENSO模擬實(shí)驗(yàn)。通過與ZC模式的對比, 發(fā)現(xiàn)在150 a的自由積分實(shí)驗(yàn)中, 帶有模式傾向誤差的F-ZC模式均能進(jìn)一步改善對ENSO事件的模擬效果, 在振幅、周期、相位鎖定和主要空間形態(tài)等方面, 模擬與觀測的相關(guān)性均有所提高。

        本研究提出的這種新的估計(jì)模式傾向誤差的方法, 計(jì)算簡便高效, 可便捷地應(yīng)用于其他模式中, 具有較高的推廣應(yīng)用價(jià)值。但不同模式中存在的傾向誤差也是不同的, 盡管本方法在ZC模式中得到了較好的應(yīng)用, 但在其他更為復(fù)雜的耦合模式中的應(yīng)用效果如何, 仍需要進(jìn)一步的工作來驗(yàn)證, 這也是本研究接下來的一個(gè)主要工作方向。而模式傾向誤差估計(jì)方法存在的一些不足也有待改進(jìn), 比如參數(shù)估計(jì)中常常遇到的濾波發(fā)散問題。此外, 模擬的最終目的仍是為了提高預(yù)報(bào)效果, 該方法是否能用于提高ENSO的預(yù)報(bào)技巧, 仍有待檢驗(yàn)。

        王琳, 張燦影, 於維櫻, 等, 2019. 厄爾尼諾-南方濤動(dòng)(ENSO)研究的戰(zhàn)略部署與研究熱點(diǎn)[J]. 世界科技研究與發(fā)展, 41(1): 32-43.

        中華人民共和國國家質(zhì)量監(jiān)督檢驗(yàn)檢疫總局, 中國國家標(biāo)準(zhǔn)化管理委員會(huì), 2017. 厄爾尼諾/拉尼娜事件判別方法: GB/T 33666-2017[S]. 北京: 中國標(biāo)準(zhǔn)出版社: 3.

        連濤, 陳大可, 唐佑民, 2017. 2014~2016厄爾尼諾事件的機(jī)制分析[J]. 中國科學(xué): 地球科學(xué), 47(9): 1014-1026.

        吳新榮, 2013. 耦合氣候模式中的參數(shù)估計(jì)研究[D]. 廣州: 中國科學(xué)院大學(xué)(南海海洋研究所): 7-12.

        沈浙奇, 唐佑民, 高艷秋, 2016. 集合資料同化方法的理論框架及其在海洋資料同化的研究展望[J]. 海洋學(xué)報(bào), 38(3): 1-14.

        岳彩軍, 陸維松, 李清泉, 等, 2004. Zebiak-Cane海氣耦合模式研究進(jìn)展[J]. 熱帶氣象學(xué)報(bào), 20(6): 723-730.

        夏詠華, 2001. “厄爾尼諾”事件及其對氣候和航海活動(dòng)的影響[D]. 大連: 大連海事大學(xué): 10-13.

        徐輝, 2006. Zebiak-Cane ENSO預(yù)報(bào)模式的可預(yù)報(bào)性問題研究[D]. 北京: 中國科學(xué)院研究生院(大氣物理研究所): 10-13.

        穆穆, 任宏利, 2017. 2014~2016年超強(qiáng)厄爾尼諾事件研究及其預(yù)測給予我們的啟示[J]. 中國科學(xué): 地球科學(xué), 47(9): 993-995.

        BARKMEIJER J, IVERSEN T, PALMER T N, 2003. Forcing singular vectors and other sensitive model structures [J]. Quarterly Journal of the Royal Meteorological Society, 129(592): 2401-2423.

        CANE M A, ZEBIAK S E, 1985. A theory for El Ni?o and the southern oscillation [J]. Science, 228(4073): 1085-1087.

        CANE M A, ZEBIAK S E, DOLAN S C, 1986. Experimental forecasts of El Ni?o [J]. Nature, 321(6073): 827-832.

        CHEN D K, CANE M A, 2008. El Ni?o prediction and predictability [J]. Journal of Computational Physics, 227(7): 3625-3640.

        CHEN D K, CANE M A, KAPLAN A,, 2004. Predictability of El Ni?o over the past 148 years [J]. Nature, 428(6984): 733-736.

        CHEN D K, CANE M A, ZEBIAK S E,, 1998. The impact of sea level data assimilation on the Lamont model prediction of the 1997/98 El Ni?o [J]. Geophysical Research Letters, 25(15): 2837-2840.

        CHEN D K, ZEBIAK S E, BUSALACCHI A J,, 1995. An improved procedure for EI Ni?o forecasting: implications for predictability [J]. Science, 269(5231): 1699-1702.

        CHEN D K, ZEBIAK S E, CANE M A,, 1997. Initialization and predictability of a coupled ENSO forecast model [J]. Monthly Weather Review, 125(5): 773-788.

        DUAN W S, TIAN B, XU H, 2014. Simulations of two types of El Ni?o events by an optimal forcing vector approach [J]. Climate Dynamics, 43(5): 1677-1692.

        DUAN W S, ZHOU F F, 2013. Non-linear forcing singular vector of a two-dimensional quasi-geostrophic model [J]. Tellus A: Dynamic Meteorology and Oceanography, 65(1): 18452.

        GAO Y Q, LIU T, SONG X S,, 2020. An extension of LDEO5 model for ENSO ensemble predictions [J]. Climate Dynamics, 55(11/12): 2979-2991.

        GAO Y Q, TANG Y M, SONG X S,, 2021. Parameter estimation based on a local ensemble transform Kalman filter applied to El Ni?o–southern oscillation ensemble prediction [J]. Remote Sensing, 13(19): 3923.

        HU J Y, DUAN W S, 2016. Relationship between optimal precursory disturbances and optimally growing initial errors associated with ENSO events: implications to target observations for ENSO prediction [J]. Journal of Geophysical Research: Oceans, 121(5): 2901-2917.

        HUANG B Y, THORNE P W, SMITH T M,, 2016. Further exploring and quantifying uncertainties for extended reconstructed sea surface temperature (ERSST) version 4 (v4) [J]. Journal of Climate, 29(9): 3119-3142.

        HUNT B R, KOSTELICH E J, SZUNYOGH I, 2007. Efficient data assimilation for spatiotemporal chaos: a local ensemble transform Kalman filter [J]. Physica D: Nonlinear Phenomena, 230(1/2): 112-126.

        MOORE A M, KLEEMAN R, 1999. Stochastic forcing of ENSO by the intraseasonal oscillation [J]. Journal of Climate, 12(5): 1199-1220.

        MU M, DUAN W S, CHEN D K,, 2015. Target observations for improving initialization of high-impact ocean-atmospheric environmental events forecasting [J]. National Science Review, 2(2): 226-236.

        PERIGAUD C, DEWITTE B, 1996. El Ni?o-La Ni?a events simulated with Cane and Zebiak’s model and observed with satellite ordata. Part I: model data comparison [J]. Journal of Climate, 9(1): 66-84.

        QI Q Q, DUAN W S, ZHENG F,, 2017. On the “spring predictability barrier” for strong El Ni?o events as derived from an intermediate coupled model ensemble prediction system [J]. Science China Earth Sciences, 60(9): 1614-1631.

        ROADS J O, 1987. Predictability in the extended range [J]. Journal of the Atmospheric Sciences, 44(23): 3495-3527.

        RUCKSTUHL Y, JANJI? T, 2020. Combined state-parameter estimation with the LETKF for convective-scale weather forecasting [J]. Monthly Weather Review, 148(4): 1607-1628.

        TANG Y M, SHEN Z Q, GAO Y Q, 2016. An introduction to ensemble-based data assimilation method in the earth sciences [M] // LEE D B, BURG T, VOLOS C. Nonlinear Systems: Design, Analysis, Estimation and Control. New York, USA: Intech: 153-187.

        TAO L J, DUAN W S, 2019. Using a nonlinear forcing singular vector approach to reduce model error effects in ENSO forecasting [J]. Weather and Forecasting, 34(5): 1321-1342.

        WU X R, HAN G J, ZHANG S Q,, 2016. A study of the impact of parameter optimization on ENSO predictability with an intermediate coupled model [J]. Climate Dynamics, 46(3/4): 711-727.

        YEH S W, KUG J S, DEWITTE B,, 2009. El Ni?o in a changing climate [J]. Nature, 461(7263): 511-514.

        ZEBIAK S E, CANE M A, 1987. A model El Ni?o–southern oscillation [J]. Monthly Weather Review, 115(10): 2262-2278.

        ZHENG F, WANG H, ZHU J, 2009. ENSO ensemble prediction: initial error perturbations vs. model error perturbations [J]. Chinese Science Bulletin, 54(14): 2516-2523.

        A NEW ALGORITHM OF ESTIMATION FOR MODEL TENDENCY ERRORS AND THE APPLICATION IN ENSO SIMULATION

        HE Qun1, 2, GAO Yan-Qiu2, 3, TANG You-Min2, 3, 4, ZHANG Ji-Cai1

        (1. Ocean College, Zhejiang University, Zhoushan 316021, China; 2. State Key Laboratory of Satellite Ocean Environment Dynamics, Second Institute of Oceanography, Ministry of Natural Resources, Hangzhou 310012, China; 3. Southern Marine Science and Engineering Guangdong Laboratory, Zhuhai 519082, China; 4. College of Oceanography, Hohai University, Nanjing 210098, China)

        Climate models are important tools for us to understand, simulate and forecast the evolution of the climate. However, even with the current state-of-the-art coupled models, due to the inevitable systematic errors in the tendency equation of model, the model tendency error, the simulations and forecasts are still far from the true state of the atmosphere/ocean. Therefore, reducing the model tendency error is of great significance to improve the simulation and forecasting effect of the model. A novel algorithm was developed for estimating the tendency error of a model using assimilation technique with local ensemble transform Kalman filter (LETKF). The new algorithm was applied to the Zebiak-Cane (ZC) model to estimate the space-time dependent tendency error by assimilating the observed data of sea surface temperature anomaly (SSTA), and the calculated tendency error was used to correct the model, and then an integral simulation was carried out. Results reveal a high correlation between the tendency error and the simulation error of the ZC model. The corrected model improved some important characteristics of the simulation of El Ni?o-Southern Oscillation (ENSO). Overall, the new algorithm is very effective and simple computationally, shows good application value in ENSO simulation, and can be easily applied to various models, and thus shall be promoted.

        model tendency error; parameter estimation; local ensemble transform Kalman filter; Zebiak-Cane model

        P731

        10.11693/hyhz20220100009

        *國家重點(diǎn)研發(fā)計(jì)劃“高影響海氣環(huán)境事件預(yù)報(bào)模式的高分辨率海洋資料同化系統(tǒng)研發(fā)”, 2017YFA0604202號(hào); 南方海洋科學(xué)與工程廣東省實(shí)驗(yàn)室(珠海)自主科研項(xiàng)目, SML2021SP314號(hào); 自然資源部第二海洋研究所基本科研業(yè)務(wù)費(fèi)專項(xiàng)項(xiàng)目資助, JG1809號(hào)。何 群, 碩士研究生, E-mail: hequn@zju.edu.cn

        高艷秋, 助理研究員, E-mail: gaoyanqiu@sio.org.cn

        2022-01-12,

        2022-03-10

        猜你喜歡
        模態(tài)特征
        抓住特征巧觀察
        新型冠狀病毒及其流行病學(xué)特征認(rèn)識(shí)
        如何表達(dá)“特征”
        不忠誠的四個(gè)特征
        抓住特征巧觀察
        車輛CAE分析中自由模態(tài)和約束模態(tài)的應(yīng)用與對比
        國內(nèi)多模態(tài)教學(xué)研究回顧與展望
        高速顫振模型設(shè)計(jì)中顫振主要模態(tài)的判斷
        基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識(shí)別
        由單個(gè)模態(tài)構(gòu)造對稱簡支梁的抗彎剛度
        欧美肥胖老妇做爰videos| 视频精品熟女一区二区三区| 国产精品片211在线观看| 99久久亚洲精品加勒比| 国产精品美女久久久久av福利| 伊人网在线视频观看| 伊人久久大香线蕉av不变影院| 狠狠色狠狠色综合日日92| 2019年92午夜视频福利| 国产精品videossex久久发布| 亚洲天堂一区二区三区视频| 久久精品中文字幕极品| 超清精品丝袜国产自在线拍| 麻神在线观看免费观看| 亚洲视频一区二区久久久| 最新亚洲人成网站在线| 中文字幕丰满乱子无码视频| 亚洲精品中文字幕91| 在线观看免费a∨网站| 色偷偷亚洲精品一区二区| 99亚洲男女激情在线观看| 亚洲乱码中文字幕一线区 | 色吧噜噜一区二区三区| 亚洲成在人线av| 98色婷婷在线| 高清国产亚洲va精品| 国产精品免费看久久久无码| 97在线视频免费| 爱情岛论坛亚洲永久入口口| 亚洲高清美女久久av| 国产综合久久久久影院| 久久乐国产精品亚洲综合| 久久亚洲综合亚洲综合| 色综合天天网| 亚洲av成人无码久久精品老人| 免费视频成人 国产精品网站| 亚洲国产天堂一区二区三区| 色综合天天综合欧美综合| 亚洲熟女国产熟女二区三区| 伊人网视频在线观看| 黑人玩弄漂亮少妇高潮大叫|