• <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视频在线免费观看| 久久热精品热| 毛片一级片免费看久久久久| 丰满乱子伦码专区| av在线app专区| 亚洲无线观看免费| 免费高清在线观看视频在线观看| 日本猛色少妇xxxxx猛交久久| 亚洲婷婷狠狠爱综合网| 欧美高清性xxxxhd video| 赤兔流量卡办理| 国产男女超爽视频在线观看| 国产伦精品一区二区三区视频9| 国产亚洲最大av| 亚洲av免费高清在线观看| 美女视频免费永久观看网站| 97超碰精品成人国产| 全区人妻精品视频| 日韩精品有码人妻一区| 一级二级三级毛片免费看| 国产在视频线精品| 中文字幕制服av| 午夜免费男女啪啪视频观看| 免费观看在线日韩| 欧美极品一区二区三区四区| 国产精品不卡视频一区二区| 亚洲一区二区三区欧美精品| 十分钟在线观看高清视频www | 久久99精品国语久久久| 国产精品免费大片| 日韩亚洲欧美综合| 99热这里只有是精品50| 国产有黄有色有爽视频| 99久久精品热视频| 国产精品国产三级国产av玫瑰| 国产精品精品国产色婷婷| av在线播放精品| 色5月婷婷丁香| 91在线精品国自产拍蜜月| 久久久久精品久久久久真实原创| 日日啪夜夜撸| 成人毛片60女人毛片免费| 老司机影院毛片| 伊人久久精品亚洲午夜| 精品少妇黑人巨大在线播放| 美女中出高潮动态图| 纯流量卡能插随身wifi吗| 精品人妻熟女av久视频| 一本久久精品| 国产成人免费无遮挡视频| 啦啦啦中文免费视频观看日本| 天堂8中文在线网| 国产视频首页在线观看| 男女国产视频网站| 一边亲一边摸免费视频| 国产精品一及| 3wmmmm亚洲av在线观看| 国产真实伦视频高清在线观看| 欧美日韩综合久久久久久| 直男gayav资源| 精品人妻熟女av久视频| 老司机影院成人| 午夜福利视频精品| av在线播放精品| 18禁在线无遮挡免费观看视频| 嫩草影院入口| 夜夜骑夜夜射夜夜干| 女性被躁到高潮视频| 少妇熟女欧美另类| 麻豆精品久久久久久蜜桃| 国产精品99久久99久久久不卡 | 国产成人精品久久久久久| 人妻制服诱惑在线中文字幕| 色5月婷婷丁香| 一个人看视频在线观看www免费| 日韩,欧美,国产一区二区三区| 国产在线免费精品| 黄色日韩在线| 日本免费在线观看一区| 天堂中文最新版在线下载| 亚洲激情五月婷婷啪啪| 久久国产乱子免费精品| 一本—道久久a久久精品蜜桃钙片| 色哟哟·www| 99九九线精品视频在线观看视频| 91精品国产国语对白视频| 久久6这里有精品| 女人十人毛片免费观看3o分钟| 日韩三级伦理在线观看| 少妇被粗大猛烈的视频| 亚洲精品日韩av片在线观看| 久久女婷五月综合色啪小说| 日韩,欧美,国产一区二区三区| 99视频精品全部免费 在线| 黄色一级大片看看| av福利片在线观看| 男女免费视频国产| 免费高清在线观看视频在线观看| 日本vs欧美在线观看视频 | 日韩av不卡免费在线播放| 久久人妻熟女aⅴ| 国产精品99久久久久久久久| 老女人水多毛片| 99久国产av精品国产电影| 久久 成人 亚洲| 亚洲成人手机| 亚洲av免费高清在线观看| 国产视频内射| 日日啪夜夜爽| 又黄又爽又刺激的免费视频.| 五月玫瑰六月丁香| 亚洲美女黄色视频免费看| 亚洲精品久久久久久婷婷小说| 亚洲av中文av极速乱| 中文字幕av成人在线电影| 欧美精品一区二区大全| 中文欧美无线码| 一本色道久久久久久精品综合| 中国国产av一级| 亚洲高清免费不卡视频| 美女内射精品一级片tv| 一区在线观看完整版| 午夜福利在线观看免费完整高清在| 国产免费福利视频在线观看| 日韩av不卡免费在线播放| 婷婷色麻豆天堂久久| av国产免费在线观看| 91久久精品电影网| 久久国产精品男人的天堂亚洲 | 美女内射精品一级片tv| 大片免费播放器 马上看| 亚洲国产色片| 亚洲av在线观看美女高潮| 直男gayav资源| 国产高清三级在线| 国产 一区 欧美 日韩| 精品人妻一区二区三区麻豆| 赤兔流量卡办理| kizo精华| 亚洲四区av| 最近2019中文字幕mv第一页| 亚洲不卡免费看| 国产免费又黄又爽又色| 十八禁网站网址无遮挡 | 老师上课跳d突然被开到最大视频| av黄色大香蕉| 国产毛片在线视频| 欧美区成人在线视频| 成人黄色视频免费在线看| 成人特级av手机在线观看| 韩国av在线不卡| 亚洲性久久影院| 少妇丰满av| 一区二区三区免费毛片| 亚洲av中文字字幕乱码综合| 街头女战士在线观看网站| 免费观看无遮挡的男女| 亚洲人成网站高清观看| 亚洲精品,欧美精品| 久久6这里有精品| 精品人妻一区二区三区麻豆| 国产永久视频网站| 超碰av人人做人人爽久久| 中文字幕免费在线视频6| 人人妻人人澡人人爽人人夜夜| 亚洲综合色惰| 国产黄色视频一区二区在线观看| 伦理电影大哥的女人| 亚洲av.av天堂| av女优亚洲男人天堂| 日韩在线高清观看一区二区三区| 亚洲综合色惰| 狂野欧美激情性bbbbbb| 国产成人精品福利久久| 一级片'在线观看视频| 久久久久性生活片| 一个人看的www免费观看视频| 七月丁香在线播放| av在线播放精品| 99久久综合免费| 男女下面进入的视频免费午夜| 中文精品一卡2卡3卡4更新| 大片免费播放器 马上看| 久热久热在线精品观看| 国产久久久一区二区三区| 人人妻人人看人人澡| 亚洲av欧美aⅴ国产| 亚洲av成人精品一二三区| 在线观看美女被高潮喷水网站| 久久久久精品性色| 国产69精品久久久久777片| 久久久久久久久久久丰满| 91精品一卡2卡3卡4卡| 精品熟女少妇av免费看| 亚洲精品日韩av片在线观看| av一本久久久久| 日韩亚洲欧美综合| 九草在线视频观看| 欧美少妇被猛烈插入视频| 91精品国产九色| 欧美日韩综合久久久久久| 免费久久久久久久精品成人欧美视频 | 国产精品一区二区三区四区免费观看| 一级毛片aaaaaa免费看小| 人妻一区二区av| 色网站视频免费| 国产 一区精品| 伦理电影免费视频| 简卡轻食公司| 免费黄频网站在线观看国产| 深夜a级毛片| av国产久精品久网站免费入址| 免费人成在线观看视频色| 国产真实伦视频高清在线观看| videossex国产| 亚洲av中文字字幕乱码综合| 久久久久久伊人网av| 97精品久久久久久久久久精品| 岛国毛片在线播放| 国产精品蜜桃在线观看| 成人毛片60女人毛片免费| 黄色怎么调成土黄色| 久久久a久久爽久久v久久| 蜜臀久久99精品久久宅男| 精品久久久久久久末码| 欧美老熟妇乱子伦牲交| 欧美精品国产亚洲| 国产精品蜜桃在线观看| 午夜日本视频在线| 久久国内精品自在自线图片| 国产精品人妻久久久久久| 日本wwww免费看| 在线免费十八禁| 日韩,欧美,国产一区二区三区| 丰满人妻一区二区三区视频av| 国产中年淑女户外野战色| 黄色日韩在线| 校园人妻丝袜中文字幕| 日韩欧美一区视频在线观看 | 精品一区在线观看国产| 丝瓜视频免费看黄片| 国产乱人视频| av不卡在线播放| 尾随美女入室| 色吧在线观看| 成年av动漫网址| 成人毛片a级毛片在线播放| 校园人妻丝袜中文字幕| 国产在线男女| 国产精品国产三级国产av玫瑰| 国产亚洲一区二区精品| 永久免费av网站大全| 亚洲精品国产成人久久av| 一本久久精品| 国产精品秋霞免费鲁丝片| 男女边吃奶边做爰视频| 成人无遮挡网站| 久热久热在线精品观看| 亚洲国产精品999| av网站免费在线观看视频| 国产亚洲91精品色在线| 久久久色成人| 欧美一级a爱片免费观看看| 国产免费福利视频在线观看| 日韩三级伦理在线观看| 在线观看免费高清a一片| 久久这里有精品视频免费| 在线观看一区二区三区| 国产69精品久久久久777片| 永久网站在线| 久久午夜福利片| 亚洲精华国产精华液的使用体验| 久久久国产一区二区| 女的被弄到高潮叫床怎么办| 精品人妻偷拍中文字幕| 天天躁日日操中文字幕| 亚洲,欧美,日韩| 国产精品人妻久久久影院| 男人狂女人下面高潮的视频| 午夜福利在线在线| 麻豆成人午夜福利视频| 免费观看av网站的网址| 久久午夜福利片| 亚洲无线观看免费| 91久久精品国产一区二区三区| 精品亚洲成国产av| 精品久久久久久久久亚洲| 黄色视频在线播放观看不卡| 欧美三级亚洲精品| 成人无遮挡网站| 十八禁网站网址无遮挡 | 在线观看美女被高潮喷水网站| 男人添女人高潮全过程视频| 十分钟在线观看高清视频www | 亚洲精品一区蜜桃| 少妇被粗大猛烈的视频| 亚洲精品456在线播放app| 国产爱豆传媒在线观看| 最近手机中文字幕大全| 国产精品一区二区在线观看99| 美女xxoo啪啪120秒动态图| 中文字幕久久专区| 国产精品国产三级国产av玫瑰| 香蕉精品网在线| 国产永久视频网站| 夜夜骑夜夜射夜夜干| 黄色一级大片看看| 久久青草综合色| 午夜免费男女啪啪视频观看| 精品午夜福利在线看| 亚洲精品国产av蜜桃| 国产黄色免费在线视频| 最近2019中文字幕mv第一页| 成人毛片a级毛片在线播放| 日韩欧美精品免费久久| 夫妻性生交免费视频一级片| 高清不卡的av网站| 尾随美女入室| www.色视频.com| av卡一久久| 亚洲国产欧美在线一区| 国产精品久久久久久av不卡| 中文字幕亚洲精品专区| 久久女婷五月综合色啪小说| 日韩三级伦理在线观看| av卡一久久| 亚洲精品一区蜜桃| 国产高清不卡午夜福利| 美女国产视频在线观看| 亚洲综合精品二区| 少妇精品久久久久久久| 成人毛片60女人毛片免费| 欧美少妇被猛烈插入视频| 国产无遮挡羞羞视频在线观看| 黄色一级大片看看| 久久久久久久久大av| 国产 精品1| 不卡视频在线观看欧美| 亚洲综合色惰| 乱系列少妇在线播放| 一个人看视频在线观看www免费| 一区二区av电影网| 春色校园在线视频观看| 久久久久久人妻| 麻豆成人午夜福利视频| 欧美一级a爱片免费观看看| 一级毛片黄色毛片免费观看视频| 啦啦啦在线观看免费高清www| a 毛片基地| 91狼人影院| 精品国产三级普通话版| 精品亚洲成国产av| 噜噜噜噜噜久久久久久91| 亚洲国产最新在线播放| 欧美 日韩 精品 国产| 一个人看的www免费观看视频| freevideosex欧美| 国产精品女同一区二区软件| 大香蕉97超碰在线| 成年免费大片在线观看| 青春草视频在线免费观看| 寂寞人妻少妇视频99o| 亚洲欧洲日产国产| 亚洲婷婷狠狠爱综合网| 高清毛片免费看| 国产有黄有色有爽视频| 亚洲国产欧美人成| 纯流量卡能插随身wifi吗| 亚洲欧美一区二区三区黑人 | 直男gayav资源| av视频免费观看在线观看| 亚洲精品乱码久久久久久按摩| 狂野欧美激情性xxxx在线观看| 国产一级毛片在线| 亚洲在久久综合| 日日啪夜夜撸| 日韩强制内射视频| 国产 精品1| 在线播放无遮挡| 一区在线观看完整版| 亚洲欧美清纯卡通| 国产成人午夜福利电影在线观看| 成人国产麻豆网| 国产中年淑女户外野战色| 免费av不卡在线播放| 亚洲欧洲日产国产| 午夜精品国产一区二区电影| 日韩,欧美,国产一区二区三区| 成人毛片a级毛片在线播放| 各种免费的搞黄视频| 男女国产视频网站| 亚洲欧美一区二区三区国产| 国产在线男女| 欧美亚洲 丝袜 人妻 在线| 久久久久人妻精品一区果冻| 国产亚洲91精品色在线| 男人和女人高潮做爰伦理| 亚洲国产精品999| 搡女人真爽免费视频火全软件| 国产成人免费无遮挡视频| 午夜精品国产一区二区电影| av.在线天堂| 国产精品av视频在线免费观看| 成人漫画全彩无遮挡| 国产精品久久久久久精品古装| 欧美日本视频| 欧美日韩在线观看h| 观看美女的网站| 国产91av在线免费观看| 欧美日韩综合久久久久久| 国产精品三级大全| 欧美成人精品欧美一级黄| 日本黄大片高清| 久久 成人 亚洲| 精品亚洲成国产av| 在线观看免费日韩欧美大片 | 国国产精品蜜臀av免费| 久久久久久人妻| 国产免费福利视频在线观看| 亚洲精品国产成人久久av| 嘟嘟电影网在线观看| 国产av精品麻豆| 久久久久视频综合| 啦啦啦视频在线资源免费观看| 少妇熟女欧美另类| 国产伦理片在线播放av一区| 久久久久精品久久久久真实原创| 国产黄色视频一区二区在线观看| 亚洲在久久综合| 黑人猛操日本美女一级片| 亚洲人与动物交配视频| 91狼人影院| 五月伊人婷婷丁香| 欧美国产精品一级二级三级 | 偷拍熟女少妇极品色| 啦啦啦啦在线视频资源| 国产熟女欧美一区二区| 最近中文字幕高清免费大全6| 精品一品国产午夜福利视频| 免费播放大片免费观看视频在线观看| 免费黄色在线免费观看| 国产av精品麻豆| 舔av片在线| 精品人妻一区二区三区麻豆| 国产高潮美女av| av一本久久久久| 国产一级毛片在线| 女人久久www免费人成看片| 日韩中文字幕视频在线看片 | 国产欧美日韩精品一区二区| 熟女电影av网| 一级毛片 在线播放| 99热6这里只有精品| 伦精品一区二区三区| 少妇的逼好多水| 最近最新中文字幕大全电影3| 免费播放大片免费观看视频在线观看| 亚洲精品第二区| 男女免费视频国产| 欧美成人一区二区免费高清观看| 伊人久久精品亚洲午夜| 亚洲精品乱久久久久久| 免费黄网站久久成人精品| av国产免费在线观看| 91午夜精品亚洲一区二区三区| 欧美xxxx黑人xx丫x性爽| 久久鲁丝午夜福利片| 久久人人爽av亚洲精品天堂 | 精品一区二区免费观看| 日本黄色片子视频| 亚洲精品自拍成人| 夜夜骑夜夜射夜夜干| 亚洲最大成人中文| 国产精品女同一区二区软件| 精品亚洲乱码少妇综合久久| 久久久久久久亚洲中文字幕| 欧美精品一区二区免费开放| 少妇人妻精品综合一区二区| 国产免费一级a男人的天堂| 天天躁日日操中文字幕| 国产成人freesex在线| 中国国产av一级| 国产色婷婷99| 草草在线视频免费看| 免费在线观看成人毛片| 人妻 亚洲 视频| 最近中文字幕2019免费版| 天堂8中文在线网| 国产精品麻豆人妻色哟哟久久| 91aial.com中文字幕在线观看| 国产淫片久久久久久久久| 女性生殖器流出的白浆| 亚洲av电影在线观看一区二区三区| 丝袜喷水一区| 亚洲成人中文字幕在线播放| 狂野欧美激情性bbbbbb| 99精国产麻豆久久婷婷| 亚洲aⅴ乱码一区二区在线播放| 国产免费又黄又爽又色| 久久99热这里只有精品18| 国产69精品久久久久777片| 国产爱豆传媒在线观看| 春色校园在线视频观看| 韩国高清视频一区二区三区| 久久精品国产a三级三级三级| 国产精品无大码| 99热6这里只有精品| 久久人妻熟女aⅴ| 精品久久久精品久久久| 久久女婷五月综合色啪小说| 日本爱情动作片www.在线观看| 久久久欧美国产精品| 国产一区亚洲一区在线观看| 国产一区亚洲一区在线观看| 少妇 在线观看| 日韩欧美 国产精品| 日本一二三区视频观看| 久久久久久久久久人人人人人人| 国产精品精品国产色婷婷| 精品久久久噜噜| 亚洲精品色激情综合| 97在线人人人人妻| 这个男人来自地球电影免费观看 | 久久精品国产自在天天线| 简卡轻食公司| 人妻夜夜爽99麻豆av| 97热精品久久久久久| .国产精品久久| 欧美成人一区二区免费高清观看| av天堂中文字幕网| 久久国产亚洲av麻豆专区| 国产欧美亚洲国产| 啦啦啦视频在线资源免费观看| 亚洲国产日韩一区二区| 亚洲精品久久午夜乱码| 亚洲,欧美,日韩| 国产精品人妻久久久影院| 亚洲av在线观看美女高潮| 国产黄频视频在线观看| 黄色视频在线播放观看不卡| 国产免费福利视频在线观看| 高清不卡的av网站| 秋霞伦理黄片| 日本vs欧美在线观看视频 | 丰满少妇做爰视频| 性高湖久久久久久久久免费观看| 黄色怎么调成土黄色| 黄色一级大片看看| 最近中文字幕高清免费大全6| 99国产精品免费福利视频| 日本爱情动作片www.在线观看| 久久鲁丝午夜福利片| 卡戴珊不雅视频在线播放| 亚洲av二区三区四区| 大片免费播放器 马上看| 欧美区成人在线视频| 久久国内精品自在自线图片| 国产免费一级a男人的天堂| 久久久久国产精品人妻一区二区| 永久网站在线| 人妻少妇偷人精品九色| 男女啪啪激烈高潮av片| 18禁裸乳无遮挡免费网站照片| 高清午夜精品一区二区三区| 欧美成人a在线观看| 国产探花极品一区二区| 国产精品秋霞免费鲁丝片| 99热6这里只有精品| av免费观看日本| 中文在线观看免费www的网站| 久久久国产一区二区| 亚洲性久久影院| 精品久久久久久久末码| 永久免费av网站大全| 天堂俺去俺来也www色官网| 国产一区二区在线观看日韩| 免费黄频网站在线观看国产| 久久久精品94久久精品| 国产探花极品一区二区| 18禁裸乳无遮挡免费网站照片| 精品人妻偷拍中文字幕| 美女视频免费永久观看网站| 亚洲综合色惰| 成人午夜精彩视频在线观看| av国产精品久久久久影院| 一级毛片黄色毛片免费观看视频| 99久国产av精品国产电影| 欧美精品一区二区大全| 国内精品宾馆在线| 日韩欧美 国产精品| 国产精品熟女久久久久浪| 永久网站在线| 女性被躁到高潮视频| 国产乱人视频| 久久久久性生活片| 成人黄色视频免费在线看| 高清毛片免费看| 日韩中文字幕视频在线看片 | 一级毛片 在线播放| 一级毛片aaaaaa免费看小| 偷拍熟女少妇极品色| 网址你懂的国产日韩在线| 国产精品久久久久成人av| 22中文网久久字幕| 久久精品久久久久久久性| 菩萨蛮人人尽说江南好唐韦庄| 久久国产精品男人的天堂亚洲 | 国产视频首页在线观看|