• <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
    久久中文字幕一级| 在线视频色国产色| 亚洲天堂国产精品一区在线| 成年人黄色毛片网站| 中文字幕久久专区| 日韩大码丰满熟妇| 长腿黑丝高跟| 18禁黄网站禁片免费观看直播| 国产亚洲精品一区二区www| 大型av网站在线播放| ponron亚洲| 午夜激情福利司机影院| 亚洲欧美日韩无卡精品| 在线播放国产精品三级| 一本久久中文字幕| 久久久国产精品麻豆| 啦啦啦免费观看视频1| 亚洲专区国产一区二区| 国产精品美女特级片免费视频播放器 | 黄色 视频免费看| 成熟少妇高潮喷水视频| 免费在线观看黄色视频的| 女生性感内裤真人,穿戴方法视频| 亚洲国产日韩欧美精品在线观看 | 一区二区三区国产精品乱码| 特大巨黑吊av在线直播| 最近在线观看免费完整版| 欧美在线黄色| 亚洲午夜理论影院| 亚洲欧洲精品一区二区精品久久久| 男插女下体视频免费在线播放| 免费av毛片视频| 国产午夜福利久久久久久| 最近最新中文字幕大全电影3| 国产午夜精品久久久久久| 日韩免费av在线播放| 久热爱精品视频在线9| 久久 成人 亚洲| 黑人欧美特级aaaaaa片| 又大又爽又粗| 黄色视频,在线免费观看| 人人妻,人人澡人人爽秒播| 麻豆一二三区av精品| 亚洲美女视频黄频| 国产黄片美女视频| 久久精品亚洲精品国产色婷小说| 天堂av国产一区二区熟女人妻 | 在线播放国产精品三级| 母亲3免费完整高清在线观看| 香蕉av资源在线| 亚洲全国av大片| 成人国产综合亚洲| 婷婷丁香在线五月| 国产成人影院久久av| 99re在线观看精品视频| 久久 成人 亚洲| 两个人看的免费小视频| 两人在一起打扑克的视频| 18禁裸乳无遮挡免费网站照片| 999久久久精品免费观看国产| 两个人的视频大全免费| 国产av麻豆久久久久久久| 午夜福利成人在线免费观看| 国产麻豆成人av免费视频| 狠狠狠狠99中文字幕| aaaaa片日本免费| 黄色片一级片一级黄色片| 99国产精品一区二区蜜桃av| av有码第一页| 欧美黑人欧美精品刺激| 久久中文看片网| 亚洲全国av大片| 999久久久精品免费观看国产| 国内精品一区二区在线观看| 亚洲成a人片在线一区二区| 免费观看人在逋| 叶爱在线成人免费视频播放| 亚洲精品美女久久久久99蜜臀| 国产主播在线观看一区二区| 国产免费av片在线观看野外av| 巨乳人妻的诱惑在线观看| 男女午夜视频在线观看| 午夜视频精品福利| 日日夜夜操网爽| 国产麻豆成人av免费视频| 久久久久久久久免费视频了| 欧美中文日本在线观看视频| 亚洲成a人片在线一区二区| 国产av又大| 亚洲七黄色美女视频| 一级片免费观看大全| 亚洲国产精品成人综合色| 国产精品久久久人人做人人爽| 在线观看舔阴道视频| 非洲黑人性xxxx精品又粗又长| 制服诱惑二区| 欧美在线一区亚洲| 两性夫妻黄色片| av视频在线观看入口| 波多野结衣高清作品| 在线永久观看黄色视频| av天堂在线播放| 久久婷婷成人综合色麻豆| 亚洲专区中文字幕在线| 日本黄色视频三级网站网址| 国产视频一区二区在线看| 色哟哟哟哟哟哟| 久热爱精品视频在线9| 午夜两性在线视频| 欧美成人午夜精品| 三级国产精品欧美在线观看 | 波多野结衣巨乳人妻| av在线播放免费不卡| 亚洲av片天天在线观看| 淫妇啪啪啪对白视频| 国产片内射在线| xxxwww97欧美| 欧美日韩一级在线毛片| 三级毛片av免费| 中文资源天堂在线| 久久午夜综合久久蜜桃| 国产人伦9x9x在线观看| 中文字幕久久专区| 麻豆一二三区av精品| 日韩欧美在线乱码| 亚洲av五月六月丁香网| 亚洲专区字幕在线| 曰老女人黄片| 午夜精品一区二区三区免费看| 怎么达到女性高潮| 国内精品一区二区在线观看| 12—13女人毛片做爰片一| 国产亚洲av高清不卡| 亚洲av片天天在线观看| 三级男女做爰猛烈吃奶摸视频| 中文亚洲av片在线观看爽| 精品国产超薄肉色丝袜足j| 成人av一区二区三区在线看| 少妇人妻一区二区三区视频| 亚洲真实伦在线观看| 亚洲国产精品sss在线观看| 午夜a级毛片| 久99久视频精品免费| 国产又黄又爽又无遮挡在线| 手机成人av网站| 久久久国产欧美日韩av| 亚洲人成77777在线视频| 亚洲人成电影免费在线| 亚洲在线自拍视频| 夜夜看夜夜爽夜夜摸| 18禁国产床啪视频网站| 日韩高清综合在线| 欧洲精品卡2卡3卡4卡5卡区| 国内精品一区二区在线观看| 老熟妇仑乱视频hdxx| 我的老师免费观看完整版| 一区二区三区高清视频在线| 成人国产一区最新在线观看| 国产av在哪里看| 午夜福利在线在线| 人人妻人人澡欧美一区二区| 精品人妻1区二区| 日本 av在线| 男女床上黄色一级片免费看| 视频区欧美日本亚洲| 最近视频中文字幕2019在线8| 国产精品久久久久久人妻精品电影| 18美女黄网站色大片免费观看| 一级作爱视频免费观看| 午夜福利成人在线免费观看| 欧美丝袜亚洲另类 | 久久久久国产精品人妻aⅴ院| 精品一区二区三区视频在线观看免费| 精华霜和精华液先用哪个| 日本 av在线| 天天躁狠狠躁夜夜躁狠狠躁| avwww免费| 免费av毛片视频| 久久久久久国产a免费观看| 午夜亚洲福利在线播放| 国产精品 国内视频| 午夜激情福利司机影院| 夜夜夜夜夜久久久久| 精品欧美一区二区三区在线| √禁漫天堂资源中文www| 久久国产精品人妻蜜桃| 精品久久久久久久久久免费视频| 一级毛片女人18水好多| 天堂av国产一区二区熟女人妻 | 夜夜看夜夜爽夜夜摸| 大型黄色视频在线免费观看| 久久精品国产亚洲av香蕉五月| 日韩三级视频一区二区三区| 午夜亚洲福利在线播放| 亚洲国产欧美一区二区综合| 香蕉国产在线看| 色av中文字幕| 国产精品野战在线观看| 久久久久免费精品人妻一区二区| 不卡av一区二区三区| 黄色 视频免费看| 巨乳人妻的诱惑在线观看| 亚洲片人在线观看| 国产精品九九99| 亚洲av第一区精品v没综合| 午夜免费激情av| 操出白浆在线播放| 九色国产91popny在线| 国产1区2区3区精品| 免费无遮挡裸体视频| 成年版毛片免费区| 一区二区三区激情视频| 少妇粗大呻吟视频| 色噜噜av男人的天堂激情| 久久久水蜜桃国产精品网| 精品国产乱码久久久久久男人| 免费搜索国产男女视频| 欧美中文日本在线观看视频| 又爽又黄无遮挡网站| 欧美久久黑人一区二区| 国产伦在线观看视频一区| 国产精品免费视频内射| 成人av一区二区三区在线看| 日韩成人在线观看一区二区三区| 99久久精品热视频| 久久精品影院6| 国产亚洲精品久久久久5区| 1024视频免费在线观看| 精品熟女少妇八av免费久了| 精品午夜福利视频在线观看一区| 天堂av国产一区二区熟女人妻 | 精品不卡国产一区二区三区| 男女午夜视频在线观看| 村上凉子中文字幕在线| 床上黄色一级片| or卡值多少钱| 18禁裸乳无遮挡免费网站照片| aaaaa片日本免费| 亚洲精品av麻豆狂野| 欧美成人午夜精品| 亚洲乱码一区二区免费版| 国产蜜桃级精品一区二区三区| 性色av乱码一区二区三区2| 亚洲国产欧美一区二区综合| e午夜精品久久久久久久| av超薄肉色丝袜交足视频| 真人一进一出gif抽搐免费| 老司机午夜福利在线观看视频| 日日爽夜夜爽网站| 久久国产精品影院| av国产免费在线观看| 亚洲中文av在线| 十八禁人妻一区二区| 午夜激情福利司机影院| 国产人伦9x9x在线观看| 久久精品国产综合久久久| 亚洲专区国产一区二区| 日韩成人在线观看一区二区三区| 亚洲无线在线观看| 成年人黄色毛片网站| 国产激情久久老熟女| 亚洲国产精品合色在线| 亚洲全国av大片| 欧美黄色片欧美黄色片| 亚洲国产精品成人综合色| 国产日本99.免费观看| 人妻丰满熟妇av一区二区三区| 一本综合久久免费| av免费在线观看网站| 美女午夜性视频免费| 叶爱在线成人免费视频播放| 黄色片一级片一级黄色片| 91老司机精品| 手机成人av网站| 男女午夜视频在线观看| 久久久精品欧美日韩精品| 亚洲七黄色美女视频| 精品一区二区三区四区五区乱码| 国产片内射在线| 欧美 亚洲 国产 日韩一| ponron亚洲| 国产91精品成人一区二区三区| 精品国产乱码久久久久久男人| 久久九九热精品免费| 欧美成人免费av一区二区三区| 波多野结衣巨乳人妻| 亚洲精品中文字幕在线视频| 少妇被粗大的猛进出69影院| 日本免费a在线| 一卡2卡三卡四卡精品乱码亚洲| 亚洲va日本ⅴa欧美va伊人久久| 欧美+亚洲+日韩+国产| 国产成年人精品一区二区| 欧美zozozo另类| 欧美黄色淫秽网站| 男男h啪啪无遮挡| 亚洲电影在线观看av| 日本精品一区二区三区蜜桃| 桃红色精品国产亚洲av| 精品无人区乱码1区二区| 成人一区二区视频在线观看| 身体一侧抽搐| 深夜精品福利| 亚洲色图 男人天堂 中文字幕| 最新美女视频免费是黄的| 麻豆成人av在线观看| 亚洲熟妇中文字幕五十中出| 丝袜人妻中文字幕| 国产高清videossex| 一级黄色大片毛片| 国产精品久久久久久精品电影| 国产精品久久视频播放| 亚洲av成人精品一区久久| 国产三级中文精品| 久久人妻福利社区极品人妻图片| xxxwww97欧美| 国产午夜精品论理片| 无遮挡黄片免费观看| 51午夜福利影视在线观看| 老司机午夜福利在线观看视频| 韩国av一区二区三区四区| 啦啦啦免费观看视频1| 国产精品av久久久久免费| 美女 人体艺术 gogo| 99久久99久久久精品蜜桃| 欧美又色又爽又黄视频| 亚洲熟妇熟女久久| 免费人成视频x8x8入口观看| 午夜激情福利司机影院| 欧美日本视频| 国产av在哪里看| 在线观看美女被高潮喷水网站 | 日本 av在线| 亚洲成av人片免费观看| 精品午夜福利视频在线观看一区| 超碰成人久久| 久久人人精品亚洲av| 久久久国产精品麻豆| 精品熟女少妇八av免费久了| 怎么达到女性高潮| 亚洲 欧美 日韩 在线 免费| 精品国产美女av久久久久小说| 亚洲七黄色美女视频| 可以免费在线观看a视频的电影网站| 黑人操中国人逼视频| 亚洲国产精品999在线| 免费看美女性在线毛片视频| 欧美日韩亚洲国产一区二区在线观看| 国产精品国产高清国产av| 国产一区二区在线观看日韩 | 久久午夜综合久久蜜桃| 亚洲中文字幕日韩| 免费av毛片视频| 国产亚洲av嫩草精品影院| 色综合婷婷激情| 日韩成人在线观看一区二区三区| 香蕉丝袜av| 亚洲国产中文字幕在线视频| 老司机在亚洲福利影院| 日本五十路高清| 一二三四社区在线视频社区8| 美女 人体艺术 gogo| 99国产精品一区二区三区| 日韩精品免费视频一区二区三区| 禁无遮挡网站| 亚洲 欧美一区二区三区| 人人妻,人人澡人人爽秒播| 我要搜黄色片| 蜜桃久久精品国产亚洲av| 中文亚洲av片在线观看爽| 老汉色av国产亚洲站长工具| 国产成人av教育| 超碰成人久久| 1024香蕉在线观看| 香蕉国产在线看| 在线国产一区二区在线| 美女午夜性视频免费| 99精品久久久久人妻精品| 亚洲av五月六月丁香网| 国产精品精品国产色婷婷| 欧洲精品卡2卡3卡4卡5卡区| 国产精品亚洲一级av第二区| 变态另类丝袜制服| xxxwww97欧美| 国产黄片美女视频| 九色成人免费人妻av| 亚洲欧美激情综合另类| 一本大道久久a久久精品| 18禁黄网站禁片免费观看直播| cao死你这个sao货| 久久伊人香网站| 亚洲电影在线观看av| 色综合站精品国产| 最新在线观看一区二区三区| 日韩国内少妇激情av| 99久久精品热视频| 在线观看午夜福利视频| 欧美大码av| 999久久久国产精品视频| 一个人免费在线观看的高清视频| 少妇的丰满在线观看| 777久久人妻少妇嫩草av网站| 深夜精品福利| 国内精品久久久久久久电影| 又黄又粗又硬又大视频| 草草在线视频免费看| 亚洲五月婷婷丁香| 99国产精品99久久久久| 精品久久久久久,| 一边摸一边抽搐一进一小说| 黑人欧美特级aaaaaa片| 大型黄色视频在线免费观看| 亚洲七黄色美女视频| 日本一本二区三区精品| 欧美一区二区精品小视频在线| 国产精品久久久久久人妻精品电影| 国产精品一区二区免费欧美| 一进一出抽搐gif免费好疼| 长腿黑丝高跟| 无人区码免费观看不卡| 亚洲午夜精品一区,二区,三区| 久久婷婷人人爽人人干人人爱| 欧美人与性动交α欧美精品济南到| 欧美日韩瑟瑟在线播放| 免费在线观看完整版高清| 亚洲成a人片在线一区二区| 手机成人av网站| av有码第一页| 日韩欧美在线二视频| 两个人的视频大全免费| 久久中文字幕人妻熟女| 午夜成年电影在线免费观看| 身体一侧抽搐| 国产一区二区三区在线臀色熟女| 婷婷丁香在线五月| 精品国产超薄肉色丝袜足j| or卡值多少钱| 国产成人啪精品午夜网站| 精品国产亚洲在线| 午夜福利免费观看在线| 亚洲精品国产一区二区精华液| 999久久久精品免费观看国产| or卡值多少钱| 亚洲国产精品久久男人天堂| 老司机靠b影院| 欧美日韩亚洲综合一区二区三区_| 黑人欧美特级aaaaaa片| 男女午夜视频在线观看| av免费在线观看网站| 国产一级毛片七仙女欲春2| 国产精品自产拍在线观看55亚洲| 亚洲男人的天堂狠狠| 国产99白浆流出| 久久久精品欧美日韩精品| 色尼玛亚洲综合影院| 精品人妻1区二区| 嫁个100分男人电影在线观看| 超碰成人久久| 女生性感内裤真人,穿戴方法视频| 亚洲欧美精品综合久久99| 中文亚洲av片在线观看爽| 十八禁网站免费在线| e午夜精品久久久久久久| 久久精品91蜜桃| 国内毛片毛片毛片毛片毛片| 久久久久久久久免费视频了| 搡老岳熟女国产| 别揉我奶头~嗯~啊~动态视频| 免费无遮挡裸体视频| av有码第一页| 精品久久久久久久久久久久久| 一级a爱片免费观看的视频| 在线观看免费午夜福利视频| 亚洲一区中文字幕在线| 亚洲国产欧美人成| www.精华液| 精品电影一区二区在线| 欧美zozozo另类| 国产av麻豆久久久久久久| avwww免费| 午夜老司机福利片| 国产黄片美女视频| 国产精品一区二区三区四区久久| 可以免费在线观看a视频的电影网站| 亚洲熟女毛片儿| 国产精品爽爽va在线观看网站| 亚洲天堂国产精品一区在线| 亚洲真实伦在线观看| tocl精华| 成人国产综合亚洲| 一个人免费在线观看的高清视频| 伦理电影免费视频| 午夜影院日韩av| 国产欧美日韩一区二区精品| 老司机在亚洲福利影院| 亚洲性夜色夜夜综合| 99热这里只有是精品50| 国产精品香港三级国产av潘金莲| 国产麻豆成人av免费视频| 亚洲成人免费电影在线观看| 黄频高清免费视频| 亚洲国产日韩欧美精品在线观看 | 男人的好看免费观看在线视频 | 国内毛片毛片毛片毛片毛片| 精品一区二区三区av网在线观看| 亚洲精品久久成人aⅴ小说| 免费在线观看成人毛片| 91九色精品人成在线观看| 美女黄网站色视频| 亚洲男人的天堂狠狠| 久久国产乱子伦精品免费另类| 欧美日韩一级在线毛片| 亚洲 欧美 日韩 在线 免费| 18禁美女被吸乳视频| 99热这里只有是精品50| av在线播放免费不卡| xxxwww97欧美| 久久久久久九九精品二区国产 | 亚洲专区字幕在线| 国产亚洲欧美在线一区二区| 91大片在线观看| 午夜a级毛片| 亚洲第一欧美日韩一区二区三区| 日韩欧美国产一区二区入口| 久久久国产成人免费| 在线观看www视频免费| 国产成人精品无人区| 不卡av一区二区三区| 日日干狠狠操夜夜爽| 一本综合久久免费| 亚洲性夜色夜夜综合| 成人欧美大片| 99精品欧美一区二区三区四区| 欧美久久黑人一区二区| 欧美黄色片欧美黄色片| 中文在线观看免费www的网站 | 在线看三级毛片| 波多野结衣高清无吗| 麻豆成人午夜福利视频| av福利片在线观看| 制服丝袜大香蕉在线| 日日摸夜夜添夜夜添小说| 国产成年人精品一区二区| 国产av一区在线观看免费| 亚洲精品久久成人aⅴ小说| 久久草成人影院| 男女午夜视频在线观看| 日本一区二区免费在线视频| 国产精品综合久久久久久久免费| 国产亚洲av嫩草精品影院| 国产亚洲av高清不卡| 久久天躁狠狠躁夜夜2o2o| 国产精品1区2区在线观看.| 精品免费久久久久久久清纯| 欧美高清成人免费视频www| 久久久久久国产a免费观看| 人妻丰满熟妇av一区二区三区| 国产91精品成人一区二区三区| 亚洲国产中文字幕在线视频| 一级a爱片免费观看的视频| 人妻久久中文字幕网| 欧美av亚洲av综合av国产av| 亚洲成a人片在线一区二区| 久久欧美精品欧美久久欧美| 99久久国产精品久久久| 亚洲 国产 在线| 黄色片一级片一级黄色片| 18禁观看日本| 免费在线观看黄色视频的| 一区二区三区高清视频在线| 欧美人与性动交α欧美精品济南到| 久久香蕉精品热| 久久天躁狠狠躁夜夜2o2o| 亚洲色图av天堂| 日日干狠狠操夜夜爽| 成人国产综合亚洲| 亚洲七黄色美女视频| 日韩精品中文字幕看吧| cao死你这个sao货| 国产乱人伦免费视频| 母亲3免费完整高清在线观看| 老熟妇乱子伦视频在线观看| 亚洲国产精品成人综合色| 国产乱人伦免费视频| 国产成人影院久久av| 精品国产亚洲在线| 一本一本综合久久| 两个人的视频大全免费| 久久久国产成人免费| 午夜福利高清视频| 亚洲成人久久爱视频| 男人舔女人的私密视频| 女警被强在线播放| 男人舔女人的私密视频| 午夜激情福利司机影院| 波多野结衣高清无吗| 欧美色欧美亚洲另类二区| 亚洲国产精品sss在线观看| 亚洲欧美精品综合久久99| 精品国产美女av久久久久小说| 国产在线精品亚洲第一网站| 久久久久久久午夜电影| 熟妇人妻久久中文字幕3abv| 在线观看美女被高潮喷水网站 | 久久这里只有精品19| 午夜福利高清视频| 首页视频小说图片口味搜索| 久久亚洲精品不卡| 亚洲人成网站高清观看| 国产av不卡久久| 亚洲av片天天在线观看| 国内毛片毛片毛片毛片毛片| 白带黄色成豆腐渣| 这个男人来自地球电影免费观看| 欧美日韩黄片免|