亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        Matlab在CSAMT數(shù)據(jù)“去噪”及近場(chǎng)校正中的應(yīng)用

        2019-11-19 02:57:06胡慶輝
        物探化探計(jì)算技術(shù) 2019年5期
        關(guān)鍵詞:近場(chǎng)頻域電阻率

        胡慶輝

        (1.山東省地質(zhì)礦產(chǎn)勘查開發(fā)局 第四地質(zhì)大隊(duì),山東 濰坊 261021; 2.山東省地礦局 海岸帶地質(zhì)環(huán)境保護(hù)重點(diǎn)實(shí)驗(yàn)室,山東 濰坊 261021)

        0 引言

        可控源音頻大地電磁測(cè)深法(CSAMT),是在大地地磁測(cè)深法(MT)和音頻大地電磁測(cè)深法(AMT)的基礎(chǔ)上發(fā)展起來的一種人工源頻率域測(cè)深方法[1]。自上世紀(jì)70年代引入國(guó)內(nèi),在油氣勘查、金屬礦普查、水文及環(huán)境地質(zhì)等多個(gè)領(lǐng)域顯示出廣闊的應(yīng)用前景。隨著經(jīng)濟(jì)社會(huì)發(fā)展,CSAMT野外數(shù)據(jù)采集所面臨的“電磁環(huán)境”越發(fā)復(fù)雜,勘探信號(hào)不可避免地受到“污染”,導(dǎo)致最終反演數(shù)據(jù)質(zhì)量不高,加大了資料解譯的難度甚至造成失誤。此外,受限于CSAMT方法的原理及假設(shè)前提,在中、低頻率上電磁場(chǎng)進(jìn)入“近區(qū)”或者“過渡區(qū)”,卡尼亞視電阻率將發(fā)生畸變,無法實(shí)現(xiàn)勘查目的[2]。底青云等[3]認(rèn)為CSAMT資料處理的質(zhì)量依賴于近場(chǎng)校正的質(zhì)量。

        國(guó)內(nèi)學(xué)者針對(duì)CSAMT方法的電磁干擾壓制及近場(chǎng)校正技術(shù)進(jìn)行了相關(guān)研究。偽隨機(jī)編碼技術(shù)在供電環(huán)節(jié)的使用能夠有效提高野外觀測(cè)數(shù)據(jù)信噪比[4-5],而雷達(dá)等采用信息熵和有理函數(shù)濾波相結(jié)合的方式對(duì)干擾數(shù)據(jù)進(jìn)行處理也取得到比較好的效果[6]。通常的近場(chǎng)校正的方法是“三角形校正法”[3,7-9],而采用“全頻域視電阻率法”進(jìn)行近場(chǎng)校正得到越來越廣泛的關(guān)注和應(yīng)用[10-11]。

        用戶界面(或接口)是指人與機(jī)器(或程序)之間交互作用的工具和方法[12]。筆者采用Matlab GUI編程,實(shí)現(xiàn)友好的可視化交互式操作界面,剔除“畸變”頻點(diǎn)數(shù)據(jù),實(shí)現(xiàn)CSAMT數(shù)據(jù)的“去噪”處理。采用數(shù)值計(jì)算功能,編寫了基于全頻域視電阻率計(jì)算的程序以實(shí)現(xiàn)近場(chǎng)校正。程序的應(yīng)用提高了CSAMT數(shù)據(jù)的處理效率和效果,保證后續(xù)處理的數(shù)據(jù)質(zhì)量,在實(shí)際生產(chǎn)中得到應(yīng)用并取得了很好的效果。

        圖1 CSAMT野外實(shí)測(cè)干擾數(shù)據(jù)典型特征Fig.1 Typical feature of interfering survey data of CSAMT(a)卡尼亞視電阻率曲線;(b)阻抗相位曲線;(c)電場(chǎng)分量幅值曲線;(d)磁場(chǎng)分量幅值曲線

        1 CSAMT干擾及“近場(chǎng)”數(shù)據(jù)特征

        CSAMT野外數(shù)據(jù)環(huán)節(jié)的電磁干擾主要來自天然電磁噪聲和人文電磁噪聲[13],其中天然電磁噪聲主要來自雷電活動(dòng),其引起的電磁場(chǎng)在CSAMT常用工作頻段內(nèi)能量較高。這類噪聲一般可采用多次疊加進(jìn)行壓制。人文電磁噪聲主要包括工業(yè)用電、大功率發(fā)射臺(tái)以及隨機(jī)的電磁脈沖干擾(如動(dòng)力電開關(guān)轉(zhuǎn)換、鐵磁性物體靠近磁探頭等)。固定頻率的電磁干擾可以通過設(shè)置陷波器進(jìn)行壓制,而能量較高、隨機(jī)的強(qiáng)干擾還不能完全克服,需要野外施工人員正確判斷電磁干擾特點(diǎn),采取必要的規(guī)避措施,盡量避開干擾源或者錯(cuò)開干擾時(shí)段,以便采集到相對(duì)高質(zhì)量的野外數(shù)據(jù)。

        隨機(jī)電磁干擾在CSAMT實(shí)測(cè)數(shù)據(jù)中主要表現(xiàn)為電場(chǎng)及磁場(chǎng)分量復(fù)測(cè)離差增大,單頻點(diǎn)數(shù)據(jù)離差加大,數(shù)據(jù)一致性較差,阻抗相位偏差較大,甚至出現(xiàn)反相,計(jì)算得出的卡尼亞電阻率出現(xiàn)“飛點(diǎn)”,曲線形態(tài)不能反映真實(shí)大地電阻率的變化特征(圖1)。

        基于上述兩個(gè)方面的問題,筆者利用Matlab數(shù)值計(jì)算及GUI功能,設(shè)計(jì)了針對(duì)野外數(shù)據(jù)預(yù)處理和近場(chǎng)校正的相關(guān)程序,在CSAMT數(shù)據(jù)“去噪”預(yù)處理及近場(chǎng)校正等方面進(jìn)行了探索,取得較好的應(yīng)用效果。

        2 CSAMT“去噪”及“近場(chǎng)校正”方法

        現(xiàn)階段,可控源音頻大地電磁方法(CSAMT)“去噪”多采用濾波方法,在干擾比較小的情況下,一般能夠得到比較好的效果。但隨著干擾的加強(qiáng),濾波方法往往難以取得理想的電磁噪聲壓制效果。

        本次“去噪”處理思路比較直觀,即從多次觀測(cè)數(shù)據(jù)中挑選“合理”數(shù)據(jù),舍棄畸變數(shù)據(jù),從而達(dá)到“去偽存真”的目的。所謂“合理”,即從電極布設(shè)及采集時(shí)就保證觀測(cè)質(zhì)量,合理避開干擾源或干擾地段,錯(cuò)開干擾時(shí)段進(jìn)行數(shù)據(jù)采集,使多數(shù)觀測(cè)數(shù)據(jù)為真實(shí)地電反映,其數(shù)據(jù)分布應(yīng)具有“集中分布”特征,重復(fù)性較好。而少數(shù)觀測(cè)值由于隨機(jī)干擾產(chǎn)生“畸變”,相關(guān)參數(shù)如相位、電場(chǎng)分量幅值、磁場(chǎng)分量幅值等偏離合理區(qū)間,出現(xiàn)明顯的“飛點(diǎn)”特征,通過友好界面對(duì)畸變數(shù)據(jù)予以剔除,從而實(shí)現(xiàn)“去噪”處理。

        圖2 CSAMT數(shù)據(jù)預(yù)處理GUI界面Fig.2 Matlab GUI interface for data pre-processing of CSAMT

        數(shù)據(jù)處理通過鼠標(biāo)在界面上框選畸變數(shù)據(jù)實(shí)現(xiàn)便捷操作(圖2)。用戶界面分為三部分:①工具欄;②繪圖區(qū);③信息提示區(qū)。工具區(qū)包括數(shù)據(jù)加載、重置、顯示參數(shù)選項(xiàng)、數(shù)據(jù)操作、數(shù)據(jù)保存等功能按鍵或下拉菜單。繪圖區(qū)則顯示相應(yīng)選項(xiàng)的散點(diǎn)圖、擬合曲線圖或者擬斷面等值線圖。界面下方為信息提示區(qū),顯示部分文件信息以及當(dāng)前點(diǎn)號(hào)等。

        加載CSAMT數(shù)據(jù)完成后,繪圖區(qū)默認(rèn)顯示首個(gè)測(cè)點(diǎn)的卡尼亞電阻率數(shù)據(jù),綠色圓點(diǎn)為某個(gè)頻率重復(fù)觀測(cè)數(shù)據(jù)點(diǎn),藍(lán)色曲線為其樣條插值擬合曲線??梢钥闯觯陔姶鸥蓴_條件下采集的數(shù)據(jù)存在明顯的離散特征,即重復(fù)觀測(cè)數(shù)據(jù)一致性較差,“飛點(diǎn)”較多。同時(shí),每個(gè)頻率的觀測(cè)數(shù)據(jù)點(diǎn)有重復(fù)性較好或相對(duì)較為集中的特征,我們認(rèn)為這樣的數(shù)據(jù)為“合理”數(shù)據(jù)。用戶使用鼠標(biāo)選擇工具欄“刪除”按鍵,在繪圖區(qū)框選需要剔除的數(shù)據(jù)點(diǎn)(圖2中紅色框及紅色點(diǎn)),它們將不參與后續(xù)的數(shù)據(jù)處理。利用下拉菜單選擇“阻抗相位”、“電場(chǎng)分量幅值”及“磁場(chǎng)分量幅值”三個(gè)選項(xiàng)在相應(yīng)數(shù)據(jù)散點(diǎn)上完成類似操作,實(shí)現(xiàn)多個(gè)參數(shù)對(duì)同一個(gè)測(cè)點(diǎn)的“去噪”處理。當(dāng)該測(cè)點(diǎn)數(shù)據(jù)選擇完畢后,通過“保存”按鍵實(shí)現(xiàn)AVG數(shù)據(jù)保存,為后續(xù)處理提供高信噪比數(shù)據(jù)。

        CSAMT數(shù)據(jù)近場(chǎng)校正是在上述AVG數(shù)據(jù)的基礎(chǔ)上,利用“等效電阻率全頻域視電阻率法”,采用迭代擬合方式計(jì)算等效電阻率。

        通常采用赤道偶極方式進(jìn)行CSAMT數(shù)據(jù)采集時(shí),采集參數(shù)為電場(chǎng)和磁場(chǎng)的水平分量Ex、Hy,卡尼亞視電阻率可以表示為式(1)。

        (1)

        (2)

        利用ρs是P的單調(diào)下降函數(shù),通過實(shí)測(cè)ρs確定P值,并通過式(3)計(jì)算大地電阻率。

        (3)

        Matlab數(shù)值計(jì)算的實(shí)現(xiàn)過程中,采用迭代方式求取F(P)=0的根。具體思路為:給定取值范圍為[a,b] ,給P賦初值c=(a+b)/2,判斷F(a).F(c)符號(hào),根據(jù)函數(shù)的單調(diào)下降性質(zhì),采用二分法逐步逼近零點(diǎn),從而求取P值進(jìn)而得到大地電阻率ρ。

        近場(chǎng)校正處理通過Matlab GUI編程實(shí)現(xiàn)友好用戶界面進(jìn)行操作。通過主界面“近場(chǎng)及靜校正”按鍵,可以調(diào)用NFC(近場(chǎng)校正)模塊(圖3)。NFC界面包含工具欄及繪圖區(qū)兩部分,工具欄主要包括數(shù)據(jù)加載、數(shù)據(jù)列表、數(shù)據(jù)編輯、近場(chǎng)校正及數(shù)據(jù)保存按鍵。在加載數(shù)據(jù)環(huán)節(jié)需要同時(shí)加載CSAMT實(shí)測(cè)數(shù)據(jù)(*.AVG)、觀測(cè)點(diǎn)坐標(biāo)文件(*.stn),用于計(jì)算每個(gè)測(cè)點(diǎn)的收發(fā)距及相對(duì)角度,計(jì)算結(jié)果將在相應(yīng)繪圖區(qū)顯示。繪圖區(qū)顯示測(cè)點(diǎn)卡尼亞視電阻率曲線和阻抗相位曲線,用戶可以通過進(jìn)一步的編輯功能對(duì)曲線進(jìn)行處理,調(diào)整曲線形態(tài),使其更符合“非波區(qū)”曲線特征。計(jì)算完成后則顯示近場(chǎng)校正前后曲線對(duì)比情況,使處理人員能夠初步了解校正效果。

        圖3 CSAMT數(shù)據(jù)近場(chǎng)校正GUI界面Fig.3 Matlab GUI interface for near-field correction of CSAMT

        3 近場(chǎng)數(shù)據(jù)正演模擬計(jì)算

        3.1 均勻半空間地電模型

        均勻半空間地電模型電阻率為1 000 Ω·m,正演計(jì)算收發(fā)距r=6 km及r=14 km兩種情況下的卡尼亞視電阻率特征(圖4)??梢钥闯?,高頻段視電阻率能夠正確反映模型電阻率特征,當(dāng)進(jìn)入到中低頻率時(shí),曲線發(fā)生畸變,近45°上揚(yáng),嚴(yán)重偏離模型電阻率。采用“等效電阻率全頻域視電阻率法”校正后,全頻域視電阻率曲線已經(jīng)基本水平,較準(zhǔn)確反映了模型電阻率數(shù)值,說明利用全頻域視電阻率法進(jìn)行CSAMT近場(chǎng)校正效果非常理想。

        圖4 均勻半空間CSAMT正演數(shù)據(jù)近場(chǎng)校正處理曲線對(duì)比Fig.4 Near-field correction curve contrast of CSAMT forward data in the homogeneous half-space model

        圖5 層狀介質(zhì)CSAMT正演數(shù)據(jù)近場(chǎng)校正處理曲線Fig.5 Near-field correction curve of CSAMT forward data in the layered model(a)H型測(cè)深曲線;(b)A型測(cè)深曲線;(c)K型測(cè)深曲線;(d)Q型測(cè)深曲線

        圖6 某測(cè)線CSAMT 1900點(diǎn)去噪及近場(chǎng)校正處理效果Fig.6 Data denoise and near-field correction process results of survey point 1900(a)實(shí)測(cè)卡尼亞電阻率數(shù)據(jù)散點(diǎn)及平均曲線;(b)實(shí)測(cè)阻抗相位數(shù)據(jù)散點(diǎn)及平均曲線;(c)“去噪”及近場(chǎng)校正后曲線對(duì)比;(d)“去噪”后阻抗相位數(shù)據(jù)點(diǎn)及平均曲線

        圖7 CSAMT測(cè)量L4線視電阻率斷面圖(預(yù)處理前)Fig.7 Apparent resistivity section of CSAMT survey of line 4(before preprocessing)

        3.2 層狀地電模型

        層狀地電模型采用三層模型,分別設(shè)置了H型、A型、K型及Q型四種測(cè)深曲線類型(圖5)。在中低頻段,各地電模型正演曲線進(jìn)入“近區(qū)”,曲線形態(tài)發(fā)生畸變。采用“全頻域視電阻率法”校正后,測(cè)深曲線形態(tài)基本反映了真實(shí)的地電模型特征。

        同時(shí)也應(yīng)該注意到,在均勻半空間模擬中,過渡區(qū)出現(xiàn)局部“下凹”,校正值偏離模型電阻率。對(duì)于層狀模型,尤其是K型測(cè)深曲線,中間高阻層雖然趨勢(shì)較明顯,但其值在校正后明顯偏低。對(duì)于過渡區(qū)“校正不足”現(xiàn)象,筆者也做了探索性工作,主要思路是在中低頻段數(shù)據(jù)近場(chǎng)校正時(shí)增加修正系數(shù)。根據(jù)大量正演模型試算數(shù)據(jù),對(duì)近區(qū)和過渡區(qū)近場(chǎng)校正修正系數(shù)采用多項(xiàng)式擬合。修正系數(shù)公式見式(4)及式(5)。

        近區(qū)修正系數(shù):

        K(f)=-0.0071f3+0.0773f2-

        0.2017f+1.0634

        (4)

        過渡區(qū)修正系數(shù):

        K(f)=-0.0018f3+0.0111f2-

        0.0164f+1.6322

        (5)

        其中,f為校正點(diǎn)頻率。在進(jìn)行近場(chǎng)校正時(shí),合理選擇近區(qū)及過渡區(qū)頻率,并在校正后數(shù)值基礎(chǔ)上進(jìn)行修正。該方法在正演數(shù)據(jù)中得到了很好的校正效果,過渡區(qū)“下凹”現(xiàn)象得到改善,近區(qū)數(shù)據(jù)也得到適當(dāng)調(diào)整,在野外實(shí)測(cè)數(shù)據(jù)中的應(yīng)用效果有待進(jìn)一步驗(yàn)證。

        圖8 CSAMT測(cè)量L4線視電阻率斷面圖(預(yù)處理后)Fig.8 Apparent resistivity section of CSAMT survey of line 4(after preprocessing)

        4 處理實(shí)例及效果分析

        某高速公路連接線擬建路線穿越煤炭采空區(qū),為確保擬建項(xiàng)目的安全,在工作區(qū)內(nèi)開展CSAMT測(cè)量工作,查清采空區(qū)位置、分布范圍等特征。觀測(cè)裝置采用赤道偶極標(biāo)量方式,采集參數(shù)為Ex和Hy,發(fā)射偶極AB=1 200 m,測(cè)量偶極MN=50 m,最大收發(fā)距為6.4 km,測(cè)量頻段為1 Hz~8 192 Hz。由于工作場(chǎng)地為城市近郊,村莊廠房密集且有高速公路和多條高壓線穿過,電磁干擾非常嚴(yán)重。在施工過程中采取多種措施保證野外采集數(shù)據(jù)質(zhì)量,如加大供電電流、增加重復(fù)觀測(cè)次數(shù)、合理避開干擾時(shí)段等。

        以某測(cè)線1900點(diǎn)的數(shù)據(jù)散點(diǎn)及卡尼亞電阻率和阻抗相位曲線為例(圖6)。圖6(a)、圖6(b)分別為卡尼亞電阻率及阻抗相位實(shí)測(cè)數(shù)據(jù)散點(diǎn)及平均值曲線,由圖6(a)、圖6(b)可以看出,各個(gè)頻段的數(shù)據(jù)點(diǎn)均存在重復(fù)性差、離散度高等特點(diǎn),說明該測(cè)點(diǎn)在采集過程中受到較強(qiáng)烈電磁干擾,造成數(shù)據(jù)信噪比下降。中低頻段進(jìn)入“近區(qū)”,視電阻率曲線以近45°上揚(yáng),且在32 Hz處出現(xiàn)局部極小值,應(yīng)為過渡區(qū)反映。阻抗相位曲線也顯示了近區(qū)及過渡區(qū)信息。

        圖6(c)、圖6(d)為數(shù)據(jù)處理后數(shù)據(jù)散點(diǎn)及曲線特征。利用預(yù)處理程序,剔除離散度大的點(diǎn),保留重復(fù)性好、相對(duì)集中分布的點(diǎn)作為“合理”的實(shí)測(cè)數(shù)據(jù),曲線形態(tài)得到適當(dāng)調(diào)整,數(shù)據(jù)信噪比得到提高。采用處理后的數(shù)據(jù)進(jìn)行近場(chǎng)校正(圖6(c)紅色曲線),選取1 Hz~8 Hz為近區(qū)頻段,16 Hz及32 Hz為過渡區(qū)頻段,增加相應(yīng)的修正系數(shù),采用“全頻域視電阻率法”進(jìn)行計(jì)算,可以看出曲線形態(tài)發(fā)生較大變化,低頻段呈45°上揚(yáng)段斜率變小,過渡區(qū)電阻率抬升,高頻段視電阻率沒有發(fā)生變化,曲線整體形態(tài)更加合理,便于反演解釋。

        選取L4線進(jìn)行平滑濾波二維反演,對(duì)處理前后的數(shù)據(jù)進(jìn)行對(duì)比,以檢驗(yàn)該程序的處理效果。L4線東西向布設(shè),穿過煤炭采區(qū),沿線多為廠區(qū),并與高壓輸電線和鐵路相交,電磁干擾較強(qiáng)烈。采用原始數(shù)據(jù)進(jìn)行反演(圖7),視電阻率等值線“下拉”現(xiàn)象明顯,剖面后段約-400 m深度近水平展布的低阻等值線出現(xiàn)波狀起伏且埋深偏淺,新太古代泰山巖群片麻巖地層起伏形態(tài)也發(fā)生改變。尤為明顯的是石24井位置埋深-200 m處出現(xiàn)團(tuán)塊狀高阻異常區(qū),與鉆井揭示巖性不符,出現(xiàn)“假異常”。

        經(jīng)過“去噪”預(yù)處理及近場(chǎng)校正,斷面等值線形態(tài)更加合理,基本反映了真實(shí)的地質(zhì)信息,異常形態(tài)符合地質(zhì)規(guī)律,地電層劃分較明顯,其形態(tài)與已知地質(zhì)斷面基本一致。淺部為不連續(xù)的高阻特征,推測(cè)為青山群八畝地組凝灰質(zhì)砂礫巖電性反映,局部低阻異常推測(cè)為斷裂構(gòu)造發(fā)育部位。中部存在較連續(xù)水平展布的低阻異常帶,推測(cè)為淄博群坊子組含煤地層,低阻異常為采空區(qū)充水所致。深部則以高阻為主要特征,推測(cè)為新太古代泰山巖群片麻巖地層,局部低阻凹槽推測(cè)為斷裂構(gòu)造發(fā)育部位。后經(jīng)鉆探驗(yàn)證,煤炭采空區(qū)實(shí)際位置及埋深與物探解譯結(jié)果基本一致,L4線驗(yàn)證孔鉆遇采空區(qū)頂板深度與解譯深度誤差小于10%。

        5 結(jié)論

        多次的項(xiàng)目實(shí)踐表明,利用Matlab編程實(shí)現(xiàn)的CSAMT數(shù)據(jù)處理程序可以有效提高處理效率和效果,已經(jīng)成為GDP32多功能電法工作站CSAMT數(shù)據(jù)處理的有力工具。

        1)隨著數(shù)據(jù)采集環(huán)境的變化,尤其當(dāng)前城市地質(zhì)勘查活動(dòng)越來越多,CSAMT方法遇到的“電磁干擾”也會(huì)越來越嚴(yán)重,甚至造成數(shù)據(jù)無法利用。所以,為了實(shí)現(xiàn)勘查目的和任務(wù),既要完善野外采集環(huán)節(jié)的“抗噪”措施,也要通過室內(nèi)適當(dāng)?shù)摹叭ピ搿碧幚?,以提高?shù)據(jù)信噪比,克服“假異?!?,保證后續(xù)處理質(zhì)量。

        2)近場(chǎng)校正可以有效提高CSAMT低頻段數(shù)據(jù)的利用程度,加大勘查深度,改善地電斷面形態(tài),提高物探解譯的準(zhǔn)確性和精度。正演地電模型通過“全頻域視電阻率法”校正,得到了很好的效果。過渡區(qū)有“校正不足”現(xiàn)象,通過數(shù)據(jù)擬合得到的修正公式與頻率選擇有關(guān)系,在實(shí)際應(yīng)用的效果有待進(jìn)一步檢驗(yàn)。

        3)Matlab作為強(qiáng)大的數(shù)值計(jì)算工具,可以在多個(gè)領(lǐng)域發(fā)揮其顯著的作用,其豐富的函數(shù)庫(kù)為數(shù)據(jù)處理功能的實(shí)現(xiàn)提供了有效支撐,并可以實(shí)現(xiàn)模塊擴(kuò)展和多語言聯(lián)合編程。

        猜你喜歡
        近場(chǎng)頻域電阻率
        超大規(guī)模智能反射面輔助的近場(chǎng)移動(dòng)通信研究
        基于反射型超表面的近場(chǎng)聚焦研究
        一種基于PDV的近場(chǎng)沖擊波高壓測(cè)量技術(shù)
        頻域稀疏毫米波人體安檢成像處理和快速成像稀疏陣列設(shè)計(jì)
        三維電阻率成像與高聚物注漿在水閘加固中的應(yīng)用
        近場(chǎng)RCS測(cè)量不確定度分析
        基于改進(jìn)Radon-Wigner變換的目標(biāo)和拖曳式誘餌頻域分離
        一種基于頻域的QPSK窄帶干擾抑制算法
        隨鉆電阻率測(cè)井的固定探測(cè)深度合成方法
        基于頻域伸縮的改進(jìn)DFT算法
        给我播放的视频在线观看| 东京热加勒比视频一区| 久久久精品人妻一区二区三区免费| 亚洲中文字幕精品久久a| 一区二区三区在线视频观看| 国产精品国产三级久久| 肥老熟女性强欲五十路| 亚洲国产精品久久艾草| 中文字幕色av一区二区三区| 最近2019年好看中文字幕视频| 亚洲av中文无码乱人伦在线播放 | 中文字幕一区二区三区四区| 青青草国产手机观看视频| 色综合久久网| 欧美日韩在线视频一区| 日韩放荡少妇无码视频| 久久夜色精品国产噜噜av| 国模无码视频一区| 亚洲AV无码精品色欲av| 一本大道久久精品一本大道久久| 国产av一区仑乱久久精品| 亚洲午夜精品第一区二区| 亚洲youwu永久无码精品| 国产亚洲精品成人aa片新蒲金| 青青草97国产精品免费观看| 国产中文字幕乱码在线| 国产真实乱XXXⅩ视频| 亚洲一区二区三区1区2区| 精品人妻av区乱码色片| 午夜福利av无码一区二区| 亚洲av成人一区二区三区av| 亚洲嫩草影院久久精品| 最大色网男人的av天堂| 久久精品国产精品亚洲艾| 中文字字幕在线中文乱码解| 免费不卡在线观看av| 久久久久久久综合狠狠综合| 99久久精品国产亚洲av天| 精品久久一区二区av| 亚洲天堂av在线免费观看| 亚洲欧美v国产一区二区|