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

    基于希爾伯特-黃變換的潮汐分析方法研究

    2020-05-11 09:37:16張佳麗張學(xué)峰張安民
    關(guān)鍵詞:模態(tài)分析方法

    張?亮,張佳麗,張學(xué)峰,張安民

    基于希爾伯特-黃變換的潮汐分析方法研究

    張?亮,張佳麗,張學(xué)峰,張安民

    (天津大學(xué)海洋科學(xué)與技術(shù)學(xué)院,天津 300072)

    經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)是希爾伯特-黃變換(HHT)的核心組成部分.在對傳統(tǒng)HHT方法研究與分析的基礎(chǔ)上,發(fā)展了改進(jìn)的集合經(jīng)驗(yàn)?zāi)B(tài)分解的方法(MEEMD)對海洋潮汐進(jìn)行分析.MEEMD方法根據(jù)水位數(shù)據(jù)序列的特點(diǎn),通過添加成對白噪聲和計(jì)算排列熵的方式來剔除分解過程中出現(xiàn)的虛假本征模態(tài)函數(shù)(IMF),從而抑制了模態(tài)混疊現(xiàn)象.?dāng)?shù)值模擬和實(shí)測水位數(shù)據(jù)實(shí)驗(yàn)表明,MEEMD方法可以在有效抑制模態(tài)混疊現(xiàn)象的基礎(chǔ)上,對含有噪聲的原始水位序列實(shí)現(xiàn)自適應(yīng)濾波.MEEMD方法分解得到的高頻IMF1和IMF2分量占原始水位數(shù)據(jù)總能量的90%以上,IMF1和IMF2分量的平均周期和平均振幅分別與調(diào)和分析得到的半日分潮與全日分潮結(jié)果基本一致.通過分析認(rèn)為,IMF1和IMF2之和是主要半日分潮、全日分潮以及淺海分潮的疊加.傳統(tǒng)的潮汐調(diào)和分析方法有賴于長期的水位記錄,且需要水位數(shù)據(jù)的一些先驗(yàn)知識,涉及到的天文知識與數(shù)學(xué)計(jì)算也較為復(fù)雜.準(zhǔn)調(diào)和分析可以處理短期的水位數(shù)據(jù),但仍然涉及天文知識與復(fù)雜數(shù)學(xué)計(jì)算,對數(shù)據(jù)的完整性也要求較高.當(dāng)收集到的水位數(shù)據(jù)因?yàn)榇嬖谠肼暬蛘咝畔⑷笔В蛔阋赃M(jìn)行傳統(tǒng)的調(diào)和或準(zhǔn)調(diào)和分析時(shí),可以使用HHT方法進(jìn)行有效的補(bǔ)充分析,以獲得主要半日分潮、全日分潮和余水位的幅值信息,因此,該方法可以作為一種很好的潮汐分析的補(bǔ)充方法.

    改進(jìn)的集合經(jīng)驗(yàn)?zāi)B(tài)分解;希爾伯特-黃變換;濾波;分潮

    海洋潮汐是海水受到月球和太陽的吸引力作用而產(chǎn)生的一種周期性升降運(yùn)動(dòng).產(chǎn)生潮汐現(xiàn)象的主要原因是由于地球上各點(diǎn)距離月球的相對位置不同,所受到的引力也不同,這種不同導(dǎo)致地球上的海水做相對運(yùn)動(dòng),引起海面升降.而實(shí)際的潮汐受到陸地、海底地形起伏變化,海底摩擦以及地球自轉(zhuǎn)等影響是十分復(fù)雜的,呈現(xiàn)出非平穩(wěn)、非線性的特征.傳統(tǒng)的潮汐分析方法有調(diào)和分析和準(zhǔn)調(diào)和分析.潮汐的調(diào)和分析方法是建立在分潮的概念上,將潮汐看成是以不同頻率傳播的各種分潮疊加產(chǎn)生的現(xiàn)象,目的是求出各個(gè)分潮的振幅和遲角[1-2].潮汐調(diào)和分析根據(jù)過去一段時(shí)間的潮汐觀測資料,計(jì)算各個(gè)分潮的調(diào)和常數(shù),用于潮汐預(yù)報(bào)以及計(jì)算各種基準(zhǔn)面等,其中涉及天文知識與數(shù)學(xué)計(jì)算,整個(gè)過程較為復(fù)雜.潮汐準(zhǔn)調(diào)和分析適用于手工計(jì)算,通過對一天或幾天的短期觀測數(shù)據(jù)計(jì)算獲得準(zhǔn)調(diào)和分潮.準(zhǔn)調(diào)和分潮是以一個(gè)分潮為主,其余分潮的調(diào)和常數(shù)與這個(gè)主要分潮建立一定關(guān)系后合并為一項(xiàng),此項(xiàng)的振幅和角速率不再是常量,而是隨時(shí)間緩慢變化[3].希爾伯特-黃變換(Hilbert-Huang transform,HHT)是Haung等于1998年提出的一種分析非線性、非平穩(wěn)信號的方法[4-6]. HHT作為一種新的信號分析理論,其關(guān)鍵部分為經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition,EMD).大部分通過EMD分解得到的固有模態(tài)函數(shù)(instrinsic mode function,IMF)可以反映真實(shí)的物理過程,同時(shí)也有大量實(shí)驗(yàn)表明,EMD方法是目前提取信號趨勢或均值的最好方法.在海洋科學(xué)研究中,EMD方法已應(yīng)用于海平面變化趨勢、風(fēng)暴潮和海溫異常預(yù)測等方面的研究[7-10].Cheng等[11]使用EMD的方法對潮汐的振幅進(jìn)行分析得到海平面上升和氣候變化的關(guān)系;Pan等[12]應(yīng)用EMD方法對具有非平穩(wěn)、非線性特征的河流潮汐進(jìn)行分析并驗(yàn)證了方法的有效性.但EMD方法本身存在一些問題,其中一個(gè)主要問題是模態(tài)混疊,針對這個(gè)問題,很多學(xué)者提出了包括集合經(jīng)驗(yàn)?zāi)B(tài)分解(ensemble empirical mode decomposition,EEMD)[13]、補(bǔ)充的EEMD(comple-mentary EEMD,CEEMD)[14]、改進(jìn)的EEMD(modi-fied EEMD,MEEMD)等解決方法[15]. 另外,黃大吉等[16]針對HHT存在的端點(diǎn)處理問題,提出了包絡(luò)的極值延拓法和鏡像閉合延拓法,兩種方法都能很好地處理端點(diǎn)問題,使得EMD分解得到合理的各個(gè)IMF模態(tài).在前述文獻(xiàn)研究的基礎(chǔ)上,本文提出了采用基于MEEMD的希爾伯特-黃變換方法對潮汐信號進(jìn)行分析,有效抑制了模態(tài)混疊問題并提高了頻率分辨率,同時(shí)該方法還具有很好的自適應(yīng)濾波效果,是潮汐分析的一種有效方法.

    1?HHT方法簡介

    1.1?MEEMD方法

    EMD方法中出現(xiàn)的模態(tài)混疊現(xiàn)象是指在一個(gè)IMF中包含差異極大的特征時(shí)間尺度,或者相近的特征時(shí)間尺度分布在不同的IMF中,造成波形混疊,難以辨認(rèn)[17-18].為了抑制模態(tài)混疊現(xiàn)象,提出了EEMD方法和CEEMD方法.CEEMD方法是EEMD方法的改進(jìn),該方法通過在原始信號中添加成對幅值有限且均值為零的白噪聲,再通過對分解結(jié)果進(jìn)行總體平均來消除白噪聲的影響,進(jìn)而達(dá)到抑制模態(tài)混疊現(xiàn)象的目的.但是分解結(jié)果易受到添加白噪聲參數(shù)的影響,若選擇不當(dāng)會(huì)產(chǎn)生虛假IMF分量并破壞分解的自適應(yīng)性.MEEMD方法為了彌補(bǔ)上述CEEMD方法的不足,提出了用計(jì)算排列熵(permutation entropy,PE)[19]的方法來剔除信號中的虛假分量,篩選出主要IMF分量.排列熵是一種能有效度量信號復(fù)雜程度的物理量,信號越復(fù)雜,排列熵值越高,反之,則熵值越低[20].而含有噪聲的信號往往具有較高的隨機(jī)性,信號復(fù)雜程度高,排列熵值相比于不含噪聲的信號更高.

    MEEMD算法具體步驟如下所述.

    式中為待分析信號()的采樣點(diǎn)個(gè)數(shù).

    通過MEEMD方法,將原始數(shù)據(jù)序列分解為個(gè)IMF分量及余項(xiàng)之和,余項(xiàng)代表數(shù)據(jù)的平均趨勢.

    1.2?希爾伯特變換

    希爾伯特變換是在MEEMD分解的基礎(chǔ)上進(jìn)行的,利用希爾伯特變換求解每階IMF的瞬時(shí)頻率,從而得到信號的振幅-時(shí)間-頻率表示,即希爾伯特譜.通過計(jì)算希爾伯特譜對整個(gè)周期的積分得到希爾伯特邊際譜.

    將式(5)分解得到的各IMF分量進(jìn)行希爾伯特變換,得到

    經(jīng)過計(jì)算獲得信號的希爾伯特譜可表示為

    希爾伯特譜是瞬時(shí)振幅在頻率-時(shí)間平面上的分布,具有良好的時(shí)頻聚集性[5].

    式中為周期.邊際譜描述了信號在整個(gè)時(shí)頻內(nèi),同一頻率幅值或能量疊加值在頻域的分布情況.

    瞬時(shí)能量表示每個(gè)時(shí)刻在整個(gè)頻域所積累的能量.

    2?數(shù)值模擬分析與方法有效性驗(yàn)證

    2.1?數(shù)值模擬水位序列的MEEMD與EMD分解

    為了驗(yàn)證HHT方法的有效性,采用數(shù)值模擬的方式來進(jìn)行驗(yàn)證.潮汐可以看作是多個(gè)分潮的疊加,調(diào)和分析方法把潮位看成是多個(gè)分潮余弦振動(dòng)之和,根據(jù)最小二乘或波譜分析原理,由實(shí)測的潮位數(shù)據(jù)計(jì)算出各分潮的平均振幅和遲角,其中分潮的平均振幅和遲角為分潮的潮汐調(diào)和常數(shù)[2].根據(jù)潮汐調(diào)和分析理論,不考慮非天文潮位,設(shè)潮汐為個(gè)分潮疊加,得到的潮高表達(dá)式為

    準(zhǔn)調(diào)和分析是通過引入不同分潮之間的關(guān)系,將不同分潮合并為一項(xiàng),獲得振幅和角速率隨時(shí)間緩慢變化的準(zhǔn)調(diào)和分潮[3].這樣,潮位的準(zhǔn)調(diào)和分潮表達(dá)式可以表示為

    式中:0代表平均海面高度加上長周期分潮;和是和天文情況有關(guān)的變量,分別叫作準(zhǔn)調(diào)和分潮的振幅系數(shù)和遲角訂正;代表準(zhǔn)調(diào)和分潮.

    這里根據(jù)資料來源處的驗(yàn)潮站數(shù)據(jù),采用固定的振幅和遲角進(jìn)行潮汐模擬計(jì)算,選取M2、S2、K1和Mf 四個(gè)分潮,其中M2分潮振幅為10m,其余分潮采用相對振幅.輸入?yún)?shù)如表1所示.

    由于水位觀測序列信息中除了水位信息外還包含了觀測噪聲等干擾,例如涌浪噪聲會(huì)使水位序列曲線受到高頻振動(dòng)的影響,在水位觀測時(shí)刻有強(qiáng)烈的上下振動(dòng),所以在數(shù)值模擬的水位序列中加入了功率為0.05的高斯白噪聲來模擬噪聲等因素的干擾.

    對上述水位模擬信號分別采用EMD與MEEMD方法進(jìn)行分解,得到的IMF分量個(gè)數(shù)分別為8和5,模擬水位序列和分解結(jié)果如圖1~圖3所示.

    表1?水位序列數(shù)值模擬參數(shù)

    Tab.1 Numerically simulated parameters of water-level series

    從圖2中可以看出EMD分解結(jié)果受噪聲影響很大,出現(xiàn)了明顯的模態(tài)混疊現(xiàn)象.高頻分量IMF1將加入的白噪聲作為一個(gè)模態(tài)與模擬水位序列中的高頻成分疊加分解出來,成為一個(gè)不具備物理意義的虛假分量;而IMF2和IMF3也出現(xiàn)了明顯的模態(tài)混疊現(xiàn)象,其分解結(jié)果也不具備任何物理意義;IMF4~I(xiàn)MF8沒有明顯的模態(tài)混疊現(xiàn)象,但是也沒有周期性規(guī)律.因此,從圖中觀察可知,EMD方法的分解效果易受到噪聲等隨機(jī)因素的干擾,產(chǎn)生模態(tài)混疊現(xiàn)象與虛假分量,使得分解結(jié)果不再具備物理意義,自適應(yīng)性較差.對于有噪聲的數(shù)據(jù)序列,在進(jìn)行EMD計(jì)算之前還需先進(jìn)行降噪等預(yù)處理.從圖3中可以看出,MEEMD方法分解得到的各分量沒有出現(xiàn)模態(tài)混疊現(xiàn)象,特別是IMF1和IMF2分量的周期性規(guī)律明顯,IMF4和IMF5也有周期性特征.從IMF分量個(gè)數(shù)上來看,MEEMD方法得到了5個(gè)IMF分量,相對于EMD方法減少了3個(gè)IMF分量,有效減少了虛假分量.其中只有IMF3不對應(yīng)原始模擬水位序列的任何一個(gè)分潮成分,且沒有明顯的周期特征,是一個(gè)虛假分量.通過分析可知,MEEMD方法相對于EMD方法有效抑制了模態(tài)混疊現(xiàn)象,減少了虛假分量.除此之外,該方法可以實(shí)現(xiàn)在分解過程中的自適應(yīng)濾波,有效解決了噪聲等因素的干擾.原始的水位模擬序列有4個(gè)分潮,MEEMD方法和EMD方法得到的IMF分量均多于4個(gè),因此要對分解結(jié)果做希爾伯特譜分析.

    希爾伯特譜具有很好的時(shí)頻聚集性,可以在圖4中直觀地看出IMF分量的時(shí)頻特性,頻率的單位為每小時(shí)循環(huán)數(shù),表示為c/h.從圖4(a)MEEMD得到的Hilbert-Huang譜中可以看到2條清晰的頻率線,這說明分解得到的其中2個(gè)高頻分量頻率集中,分解效果較好.從圖4(b)的邊際譜中可以看到3處頻率集中的位置,其中有1處位于低頻段,這說明分解有效,獲得了低頻分量,但所占能量比例較?。畧D4(c)中EMD得到的Hilbert-Huang譜中頻率線模糊,這說明分量的頻率分散,從圖4(d)的邊際譜中觀察可知分解結(jié)果頻率分散,有較多虛假分量存在,分解效果較差.

    上述數(shù)值模擬的結(jié)果初步表明,MEEMD方法能夠?qū)性肼暤乃粩?shù)據(jù)進(jìn)行自適應(yīng)濾波并有效抑制模態(tài)混疊現(xiàn)象,分解效果較好.但是該方法還具有局限性,不能將頻率十分相近的分潮分離出來,而是疊加在了一起,如IMF1是半日分潮M2和S2的疊加,分解結(jié)果對于長周期分潮也沒有很好的體現(xiàn).因此,基于MEEMD的方法可以用于提取水位序列中的半日潮族與全日潮族;對于淺海海域,頻率較高的淺海分潮也包含在高頻IMF1分量中被提取出來,這樣得到的IMF1和IMF2之和除了全日分潮和半日分潮成分,還有淺海分潮成分.

    2.2?數(shù)值模擬水位序列的能量譜

    根據(jù)式(12)計(jì)算希爾伯特能量,可以得到水位序列的能量在各IMF分量的分布情況,如圖5所示.從圖中可以看出,高頻的IMF1和IMF2分量占總能量的絕大部分,經(jīng)計(jì)算可知其能量值占總能量的96.39%,是主要的波動(dòng)部分.

    (d)EMD分解結(jié)果的邊際譜

    圖5?IMF分量能量分布

    3?實(shí)驗(yàn)分析

    3.1?資料來源

    為驗(yàn)證HHT方法的適用性,本文對實(shí)測的水位記錄進(jìn)行了實(shí)驗(yàn)與分析.資料來源于美國國家海洋和大氣管理局(NOAA),本研究選取了位于美國南灘編號為9435380驗(yàn)潮站水位記錄中,2016—2017年的逐時(shí)水位記錄作為實(shí)測分析數(shù)據(jù).NOAA提供的資料顯示,該驗(yàn)潮站建于1967年1月28日,位于北緯44°37.5′N,西經(jīng)124°2.7′W,時(shí)間子午線為120°W,起算基準(zhǔn)面為平均低低潮面(MLLW).將資料數(shù)據(jù)進(jìn)行MEEMD分解,對分解得到的各IMF分量進(jìn)行希爾伯特變換后得到希爾伯特譜、邊際譜與能量譜,并計(jì)算各IMF分量的平均振幅與平均周期.將基于HHT方法分析得到的結(jié)果與NOAA提供的該驗(yàn)潮站潮汐調(diào)和分析數(shù)據(jù)以及經(jīng)過準(zhǔn)調(diào)和分析獲得的結(jié)果進(jìn)行對比,驗(yàn)證了方法的適用性.

    3.2?實(shí)測水位記錄的分析結(jié)果

    圖6是2016—2017年的長期水位序列經(jīng)MEEMD分解后得到的9個(gè)IMF分量與余項(xiàng).從圖中可以看出高頻IMF分量的分解效果較好,同一IMF分量中沒有出現(xiàn)模態(tài)混疊現(xiàn)象.其中IMF1和IMF2兩個(gè)高頻分量有明顯的周期性波動(dòng).因?yàn)槌毕牟▌?dòng)是由多個(gè)分潮疊加而成的復(fù)雜振動(dòng),并且原始數(shù)據(jù)序列中還帶有噪聲等因素的干擾,因此分解的結(jié)果中IMF分量較多,需要進(jìn)一步的分析.

    為了更加清晰地了解各IMF分量的時(shí)頻特性,對IMF1~I(xiàn)MF9進(jìn)行Hilbert譜分析,得到Hilbert-Huang譜以及邊際譜,如圖7所示.圖7(a)是Hilbert-Huang譜,反映了水位序列的時(shí)間-頻率特性.從圖7(a)中可以發(fā)現(xiàn)兩條清晰的頻率線,分別在頻率為0.08c/h與0.04c/h附近,這表示其中兩個(gè)主要高頻IMF分量的頻率聚集性良好,具有明顯的周期性;在低頻0~0.02c/h之間也有頻率聚集現(xiàn)象,這表明在低頻范圍內(nèi)IMF分量的頻率也較為集中,但與兩處高頻聚集區(qū)相比并不明顯.圖7(b)是邊際譜,反映了水位序列的幅值-頻率特性.與Hilbert-Huang譜相比,在低頻0~0.02c/h之間的頻率聚集現(xiàn)象更為明顯.這部分頻率可能有長周期分潮與之對應(yīng),但是相對于高頻分量所占比例較?。ㄟ^實(shí)驗(yàn)分析可知,實(shí)測數(shù)據(jù)的分析結(jié)果與數(shù)值模擬的結(jié)果相符合.

    除了1年以上的長期觀測數(shù)據(jù),還需要對短期的水位數(shù)據(jù)進(jìn)行分析.圖8(a)和圖8(b)分別為選取了2017年中某3d和7d的水位數(shù)據(jù)分解結(jié)果.短期觀測的水位數(shù)據(jù)得到較少的IMF分量,3d和7d兩組水位序列均分解出3個(gè)IMF分量.

    圖6?2016—2017年原始水位序列與IMF1~I(xiàn)MF9及余項(xiàng)

    圖7?2016—2017年水位序列Hilbert-Huang譜與邊際譜

    從圖8中可以看出,短期水位序列的分解結(jié)果沒有出現(xiàn)明顯的模態(tài)混疊現(xiàn)象,各IMF分量均具有周期波動(dòng)特性.將選取的3d水位序列進(jìn)行時(shí)頻分析,得到圖9(a)和圖9(b)的Hilbert-Huang譜和邊際譜圖.從圖中可以看出,短期水位序列也具有類似于長期水位序列的時(shí)間-頻率特性,頻率主要在0.08c/h和0.04c/h附近聚集,且振動(dòng)能量在0.08c/h附近更為集中;在低頻和更高頻的區(qū)域幾乎沒有能量分布.

    圖8?短期原始水位序列及分解結(jié)果

    圖9?3 d水位序列Hilbert-Huang譜與邊際譜

    圖10是歸一化的各IMF分量的能量以及各IMF分量與原始水位序列相關(guān)系數(shù)的對比分析.其中相關(guān)系數(shù)表征IMF分量與原始序列的相似程度,相關(guān)系數(shù)較大的IMF分量與原始序列相似程度高,帶有更多能量.從圖中的對比觀察可知,各IMF的能量分布與相關(guān)系數(shù)分布趨勢一致,高頻IMF分量占據(jù)原始水位序列的絕大部分振動(dòng)能量,經(jīng)計(jì)算能量占有率為94.98%,這與譜分析的結(jié)果一致.短期水位序列的分析結(jié)果也是如此,高頻IMF1和IMF2分量幾乎占據(jù)了振動(dòng)的全部能量.通過實(shí)驗(yàn)分析可知,不管是長期還是短期實(shí)測水位數(shù)據(jù),分析結(jié)果均與數(shù)值模擬的結(jié)果相符合.

    圖10?歸一化能量與相關(guān)系數(shù)對比

    3.3?HHT分析方法與潮汐調(diào)和分析結(jié)果對比

    NOAA提供的資料顯示,該驗(yàn)潮站處由潮汐調(diào)和分析得到共37個(gè)分潮,其中主要半日分潮、全日分潮和淺水分潮的調(diào)和常數(shù)與周期如表2所示,振幅高度以當(dāng)?shù)亻L期平均海面(MSL)為零點(diǎn)起算.

    表2?驗(yàn)潮站分潮參數(shù)

    Tab.2?Tidal parameters of tide gauge stations

    為了與調(diào)和分析的結(jié)果進(jìn)行對比,計(jì)算出2016—2017年的長期水位序列各IMF的平均周期和平均振幅.根據(jù)式(8)與式(9)計(jì)算得到各IMF的瞬時(shí)幅值與瞬時(shí)頻率序列后,再分別求這兩個(gè)序列的均值,可得平均振幅和平均周期如表3所示.

    通過對比表2和表3的分析結(jié)果可知,IMF1和IMF2分量的平均周期與平均振幅分別與半日分潮和全日分潮接近.IMF1分量的平均周期與半日分潮周期一致,特別是與M2分潮的周期最為接近,平均振幅略大于M2分潮的振幅,這說明該分量包含的不是單一分潮,而是多個(gè)半日分潮的疊加,其中可能還包含頻率更高但是比例很小的淺海分潮成分,因?yàn)镸2分潮占比最大,所以更多體現(xiàn)出M2分潮的性質(zhì).IMF2分量的平均周期與全日分潮周期一致,特別是與K1分潮的周期最為接近,但是平均振幅卻小于K1分潮的振幅,這說明該分量中部分成分被分解到了其他模態(tài)中.因?yàn)樵诜纸膺^程中各IMF分量之間的正交性不能完全保證,MEEMD方法雖然可以有效抑制模態(tài)混疊現(xiàn)象卻不能完整解決,在分解時(shí)會(huì)存在能量損失.

    表3?各IMF分量的平均振幅與平均周期

    Tab.3?Average amplitude and period of IMF components

    3.4?HHT分析方法與潮汐準(zhǔn)調(diào)和分析結(jié)果對比

    根據(jù)潮汐的準(zhǔn)調(diào)和分析方法,分別選取3組長度為3d的水位記錄,進(jìn)行準(zhǔn)調(diào)和分析,得到K1、O1、M2、S2、M4、MS4準(zhǔn)調(diào)和分潮的振幅(單位:m),如表4所示.

    表4?準(zhǔn)調(diào)和分潮振幅

    Tab.4?Quasi-harmonic tidal amplitude m

    對3組實(shí)驗(yàn)數(shù)據(jù)進(jìn)行HHT分析,得到各IMF分量的平均周期和平均振幅如表5所示.

    表5?各次實(shí)驗(yàn)IMF分量的平均振幅和平均周期

    Tab.5 Average amplitude and period of IMF components in each experiment

    通過表5中的數(shù)據(jù)可知,實(shí)驗(yàn)1和實(shí)驗(yàn)3分別獲得了3個(gè)IMF分量,而實(shí)驗(yàn)2獲得了2個(gè)IMF分量,但是3次實(shí)驗(yàn)獲得的IMF1和IMF2分量的平均振幅和平均周期十分相近,平均振幅可分別對應(yīng)于表4中準(zhǔn)調(diào)和分潮結(jié)果中M2和K1分潮的振幅.而通過表3和表5的對比可知,在平均周期上IMF1和IMF2分量的平均周期分別介于M2、S2和K1、O1之間,所以可以認(rèn)為這兩個(gè)分量分別為半日準(zhǔn)調(diào)和分潮和全日準(zhǔn)調(diào)和分潮的疊加;因?yàn)樵谡穹螹2和K1分潮占主導(dǎo)地位,且在分解過程中存在能量損失,因此IMF1和IMF2分量分別主要對應(yīng)于M2和K1準(zhǔn)調(diào)和分潮.IMF3分量的平均振幅和平均周期與表4中的準(zhǔn)調(diào)和分析結(jié)果沒有明顯的對應(yīng)關(guān)系,且實(shí)驗(yàn)1和實(shí)驗(yàn)3的結(jié)果差別很大,所以IMF3分量不具備明顯的分潮特性,可能為其他長周期分潮的疊加.因?yàn)槊看芜x取的實(shí)驗(yàn)數(shù)據(jù)長度只有3d,對應(yīng)于一個(gè)朔望月(29.5306d)的不同時(shí)期,所以獲得的平均振幅結(jié)果會(huì)有一定差異.例如當(dāng)朔或望時(shí)為大潮期,獲得的分量平均振幅較大;而當(dāng)上、下弦月時(shí)為小潮期,獲得的分量平均振幅較小,表3和表4的振幅分析結(jié)果與此相一致.3d水位數(shù)據(jù)不能很好地體現(xiàn)長周期分潮的性質(zhì),因此雖然IMF3具有明顯的波動(dòng)規(guī)律,但選取不同時(shí)間的水位數(shù)據(jù)得到的IMF3分量的平均周期結(jié)果可能不同.準(zhǔn)調(diào)和分潮獲得了淺水分潮成分,在HHT分析結(jié)果中并沒有體現(xiàn),可能包含在IMF1分量中或是損失于分解過程中.

    4?結(jié)?語

    希爾伯特-黃變換(HHT)是一種新的處理非平穩(wěn)、非線性信號或數(shù)據(jù)的方法,本文采用基于改進(jìn)的總體平均經(jīng)驗(yàn)?zāi)B(tài)分解(MEEMD)的HHT方法對原始水位數(shù)據(jù)進(jìn)行分析.通過MEEMD方法,把具有潮汐波動(dòng)與其他噪聲干擾的復(fù)雜驗(yàn)潮站水位數(shù)據(jù)記錄,分解為有限數(shù)目且具有一定物理意義的本征模態(tài)函數(shù)(IMF).通過希爾伯特變換計(jì)算各IMF的瞬時(shí)頻率、幅值等得到數(shù)據(jù)序列的時(shí)頻表示與能量分布情況.通過數(shù)值模擬與實(shí)測數(shù)據(jù)實(shí)驗(yàn)分析,HHT方法可以對原始數(shù)據(jù)進(jìn)行自適應(yīng)濾波分解,得到的高頻IMF分量的時(shí)頻聚集性良好,但是包括長周期分潮在內(nèi)的低頻分量的分解效果并不理想.將分解結(jié)果與傳統(tǒng)的潮汐調(diào)和分析得到的分潮進(jìn)行對比可知,IMF1和IMF2分量分別與主要半日分潮和全日分潮一致,其和可以認(rèn)為是主要半日分潮、全日分潮以及高頻淺海分潮的疊加.雖然MEEMD方法有效抑制了模態(tài)混疊現(xiàn)象,但是對于結(jié)構(gòu)復(fù)雜的信號或數(shù)據(jù)仍有模態(tài)混疊現(xiàn)象存在,需要對各分量做進(jìn)一步細(xì)致的分析.對于3~7d短期的實(shí)測數(shù)據(jù),HHT方法的分解結(jié)果沒有出現(xiàn)明顯的模態(tài)混疊現(xiàn)象,高頻IMF1和IMF2分量的平均振幅分別與準(zhǔn)調(diào)和分析獲得的M2和K1相一致,且分析結(jié)果具有較高精度.雖然本研究沒有使用HHT方法對淺水分潮進(jìn)行進(jìn)一步的分析,但該方法相比于傳統(tǒng)調(diào)和分析方法來說,仍然具有一定優(yōu)勢,尤其當(dāng)收集到的水位數(shù)據(jù)沒有經(jīng)過預(yù)處理(如含有涌浪等高頻噪聲)或水位數(shù)據(jù)殘缺不全、數(shù)據(jù)來源信息不明而不足以進(jìn)行傳統(tǒng)的調(diào)和或準(zhǔn)調(diào)和分析時(shí),可以使用HHT方法進(jìn)行有效的補(bǔ)充分析,以獲得主要半日分潮、全日分潮和余水位的幅值信息.該方法目前還不能夠用于精確地預(yù)報(bào)潮位,但是在采用該方法提取有效水位特征信息基礎(chǔ)上,可以借助神經(jīng)網(wǎng)絡(luò)等方法實(shí)現(xiàn)對潮位的預(yù)測和預(yù)報(bào),這部分內(nèi)容需要做進(jìn)一步的研究.

    [1] 張鳳燁,魏澤勛,王新怡,等. 潮汐調(diào)和分析方法的探討[J]. 海洋科學(xué),2011,35(6):68-75.

    Zhang Fengye,Wei Zexun,Wang Xinyi,et al. Tidal harmonic analysis[J]. Marine Sciences,2011,35(6):68-75(in Chinese).

    [2] 宋艷朋. 潮汐調(diào)和分析預(yù)報(bào)與基準(zhǔn)面計(jì)算軟件實(shí)現(xiàn)及南海應(yīng)用研究[D]. 青島:山東科技大學(xué),2017.

    Song Yanpeng. The Implementation of Tide Harmonic Analysis Forecasting and Datum Calculation Software and Application in South China Sea[D]. Qingdao:Shandong University Of Science and Technology,2017(in Chinese).

    [3] 方國洪. 潮汐和潮流的分析和預(yù)報(bào)[M]. 北京:海洋出版社,1986.

    Fang Guohong. Analysis and Forecast of Tides and Tidal [M]. Beijing:China Ocean Press,1986(in Chinese).

    [4] Huang N E,Shen Z,Long S R,et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]. Proceedings:Mathematical,Physical and Engineering Sciences,1998,454(1971):903-995.

    [5] Huang N E,Shen Z,Long S R. A new view of nonlinear water waves:The Hilbert spectrum[J]. Ann Rrev Fluid Mech,1999,31(1):417-457.

    [6] 黃誠惕. 希爾伯特-黃變換及其應(yīng)用研究[D]. 成都:西南交通大學(xué),2006.

    Huang Chengti. Study on Hilbert-Huang Transform and Its Application[D]. Chengdu:Southwest Jiaotong University,2006(in Chinese)

    [7] Chen X,F(xiàn)eng Y,Huang N E. Global sea level trend during 1993—2012[J]. Global and Planetary Change,2014,112:26-32.

    [8] Chambers D P. Evaluation of empirical mode decomposition for quantifying multi-decadal variations in sea level records[J]. Nonlinear Processes in Geophysics,2015,1(2):157-166.

    [9] 楊金湘,袁方超,李郅明,等. EMD方法在中國沿岸風(fēng)暴潮增水分析中的應(yīng)用[J]. 海洋通報(bào),2016,35(2):157-162,169.

    Yang Jinxiang,Yuan Fangchao,Li Zhiming,et al. Application of the EMD method in the storm surges at the coast of China[J]. Marine Science Bulletin,2016,35(2):157-162,169(in Chinese).

    [10] 李嘉康,李其杰,趙?穎,等. 基于CEEMD-BP神經(jīng)網(wǎng)絡(luò)的海溫異常預(yù)測研究[J]. 數(shù)學(xué)的實(shí)踐與認(rèn)識,2017,47(24):163-171.

    Li Jiakang,Li Qijie,Zhao Ying,et al. Sea surface temperature anomaly forecast research based on CEEMD-BP neural network[J]. Mathematics in Practice and Theory,2017,47(24):163-171(in Chinese).

    [11] Cheng Y,Ezer T,Atkinson L P,et al. Analysis of tidal amplitude changes using the EMD method[J]. Continental Shelf Research,2017,148:44-52.

    [12] Pan H,Guo Z,Wang Y,et al. Application of the EMD method to river tides[J]. Journal of Atmospheric & Oceanic Technology,2018,35:809-819.

    [13] Wu Z H,Huang N E. Ensemble empirical mode decomposition:A noise-assisted data analysis method[J]. Advances in Adaptive Data Analysis,2009,1(1):1-41.

    [14] Yeh J R,Shinh J S,Huang N E. Complementary ensemble empirical mode decomposition:A novel noise enhanced data analysis method[J]. Advances in Adaptive Data Analysis,2010,2(2):135-156.

    [15] 鄭近德,程軍圣,楊?宇. 改進(jìn)的EEMD算法及其應(yīng)用研究[J]. 振動(dòng)與沖擊,2013,32(21):21-26,46.

    Zheng Jinde,Cheng Junsheng,Yang Yu. Modified EEMD algorithm and its applications[J]. Journal of Vi-

    bration and Shock,2013,32(21):21-26,46(in Chinese).

    [16] 黃大吉,趙進(jìn)平,蘇紀(jì)蘭. 希爾伯特-黃變換的端點(diǎn)延拓[J]. 海洋學(xué)報(bào),2003,25(1):1-11.

    Huang Daji,Zhao Jinping,Su Jilan. Practical implementation of the Hilbert-Huang transform algorithm[J]. Acta Oceanologica Sinica,2003,25(1):1-11(in Chinese).

    [17] 胡愛軍,孫敬敬,向?玲. 經(jīng)驗(yàn)?zāi)B(tài)分解中的模態(tài)混疊問題[J]. 振動(dòng)、測試與診斷,2011,31(4):429-434,532-533.

    Hu Aijun,Sun Jingjing,Xiang Ling. Mode mixing in empirical mode decomposition[J]. Journal of Vibration,Measurement & Diagnosis,2011,31(4):429-434,532-533(in Chinese).

    [18] 練繼建,榮欽彪,董霄峰,等. 抑制模態(tài)混疊的HHT結(jié)構(gòu)模態(tài)參數(shù)識別方法研究[J]. 振動(dòng)與沖擊,2018,37(18):1-8.

    Lian Jijian,Rong Qinbiao,Dong Xiaofeng,et al. Structural model parameter identification method based on an improved HHT for suppressing mode mixing[J]. Journal of Vibration and Shock,2018,37(18):1-8(in Chinese).

    [19] Bandt C,Pompe B. Permutation entropy:A natural complexity measure for time series[J]. Physical Review Letters,The American Physiological Society,2002,88(17):174102-1-174102-4.

    [20] 陳?哲,李亞安. 基于多尺度排列熵的艦船輻射噪聲復(fù)雜度特征提取研究[J]. 振動(dòng)與沖擊,2019,38(12):225-230.

    Chen Zhe,Li Ya’an. Research on complexity feature extraction of ship radiated noise based on multi-scale permutation entropy[J]. Journal of Vibration and Shock,2019,38(12):225-230(in Chinese).

    Study of the Analysis of Tides Based on the HHT Method

    Zhang Liang,Zhang Jiali,Zhang Xuefeng,Zhang Anmin

    (School of Marine Science and Technology,Tianjin University,Tianjin 300072,China)

    Empirical mode decomposition(EMD)is the core component of the Hilbert-Huang transform(HHT).Based on research and analysis of the traditional HHT method,we developed a method known as modified ensemble empirical mode decomposition(MEEMD)to analyze ocean tides.With respect to the characteristics of a water-level data series,we used MEEMD method to eliminate the false intrinsic mode function(IMF)in the decomposition process by adding paired white noise and calculating the permutation entropy,thus suppressing the occurrence of mode aliasing.A comparison of the numerical simulation results with measured water-level data shows that MEEMD method can effectively suppress modal aliasing and realize adaptive filtering for original water-level sequences with noise.In this study,the high-frequency IMF1and IMF2components decomposed by MEEMD method accounted for more than 90% of the total energy of the original water-level data.The average period and amplitude of the IMF1and IMF2components were basically consistent with the results of the semidiurnal-constituent and diurnal-constituent tidal components obtained by harmonic analysis,respectively.The results show that the sum of IMF1and IMF2is the superposition of the main semidiurnal,diurnal and shallow water tidal components,respectively. Traditional tidal harmonic analysis methods rely on long-term water-level records and require some prior knowledge of the water-level data.The astronomical knowledge and mathematical calculations involved are also more complex.Quasi-harmonic analysis can process short-term water-level data,but still requires astronomical knowledge and complex mathematical calculations,as well as high data integrity.When the collected water-level data are inadequate for traditional harmonic or quasi-harmonic analysis due to noise or lack of information,the HHT method can be used for effective supplementary analysis to obtain the magnitude information of the main semidiurnal constituent,diurnal constituent,and residual water levels.Therefore,MEEMD-HHT method can be used as a good supplementary method for tidal analysis.

    modified ensemble empirical mode decomposition(MEEMD);Hilbert-Huang transform(HHT);filtering;tidal constituent

    P731.23

    A

    0493-2137(2020)07-0725-11

    10.11784/tdxbz201908071

    2019-08-30;

    2019-10-07.

    張?亮(1983—??),男,博士,講師,liangzhang@tju.edu.cn.

    張學(xué)峰,xfz_nmdis@126.com.

    國家自然科學(xué)基金資助項(xiàng)目(41876014,11801402).

    Supported by the National Natural Science Foundation of China(No. 41876014,No. 11801402).

    (責(zé)任編輯:金順愛)

    猜你喜歡
    模態(tài)分析方法
    隱蔽失效適航要求符合性驗(yàn)證分析
    電力系統(tǒng)不平衡分析
    電子制作(2018年18期)2018-11-14 01:48:24
    電力系統(tǒng)及其自動(dòng)化發(fā)展趨勢分析
    可能是方法不對
    用對方法才能瘦
    Coco薇(2016年2期)2016-03-22 02:42:52
    國內(nèi)多模態(tài)教學(xué)研究回顧與展望
    四大方法 教你不再“坐以待病”!
    Coco薇(2015年1期)2015-08-13 02:47:34
    捕魚
    基于HHT和Prony算法的電力系統(tǒng)低頻振蕩模態(tài)識別
    由單個(gè)模態(tài)構(gòu)造對稱簡支梁的抗彎剛度
    精品人妻1区二区| 侵犯人妻中文字幕一二三四区| 色婷婷久久久亚洲欧美| 桃红色精品国产亚洲av| 国产91精品成人一区二区三区 | h视频一区二区三区| 国产成人免费无遮挡视频| 国产熟女午夜一区二区三区| 欧美国产精品va在线观看不卡| 老司机福利观看| 99香蕉大伊视频| 波多野结衣一区麻豆| 久久久水蜜桃国产精品网| 日韩熟女老妇一区二区性免费视频| 国产精品99久久99久久久不卡| 搡老熟女国产l中国老女人| 免费人妻精品一区二区三区视频| 麻豆av在线久日| 久久精品成人免费网站| 国产精品 国内视频| 99国产精品一区二区蜜桃av | 国产成人系列免费观看| 自线自在国产av| 黄色 视频免费看| 每晚都被弄得嗷嗷叫到高潮| 黄色 视频免费看| 日韩中文字幕欧美一区二区| 欧美日韩黄片免| 亚洲精品中文字幕在线视频| 高潮久久久久久久久久久不卡| 一本—道久久a久久精品蜜桃钙片| 国产成人免费观看mmmm| 99久久国产精品久久久| 极品人妻少妇av视频| 日韩一卡2卡3卡4卡2021年| 一区在线观看完整版| 视频在线观看一区二区三区| 久久精品91无色码中文字幕| 最近最新免费中文字幕在线| 99国产精品一区二区蜜桃av | 精品久久久久久久毛片微露脸| 人人妻人人爽人人添夜夜欢视频| 中文字幕人妻丝袜一区二区| 亚洲久久久国产精品| 91精品三级在线观看| 亚洲视频免费观看视频| 国产伦理片在线播放av一区| 精品国产一区二区三区四区第35| 国产在线观看jvid| av欧美777| 亚洲成人国产一区在线观看| 久久天堂一区二区三区四区| 欧美黑人精品巨大| 亚洲精品一二三| 国产不卡av网站在线观看| 亚洲视频免费观看视频| 午夜福利一区二区在线看| 麻豆乱淫一区二区| 老熟妇乱子伦视频在线观看| 国产精品1区2区在线观看. | 亚洲成国产人片在线观看| 777米奇影视久久| 性少妇av在线| 视频区图区小说| 精品久久蜜臀av无| 多毛熟女@视频| 正在播放国产对白刺激| 正在播放国产对白刺激| 我的亚洲天堂| 男男h啪啪无遮挡| 9191精品国产免费久久| 黑人操中国人逼视频| 五月天丁香电影| 亚洲色图av天堂| 麻豆乱淫一区二区| 18禁观看日本| 日韩制服丝袜自拍偷拍| 色在线成人网| 欧美变态另类bdsm刘玥| 18禁国产床啪视频网站| 日韩三级视频一区二区三区| 国产在线免费精品| 国产精品一区二区在线观看99| 99在线人妻在线中文字幕 | √禁漫天堂资源中文www| 免费久久久久久久精品成人欧美视频| 国产亚洲午夜精品一区二区久久| 伦理电影免费视频| 国产精品免费大片| 可以免费在线观看a视频的电影网站| 2018国产大陆天天弄谢| 老司机午夜福利在线观看视频 | 9191精品国产免费久久| 国产精品久久久久成人av| 国产亚洲精品一区二区www | 两性午夜刺激爽爽歪歪视频在线观看 | 国产不卡一卡二| 搡老岳熟女国产| 午夜福利乱码中文字幕| av欧美777| 在线永久观看黄色视频| 国产精品免费一区二区三区在线 | 亚洲色图av天堂| 亚洲成人免费电影在线观看| 午夜视频精品福利| 丁香六月天网| 国产淫语在线视频| 精品人妻1区二区| 80岁老熟妇乱子伦牲交| 天堂动漫精品| 男女下面插进去视频免费观看| 日本av免费视频播放| 国产黄频视频在线观看| 电影成人av| 成人黄色视频免费在线看| 国产精品秋霞免费鲁丝片| 欧美日韩精品网址| 一本久久精品| 国产视频一区二区在线看| 国产亚洲av高清不卡| 韩国精品一区二区三区| 男女无遮挡免费网站观看| 国产在线视频一区二区| 日本a在线网址| 桃红色精品国产亚洲av| 国产一区二区在线观看av| 国产人伦9x9x在线观看| 亚洲成人手机| 757午夜福利合集在线观看| 亚洲欧美激情在线| 无遮挡黄片免费观看| 亚洲,欧美精品.| 90打野战视频偷拍视频| 这个男人来自地球电影免费观看| 热99re8久久精品国产| 国产精品久久电影中文字幕 | 三级毛片av免费| 正在播放国产对白刺激| 丝袜美足系列| 黄网站色视频无遮挡免费观看| 99re在线观看精品视频| 亚洲人成77777在线视频| 黄片大片在线免费观看| 在线av久久热| 色视频在线一区二区三区| 午夜福利影视在线免费观看| 侵犯人妻中文字幕一二三四区| 18禁观看日本| 一本一本久久a久久精品综合妖精| 午夜福利免费观看在线| 久久天躁狠狠躁夜夜2o2o| 午夜老司机福利片| 最新美女视频免费是黄的| 天天躁日日躁夜夜躁夜夜| 999久久久精品免费观看国产| 丝袜喷水一区| 人人妻人人澡人人爽人人夜夜| a级毛片在线看网站| 国产97色在线日韩免费| 脱女人内裤的视频| 九色亚洲精品在线播放| videos熟女内射| 91精品三级在线观看| 视频在线观看一区二区三区| 日韩中文字幕视频在线看片| 国产成人av激情在线播放| 亚洲欧美精品综合一区二区三区| 无遮挡黄片免费观看| 日韩三级视频一区二区三区| 久久久久久久精品吃奶| 中文字幕色久视频| 亚洲 欧美一区二区三区| 少妇 在线观看| av视频免费观看在线观看| 免费观看av网站的网址| av福利片在线| 高清在线国产一区| 国产在视频线精品| 每晚都被弄得嗷嗷叫到高潮| 久久精品国产亚洲av香蕉五月 | 精品亚洲乱码少妇综合久久| 久久av网站| 久久亚洲真实| 欧美人与性动交α欧美精品济南到| 国产成人免费观看mmmm| 免费在线观看日本一区| 下体分泌物呈黄色| 亚洲少妇的诱惑av| 亚洲中文av在线| 久久九九热精品免费| 国产精品 欧美亚洲| av电影中文网址| 黄色 视频免费看| 午夜福利一区二区在线看| 国产男女超爽视频在线观看| 如日韩欧美国产精品一区二区三区| 国产97色在线日韩免费| 免费一级毛片在线播放高清视频 | 狠狠精品人妻久久久久久综合| a级片在线免费高清观看视频| 国产一区二区激情短视频| 久久久久久久久久久久大奶| 视频区图区小说| 黄片播放在线免费| 亚洲精品一卡2卡三卡4卡5卡| 免费少妇av软件| 国产亚洲精品久久久久5区| 亚洲国产av新网站| 久久热在线av| 麻豆成人av在线观看| 国产又爽黄色视频| 国产精品一区二区在线观看99| 免费观看a级毛片全部| 国产aⅴ精品一区二区三区波| 丰满少妇做爰视频| 亚洲国产成人一精品久久久| 久久中文看片网| 亚洲精华国产精华精| 国产97色在线日韩免费| 亚洲av片天天在线观看| 飞空精品影院首页| 人妻一区二区av| 国产黄色免费在线视频| 国产亚洲欧美精品永久| 欧美人与性动交α欧美软件| 黑人巨大精品欧美一区二区蜜桃| 成年人免费黄色播放视频| 汤姆久久久久久久影院中文字幕| 亚洲综合色网址| 女人爽到高潮嗷嗷叫在线视频| 自线自在国产av| 国产在线免费精品| 欧美日韩福利视频一区二区| 国精品久久久久久国模美| 每晚都被弄得嗷嗷叫到高潮| 日韩三级视频一区二区三区| 性少妇av在线| 老熟妇乱子伦视频在线观看| 日韩欧美一区视频在线观看| 亚洲精品自拍成人| 最近最新中文字幕大全电影3 | 18禁美女被吸乳视频| 欧美乱妇无乱码| 国产精品免费大片| 老司机在亚洲福利影院| 国产精品香港三级国产av潘金莲| 18禁美女被吸乳视频| 两人在一起打扑克的视频| 亚洲中文日韩欧美视频| 老司机靠b影院| www.999成人在线观看| 美女主播在线视频| 黄色怎么调成土黄色| 妹子高潮喷水视频| 国产欧美日韩精品亚洲av| 男人操女人黄网站| 久久影院123| 日韩一区二区三区影片| 一区二区日韩欧美中文字幕| 久久久久久久精品吃奶| 久久久久视频综合| 黑人操中国人逼视频| 伊人久久大香线蕉亚洲五| 性高湖久久久久久久久免费观看| 久9热在线精品视频| 淫妇啪啪啪对白视频| 老汉色∧v一级毛片| 一二三四社区在线视频社区8| 久久九九热精品免费| 人人妻人人添人人爽欧美一区卜| 少妇 在线观看| av视频免费观看在线观看| a级毛片在线看网站| 色尼玛亚洲综合影院| 中文亚洲av片在线观看爽 | 精品人妻在线不人妻| 欧美精品一区二区大全| 国产激情久久老熟女| 欧美精品亚洲一区二区| 少妇被粗大的猛进出69影院| av超薄肉色丝袜交足视频| 黑人猛操日本美女一级片| 国产亚洲一区二区精品| www.999成人在线观看| 久久久国产成人免费| 美女福利国产在线| av线在线观看网站| 国产精品自产拍在线观看55亚洲 | 亚洲国产av新网站| av不卡在线播放| 亚洲国产毛片av蜜桃av| 91av网站免费观看| 国产精品免费大片| 亚洲欧美激情在线| 捣出白浆h1v1| 视频区图区小说| 精品国产国语对白av| 丁香六月天网| 黑丝袜美女国产一区| 久久久久国产一级毛片高清牌| 热99国产精品久久久久久7| 欧美在线黄色| 视频在线观看一区二区三区| 两个人看的免费小视频| 亚洲精品国产区一区二| 国产成人一区二区三区免费视频网站| 黄片播放在线免费| 久久久久久久大尺度免费视频| 黄片小视频在线播放| 制服诱惑二区| 亚洲中文日韩欧美视频| 午夜福利视频在线观看免费| 狠狠狠狠99中文字幕| 在线播放国产精品三级| 精品国产亚洲在线| 老汉色av国产亚洲站长工具| 国产国语露脸激情在线看| 亚洲国产欧美网| 国产在线观看jvid| 久久天躁狠狠躁夜夜2o2o| 黑人欧美特级aaaaaa片| 久久婷婷成人综合色麻豆| 美女高潮到喷水免费观看| 国产区一区二久久| 97人妻天天添夜夜摸| 精品少妇黑人巨大在线播放| 日本wwww免费看| 老司机深夜福利视频在线观看| 国产精品麻豆人妻色哟哟久久| cao死你这个sao货| 黄色丝袜av网址大全| 婷婷成人精品国产| 天天添夜夜摸| 国产单亲对白刺激| 在线观看66精品国产| 人妻 亚洲 视频| 欧美精品一区二区免费开放| 人人澡人人妻人| 成人av一区二区三区在线看| 久久国产精品大桥未久av| 777米奇影视久久| 777久久人妻少妇嫩草av网站| 两个人免费观看高清视频| 亚洲熟女精品中文字幕| 精品久久久精品久久久| 欧美黑人精品巨大| 制服诱惑二区| 99久久99久久久精品蜜桃| 国产在线观看jvid| 一级,二级,三级黄色视频| 窝窝影院91人妻| videos熟女内射| 亚洲人成电影观看| 一级毛片精品| 亚洲少妇的诱惑av| 日本精品一区二区三区蜜桃| 欧美日韩一级在线毛片| 无限看片的www在线观看| 久久精品国产亚洲av香蕉五月 | 亚洲国产毛片av蜜桃av| 一级片免费观看大全| 久久久久国产一级毛片高清牌| 亚洲国产看品久久| 人妻一区二区av| 超色免费av| 在线观看免费高清a一片| 国产精品麻豆人妻色哟哟久久| 成人手机av| 色综合婷婷激情| 日韩成人在线观看一区二区三区| 精品福利永久在线观看| 久久青草综合色| 久久亚洲精品不卡| 99精品久久久久人妻精品| 超碰成人久久| 欧美黄色片欧美黄色片| 一区二区av电影网| 成人特级黄色片久久久久久久 | 男女之事视频高清在线观看| 成人18禁在线播放| 黑人巨大精品欧美一区二区蜜桃| 夜夜爽天天搞| 国产1区2区3区精品| 在线观看一区二区三区激情| 亚洲精华国产精华精| 日日爽夜夜爽网站| 国产精品麻豆人妻色哟哟久久| 麻豆av在线久日| 精品少妇内射三级| 午夜福利欧美成人| 国产精品美女特级片免费视频播放器 | 久久精品aⅴ一区二区三区四区| 免费观看a级毛片全部| 我要看黄色一级片免费的| 亚洲久久久国产精品| e午夜精品久久久久久久| 欧美精品亚洲一区二区| 啦啦啦 在线观看视频| 国产激情久久老熟女| 亚洲专区中文字幕在线| 麻豆av在线久日| 中文字幕高清在线视频| 少妇被粗大的猛进出69影院| 国产成人av教育| 一级黄色大片毛片| 午夜福利在线免费观看网站| 亚洲国产中文字幕在线视频| 久久久久久人人人人人| 一级片免费观看大全| 50天的宝宝边吃奶边哭怎么回事| 精品一区二区三区视频在线观看免费 | 国产精品自产拍在线观看55亚洲 | 窝窝影院91人妻| 久久av网站| 婷婷丁香在线五月| 91字幕亚洲| 最新的欧美精品一区二区| 国产亚洲精品一区二区www | 日本黄色日本黄色录像| 狂野欧美激情性xxxx| 欧美日韩黄片免| 啦啦啦 在线观看视频| 欧美日韩一级在线毛片| 国产精品1区2区在线观看. | 淫妇啪啪啪对白视频| 日韩成人在线观看一区二区三区| 国产精品成人在线| 免费在线观看视频国产中文字幕亚洲| 中文字幕人妻丝袜一区二区| 夜夜夜夜夜久久久久| 国产成人欧美| 99国产精品一区二区蜜桃av | 亚洲精品av麻豆狂野| av在线播放免费不卡| 日本一区二区免费在线视频| 欧美日韩精品网址| 午夜福利在线免费观看网站| 久久久久久久久免费视频了| 国产在线观看jvid| 午夜精品久久久久久毛片777| 黄色视频在线播放观看不卡| 亚洲av日韩精品久久久久久密| 一级毛片电影观看| 97在线人人人人妻| 一区在线观看完整版| 老熟妇乱子伦视频在线观看| 90打野战视频偷拍视频| 丰满人妻熟妇乱又伦精品不卡| 嫩草影视91久久| 欧美精品一区二区免费开放| 国产精品熟女久久久久浪| www.999成人在线观看| svipshipincom国产片| 精品人妻1区二区| 国产精品成人在线| 9热在线视频观看99| 91av网站免费观看| 王馨瑶露胸无遮挡在线观看| 黄片小视频在线播放| 欧美日韩亚洲国产一区二区在线观看 | 午夜福利视频精品| 大片免费播放器 马上看| 国产一卡二卡三卡精品| 日韩中文字幕视频在线看片| 免费在线观看影片大全网站| 国产精品 国内视频| 超碰成人久久| 成人亚洲精品一区在线观看| 国产精品 欧美亚洲| 男女午夜视频在线观看| 欧美乱妇无乱码| 一边摸一边抽搐一进一出视频| 大陆偷拍与自拍| 91大片在线观看| 99re在线观看精品视频| 亚洲国产av影院在线观看| 欧美激情久久久久久爽电影 | 日韩人妻精品一区2区三区| 99九九在线精品视频| 亚洲精品粉嫩美女一区| 久久午夜综合久久蜜桃| h视频一区二区三区| 亚洲第一欧美日韩一区二区三区 | 极品教师在线免费播放| 亚洲 欧美一区二区三区| 91麻豆精品激情在线观看国产 | svipshipincom国产片| 丁香欧美五月| 国产成人一区二区三区免费视频网站| 久久久精品区二区三区| 日韩人妻精品一区2区三区| av不卡在线播放| 中文字幕另类日韩欧美亚洲嫩草| 久久 成人 亚洲| 天堂动漫精品| 亚洲精品国产色婷婷电影| 亚洲第一av免费看| 国产一区二区在线观看av| 午夜视频精品福利| 日韩大码丰满熟妇| 97人妻天天添夜夜摸| 香蕉丝袜av| 亚洲国产欧美网| 国产欧美亚洲国产| 亚洲情色 制服丝袜| 国产91精品成人一区二区三区 | 国产xxxxx性猛交| 亚洲人成伊人成综合网2020| 免费在线观看影片大全网站| 欧美亚洲 丝袜 人妻 在线| 国产精品免费视频内射| 亚洲精品粉嫩美女一区| 我的亚洲天堂| 美女午夜性视频免费| 成人国产av品久久久| 久久久国产精品麻豆| 免费在线观看日本一区| 欧美黄色片欧美黄色片| 在线播放国产精品三级| avwww免费| av网站在线播放免费| 色在线成人网| 最新在线观看一区二区三区| 十八禁高潮呻吟视频| 国产麻豆69| 极品少妇高潮喷水抽搐| 久久久国产成人免费| 悠悠久久av| 精品熟女少妇八av免费久了| 不卡一级毛片| 国产色视频综合| 久久久精品区二区三区| 汤姆久久久久久久影院中文字幕| 亚洲欧美精品综合一区二区三区| 考比视频在线观看| 久久免费观看电影| 狂野欧美激情性xxxx| 午夜免费鲁丝| 免费在线观看视频国产中文字幕亚洲| 亚洲三区欧美一区| 国产精品久久久av美女十八| 亚洲成人国产一区在线观看| 性高湖久久久久久久久免费观看| 嫁个100分男人电影在线观看| 国产单亲对白刺激| 一级毛片女人18水好多| 日韩欧美三级三区| 午夜激情久久久久久久| 欧美在线黄色| 夫妻午夜视频| 久久中文看片网| 怎么达到女性高潮| 天堂动漫精品| 国精品久久久久久国模美| 成年人免费黄色播放视频| 啦啦啦在线免费观看视频4| www.999成人在线观看| 日日爽夜夜爽网站| 国产亚洲欧美在线一区二区| 热re99久久国产66热| 丝袜美腿诱惑在线| 老司机深夜福利视频在线观看| 狠狠婷婷综合久久久久久88av| 我要看黄色一级片免费的| 免费观看人在逋| 99国产精品一区二区三区| 男女下面插进去视频免费观看| a级毛片黄视频| 色精品久久人妻99蜜桃| 激情视频va一区二区三区| 中文欧美无线码| 少妇粗大呻吟视频| 国产国语露脸激情在线看| 十八禁高潮呻吟视频| 久久精品国产99精品国产亚洲性色 | 欧美大码av| 91精品三级在线观看| 中文字幕高清在线视频| 成在线人永久免费视频| 国产精品熟女久久久久浪| 老鸭窝网址在线观看| 亚洲成人免费电影在线观看| 日本av免费视频播放| 老熟妇仑乱视频hdxx| 精品国产一区二区久久| 十八禁高潮呻吟视频| 亚洲第一欧美日韩一区二区三区 | 久久精品91无色码中文字幕| 亚洲三区欧美一区| 欧美人与性动交α欧美精品济南到| 午夜精品国产一区二区电影| 国产欧美日韩一区二区精品| 丝瓜视频免费看黄片| 99热网站在线观看| 亚洲国产看品久久| 男女高潮啪啪啪动态图| 国产精品二区激情视频| 9191精品国产免费久久| 免费观看av网站的网址| 激情在线观看视频在线高清 | 69精品国产乱码久久久| www.自偷自拍.com| 老司机福利观看| a级片在线免费高清观看视频| 久久影院123| 人人妻人人澡人人爽人人夜夜| 亚洲自偷自拍图片 自拍| 亚洲黑人精品在线| 999精品在线视频| 99久久精品国产亚洲精品| av欧美777| 久久久久久久久免费视频了| 国产亚洲欧美在线一区二区| 国产精品久久久久成人av| 久久免费观看电影|