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

    一種新的模式傾向誤差估計算法及其在ENSO模擬中的應用*

    2022-09-21 02:37:10高艷秋唐佑民張繼才
    海洋與湖沼 2022年5期
    關鍵詞:估計值模態(tài)觀測

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

    一種新的模式傾向誤差估計算法及其在ENSO模擬中的應用*

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

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

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

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

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

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

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

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

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

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

    1 模式與方法

    1.1 Zebiak-Cane模式

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

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

    我們在ZC模式的SSTA訂正方程中引入代表模式傾向誤差項的參數(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ū)域一致, 時間長度選取1950年1月至2020年12月, 共852個月。需要對選取的SSTA數(shù)據(jù)進行預處理: 為了與ZC模式網(wǎng)格分辨率一致, 將數(shù)據(jù)插值到5.625°×2.000°的網(wǎng)格分辨率, 然后使用100個符合均值為0, 方差為0.1 °C(Huang, 2016)的高斯分布的隨機數(shù)對其進行擾動。

    1.3 LETKF框架下的傾向誤差估計算法

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

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

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

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

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

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

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

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

    2 結(jié)果與分析

    2.1 模式傾向誤差的計算及分析

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

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

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

    注: OSE: only state estimate, 僅進行狀態(tài)估計實驗

    表1 ZC模式同化實驗效果對比

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    圖10 觀測(a, b)、ZC模式(c, d)和F-ZC模式(e, f)模擬的SSTA值EOF分解的前兩個主要特征模態(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同化技術對模式傾向誤差進行估計, 發(fā)展出一種新的模式傾向誤差的計算方法。在這基礎上, 基于1950~2020年的SSTA觀測數(shù)據(jù), 計算得到了對應于1950~2011年的模式傾向誤差估計值。通過相關性分析和EOF分析發(fā)現(xiàn), 模式傾向誤差估計值和模擬誤差之間具有高度相關關系。進一步地將傾向誤差估計項加入ZC模式, 重復開展了多次ENSO模擬實驗。通過與ZC模式的對比, 發(fā)現(xiàn)在150 a的自由積分實驗中, 帶有模式傾向誤差的F-ZC模式均能進一步改善對ENSO事件的模擬效果, 在振幅、周期、相位鎖定和主要空間形態(tài)等方面, 模擬與觀測的相關性均有所提高。

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

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

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

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

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

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

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

    夏詠華, 2001. “厄爾尼諾”事件及其對氣候和航?;顒拥挠绊慬D]. 大連: 大連海事大學: 10-13.

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

    穆穆, 任宏利, 2017. 2014~2016年超強厄爾尼諾事件研究及其預測給予我們的啟示[J]. 中國科學: 地球科學, 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

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

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

    2022-01-12,

    2022-03-10

    猜你喜歡
    估計值模態(tài)觀測
    觀測到恒星死亡瞬間
    軍事文摘(2023年18期)2023-11-03 09:45:42
    一道樣本的數(shù)字特征與頻率分布直方圖的交匯問題
    統(tǒng)計信息
    2018年4月世界粗鋼產(chǎn)量表(續(xù))萬噸
    天測與測地VLBI 測地站周圍地形觀測遮掩的討論
    可觀測宇宙
    太空探索(2016年7期)2016-07-10 12:10:15
    國內(nèi)多模態(tài)教學研究回顧與展望
    高分辨率對地觀測系統(tǒng)
    太空探索(2015年8期)2015-07-18 11:04:44
    基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識別
    由單個模態(tài)構(gòu)造對稱簡支梁的抗彎剛度
    計算物理(2014年2期)2014-03-11 17:01:39
    免费在线观看完整版高清| 精品一区二区三区视频在线| av免费在线看不卡| 国国产精品蜜臀av免费| 成人午夜精彩视频在线观看| 全区人妻精品视频| 一二三四中文在线观看免费高清| 26uuu在线亚洲综合色| 在线天堂最新版资源| 人妻 亚洲 视频| 亚洲国产欧美日韩在线播放| 大香蕉97超碰在线| 91在线精品国自产拍蜜月| 亚洲av中文av极速乱| 中文字幕另类日韩欧美亚洲嫩草| 日本与韩国留学比较| 久久这里只有精品19| 九色成人免费人妻av| 狂野欧美激情性xxxx在线观看| 久久久久久久精品精品| 国内精品宾馆在线| 国产精品国产三级国产专区5o| 中文字幕亚洲精品专区| 国产精品人妻久久久久久| 街头女战士在线观看网站| 免费女性裸体啪啪无遮挡网站| 最新的欧美精品一区二区| 又大又黄又爽视频免费| 777米奇影视久久| 天堂8中文在线网| av电影中文网址| 亚洲精品成人av观看孕妇| 97精品久久久久久久久久精品| 久久久久久久国产电影| 国产精品人妻久久久久久| 国产高清不卡午夜福利| 超色免费av| 少妇人妻久久综合中文| 亚洲国产精品一区三区| a级毛片黄视频| 国产精品蜜桃在线观看| 国产深夜福利视频在线观看| 少妇的逼水好多| 国产精品.久久久| 午夜福利网站1000一区二区三区| 人人妻人人澡人人看| 国产成人免费观看mmmm| 天堂8中文在线网| 欧美 亚洲 国产 日韩一| 热re99久久精品国产66热6| 丝瓜视频免费看黄片| 午夜视频国产福利| 日本-黄色视频高清免费观看| 日本与韩国留学比较| 女人精品久久久久毛片| 日日爽夜夜爽网站| 秋霞在线观看毛片| 男人添女人高潮全过程视频| 黄色配什么色好看| 一区二区三区精品91| 亚洲成国产人片在线观看| 9色porny在线观看| 高清欧美精品videossex| 欧美 亚洲 国产 日韩一| 香蕉精品网在线| 亚洲欧美日韩卡通动漫| 国产精品久久久久久久电影| 五月伊人婷婷丁香| 中国国产av一级| 2021少妇久久久久久久久久久| 22中文网久久字幕| 欧美性感艳星| 三上悠亚av全集在线观看| 不卡视频在线观看欧美| 成年美女黄网站色视频大全免费| 99热国产这里只有精品6| 午夜激情久久久久久久| 国语对白做爰xxxⅹ性视频网站| 侵犯人妻中文字幕一二三四区| 国产乱人偷精品视频| 亚洲精品一二三| 亚洲成人av在线免费| 精品酒店卫生间| 国产免费福利视频在线观看| 男的添女的下面高潮视频| 高清黄色对白视频在线免费看| 五月伊人婷婷丁香| 亚洲欧美日韩卡通动漫| 国产黄色视频一区二区在线观看| 亚洲图色成人| 欧美日本中文国产一区发布| 丝袜人妻中文字幕| 午夜久久久在线观看| 毛片一级片免费看久久久久| 亚洲精品中文字幕在线视频| 国产av一区二区精品久久| 成人亚洲精品一区在线观看| 亚洲经典国产精华液单| 欧美性感艳星| 久久国产精品大桥未久av| 青春草亚洲视频在线观看| 亚洲国产精品一区二区三区在线| 777米奇影视久久| 欧美精品亚洲一区二区| 精品人妻熟女毛片av久久网站| av免费观看日本| 一二三四中文在线观看免费高清| 久久狼人影院| 制服诱惑二区| 2021少妇久久久久久久久久久| 精品一品国产午夜福利视频| 国产精品一国产av| 美国免费a级毛片| 日韩电影二区| 日本色播在线视频| 我要看黄色一级片免费的| 久久国产精品男人的天堂亚洲 | 在线观看三级黄色| 亚洲精品第二区| 99国产精品免费福利视频| 少妇被粗大猛烈的视频| 一级毛片我不卡| 国产又爽黄色视频| 中文欧美无线码| 少妇猛男粗大的猛烈进出视频| 久久精品国产鲁丝片午夜精品| 熟女电影av网| 少妇的逼好多水| 侵犯人妻中文字幕一二三四区| 不卡视频在线观看欧美| 亚洲精品色激情综合| 成年人午夜在线观看视频| 精品人妻在线不人妻| 免费人妻精品一区二区三区视频| 精品卡一卡二卡四卡免费| 国产高清三级在线| 国产精品一国产av| 国产精品嫩草影院av在线观看| 日韩欧美一区视频在线观看| 国产黄色免费在线视频| 亚洲婷婷狠狠爱综合网| 欧美变态另类bdsm刘玥| 亚洲欧洲国产日韩| 欧美日韩av久久| 国产毛片在线视频| 精品午夜福利在线看| 国产黄色视频一区二区在线观看| 久久女婷五月综合色啪小说| 国产片内射在线| 国产成人欧美| 亚洲av日韩在线播放| av一本久久久久| 晚上一个人看的免费电影| 精品国产露脸久久av麻豆| 夫妻午夜视频| 天天躁夜夜躁狠狠躁躁| 视频在线观看一区二区三区| 免费观看性生交大片5| 老女人水多毛片| 1024视频免费在线观看| 9191精品国产免费久久| 国产精品国产av在线观看| 亚洲天堂av无毛| videossex国产| 精品一区二区三卡| 精品国产乱码久久久久久小说| 午夜老司机福利剧场| 日韩精品免费视频一区二区三区 | 在现免费观看毛片| 国产毛片在线视频| 国产国拍精品亚洲av在线观看| 黄网站色视频无遮挡免费观看| 大香蕉久久网| 九色成人免费人妻av| 欧美亚洲日本最大视频资源| 国产成人精品无人区| 天堂8中文在线网| 亚洲图色成人| 街头女战士在线观看网站| 99九九在线精品视频| 在线观看免费视频网站a站| 亚洲成人手机| 日韩视频在线欧美| 国产极品粉嫩免费观看在线| 国产精品一区www在线观看| 亚洲,欧美,日韩| 亚洲伊人久久精品综合| 激情视频va一区二区三区| 人体艺术视频欧美日本| 人人妻人人添人人爽欧美一区卜| 宅男免费午夜| 日韩大片免费观看网站| 久久久久久人人人人人| 日韩av在线免费看完整版不卡| 免费观看性生交大片5| 韩国高清视频一区二区三区| 国产成人精品久久久久久| 搡老乐熟女国产| 久久人人97超碰香蕉20202| 精品福利永久在线观看| 好男人视频免费观看在线| 一区二区日韩欧美中文字幕 | 亚洲精品第二区| 看十八女毛片水多多多| 少妇的逼好多水| 女的被弄到高潮叫床怎么办| 伊人亚洲综合成人网| 亚洲精品av麻豆狂野| 免费少妇av软件| 亚洲美女黄色视频免费看| videossex国产| 蜜桃在线观看..| 日本wwww免费看| 欧美国产精品va在线观看不卡| 免费看光身美女| 日韩三级伦理在线观看| 18禁观看日本| 亚洲美女黄色视频免费看| 我要看黄色一级片免费的| 99re6热这里在线精品视频| 日韩三级伦理在线观看| 高清毛片免费看| 国产av精品麻豆| 欧美日韩成人在线一区二区| 亚洲精品,欧美精品| 久久精品aⅴ一区二区三区四区 | 国产精品久久久久久久久免| 精品人妻偷拍中文字幕| 久久综合国产亚洲精品| 两个人看的免费小视频| 日韩一本色道免费dvd| 日本与韩国留学比较| 午夜久久久在线观看| 欧美人与性动交α欧美软件 | 寂寞人妻少妇视频99o| 少妇人妻久久综合中文| 卡戴珊不雅视频在线播放| 亚洲精品中文字幕在线视频| 成人综合一区亚洲| 亚洲经典国产精华液单| 少妇人妻 视频| 国产成人免费无遮挡视频| 午夜精品国产一区二区电影| 午夜免费男女啪啪视频观看| av片东京热男人的天堂| 中文字幕制服av| 母亲3免费完整高清在线观看 | 亚洲欧美中文字幕日韩二区| 中文字幕人妻丝袜制服| 啦啦啦中文免费视频观看日本| 久久这里有精品视频免费| 日本vs欧美在线观看视频| 妹子高潮喷水视频| 日韩av在线免费看完整版不卡| 日韩欧美精品免费久久| 精品酒店卫生间| 中文字幕最新亚洲高清| 精品视频人人做人人爽| 亚洲综合色惰| 国产爽快片一区二区三区| 91午夜精品亚洲一区二区三区| 精品国产国语对白av| 在线精品无人区一区二区三| 国产视频首页在线观看| 精品人妻在线不人妻| 我的女老师完整版在线观看| 男人舔女人的私密视频| av福利片在线| 欧美成人精品欧美一级黄| 视频在线观看一区二区三区| 水蜜桃什么品种好| 人成视频在线观看免费观看| 春色校园在线视频观看| 欧美老熟妇乱子伦牲交| √禁漫天堂资源中文www| 久久精品国产亚洲av天美| 欧美日韩成人在线一区二区| 久久人人爽人人片av| 免费观看在线日韩| 99九九在线精品视频| 免费看光身美女| 国产亚洲最大av| 少妇精品久久久久久久| 中国三级夫妇交换| 内地一区二区视频在线| 哪个播放器可以免费观看大片| 人妻少妇偷人精品九色| 国产女主播在线喷水免费视频网站| av卡一久久| 午夜福利,免费看| 国产亚洲午夜精品一区二区久久| 国产亚洲av片在线观看秒播厂| 26uuu在线亚洲综合色| 国产69精品久久久久777片| 午夜福利网站1000一区二区三区| 国产精品嫩草影院av在线观看| 亚洲av电影在线进入| av在线播放精品| 亚洲国产精品999| 日韩av免费高清视频| 亚洲少妇的诱惑av| 色婷婷av一区二区三区视频| 亚洲精品一区蜜桃| 搡女人真爽免费视频火全软件| 国产一级毛片在线| 成人手机av| 亚洲五月色婷婷综合| 午夜福利网站1000一区二区三区| 精品一区二区三区视频在线| 午夜福利网站1000一区二区三区| 精品一区二区三区视频在线| 交换朋友夫妻互换小说| 欧美日韩av久久| 亚洲欧美一区二区三区黑人 | 精品亚洲成国产av| 亚洲伊人久久精品综合| 精品人妻一区二区三区麻豆| 国产精品秋霞免费鲁丝片| 久久久久久人人人人人| 久久精品久久久久久噜噜老黄| 亚洲av综合色区一区| 免费高清在线观看视频在线观看| 99热6这里只有精品| 纵有疾风起免费观看全集完整版| a 毛片基地| 免费观看av网站的网址| 婷婷色综合www| 成人亚洲欧美一区二区av| 黄色毛片三级朝国网站| 成人亚洲精品一区在线观看| 18禁在线无遮挡免费观看视频| 精品酒店卫生间| 久久青草综合色| 亚洲成av片中文字幕在线观看 | 成人无遮挡网站| 欧美激情极品国产一区二区三区 | 久久人妻熟女aⅴ| av卡一久久| 亚洲中文av在线| 亚洲色图综合在线观看| 在线天堂中文资源库| 亚洲精品久久午夜乱码| 国产在线免费精品| 嫩草影院入口| 国产av码专区亚洲av| 嫩草影院入口| 精品国产露脸久久av麻豆| 亚洲一区二区三区欧美精品| 人妻系列 视频| 纵有疾风起免费观看全集完整版| 人妻少妇偷人精品九色| 亚洲av中文av极速乱| 最近的中文字幕免费完整| 欧美精品一区二区大全| 国产 精品1| av网站免费在线观看视频| www日本在线高清视频| 国产av码专区亚洲av| 精品久久久精品久久久| 亚洲欧美一区二区三区黑人 | 精品亚洲成a人片在线观看| www.色视频.com| 自线自在国产av| 免费在线观看黄色视频的| 纯流量卡能插随身wifi吗| 精品国产乱码久久久久久小说| 欧美日韩亚洲高清精品| 久久这里只有精品19| 国产在线一区二区三区精| 人人妻人人澡人人看| 免费黄网站久久成人精品| 五月伊人婷婷丁香| 亚洲av在线观看美女高潮| 亚洲国产av新网站| av在线老鸭窝| 爱豆传媒免费全集在线观看| 亚洲丝袜综合中文字幕| 激情五月婷婷亚洲| 欧美+日韩+精品| 久久人人爽av亚洲精品天堂| 久久精品国产亚洲av涩爱| 最近中文字幕高清免费大全6| 久久久久国产网址| 久久韩国三级中文字幕| 香蕉丝袜av| 亚洲欧美清纯卡通| 交换朋友夫妻互换小说| 国产日韩欧美视频二区| 国产一区二区三区av在线| 欧美3d第一页| 人人妻人人澡人人看| 天天躁夜夜躁狠狠躁躁| 亚洲精品日韩在线中文字幕| 亚洲高清免费不卡视频| av又黄又爽大尺度在线免费看| 久久精品久久精品一区二区三区| 另类精品久久| 亚洲经典国产精华液单| 18禁在线无遮挡免费观看视频| 成人亚洲欧美一区二区av| 男女国产视频网站| 春色校园在线视频观看| 国产综合精华液| 国产亚洲一区二区精品| 午夜福利网站1000一区二区三区| 少妇 在线观看| 久久99蜜桃精品久久| 丝袜在线中文字幕| 少妇人妻久久综合中文| 亚洲国产看品久久| www.熟女人妻精品国产 | 各种免费的搞黄视频| 日韩视频在线欧美| 中文字幕最新亚洲高清| 国产欧美亚洲国产| 国产在视频线精品| 中文字幕另类日韩欧美亚洲嫩草| 久久ye,这里只有精品| av女优亚洲男人天堂| 高清视频免费观看一区二区| 日韩中字成人| 一边摸一边做爽爽视频免费| 国产精品久久久久久久电影| 久久青草综合色| 热99久久久久精品小说推荐| 日韩中文字幕视频在线看片| 亚洲图色成人| 午夜av观看不卡| 熟女人妻精品中文字幕| av.在线天堂| 丝袜美足系列| 免费看光身美女| 亚洲欧洲精品一区二区精品久久久 | 秋霞伦理黄片| 青春草亚洲视频在线观看| 国产在线视频一区二区| 美女内射精品一级片tv| 精品熟女少妇av免费看| 97人妻天天添夜夜摸| 一区二区三区精品91| 午夜精品国产一区二区电影| 少妇的丰满在线观看| 久久精品aⅴ一区二区三区四区 | 高清欧美精品videossex| 美女大奶头黄色视频| 99热全是精品| 成人国产麻豆网| 国产一区有黄有色的免费视频| 多毛熟女@视频| 亚洲成av片中文字幕在线观看 | 黄色视频在线播放观看不卡| 最新的欧美精品一区二区| 国产精品蜜桃在线观看| 亚洲精品乱久久久久久| 蜜桃在线观看..| 日日爽夜夜爽网站| 亚洲欧美一区二区三区黑人 | 高清视频免费观看一区二区| 久久毛片免费看一区二区三区| 一边亲一边摸免费视频| 国产午夜精品一二区理论片| 高清不卡的av网站| 日本免费在线观看一区| 51国产日韩欧美| 制服诱惑二区| 精品亚洲成国产av| 精品一区在线观看国产| 女性生殖器流出的白浆| 日本欧美国产在线视频| 亚洲精品久久久久久婷婷小说| 久久久久久人人人人人| 日韩中文字幕视频在线看片| 熟女av电影| 亚洲精品av麻豆狂野| 99国产综合亚洲精品| 亚洲国产精品成人久久小说| 亚洲天堂av无毛| 91成人精品电影| 国产男女内射视频| 亚洲精品美女久久av网站| 亚洲精品国产av成人精品| 91午夜精品亚洲一区二区三区| 如何舔出高潮| 大码成人一级视频| 亚洲成av片中文字幕在线观看 | 最近2019中文字幕mv第一页| 最黄视频免费看| 精品亚洲成国产av| 中国国产av一级| 国产精品一区二区在线观看99| 亚洲国产精品成人久久小说| 男人操女人黄网站| 在线观看人妻少妇| 成人漫画全彩无遮挡| 国产精品一区二区在线观看99| 曰老女人黄片| 最黄视频免费看| 国产又爽黄色视频| 天美传媒精品一区二区| 久久人妻熟女aⅴ| 国产精品一二三区在线看| 精品久久久久久电影网| 另类亚洲欧美激情| 好男人视频免费观看在线| 男女下面插进去视频免费观看 | 日韩三级伦理在线观看| 人妻一区二区av| 美国免费a级毛片| 亚洲精品久久久久久婷婷小说| 中文字幕免费在线视频6| 大片电影免费在线观看免费| 日韩,欧美,国产一区二区三区| 亚洲四区av| 91aial.com中文字幕在线观看| 高清在线视频一区二区三区| 男女高潮啪啪啪动态图| 欧美最新免费一区二区三区| 国产成人aa在线观看| 美女主播在线视频| 精品视频人人做人人爽| 一二三四中文在线观看免费高清| 中文天堂在线官网| 精品99又大又爽又粗少妇毛片| 少妇的丰满在线观看| 人人妻人人爽人人添夜夜欢视频| 欧美人与性动交α欧美软件 | 91成人精品电影| 久久午夜综合久久蜜桃| 免费高清在线观看日韩| 午夜久久久在线观看| 国产在线视频一区二区| 精品一区二区三区视频在线| av福利片在线| 日韩一区二区视频免费看| 国产成人av激情在线播放| 国产有黄有色有爽视频| freevideosex欧美| 韩国精品一区二区三区 | a级片在线免费高清观看视频| 亚洲欧美日韩卡通动漫| 国产永久视频网站| av线在线观看网站| 一级毛片电影观看| 亚洲成人av在线免费| 人妻系列 视频| 成人国产av品久久久| 久久久久久人人人人人| 九色成人免费人妻av| 少妇高潮的动态图| 一二三四在线观看免费中文在 | 欧美成人精品欧美一级黄| 18禁裸乳无遮挡动漫免费视频| 人人妻人人爽人人添夜夜欢视频| 99久国产av精品国产电影| 日本色播在线视频| 成年动漫av网址| 国产精品久久久久久av不卡| 18在线观看网站| 十八禁高潮呻吟视频| 男女国产视频网站| 九九在线视频观看精品| 九色亚洲精品在线播放| 亚洲av中文av极速乱| 国产精品久久久久久久久免| 久久这里只有精品19| 青春草亚洲视频在线观看| 午夜影院在线不卡| 有码 亚洲区| 久久久欧美国产精品| 免费大片18禁| 岛国毛片在线播放| 99久久综合免费| 天堂中文最新版在线下载| 欧美精品亚洲一区二区| 久久99蜜桃精品久久| 中文字幕亚洲精品专区| 国产成人91sexporn| 国产色婷婷99| 激情五月婷婷亚洲| 日韩成人伦理影院| 91在线精品国自产拍蜜月| 国产精品.久久久| 久久精品国产综合久久久 | 搡老乐熟女国产| 18+在线观看网站| 日日爽夜夜爽网站| 肉色欧美久久久久久久蜜桃| 国产精品99久久99久久久不卡 | 国产免费又黄又爽又色| 激情五月婷婷亚洲| 捣出白浆h1v1| videossex国产| 国产极品天堂在线| 日本欧美视频一区| 黄色毛片三级朝国网站| 日产精品乱码卡一卡2卡三| 亚洲第一av免费看| 嫩草影院入口| 日本av免费视频播放| 久久狼人影院| 丰满少妇做爰视频| 观看av在线不卡| 久热这里只有精品99| 欧美日韩国产mv在线观看视频| 国产一区二区三区av在线| 一区二区日韩欧美中文字幕 | 日韩熟女老妇一区二区性免费视频| 韩国av在线不卡| 又粗又硬又长又爽又黄的视频| 伦理电影大哥的女人| 国产亚洲一区二区精品| 日本猛色少妇xxxxx猛交久久| 免费av不卡在线播放| 999精品在线视频| √禁漫天堂资源中文www|