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

    基于結(jié)構(gòu)自適應(yīng)中值濾波器的隨機噪聲衰減方法

    2012-12-18 05:28:48高靜懷陳文超朱振宇
    地球物理學報 2012年5期
    關(guān)鍵詞:同相軸置信度剖面

    王 偉,高靜懷,陳文超*,朱振宇

    1 西安交通大學電子與信息工程學院波動與信息研究所,西安 710049

    2 中海石油研究總院,北京 100027

    基于結(jié)構(gòu)自適應(yīng)中值濾波器的隨機噪聲衰減方法

    王 偉1,高靜懷1,陳文超1*,朱振宇2

    1 西安交通大學電子與信息工程學院波動與信息研究所,西安 710049

    2 中海石油研究總院,北京 100027

    本文提出一種保護斷層、裂縫等地層邊緣特征的結(jié)構(gòu)自適應(yīng)中值濾波器,用于衰減地震資料中的隨機噪聲.基于地震反射同相軸局部呈線型結(jié)構(gòu)的假設(shè),采用梯度結(jié)構(gòu)張量估計地層傾向,分析地層結(jié)構(gòu)的規(guī)則程度,在此基礎(chǔ)上引入地震剖面中線型和橫向不連續(xù)性兩種結(jié)構(gòu)特征的置信度量.結(jié)構(gòu)自適應(yīng)中值濾波器根據(jù)這兩種置信度量調(diào)整濾波器窗函數(shù)的尺度和形狀,根據(jù)地層傾角調(diào)整濾波器窗函數(shù)的方向,從而使得濾波操作窗能夠最佳匹配信號的局部結(jié)構(gòu)特征.將本文方法用于合成和實際數(shù)據(jù)的處理,并與兩種常用中值濾波方法進行對比,結(jié)果表明,該方法能夠更好地解決地震剖面的隨機噪聲衰減和有效信號保真的問題,在增強反射同相軸的橫向一致性的同時有效保持了剖面內(nèi)的地層邊緣和細節(jié)特征,顯著改善了地震資料的品質(zhì).

    中值濾波,自適應(yīng)濾波,噪聲衰減,構(gòu)造保護,梯度結(jié)構(gòu)張量

    1 引 言

    在石油天然氣勘探開發(fā)中,斷層、通道和裂縫等結(jié)構(gòu)特征以及河道砂體等沉積特征是發(fā)現(xiàn)和描述油藏的基礎(chǔ),分辨和分析地震資料中這些結(jié)構(gòu)特征具有重要的意義.然而,由于這些區(qū)域地質(zhì)形態(tài)的復雜性,使得地震數(shù)據(jù)中這些區(qū)域易受噪聲的干擾,無論采用人工解釋還是相干和邊緣檢測等自動檢測算法來描述和刻畫這些區(qū)域都變得較為困難.傳統(tǒng)的地震資料去噪方法,如F-K域傾向濾波、f-x域預測濾波和多道相干濾波等,多是利用地震反射同相軸的橫向一致性,易造成對傾斜同相軸振幅的壓制,也會模糊小的斷層和裂縫等不連續(xù)結(jié)構(gòu),甚至會引起大的斷層兩側(cè)同相軸的錯誤連接.如何能夠較好地平衡地震噪聲衰減和有效信號保護的關(guān)系受到了廣泛的關(guān)注[1-6],而能夠保護和增強斷裂等邊緣結(jié)構(gòu)的疊后噪聲衰減技術(shù)尤其受到地震解釋工作者的重視[7-9].

    Luo等[10]基于次序統(tǒng)計的思想給出一種地震數(shù)據(jù)保邊濾波(EPS,Edge-Preserving Smoothing)算法,該方法采用Kuwahara多窗分析技術(shù),計算當前分析點周圍各個子窗內(nèi)數(shù)據(jù)的均值和方差,將最小方差對應(yīng)的子窗的均值作為當前點的濾波輸出.AlBinHassan等[11]發(fā)展了這種保邊濾波思想,給出了三維情況下的子窗剖分方法,將EPS算法用于三維地震數(shù)據(jù)的保邊濾波處理.Lu等[12]基于多項式擬合的思想推廣了一維的EPS算法,該方法通過計算包含當前分析點的一系列平移子窗內(nèi)信號對給定階次的多項式擬合誤差,選擇出最小擬合誤差對應(yīng)的子窗的信號估計值作為當前點的濾波輸出.楊培杰等[13]也基于這種多窗分析思想給出一種方向性邊界保持斷層增強技術(shù),使得相干體中的斷層圖像得到更加清晰顯示.Liu等[14-15]應(yīng)用圖像處理中的二維多級中值濾波器消除地震隨機噪聲,并分析了濾波器尺度參數(shù)選擇對噪聲衰減和信號細節(jié)保護的影響.此外,F(xiàn)ehmers等[16]引入微分方程濾波法用于地震資料的保邊濾波處理,該方法采用梯度結(jié)構(gòu)張量估計地層方向,度量地層的連續(xù)性,以此約束擴散濾波方程的濾波方向和濾波程度.Lavialle等[17]結(jié)合線型和面型擴散濾波模型給出專門針對保護和增強斷層構(gòu)造的斷層保持擴散濾波器.孫夕平等[18]和王旭松等[19]通過在擴散方程中引入二階導數(shù)項、張爾華等[20]通過修正擴散方程中的保邊約束項進一步改善非線性各向異性擴散濾波器的邊緣保護性能.然而,基于各向異性擴散方程的濾波方法以原始地震數(shù)據(jù)作為初始條件,通過擴散時間的迭代推移形成濾波作用,由于不好估計出獲得最佳濾波效果的擴散時間,導致不能準確預測其濾波性態(tài)[21].

    經(jīng)典的二維中值濾波法具有優(yōu)異的噪聲衰減性能,尤其適用于消除非平穩(wěn)信號中的峰值噪聲.但是,由于在濾波過程中其濾波窗口保持恒定,且濾波操作沒有方向特性,使得采用中值濾波法處理具有層狀結(jié)構(gòu)的地震信號時,必然會引起有用信號的損傷和不連續(xù)結(jié)構(gòu)特征的模糊.一些改進的中值濾波方法如二維多級中值濾波器[14-15],雖然有較強的信號細節(jié)保護能力,但是其對噪聲的衰減能力不夠.本文結(jié)合微分方程濾波法中利用梯度結(jié)構(gòu)張量刻畫地震反射層局部結(jié)構(gòu)特征的方法,定義了帶有能量變化強度信息的橫向不連續(xù)性置信度量,結(jié)合線型置信度量和橫向不連續(xù)性置信度量,給出一種二維結(jié)構(gòu)自適應(yīng)中值濾波方法,使中值濾波操作窗函數(shù)根據(jù)信號局部結(jié)構(gòu)特征自適應(yīng)地作方向延伸和尺度伸縮,從而既能最大限度地消除隨機噪聲,又能有效保護有用信號和地層邊緣等細節(jié)特征.

    2 結(jié)構(gòu)自適應(yīng)中值濾波原理

    2.1 地層傾向估計及橫向不連續(xù)性度量

    有多種方法可用于地層結(jié)構(gòu)的規(guī)則性分析[16,22-25],梯度結(jié)構(gòu)張量法是其中一種簡單而又有效的途徑.記二維地震剖面為u(x,t),則該二維剖面對應(yīng)的梯度結(jié)構(gòu)張量由下式描述

    其中,和T分別為梯度算子和轉(zhuǎn)置算子,?表示卷積運算,Gρ表示尺度為ρ的二維高斯函數(shù)Gρ(x,t)=exp(-(x2+t2)/2ρ2).在梯度結(jié)構(gòu)張量的定義中,由梯度向量Δu及其轉(zhuǎn)置向量形成的張量矩陣與尺度為ρ的高斯低通濾波器作用,確定了梯度結(jié)構(gòu)張量可度量的信號特征的最大尺度.將該對稱半正定矩陣Sρ(Δu)作矩陣特征分解

    其中,特征值μ1≥μ2≥0反應(yīng)了信號在位置(x,t)的高斯鄰域Gρ內(nèi)沿特征方向的振幅變化強度;特征向量v1和v2給出局部正交方位,v1對應(yīng)于局部振幅變化最大的方向,即信號梯度方向,而v2對應(yīng)于局部振幅變化最小的方向,即反射同相軸的取向.根據(jù)梯度結(jié)構(gòu)張量兩個特征值的取值情況,可分辨不同的信號結(jié)構(gòu)特征,文獻[26]中引入線型信號結(jié)構(gòu)的置信度量:

    該置信度量在區(qū)間[0,1]之間取值,衡量局部信號特征與線型結(jié)構(gòu)的相似程度,同時亦表征局部方位估計的準確程度.若CL→1,則表示局部信號特征趨于線型結(jié)構(gòu),估計的地層傾角的準確度越高;反之,若CL→0,則局部信號結(jié)構(gòu)背離了線型假設(shè),這可能由干擾噪聲和斷層、裂縫等橫向不連續(xù)性造成.根據(jù)線型結(jié)構(gòu)置信度量可較好地識別局部呈線型結(jié)構(gòu)的反射同相軸;但是,由于該置信度量值判斷的是信號能量的相對變化程度,容易受到噪聲的干擾,并且在無明顯信號結(jié)構(gòu)的區(qū)域內(nèi)取值沒有意義,從而很難從干擾背景中明確地識別斷層、裂縫等橫向不連續(xù)結(jié)構(gòu).

    基于上述分析,本文在線型結(jié)構(gòu)置信度量的基礎(chǔ)上,給出一種含有信號橫向能量變化強度信息的置信度量值,用于識別斷層、裂縫等橫向不連續(xù)結(jié)構(gòu).

    (4)式中(1-CL)項表示信號局部結(jié)構(gòu)特征相對線型結(jié)構(gòu)的背離程度,μ2表征在給定度量鄰域內(nèi)在均方誤差最小意義下信號沿局部一致方向的能量變化強度[24].在反射同相軸呈線型結(jié)構(gòu)的區(qū)域,由于CL→1,即(1-CL)→0,此時沿地層傾向的信號變化相對較小,即μ2→0,從而使得CI→0.在含有斷層、裂縫等橫向不連續(xù)性的區(qū)域,由于CL→0,即(1-CL)→1;且在橫向不連續(xù)區(qū)域,由于信號結(jié)構(gòu)背離了局部線型結(jié)構(gòu)假設(shè),此處信號的一致性取向不明確,因而沿局部一致性方向的信號變化強度也相對較大,即μ2?0,從而使得CI?0.因此,根據(jù)橫向不連續(xù)性置信度量CI可以很好地識別斷層和裂縫等橫向不連續(xù)結(jié)構(gòu).由于信號橫向能量變化項μ2的存在,對于噪聲能量弱于信號能量的情況,CI具有良好的抗噪能力;此外,CI也能夠判別無明顯信號結(jié)構(gòu)的區(qū)域(此時μ2→0,即CI→0),這對采用CI指導濾波過程尤為重要,可以避免一些濾波假象的產(chǎn)生.

    2.2 結(jié)構(gòu)自適應(yīng)中值濾波器

    對于二維地震剖面u(x),記x=(x,t)為采樣點的空間、時間位置,定義位置x處的結(jié)構(gòu)自適應(yīng)中值濾波操作為:

    式中,M(x,y)表示位置x處的結(jié)構(gòu)自適應(yīng)濾波窗口(見圖1a),y為當前濾波窗口包含的采樣點的空間、時間位置,(x)表示位置x處的結(jié)構(gòu)自適應(yīng)中值濾波的輸出結(jié)果.結(jié)構(gòu)自適應(yīng)濾波窗M(x,y)是濾波位置x的函數(shù),其大小和形狀取決于x處分析鄰域內(nèi)的信號結(jié)構(gòu)特征:

    其中,·表示內(nèi)積運算;n(x)和n⊥(x)分別對應(yīng)于x處的地層傾向及梯度方向,取為梯度結(jié)構(gòu)張量的兩個特征向量,即n(x)=v2(x),n⊥(x)=v1(x);σ1(x)和σ2(x)分別表示橢圓形濾波窗口的長軸和短軸,它們共同決定了結(jié)構(gòu)自適應(yīng)中值濾波器的濾波尺度以及濾波操作的方向選擇性.濾波器尺度參數(shù)的設(shè)計是影響結(jié)構(gòu)自適應(yīng)中值濾波器性態(tài)的關(guān)鍵,根據(jù)結(jié)構(gòu)自適應(yīng)濾波操作的要求,本文結(jié)合線型結(jié)構(gòu)置信度量CL和橫向不連續(xù)性置信度量CI,給出如下的濾波器尺度參數(shù)自適應(yīng)選擇策略:

    式中Rmax規(guī)定了橢圓形濾波窗口的最大尺寸,g(·)為關(guān)于CI(x)的單調(diào)減函數(shù),其取值范圍限定為(0,1],文中取g(·)為指數(shù)函數(shù)

    其中,β為針對CI(x)的閾值參數(shù),控制指數(shù)函數(shù)的衰減速度,β取值越小,指數(shù)函數(shù)的下降越快;反之,指數(shù)函數(shù)的衰減變緩.

    結(jié)構(gòu)自適應(yīng)中值濾波器根據(jù)信號的局部結(jié)構(gòu)特征具有靈活的濾波尺度調(diào)節(jié)功能和方向選擇特性,控制最大濾波尺度的σ1(x)由橫向不連續(xù)性置信度量CI(x)決定,而調(diào)整方向選擇特性的σ2(x)由線型結(jié)構(gòu)置信度量CL(x)決定.在斷層、裂縫等橫向不連續(xù)區(qū)域由于CI(x)?0,使得橢圓型窗函數(shù)的長軸σ1(x)趨短;而在其他區(qū)域隨CI(x)→0,有σ1(x)→Rmax.在反射同相軸呈線型結(jié)構(gòu)區(qū)域,由于CL(x)→1,使得橢圓型窗函數(shù)的短軸σ2(x)萎縮;而在其他區(qū)域隨CL(x)→0,有σ2(x)→σ1(x).因此,如圖1b所示,在反射同相軸局部呈線型延伸的區(qū)域,濾波窗口沿同相軸取向拉伸為狹長的各向異性窗口(位置x1區(qū)域);在斷層等橫向不連續(xù)區(qū)域,濾波窗口自適應(yīng)地萎縮為小窗以匹配此處的不連續(xù)結(jié)構(gòu)特征(位置x2區(qū)域);而在無明顯信號結(jié)構(gòu)的噪聲干擾區(qū)域,濾波窗口延展為大的近似各向同性窗口(位置x3區(qū)域).

    圖1 結(jié)構(gòu)自適應(yīng)中值濾波器窗口的構(gòu)造(a)濾波器窗口函數(shù);(b)信號結(jié)構(gòu)自適應(yīng)濾波策略.Fig.1 Filter window construction of structure-adaptive median filter(a)Filter window function;(b)Signal structure-consistent filtering mechanism.

    聯(lián)合式(7)和(8)可見,結(jié)構(gòu)自適應(yīng)中值濾波器的濾波性能受到參數(shù)Rmax和β共同約束,最大尺寸參數(shù)Rmax控制濾波器的平滑性能,閾值參數(shù)β則調(diào)整濾波器對斷層等邊緣結(jié)構(gòu)的保持程度.考慮到文中采用梯度結(jié)構(gòu)張量估計地層的傾角和規(guī)則程度,梯度結(jié)構(gòu)張量的尺度參數(shù)ρ決定了其對局部信號結(jié)構(gòu)的度量孔徑,因此為了確?;趦煞N置信度量約束的濾波過程可信,要求濾波窗口的最大尺寸Rmax不超過2ρ,即取Rmax≤2ρ.對于閾值參數(shù)β,由圖2可見,當β取值偏小時,稍大的橫向不連續(xù)性置信度量CI(x)即引起濾波窗口的尺寸減小,因而會使一些弱的斷層、裂縫等邊緣結(jié)構(gòu)得到有效保護,但同時也會減弱濾波器在橫向不連續(xù)性區(qū)域的濾波效果;當β取值偏大時,只有很大的CI(x)才會導致濾波窗口萎縮到足夠小,因而在這種情況下濾波器的噪聲衰減能力得到增強,但在一定程度上會造成弱邊緣結(jié)構(gòu)的模糊和壓制.

    圖2 結(jié)構(gòu)自適應(yīng)中值濾波窗的長軸控制函數(shù)Fig.2 Controlling functions of the main axis of the structure-adaptive median filter kernel

    2.3 保邊性態(tài)自適應(yīng)控制方法

    結(jié)構(gòu)自適應(yīng)中值濾波器的邊緣結(jié)構(gòu)保持性能受閾值參數(shù)β控制,選擇合適的β是保證濾波過程中能夠有效保持斷層、裂縫等地層邊緣的關(guān)鍵.然而,地震資料通常從淺層到深層,即使在同一層段內(nèi)信號橫向能量變化的動態(tài)范圍波動也很大,很難統(tǒng)一地選擇一個全局閾值參數(shù)β.本文給出一種區(qū)域分割的閾值參數(shù)自適應(yīng)選擇策略,即根據(jù)待處理剖面內(nèi)橫向能量變化在時間和空間上的分布情況,將剖面Ω分割為若干具有一定邊界重疊的子區(qū)域Ω=∪Ωi,在每個子區(qū)域Ωi內(nèi)取橫向不連續(xù)性置信度量CI(x)的極大值,即

    從而通過下式得到該區(qū)域內(nèi)的自適應(yīng)閾值參數(shù)β為

    式中α為根據(jù)保邊緣濾波要求而全局確定的百分比因子,thr是由剖面的整體噪聲干擾水平所確定的基底噪聲閾值.

    通過將待處理剖面分割為若干子區(qū)域,由(9)和(10)式將(8)式中對閾值參數(shù)β的絕對取值轉(zhuǎn)化為對百分比因子α的相對取值,從而能夠統(tǒng)一調(diào)整結(jié)構(gòu)自適應(yīng)中值濾波器的邊緣保持性態(tài).在實際資料處理中,并不需要精確地確定子區(qū)域的大小,一種典型的選擇是每個子區(qū)域內(nèi)包含100道,每道包含150個采樣點,即可滿足濾波過程保持斷層、裂縫等邊緣特征的要求.

    3 合成模型分析

    為了驗證結(jié)構(gòu)自適應(yīng)中值濾波器的保護邊緣去噪的有效性,本文選取主頻為20Hz的Ricker子波,采用射線追蹤方法得到一個含有斷層的復雜地層結(jié)構(gòu)模型的合成地震數(shù)據(jù)(見圖3a),在合成數(shù)據(jù)中加入5dB的帶限高斯隨機噪聲,得到圖3b的含噪數(shù)據(jù).

    圖3 含有斷層的復雜地層結(jié)構(gòu)的合成地震數(shù)據(jù)(a)原始數(shù)據(jù);(b)加噪數(shù)據(jù).Fig.3 Synthetic seismic data of complex geology model with a fault(a)Clean data;(b)Noisy data.

    采用結(jié)構(gòu)自適應(yīng)中值濾波器處理加噪的合成數(shù)據(jù),考慮到合成數(shù)據(jù)的信噪比很低,在梯度結(jié)構(gòu)張量的計算中,二維高斯低通濾波函數(shù)取較大的尺度參數(shù)ρ=4.0,結(jié)構(gòu)自適應(yīng)中值濾波器的最大濾波尺寸取為Rmax=ρ=4.0,即對應(yīng)最大尺度為9點的濾波窗口.針對圖3b中黑色長方形框所選定的典型區(qū)域,分析閾值百分比因子α對濾波結(jié)果的影響.調(diào)整α從10%開始以10%的間隔變化到200%,記錄每次濾波結(jié)果的信噪比,得到該區(qū)域的濾波結(jié)果的信噪比隨閾值百分比因子α的變化曲線,如圖4所示.從圖中可以看出,信噪比曲線隨α的增大而升高,且曲線在α=90%后升高逐漸變緩.這是由于隨著α的增加,濾波器的保邊緣性能減弱,當α>90%后,濾波過程已造成斷面處波形的損傷,在一定程度上抵消了濾波性能增強帶來的信噪比增加.因而在處理該加噪的合成數(shù)據(jù)時,選擇閾值百分比因子α=90%,此時對應(yīng)的結(jié)構(gòu)自適應(yīng)中值濾波器的各向異性窗函數(shù)的長軸σ1和短軸σ2的取值分布如圖5所示.對比圖5a和圖5b可見,結(jié)構(gòu)自適應(yīng)中值濾波器的尺度參數(shù)能夠很好地匹配合成數(shù)據(jù)的信號結(jié)構(gòu)特征,在局部呈線型結(jié)構(gòu)的同相軸區(qū)域橢圓形濾波窗口的長軸拉伸、短軸萎縮,形成有強方向選擇特性的各向異性窗口;在斷層區(qū)域長軸和短軸同時萎縮,形成小尺度分析窗口;而在無信號特征的區(qū)域,長軸和短軸均被拉伸,形成大尺度的各向同性窗口.

    圖4 不同α對應(yīng)的結(jié)構(gòu)自適應(yīng)中值濾波結(jié)果的信噪比Fig.4 SNR of filtered results by structure-adaptive median filter using differentα

    圖5 針對加噪合成數(shù)據(jù)的結(jié)構(gòu)自適應(yīng)中值濾波器的窗口參數(shù)(a)橢圓形濾波窗口的長軸σ1;(b)橢圓形濾波窗口的短軸σ2.Fig.5 Window parameters of structure-adaptive median filter corresponding to the noisy synthetic data(a)Main axesσ1of the elliptic filter window;(b)The second axesσ2of the elliptic filter window.

    加噪合成數(shù)據(jù)的結(jié)構(gòu)自適應(yīng)中值濾波結(jié)果及其得到的噪聲剖面分別如圖6a和圖6d所示.為了對比本文方法的處理效果,使用經(jīng)典的二維中值濾波器和具有信號細節(jié)保護能力的二維多級中值濾波器處理加噪合成數(shù)據(jù).取兩種中值濾波器的窗口尺寸與結(jié)構(gòu)自適應(yīng)中值濾波器的最大窗口尺寸相同,即二維中值濾波器選擇9×9的濾波窗口,二維多級中值濾波器取9點的濾波長度,兩者對應(yīng)的濾波結(jié)果以及得到的噪聲剖面分別如圖6b-c和圖6e-f所示.

    從圖6a可見,經(jīng)過結(jié)構(gòu)自適應(yīng)中值濾波后,隨機噪聲得到很大程度地衰減,反射同相軸的一致性得到增強,剖面整體變得非常干凈;從圖6d可見,結(jié)構(gòu)自適應(yīng)中值濾波器得到的噪聲剖面中只有隨機分布的干擾噪聲,無明顯的有效信號結(jié)構(gòu),斷層結(jié)構(gòu)在濾波過程中得到了有效保持.對比圖6b-c中兩種中值濾波方法的濾波結(jié)果,9×9的二維中值濾波器消除了大部分的隨機噪聲,但也造成有效信號的損傷,從其得到的噪聲剖面圖6e可見,濾波過程引起斷層信息的模糊,也造成反射同相軸的損傷,尤其對傾斜反射結(jié)構(gòu)的損傷更為明顯;9點濾波長度的二維多級中值濾波器雖然有較強的信號細節(jié)保護能力,但是從其得到的噪聲剖面圖6f可見,二維多級中值濾波器對隨機噪聲衰減程度有限.

    從上述分析、比較可以看出:在取相同濾波窗口尺度的條件下,相比二維中值濾波器和二維多級中值濾波器,結(jié)構(gòu)自適應(yīng)中值濾波器能夠根據(jù)地震剖面的信號結(jié)構(gòu)特征,自適應(yīng)地構(gòu)造與信號結(jié)構(gòu)相吻合的濾波窗口,在消除隨機噪聲和保護斷層等邊緣結(jié)構(gòu)之間找到了很好折中;結(jié)構(gòu)自適應(yīng)中值濾波器的噪聲衰減性能接近9×9的二維中值濾波器,而其對斷層及有效信號保護能力又優(yōu)于9點濾波長度的二維多級中值濾波器.

    4 實際資料算例

    圖6 加噪合成數(shù)據(jù)的濾波結(jié)果比較(a)結(jié)構(gòu)自適應(yīng)中值濾波結(jié)果;(b)二維中值濾波結(jié)果;(c)二維多級中值濾波結(jié)果;(d)結(jié)構(gòu)自適應(yīng)中值濾波器得到的噪聲剖面;(e)二維中值濾波器得到的噪聲剖面;(f)二維多級中值濾波器得到的噪聲剖面.Fig.6 Comparison of filtering results of the noisy synthetic data(a)Result obtained by structure-adaptive median filter;(b)Result obtained by 2Dmedian filter;(c)Result obtained by 2Dmultistage median filter;(d)Noise removed by structure-adaptive median filter;(e)Noise removed by 2Dmedian filter;(f)Noise removed by 2D multistage median filter.

    圖7 XX油田實際地震剖面Fig.7 Real seismic profile from XX oil field

    基于理論分析結(jié)果,進一步將結(jié)構(gòu)自適應(yīng)中值濾波器用于疊后地震剖面的濾波處理.圖7中給出某油田的疊后純波地震資料,剖面中含有豐富的不連續(xù)結(jié)構(gòu)特征,除了兩條主要的斷層外,在中央部位有一處雜亂反射區(qū)域,在右側(cè)偏下部有一處振幅異常區(qū)域.由于受到隨機噪聲和采集腳印的干擾,同相軸的連續(xù)性受到一定影響,不利于層位的橫向追蹤、不連續(xù)結(jié)構(gòu)的自動檢測和拾取等后繼的解釋和處理.運用本文的結(jié)構(gòu)自適應(yīng)中值濾波器處理該剖面,選擇100道、每道150個采樣點的子區(qū)域分割;在梯度結(jié)構(gòu)張量的計算中取高斯濾波器的尺度參數(shù)ρ=3.0;結(jié)構(gòu)自適應(yīng)中值濾波器的最大濾波尺寸取為Rmax=4.0,即對應(yīng)最大為9點的濾波窗口;為了達到較好的邊緣保持性能,取閾值百分比因子α=50%,對應(yīng)的濾波結(jié)果和得到的噪聲剖面分別如圖8a和圖8d所示.作為對比,同樣采用兩種中值濾波方法處理該剖面.由于二維中值濾波器的濾波窗口選擇得過大,對有效信號損傷嚴重,而選擇得過小,會導致噪聲衰減不夠,故此處取二維中值濾波器的窗口尺寸為5×5,而二維多級中值濾波器采用9點濾波長度,即與結(jié)構(gòu)自適應(yīng)中值濾波器的最大濾波尺寸相當,對應(yīng)的濾波結(jié)果和它們得到的噪聲剖面分別如圖8b-c、e-f所示.

    圖8 實際地震資料濾波結(jié)果比較(a)結(jié)構(gòu)自適應(yīng)中值濾波結(jié)果;(b)二維中值濾波結(jié)果;(c)二維多級中值濾波結(jié)果;(d)結(jié)構(gòu)自適應(yīng)中值濾波器得到的噪聲剖面;(e)二維中值濾波器得到的噪聲剖面;(f)二維多級中值濾波器得到的噪聲剖面.Fig.8 Comparison of filtering results of real seismic profile(a)Result obtained by structure-adaptive median filter;(b)Result obtained by 2Dmedian filter;(c)Result obtained by 2Dmultistage median filter;(d)Noise removed by structure-adaptive median filter;(e)Noise removed by 2Dmedian filter;(f)Noise removed by 2Dmultistage median filter.

    對比結(jié)構(gòu)自適應(yīng)中值濾波前后的剖面圖7和圖8a可見,結(jié)構(gòu)自適應(yīng)中值濾波過程消除了大部分的干擾噪聲,改善了反射同相軸的一致性,斷層構(gòu)造得到增強,振幅異常區(qū)域變得清晰.分析結(jié)構(gòu)自適應(yīng)中值濾波得到的噪聲剖面(圖8d)可見,淺層的干擾噪聲非常強,中深層的干擾噪聲相對較弱,剖面中除隨機噪聲外還有一些大角度的相干噪聲;對比原始地震剖面圖7,噪聲剖面中看不到明顯的有效信號結(jié)構(gòu).分析圖8b中二維中值濾波器的濾波結(jié)果,從視覺上看噪聲基本被壓制,信噪比得到很大的提高,但是信號的邊緣結(jié)構(gòu)尤其斷層、雜亂反射等被嚴重模糊,這一點可以從其得到的噪聲剖面(圖8e)得到印證.從圖8c和圖8f可見,二維多級中值濾波器僅僅沿同相軸走向有微弱的噪聲衰減作用,濾波能力非常有限.以上結(jié)果表明,結(jié)構(gòu)自適應(yīng)中值濾波器對淺層的強噪聲和深層的弱噪聲均能夠有效壓制,不僅能夠衰減隨機噪聲,也能夠消除部分相干噪聲,并且很好地保持了有效信號和斷層、振幅異常等地層邊緣和細節(jié)特征;結(jié)構(gòu)自適應(yīng)中值濾波器不僅具有二維中值濾波器優(yōu)異的噪聲衰減能力,而且又不會因保護地層邊緣和細節(jié)特征的需求而像二維多級中值濾波器那樣引起濾波能力的明顯降低.

    對原始地震剖面和三種濾波方法的處理結(jié)果進行頻譜分析,圖9給出相應(yīng)剖面的平均振幅譜.對比處理前后的頻譜發(fā)現(xiàn),三種方法處理后的高頻段隨機噪聲幅度均要低于處理前,且頻譜的主要結(jié)構(gòu)特點沒有被破壞;然而結(jié)構(gòu)自適應(yīng)中值濾波器和二維多級中值濾波器均較好地保持了主頻和有用信息頻段的振幅譜,而二維中值濾波器對主頻和有用信息頻段的振幅譜的損傷較為嚴重.以上結(jié)果說明結(jié)構(gòu)自適應(yīng)中值濾波器能夠有效壓制地震隨機噪聲,并且?guī)缀醪粨p傷有用信息頻段的頻譜.

    圖9 實際地震資料濾波結(jié)果平均振幅譜對比(a)圖7的平均振幅譜;(b)圖8(a)的平均振幅譜;(c)圖8(b)的平均振幅譜;(d)圖8(c)的平均振幅譜.Fig.9 Comparison of average spectrum of filtering results of real seismic profile(a)Average spectrum of Fig.7;(b)Average spectrum of Fig.8(a);(c)Average spectrum of Fig.8(b);(d)Average spectrum of Fig.8(c).

    5 結(jié) 論

    在實際地震資料的處理中,中值濾波器的應(yīng)用已經(jīng)比較成熟.本文針對現(xiàn)有的一些中值濾波方法在疊后地震資料隨機噪聲衰減方面存在的缺陷,提出了一種結(jié)構(gòu)自適應(yīng)中值濾波器.該方法依據(jù)梯度結(jié)構(gòu)張量對分析鄰域內(nèi)地層結(jié)構(gòu)規(guī)則程度的估計情況,引入了地震信號的兩種結(jié)構(gòu)特征的置信度量,利用它們調(diào)控結(jié)構(gòu)自適應(yīng)中值濾波器的濾波窗口的尺度和形狀,使得濾波窗口能夠最佳地匹配處理鄰域內(nèi)的地層結(jié)構(gòu)特征.通過合成模型和實際地震資料的處理結(jié)果表明,相比經(jīng)典的二維中值濾波器及二維多級中值濾波器,本文方法不僅能夠有效地衰減隨機噪聲、提高剖面的信噪比,而且能夠較好地保持有效信號和斷層、裂縫等地層邊緣和細節(jié)特征.本文方法可作為一個預處理過程,濾波結(jié)果進一步用于層位的自動追蹤、對噪聲敏感屬性的提取等后繼的解釋和處理流程.

    (References)

    [1] 李月,楊寶俊,趙雪平等.檢測地震勘探微弱同相軸的混沌振子算法.地球物理學報,2005,48(6):1428-1433.Li Y,Yang B J,Zhao X P,et al.An algorithm of chaotic vibrator to detect weak events in seismic prospecting records.Chinese J.Geophys.(in Chinese),2005,48(6):1428-1433.

    [2] 高靜懷,毛劍,滿蔚仕等.疊前地震資料噪聲衰減的小波域方法研究.地球物理學報,2006,49(4):1155-1163.Gao J H,Mao J,Man W S,et al.On the denoising method of prestack seismic data in wavelet domain.Chinese J.Geophys.(in Chinese),2006,49(4):1155-1163.

    [3] Neelamani R,Baumstein A I,Gillard D G,et al.Coherent and random noise attenuation using the curvelet transform.The Leading Edge,2008,27(2):240-248.

    [4] 劉洋,F(xiàn)omel S,劉財?shù)?高階seislet變換及其在隨機噪聲消除中的應(yīng)用.地球物理學報,2009,52(8):2142-2151.Liu Y,F(xiàn)omel S,Liu C,et al.High-order seislet transform and its application of random noise attenuation.Chinese J.Geophys.(in Chinese),2009,52(8):2142-2151.

    [5] 陸文凱.基于離散余弦變換的地震隨機噪聲壓制技術(shù).石油地球物理勘探,2011,46(2):202-206.Lu W K.Seismic random noise suppression based on the discrete cosine transform.Oil Geophysical Prospecting(in Chinese),2011,46(2):202-206.

    [6] 林紅波,李月,徐學純.壓制地震勘探隨機噪聲的分段時頻峰值濾波方法.地球物理學報,2011,54(5):1358-1366.Lin H B,Li Y,Xu X C.Segmenting time-frequency peak filtering method to attenuation of seismic random noise.Chinese J.Geophys.(in Chinese),2011,54(5):1358-1366.

    [7] Chopra S,Marfurt K J.Seismic attributes–a historical perspective.Geophysics,2005,70(5):3SO-28SO.

    [8] Chopra S,Marfurt K J.Emerging and future trends in seismic attributes.The Leading Edge,2008,27(3):298-318.

    [9] Chopra S,Misra S,Marfurt K J.Coherence and curvature attributes on preconditioned seismic data.The Leading Edge,2011,30(4)386-393.

    [10] Luo Y,Marhoon M,Al-Dossary S,et al.Acquistion Processing—Edge-preserving smoothing and applications.The Leading Edge,2002,21(2):136-158.

    [11] AlBinHassan N M,Luo Y,Al-Faraj M N.3Dedgepreserving smoothing and applications.Geophysics,2006,71(4):P5-P11.

    [12] Lu Y H,Lu W K.Edge-preserving polynomial fitting method to suppress random seismic noise.Geophysics,2009,74(4):V69-V73.

    [13] 楊培杰,穆星,張景濤.方向性邊界保持斷層增強技術(shù).地球物理學報,2010,53(12):2992-2997.Yang P J,Mu X,Zhang J T.Orientational edge preserving fault enhance.Chinese J.Gephys.(in Chinese),2010,53(12):2992-2997.

    [14] 劉財,王典,劉洋等.二維多級中值濾波技術(shù)在隨機噪聲消除中的應(yīng)用初探.石油地球物理勘探,2005,40(2):163-167.Liu C,Wang D,Liu Y,et al.Preliminary study of using 2D multi-level median filtering technique to eliminate random noises.Oil Geophysical Prospecting(in Chinese),2005,40(2):163-167.

    [15] Liu C,Liu Y,Yang B J,et al.A 2Dmultistage median filter to reduce random seismic noise.Geophysics,2006,71(5):V105-V110.

    [16] Fehmers G C,H?cker C F W.Fast structural interpretation with structure-oriented filtering.Geophysics,2003,68(4):1286-1293.

    [17] Lavialle O,Pop S,Germain C,et al.Seismic fault preserving diffusion.Journal of Applied Geophysics,2007,61(2):132-141.

    [18] 孫夕平,杜世通,湯磊.相干增強各向異性擴散濾波技術(shù).石油地球物理勘探,2004,39(6):651-655.Sun X P,Du S T,Tang L.Coherent-enhancing anisotropic diffusion filtering technique.Oil Geophysical Prospecting(in Chinese),2004,39(6):651-655.

    [19] 王旭松,楊長春.對地震圖像進行保邊濾波的非線性各向異性擴散算法.地球物理學進展,2006,21(2):452-457.Wang X S,Yang C C.An edge-preserving smoothing algorithm of seismic image using nonlinear anisotropic diffusion equation.Progress in Geophys.(in Chinese),2006,21(2):452-457.

    [20] 張爾華,王偉,高靜懷等.非線性各向異性擴散濾波器用于三維地震資料噪聲衰減與結(jié)構(gòu)特征增強.地球物理學進展,2010,25(3):866-870.Zhang E H,Wang W,Gao J H,et al.Non-linear anisotropic diffusion filtering for 3Dseismic noise removal and structure enhancement.Progress in Geophys.(in Chinese),2010,25(3):866-870.

    [21] Gómez G.Local smoothness in terms of variance.Proceedings of the BMVC,2000:815-824.

    [22] Marfurt K J,Sudhaker V,Gersztenkorn A,et al.Coherency calculations in the presence of structural dip.Geophysics,1999,64(1):104-111.

    [23] Randen T,Monsen E,Signer C,et al.Three-dimensional texture attributes for seismic data analysis.70th Annual International Meeting,SEG Expanded Abstracts,2000:668-671.

    [24] Wang Y E,Luo Y,Al-faraj M N.Inverse-vector approach for smoothing dips and azimuths.Geophysics,2008,73(6):P9-P14.

    [25] Gao D L.Latest developments in seismic texture analysis for subsurface structure,facies,and reservoir characterization:A review.Geophysics,2011,76(2):W1-W13.

    [26] Bakker P.Image structure analysis for seismic interpretation[Ph.D′s thesis].Netherlands:Technische Universiteit Delft,2002.

    Random seismic noise suppression via structure-adaptive median filter

    WANG Wei1,GAO Jing-Huai1,CHEN Wen-Chao1*,ZHU Zhen-Yu2
    1 School of Electronic and Information Engineering,Xi′an Jiaotong University,Xi′an 710049,China
    2 CNOOC Research Institute,Beijing100027,China

    This paper presents a structure-adaptive median filter for reducing random seismic noise which preserves seismic edges and details such as faults and fractures.By considering seismic horizons as line-like structures,this filtering framework relies on the computation of the Gradient Structure Tensors,which provides dips of seismic events and the regularity of local seismic structures,based on which two confidence measures are defined for different seismic structures.The structure-adaptive median filter adjusts the shape of the filter kernel according to the two confidence measures and aligns the filtering kernel along local dips to make the filtering kernel optimally matched with different local geological features.The proposed filter has been applied to both the synthetic and real data,and compared with two widely used median filtering schemes.The processing and comparison results show that the proposed structure-adaptive median filter is more suitable to balance suppressing random seismic noise and preserving signals which preserves seismic edges and details while enhancing the coherence of seismic horizons,and improves the quality of seismic images significantly.

    Median filter,Adaptive filtering,Noise suppression,Structure preserving,Gradient structure tensor

    10.6038/j.issn.0001-5733.2012.05.030

    P631

    2011-04-06,2011-06-28收修定稿

    國家自然科學基金重點項目(40730424)、國家自然科學基金面上項目(40674064)、國家油氣重大專項(2011ZX05023-005-009)資助.

    王偉,男,1983年生,現(xiàn)為西安交通大學在讀博士生,主要從事地震信號處理和地震屬性分析方面的研究.E-mail:feitianliuhuo@gmail.com

    *通訊作者陳文超,男,1970年生,副教授,2000年于西安交通大學獲博士學位,主要從事小波分析及盲信號分離在地震信號處理、解釋中的應(yīng)用研究.E-mail:wencchen@m(xù)ail.xjtu.edu.cn

    王偉,高靜懷,陳文超等.基于結(jié)構(gòu)自適應(yīng)中值濾波器的隨機噪聲衰減方法.地球物理學報,2012,55(5):1732-1741,

    10.6038/j.issn.0001-5733.2012.05.030.

    Wang W,Gao J H,Chen W C,et al.Random seismic noise suppression via structure-adaptive median filter.Chinese J.Geophys.(in Chinese),2012,55(5):1732-1741,doi:10.6038/j.issn.0001-5733.2012.05.030.

    (本文編輯 劉少華)

    猜你喜歡
    同相軸置信度剖面
    硼鋁復合材料硼含量置信度臨界安全分析研究
    三點法定交叉剖面方法
    ——工程地質(zhì)勘察中,一種做交叉剖面的新方法
    虛同相軸方法及其在陸上地震層間多次波壓制中的應(yīng)用
    正負關(guān)聯(lián)規(guī)則兩級置信度閾值設(shè)置方法
    基于曲線擬合的投棄式剖面儀電感量算法
    電子測試(2017年12期)2017-12-18 06:35:46
    一種改進的相關(guān)法自動拾取同相軸
    復雜多約束條件通航飛行垂直剖面規(guī)劃方法
    一種反射同相軸自動拾取算法
    置信度條件下軸承壽命的可靠度分析
    軸承(2015年2期)2015-07-25 03:51:04
    近年來龍門山斷裂GPS剖面變形與應(yīng)變積累分析
    地震研究(2014年3期)2014-02-27 09:30:50
    亚洲中文日韩欧美视频| 黄色欧美视频在线观看| 直男gayav资源| 欧美色欧美亚洲另类二区| 国内久久婷婷六月综合欲色啪| 尾随美女入室| 亚洲欧美日韩无卡精品| 日本爱情动作片www.在线观看 | 97人妻精品一区二区三区麻豆| 亚洲熟妇熟女久久| 国产黄色视频一区二区在线观看 | 如何舔出高潮| 91狼人影院| 久久热精品热| or卡值多少钱| 成人鲁丝片一二三区免费| 小蜜桃在线观看免费完整版高清| 蜜臀久久99精品久久宅男| 国产精品三级大全| 欧美xxxx性猛交bbbb| 成人综合一区亚洲| 在线观看av片永久免费下载| 美女 人体艺术 gogo| 男插女下体视频免费在线播放| 老司机影院成人| 亚洲第一电影网av| 国产探花极品一区二区| h日本视频在线播放| 自拍偷自拍亚洲精品老妇| 国产男人的电影天堂91| 六月丁香七月| 成年女人毛片免费观看观看9| 成人一区二区视频在线观看| 亚洲成人av在线免费| 国产三级中文精品| 亚洲成人精品中文字幕电影| 欧美日韩在线观看h| 国产一区二区在线av高清观看| 日本爱情动作片www.在线观看 | 国产探花极品一区二区| 人妻夜夜爽99麻豆av| 高清午夜精品一区二区三区 | 女人被狂操c到高潮| 直男gayav资源| 99久久精品国产国产毛片| 精品午夜福利在线看| 听说在线观看完整版免费高清| 国产高清激情床上av| 网址你懂的国产日韩在线| 人妻丰满熟妇av一区二区三区| 国产免费男女视频| 插逼视频在线观看| 国产亚洲精品久久久久久毛片| 中文亚洲av片在线观看爽| 在线观看av片永久免费下载| 国产精品99久久久久久久久| 日本黄色视频三级网站网址| 午夜福利在线在线| 国产爱豆传媒在线观看| 亚洲精品一区av在线观看| 亚洲三级黄色毛片| 亚洲国产精品国产精品| 精品午夜福利视频在线观看一区| 亚洲aⅴ乱码一区二区在线播放| 国内精品宾馆在线| 国产高清三级在线| 欧美最新免费一区二区三区| 亚洲成人中文字幕在线播放| 亚洲最大成人av| 欧美成人一区二区免费高清观看| 天堂网av新在线| 寂寞人妻少妇视频99o| 99热只有精品国产| 成人无遮挡网站| av黄色大香蕉| 成人特级黄色片久久久久久久| 丰满人妻一区二区三区视频av| av卡一久久| 国产69精品久久久久777片| 久久久久久久久久久丰满| 欧美潮喷喷水| 夜夜看夜夜爽夜夜摸| 久久久久久久久久黄片| 日本欧美国产在线视频| 国产精品不卡视频一区二区| 男人舔女人下体高潮全视频| 搡老妇女老女人老熟妇| 成人漫画全彩无遮挡| 岛国在线免费视频观看| 久久这里只有精品中国| 久久久成人免费电影| av国产免费在线观看| 九九久久精品国产亚洲av麻豆| 成年女人永久免费观看视频| 午夜a级毛片| 搡女人真爽免费视频火全软件 | 午夜老司机福利剧场| 亚洲最大成人手机在线| 99在线视频只有这里精品首页| 久久99热6这里只有精品| 热99re8久久精品国产| 精品久久久久久久久久久久久| 欧洲精品卡2卡3卡4卡5卡区| 国产aⅴ精品一区二区三区波| 国产探花在线观看一区二区| 精品福利观看| 亚洲欧美日韩无卡精品| 噜噜噜噜噜久久久久久91| 在线免费十八禁| 国产aⅴ精品一区二区三区波| 亚洲精品一卡2卡三卡4卡5卡| 日韩精品青青久久久久久| 一级毛片久久久久久久久女| 精品久久久久久成人av| 亚洲人成网站在线播放欧美日韩| 久久久久久大精品| 成熟少妇高潮喷水视频| 久久久久免费精品人妻一区二区| 99久久九九国产精品国产免费| 欧美日韩精品成人综合77777| 日韩一区二区视频免费看| 欧美+日韩+精品| 不卡一级毛片| 99久久精品热视频| 18禁在线播放成人免费| 97热精品久久久久久| 成年免费大片在线观看| 99精品在免费线老司机午夜| 亚洲av二区三区四区| 黄片wwwwww| 久久精品夜夜夜夜夜久久蜜豆| 香蕉av资源在线| 2021天堂中文幕一二区在线观| or卡值多少钱| 18禁裸乳无遮挡免费网站照片| 综合色丁香网| 免费在线观看成人毛片| 国产在线精品亚洲第一网站| 精品99又大又爽又粗少妇毛片| 国产免费男女视频| 自拍偷自拍亚洲精品老妇| 免费不卡的大黄色大毛片视频在线观看 | 一本久久中文字幕| 精品一区二区免费观看| 久久6这里有精品| 国产黄色视频一区二区在线观看 | 日本免费a在线| 在线看三级毛片| 亚洲欧美成人精品一区二区| 此物有八面人人有两片| 丰满人妻一区二区三区视频av| 在现免费观看毛片| 成人一区二区视频在线观看| 国产精品久久久久久亚洲av鲁大| a级毛片a级免费在线| 最新在线观看一区二区三区| 伦精品一区二区三区| 3wmmmm亚洲av在线观看| 性插视频无遮挡在线免费观看| 欧洲精品卡2卡3卡4卡5卡区| 日本黄色视频三级网站网址| 美女黄网站色视频| 日本欧美国产在线视频| 麻豆一二三区av精品| 亚洲真实伦在线观看| 欧美性感艳星| 国产69精品久久久久777片| 日韩国内少妇激情av| 自拍偷自拍亚洲精品老妇| 日本a在线网址| 日本欧美国产在线视频| 国产又黄又爽又无遮挡在线| 欧美+日韩+精品| 日韩一本色道免费dvd| 床上黄色一级片| 久久99热6这里只有精品| 亚洲国产精品合色在线| 少妇人妻一区二区三区视频| 综合色丁香网| 老司机福利观看| 婷婷亚洲欧美| 久久精品国产自在天天线| 久久久精品欧美日韩精品| 国产精品嫩草影院av在线观看| 国产av一区在线观看免费| 成年版毛片免费区| 免费看a级黄色片| 国产蜜桃级精品一区二区三区| 国产成人福利小说| 超碰av人人做人人爽久久| 乱系列少妇在线播放| 久久精品国产亚洲av天美| 国产精品不卡视频一区二区| 国产三级在线视频| 亚洲熟妇中文字幕五十中出| 国内精品宾馆在线| 亚洲三级黄色毛片| 久久久久久九九精品二区国产| 床上黄色一级片| 国模一区二区三区四区视频| 久久草成人影院| 伊人久久精品亚洲午夜| 3wmmmm亚洲av在线观看| 美女xxoo啪啪120秒动态图| 99热6这里只有精品| av在线老鸭窝| av.在线天堂| 日日撸夜夜添| 亚洲欧美中文字幕日韩二区| 国产在线男女| 99国产极品粉嫩在线观看| 搡老妇女老女人老熟妇| 免费黄网站久久成人精品| 在线天堂最新版资源| 久久久久久久午夜电影| 国产真实乱freesex| 波多野结衣巨乳人妻| 国产真实伦视频高清在线观看| 日日摸夜夜添夜夜爱| 亚洲天堂国产精品一区在线| 亚洲成人久久性| 97在线视频观看| 午夜激情欧美在线| 内射极品少妇av片p| 男女下面进入的视频免费午夜| 久久精品国产清高在天天线| 99久久久亚洲精品蜜臀av| 久久精品综合一区二区三区| 最近在线观看免费完整版| 人妻丰满熟妇av一区二区三区| 三级国产精品欧美在线观看| 深夜精品福利| 18禁在线无遮挡免费观看视频 | 国产精品人妻久久久久久| 91av网一区二区| 国产精品日韩av在线免费观看| 亚洲最大成人手机在线| 国产色爽女视频免费观看| 成人永久免费在线观看视频| 欧美又色又爽又黄视频| 综合色av麻豆| 国产亚洲欧美98| 三级经典国产精品| 国产午夜精品久久久久久一区二区三区 | 九色成人免费人妻av| 精品人妻熟女av久视频| 亚洲丝袜综合中文字幕| 两个人的视频大全免费| 99国产极品粉嫩在线观看| 少妇丰满av| 婷婷精品国产亚洲av在线| 乱人视频在线观看| 亚洲人成网站在线播| 小说图片视频综合网站| 精品午夜福利在线看| 一级黄片播放器| 特级一级黄色大片| 亚洲美女黄片视频| 97在线视频观看| 老司机影院成人| 看片在线看免费视频| av天堂在线播放| 女人十人毛片免费观看3o分钟| av黄色大香蕉| 日韩一本色道免费dvd| 日本爱情动作片www.在线观看 | 男女下面进入的视频免费午夜| 亚洲最大成人手机在线| 国产精品不卡视频一区二区| 一a级毛片在线观看| 一级毛片我不卡| 99久久精品一区二区三区| 免费看a级黄色片| av卡一久久| 亚洲欧美日韩高清在线视频| 中文字幕av成人在线电影| 欧美一区二区亚洲| 蜜桃久久精品国产亚洲av| av.在线天堂| 日韩强制内射视频| 91精品国产九色| 亚洲成人av在线免费| 97在线视频观看| 晚上一个人看的免费电影| 老熟妇乱子伦视频在线观看| 床上黄色一级片| 久久久国产成人精品二区| 两个人的视频大全免费| 久久久色成人| 身体一侧抽搐| 亚洲最大成人手机在线| 精品熟女少妇av免费看| 久久国内精品自在自线图片| 插逼视频在线观看| 久久久久九九精品影院| 在线观看免费视频日本深夜| 精品久久久久久久末码| 国产白丝娇喘喷水9色精品| 在线观看美女被高潮喷水网站| 国产探花在线观看一区二区| 天堂√8在线中文| 日韩欧美国产在线观看| 国产精品一及| 日本三级黄在线观看| av天堂在线播放| 国产精品美女特级片免费视频播放器| 人妻夜夜爽99麻豆av| 夜夜爽天天搞| 网址你懂的国产日韩在线| 搡老妇女老女人老熟妇| 免费看光身美女| 欧美一区二区亚洲| 特级一级黄色大片| 日本爱情动作片www.在线观看 | 成人无遮挡网站| 白带黄色成豆腐渣| 亚洲aⅴ乱码一区二区在线播放| 久久久久久久久久成人| 免费大片18禁| 欧美激情国产日韩精品一区| 一个人免费在线观看电影| 亚洲人成网站在线观看播放| 久久久久久久午夜电影| 免费人成在线观看视频色| 国产精品国产三级国产av玫瑰| 亚洲久久久久久中文字幕| 18禁黄网站禁片免费观看直播| 久久精品国产亚洲网站| 人人妻,人人澡人人爽秒播| 一进一出抽搐动态| 一进一出抽搐gif免费好疼| 久久久久九九精品影院| 精品一区二区三区视频在线观看免费| 男插女下体视频免费在线播放| 在线a可以看的网站| 成人午夜高清在线视频| 黄色配什么色好看| 卡戴珊不雅视频在线播放| 一个人看视频在线观看www免费| 国产精品日韩av在线免费观看| 联通29元200g的流量卡| 长腿黑丝高跟| 精品久久久久久久人妻蜜臀av| АⅤ资源中文在线天堂| av专区在线播放| 亚洲熟妇中文字幕五十中出| 欧美潮喷喷水| 国产精品三级大全| 在线国产一区二区在线| 黄色视频,在线免费观看| 久久久精品大字幕| 麻豆国产97在线/欧美| 亚洲欧美日韩高清在线视频| 麻豆国产97在线/欧美| 一级毛片aaaaaa免费看小| 国产欧美日韩精品亚洲av| 人妻丰满熟妇av一区二区三区| 亚洲18禁久久av| 精品一区二区三区av网在线观看| 国产精品一区二区三区四区久久| 亚洲丝袜综合中文字幕| 露出奶头的视频| 欧美高清成人免费视频www| 露出奶头的视频| 免费av观看视频| 色综合站精品国产| 国产亚洲精品综合一区在线观看| 免费看美女性在线毛片视频| 国产精品乱码一区二三区的特点| 一个人观看的视频www高清免费观看| 久久久久国产精品人妻aⅴ院| 91在线精品国自产拍蜜月| 少妇熟女aⅴ在线视频| 岛国在线免费视频观看| 波野结衣二区三区在线| 禁无遮挡网站| 午夜福利高清视频| 亚洲成a人片在线一区二区| 日本与韩国留学比较| 国产69精品久久久久777片| 国产精品无大码| 蜜桃亚洲精品一区二区三区| 国产精品av视频在线免费观看| 亚洲欧美精品综合久久99| 亚洲乱码一区二区免费版| 狂野欧美白嫩少妇大欣赏| 日韩一本色道免费dvd| 秋霞在线观看毛片| 国产真实伦视频高清在线观看| 蜜桃亚洲精品一区二区三区| 男女做爰动态图高潮gif福利片| 午夜精品国产一区二区电影 | АⅤ资源中文在线天堂| 久久久欧美国产精品| 亚洲性夜色夜夜综合| 国产亚洲精品综合一区在线观看| 精品免费久久久久久久清纯| 欧美不卡视频在线免费观看| 97热精品久久久久久| 夜夜看夜夜爽夜夜摸| 好男人在线观看高清免费视频| 人妻夜夜爽99麻豆av| 一进一出抽搐动态| 国产一区二区三区av在线 | 国产黄片美女视频| 免费观看的影片在线观看| 永久网站在线| 精品无人区乱码1区二区| 免费看a级黄色片| 国产精品久久久久久久久免| 亚洲美女视频黄频| 十八禁国产超污无遮挡网站| 97在线视频观看| 中国美白少妇内射xxxbb| 搡老妇女老女人老熟妇| 国产高清不卡午夜福利| 国内精品宾馆在线| 亚洲精品日韩在线中文字幕 | 内射极品少妇av片p| 欧美激情在线99| 一级黄色大片毛片| 最近手机中文字幕大全| 日韩人妻高清精品专区| 国产精品伦人一区二区| av在线老鸭窝| 色综合亚洲欧美另类图片| 又黄又爽又刺激的免费视频.| 亚洲18禁久久av| 真人做人爱边吃奶动态| 搡老妇女老女人老熟妇| 在线观看免费视频日本深夜| 国产人妻一区二区三区在| 久久久午夜欧美精品| 啦啦啦啦在线视频资源| 久久精品国产亚洲av香蕉五月| 又黄又爽又刺激的免费视频.| 亚洲精品粉嫩美女一区| 欧美日韩国产亚洲二区| 99久久无色码亚洲精品果冻| 不卡一级毛片| 国产视频内射| 国内揄拍国产精品人妻在线| av在线老鸭窝| 狂野欧美激情性xxxx在线观看| 色5月婷婷丁香| 日韩av在线大香蕉| 国产精品亚洲一级av第二区| 成人性生交大片免费视频hd| 少妇人妻一区二区三区视频| 中文字幕人妻熟人妻熟丝袜美| 久久人人爽人人爽人人片va| 91狼人影院| 毛片一级片免费看久久久久| 人妻夜夜爽99麻豆av| 麻豆乱淫一区二区| 嫩草影院精品99| 久久久久久久久久黄片| 久久草成人影院| 免费观看在线日韩| 97热精品久久久久久| 99热网站在线观看| 最近2019中文字幕mv第一页| 身体一侧抽搐| 亚洲乱码一区二区免费版| 国产成人aa在线观看| 亚洲av成人精品一区久久| 国产精品一区二区三区四区久久| 人妻丰满熟妇av一区二区三区| 国产老妇女一区| 乱人视频在线观看| av在线亚洲专区| 超碰av人人做人人爽久久| 久久亚洲国产成人精品v| 麻豆国产97在线/欧美| 精品人妻一区二区三区麻豆 | 久久中文看片网| 久久人妻av系列| 91麻豆精品激情在线观看国产| 色尼玛亚洲综合影院| 国产精品久久久久久精品电影| 如何舔出高潮| 久久国产乱子免费精品| 国产单亲对白刺激| 亚洲真实伦在线观看| 成人二区视频| 亚洲国产精品sss在线观看| 精品无人区乱码1区二区| 久久精品人妻少妇| 韩国av在线不卡| 国产av在哪里看| 国产亚洲欧美98| 欧美日本亚洲视频在线播放| 日韩 亚洲 欧美在线| ponron亚洲| 国产精品一区www在线观看| 国内精品美女久久久久久| 老司机影院成人| 一区二区三区高清视频在线| 国产欧美日韩一区二区精品| 国产免费一级a男人的天堂| 人妻少妇偷人精品九色| 国产精品一区二区性色av| 99在线视频只有这里精品首页| 亚洲中文字幕一区二区三区有码在线看| 久久九九热精品免费| 1000部很黄的大片| 成人性生交大片免费视频hd| 天堂影院成人在线观看| 久久精品人妻少妇| 一个人观看的视频www高清免费观看| 亚洲专区国产一区二区| 嫩草影院入口| 在线观看66精品国产| 在线观看午夜福利视频| 亚洲人成网站在线播放欧美日韩| 亚洲精品久久国产高清桃花| 国产高清不卡午夜福利| 一区二区三区免费毛片| 亚洲人成网站在线播放欧美日韩| 最近手机中文字幕大全| 国产视频一区二区在线看| 春色校园在线视频观看| 在线天堂最新版资源| 国产亚洲精品久久久com| 寂寞人妻少妇视频99o| 久久久国产成人免费| 国产乱人视频| 在线免费观看不下载黄p国产| 99久久久亚洲精品蜜臀av| 三级毛片av免费| 国产精品福利在线免费观看| 如何舔出高潮| 久久婷婷人人爽人人干人人爱| 又爽又黄a免费视频| 久久久久久久久久黄片| 国产精品美女特级片免费视频播放器| 国产精品99久久久久久久久| 欧美zozozo另类| 亚洲精品一卡2卡三卡4卡5卡| 国产午夜福利久久久久久| av女优亚洲男人天堂| 日韩精品青青久久久久久| 少妇人妻精品综合一区二区 | 91狼人影院| 国产高清有码在线观看视频| av福利片在线观看| 日本爱情动作片www.在线观看 | 干丝袜人妻中文字幕| 成人三级黄色视频| 色视频www国产| 神马国产精品三级电影在线观看| 中文字幕精品亚洲无线码一区| 国产一区二区三区av在线 | 久久人妻av系列| 国产国拍精品亚洲av在线观看| 真人做人爱边吃奶动态| 亚洲av一区综合| 国语自产精品视频在线第100页| 国模一区二区三区四区视频| 精品人妻熟女av久视频| av在线蜜桃| 日本a在线网址| 亚洲国产日韩欧美精品在线观看| 综合色丁香网| 国产白丝娇喘喷水9色精品| 99久国产av精品国产电影| 国产av在哪里看| 午夜福利成人在线免费观看| 亚洲欧美成人综合另类久久久 | 一卡2卡三卡四卡精品乱码亚洲| 成人亚洲欧美一区二区av| 国产成年人精品一区二区| 久久亚洲国产成人精品v| 欧美成人一区二区免费高清观看| 国产精品综合久久久久久久免费| 国产麻豆成人av免费视频| 亚洲欧美日韩无卡精品| 精品久久久久久久久久久久久| 国国产精品蜜臀av免费| 久久久久久久久久成人| 国产成人影院久久av| 国产三级中文精品| 一级av片app| 久久久久久久久大av| 亚洲人成网站在线观看播放| 午夜亚洲福利在线播放| 欧美不卡视频在线免费观看| 免费观看的影片在线观看| 自拍偷自拍亚洲精品老妇| 乱码一卡2卡4卡精品| 中国国产av一级| 国内精品一区二区在线观看| 亚洲四区av| 国产一区二区三区av在线 | 两性午夜刺激爽爽歪歪视频在线观看| 亚洲精品国产av成人精品 | 国产亚洲欧美98| 日日撸夜夜添| 日日摸夜夜添夜夜添av毛片| 91在线精品国自产拍蜜月| 国产免费男女视频| 亚洲欧美清纯卡通| 日本与韩国留学比较| 亚洲图色成人| 波多野结衣高清作品| 国产伦在线观看视频一区| 91久久精品国产一区二区成人| 熟妇人妻久久中文字幕3abv| 成人午夜高清在线视频| 国产午夜福利久久久久久| 亚洲一区高清亚洲精品| 免费电影在线观看免费观看| 91狼人影院| 久久韩国三级中文字幕| 亚洲中文字幕一区二区三区有码在线看|