999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

基于近似熵的斯隆數(shù)字化巡天中類(lèi)星體光變復(fù)雜性分析*

2019-10-23 01:22:40唐潔劉曉琴
物理學(xué)報(bào) 2019年14期
關(guān)鍵詞:分析方法

唐潔 劉曉琴

1)(陜西理工大學(xué)物理與電信工程學(xué)院,漢中 723001)

2)(漢中職業(yè)技術(shù)學(xué)院藥學(xué)與醫(yī)學(xué)技術(shù)系,漢中 723002)

光變是類(lèi)星體的重要觀測(cè)特征之一,類(lèi)星體在多個(gè)波段存在劇烈的光變現(xiàn)象.光變非常復(fù)雜,具有非線(xiàn)性特征.以斯隆數(shù)字化巡天(Sloan digital sky survey,SDSS)stripe 82天區(qū)中的類(lèi)星體為研究對(duì)象,利用近似熵方法分析了類(lèi)星體光變的復(fù)雜性.首先應(yīng)用模擬信號(hào)檢驗(yàn)了近似熵方法對(duì)周期序列、白噪聲序列、混沌序列和組合序列的區(qū)分能力,驗(yàn)證了近似熵方法是一種識(shí)別不同類(lèi)型時(shí)間序列的有效方法.再計(jì)算了SDSS第7次釋放數(shù)據(jù)中光譜證認(rèn)過(guò)的類(lèi)星體光變的近似熵,并分析了它們的復(fù)雜性.結(jié)果表明: SDSS類(lèi)星體光變的近似熵值最大值為0.58,類(lèi)星體光變的復(fù)雜性介于周期序列和白噪聲序列的復(fù)雜性之間,近一半樣本的復(fù)雜性與混沌序列基本一致.

1 引 言

作為活動(dòng)性最強(qiáng)活動(dòng)星系核,類(lèi)星體除了強(qiáng)發(fā)射線(xiàn)、較大紅移、高光度等觀測(cè)特征外,還具有快速、不規(guī)則、劇烈的光變現(xiàn)象[1,2],從長(zhǎng)波的射電波段到高能的伽瑪波段都觀測(cè)到較強(qiáng)的光變信號(hào).對(duì)光變的研究可以為我們提供許多重要信息,有利于更深入地了解光變產(chǎn)生的物理機(jī)制.獲得了很多觀測(cè)現(xiàn)象支持的吸積盤(pán)模型認(rèn)為[3],吸積盤(pán)結(jié)構(gòu)的改變,吸積盤(pán)自身的不穩(wěn)定性和變化的吸積率都有可能導(dǎo)致類(lèi)星體光變現(xiàn)象的發(fā)生.類(lèi)星體的光變資料主要來(lái)源于長(zhǎng)時(shí)間的多波段測(cè)光監(jiān)測(cè),目前已經(jīng)積累了大量的觀測(cè)資料,為進(jìn)一步研究類(lèi)星體光變提供了便利條件.

類(lèi)星體被發(fā)現(xiàn)后不久,其輻射流量的變化就被觀測(cè)到了.已有研究對(duì)類(lèi)星體光變現(xiàn)象從不同角度做了大量的分析.MacLeod等[4]認(rèn)為類(lèi)星體光變可以利用隨機(jī)游走模型來(lái)擬合,但Kasliwal等[5]卻認(rèn)為該模型并非適合所有的類(lèi)星體,Guo等[6]認(rèn)為隨機(jī)游走模型不能極好的擬合部分類(lèi)星體光線(xiàn)曲線(xiàn).解釋光變現(xiàn)象發(fā)生常用的吸積盤(pán)不穩(wěn)定模型也是將光變視為隨機(jī)事件,盤(pán)通過(guò)黏滯耗散和隨機(jī)因素同外部介質(zhì)產(chǎn)生作用[3,7].隨著觀測(cè)技術(shù)水平的提高和一些大型望遠(yuǎn)鏡的使用,在對(duì)一些長(zhǎng)期監(jiān)測(cè)的類(lèi)星體光變資料分析時(shí)發(fā)現(xiàn)光變呈顯出準(zhǔn)周期現(xiàn)象[8,9],如OJ 287長(zhǎng)期光變曲線(xiàn)存在約12年的長(zhǎng)周期[10,11],運(yùn)用該周期推測(cè)OJ 287下一次可能爆發(fā)的時(shí)間,其后觀測(cè)到的爆發(fā)時(shí)間與該預(yù)測(cè)時(shí)間非常吻合[12,13].隨著非線(xiàn)性動(dòng)力學(xué)理論的發(fā)展,非線(xiàn)性動(dòng)力學(xué)分析方法被引入到類(lèi)星體光變的研究中來(lái)[14-16].如果光變具有混沌特性,非線(xiàn)性機(jī)制就可以被用來(lái)解釋類(lèi)星體光變.

類(lèi)星體光變資料里隱含有線(xiàn)性和非線(xiàn)性的成分,是極其復(fù)雜的[17].將不規(guī)則的光變現(xiàn)象視為復(fù)雜的,僅是一種假設(shè),簡(jiǎn)單的定性分析缺乏說(shuō)服力.類(lèi)星體光變是不是屬于復(fù)雜系統(tǒng)? 若其存在復(fù)雜性,復(fù)雜到何種程度? 這些問(wèn)題值得進(jìn)一步深入研究.復(fù)雜性研究處于有序與混沌的邊緣,現(xiàn)在主要指的是對(duì)非線(xiàn)性現(xiàn)象的研究[18].混沌動(dòng)力學(xué)理論可以進(jìn)行復(fù)雜性分析,但混沌分析方法一般要求數(shù)據(jù)量非常大,少則上千,多則上萬(wàn),對(duì)觀測(cè)資料的質(zhì)量要求也較高,要想獲得較完整的類(lèi)星體光變的資料需要占用望遠(yuǎn)鏡較長(zhǎng)的觀測(cè)時(shí)間.盡管與過(guò)去相比,目前類(lèi)星體的監(jiān)測(cè)項(xiàng)目增多,但監(jiān)測(cè)樣本數(shù)還是較少,觀測(cè)時(shí)間的跨度也較短.

斯隆數(shù)字化巡天(Sloan digital sky survey,SDSS)項(xiàng)目提供了大量類(lèi)星體測(cè)光監(jiān)測(cè)數(shù)據(jù),SDSS第7次釋放數(shù)據(jù)(data release 7,DR7)釋放有光譜證認(rèn)過(guò)的類(lèi)星體5個(gè)波段10年觀測(cè)的數(shù)據(jù),但單個(gè)類(lèi)星體樣本光變資料數(shù)據(jù)量少.常用的、對(duì)時(shí)間序列長(zhǎng)度要求相對(duì)較低的復(fù)雜性分析方法有近似熵[19,20]、Lempel-Ziv 算法[21]、樣本熵[22]和Kolmogrov熵[23]等方法.超過(guò)100個(gè)數(shù)據(jù)量時(shí),利用近似熵方法就能獲得比較可靠的估計(jì)值,并且沒(méi)有粗粒化處理要求,這樣可以保留原始數(shù)據(jù)的有效信息[19],近似熵方法優(yōu)于后兩種方法.本文利用近似熵方法計(jì)算和度量類(lèi)星體光變的復(fù)雜性,分析類(lèi)星體光變的復(fù)雜度有助于我們更好地了解光變的復(fù)雜特征,揭示復(fù)雜光變的非線(xiàn)性動(dòng)力學(xué)特性.

2 近似熵方法

混沌動(dòng)力學(xué)理論常用Kolmogrov熵來(lái)度量動(dòng)力學(xué)系統(tǒng)的復(fù)雜性,要獲得Kolmogrov熵值,理論上要求被分析的對(duì)象數(shù)據(jù)量足夠多,且該方法抗噪能力差,對(duì)于信噪比低的時(shí)間序列較難獲得可靠的熵值,因此最好別疊加有噪聲,但實(shí)測(cè)數(shù)據(jù)很難達(dá)到這些要求.Pincus和Huang[19]提出了對(duì)數(shù)據(jù)量長(zhǎng)度要求相對(duì)較低,有較好的抗干擾、抗噪能力的近似熵 (approximate entropy,ApEn).不同于Kolmogrov熵,近似熵不需要利用大量數(shù)據(jù)重構(gòu)奇異吸引子,主要是從統(tǒng)計(jì)角度來(lái)度量時(shí)間序列產(chǎn)生新模式的概率大小,概率越大復(fù)雜性越大,近似熵值也越大,隨機(jī)成分增多; 近似熵值越小,時(shí)間序列周期成分增多.

若存在一個(gè)給定長(zhǎng)度的一維時(shí)間序列X={x(n),n=1,2,···,N},如模式維數(shù)m和相似容限r(nóng)已經(jīng)確定,構(gòu)建一個(gè)m維的向量空間:

將任意向量X(i)與其他向量X(j)的對(duì)應(yīng)元素差值最大的記為

對(duì)每個(gè)向量X(i)計(jì)算Dij小于相似容限r(nóng)的數(shù)目,將該數(shù)目與向量總數(shù)N-m的比記為將模式維數(shù)增加1,按上面的方法計(jì)算維數(shù)m+1的時(shí)間序列的近似熵Sapen可表示為[19]

3 典型時(shí)間序列的近似熵

為了檢驗(yàn)近似熵對(duì)于不同時(shí)間序列的區(qū)分能力,選取周期序列、混沌序列、白噪聲序列和組合序列,分析5種典型時(shí)間序列的近似熵,通過(guò)分析典型不同類(lèi)型時(shí)間序列在復(fù)雜度的取值范圍,可以作為類(lèi)星體光變屬于哪種類(lèi)型時(shí)間序列的參照.

5種典型時(shí)間序列如圖1所示,其中周期序列(用Sine表示)為y=sin(2πt/10);混沌序列(用Logistic表示)為x(n+1)=ωx(n)[1-x(n)],其中初始值x(0)=0.7,ω=0.75,由該方程產(chǎn)生混沌序列;白噪聲序列(用White noise表示)由均值為0、方差為1的隨機(jī)序列組成;組合序列包括混沌序列疊加周期序列(用Logistic+Sine表示)和混沌序列疊加白噪聲序列(用Logistic+White noise表示).

對(duì)于時(shí)間序列長(zhǎng)度的選取,Pincus和Huang[19]通過(guò)他們的大量研究,建議數(shù)據(jù)點(diǎn)數(shù)最好取100以上,一般也不要超過(guò)5000.我們將時(shí)間序列長(zhǎng)度都取900,模式維數(shù)m=2,相似容限r(nóng)取0.05—0.50,5種典型時(shí)間序列的近似熵分析結(jié)果見(jiàn)圖2.從近似熵值分布可以看出: 周期序列的近似熵最小,白噪聲序列的近似熵最大,混沌序列的近似熵值介于周期序列和白噪聲序列之間,混沌序列疊加周期序列比混沌序列的近似熵值略大,相似容限小于0.35時(shí),混沌序列疊加白噪聲序列比白噪聲序列的近似熵值小; 相似容限大于0.35時(shí),混沌序列疊加白噪聲序列比白噪聲序列的大.對(duì)于不同類(lèi)型的時(shí)間序列,近似熵值按周期序列、混沌序列、混沌序列疊加周期序列、混沌序列疊加白噪聲序列到白噪聲序列的順序逐漸增大,近似熵值越大復(fù)雜性也越高,周期序列、混沌序列、白噪聲序列的近似熵值差別很大,各自分別在不同的區(qū)間,白噪聲序列是混沌序列的近似熵值近5倍,白噪聲序列和混沌序列疊加白噪聲序列復(fù)雜度差別不大,白噪聲序列近似熵值接近組合序列,近似熵值都超過(guò)1.

這些結(jié)果說(shuō)明近似熵能有效將混沌序列、周期序列、白噪聲序列和不同成分組合序列區(qū)分開(kāi),復(fù)雜性能夠區(qū)分不同隨機(jī)程度的時(shí)間序列,驗(yàn)證了近似熵可以描述時(shí)間序列隨機(jī)程度,對(duì)混沌有一定的識(shí)別能力,較好地表征了非線(xiàn)性結(jié)構(gòu)的復(fù)雜性,且復(fù)雜度隨著隨機(jī)成分增加而增大.

圖1 典型時(shí)間序列Fig.1.Typical time series.

圖2 典型時(shí)間序列的近似熵Fig.2.The approximate entropy of typical time series.

4 SDSS類(lèi)星體光變的復(fù)雜性分析

SDSS項(xiàng)目使用2.5 m專(zhuān)用望遠(yuǎn)鏡,裝備了先進(jìn)的大視場(chǎng)掃描成像裝置,觀測(cè)底板的640條光纖可以同時(shí)拍攝640個(gè)目標(biāo),從2000年項(xiàng)目開(kāi)始巡天觀測(cè)后,已經(jīng)積累了許多高質(zhì)量的u,g,r,i,z波段測(cè)光數(shù)據(jù).

為了展示近似熵方法在度量SDSS類(lèi)星體光變復(fù)雜性的能力,我們先以類(lèi)星體SDSS J012228.72-001332.0為例來(lái)分析該方法的可行性.由于SDSS有很多觀測(cè)任務(wù),望遠(yuǎn)鏡觀測(cè)時(shí)間有限,不可能長(zhǎng)時(shí)間觀測(cè)某個(gè)類(lèi)星體,這些類(lèi)星體大多是間隔90—270 d才有測(cè)光數(shù)據(jù).為了便于計(jì)算近似熵,缺失值用插值方法求得,圖3給出了類(lèi)星體SDSS J012228.72-001332.0原始和插值處理后的5個(gè)波段光變曲線(xiàn).

在計(jì)算近似熵的值前,需要先選取模式維數(shù)m、相似容限r(nóng)和時(shí)間序列的長(zhǎng)度N.對(duì)于N的選取,Pincus和Huang[19]通過(guò)他們的大量研究,建議數(shù)據(jù)點(diǎn)數(shù)最好取100以上,一般也不要超過(guò)5000,本文對(duì)原觀測(cè)數(shù)據(jù)每隔10 d插值一個(gè)數(shù)據(jù),數(shù)據(jù)量完全滿(mǎn)足條件.對(duì)于模式維數(shù)m,一般選取m=2,相對(duì)于m=1,可以在時(shí)間序列的聯(lián)合概率重構(gòu)過(guò)程中,獲得更多更詳細(xì)的信息; 若取m大于2,相似容限r(nóng)要求比較大,但這樣會(huì)導(dǎo)致丟失許多信息.對(duì)于參數(shù)相似容限r(nóng),Pincus和Huang[19]認(rèn)為相似容限r(nóng)取0.10—0.25倍原始數(shù)據(jù)的標(biāo)準(zhǔn)差,可以估計(jì)出比較理想的統(tǒng)計(jì)特性.

圖4給出了類(lèi)星體SDSS J012228.72-001332.0光變資料的相似容限-近似熵圖,在相似容限為0.2前,近似熵值隨相似容限的增加而增大,相似容限r(nóng)在增為0.2倍原始數(shù)據(jù)的標(biāo)準(zhǔn)差后,5個(gè)波段光變資料的近似熵值基本保持一個(gè)比較穩(wěn)定的值.g波段近似熵最小,在0.3左右振蕩; u波段近似熵最大,最大值接近0.35,相似容限在區(qū)間0.4—0.6時(shí),近似熵值隨相似容限增加有減小的趨勢(shì).因此,在本文后面的大樣本近似熵值分析時(shí),將相似容限r(nóng)取為0.2倍原始數(shù)據(jù)的標(biāo)準(zhǔn)差.從單個(gè)樣本近似熵分析可以看出,5個(gè)波段光變資料的近似熵值相差不多,同一類(lèi)星體不同波段光變曲線(xiàn)復(fù)雜性基本相同,這可能是因?yàn)樗鼈兊墓庾兦€(xiàn)波形非常相似.

本文選取的類(lèi)星體樣本來(lái)源于SDSS DR7中的stripe 82天區(qū)光譜證認(rèn)過(guò)的類(lèi)星體,Macleod等[4]搜集整理提供了9258個(gè)類(lèi)星體測(cè)光數(shù)據(jù),為了獲得可靠的分析結(jié)果,要求每隔360 d需要至少一次的觀測(cè)數(shù)據(jù),由于2000—2003年觀測(cè)數(shù)據(jù)太少,起始時(shí)間選擇從2003年開(kāi)始,并除去測(cè)光誤差大于0.1個(gè)星等的數(shù)據(jù)點(diǎn),并且要求這些數(shù)據(jù)大致均勻分布,按這些條件再選擇觀測(cè)數(shù)據(jù)點(diǎn)數(shù)超過(guò)20的作為樣本.最后,我們得到i,r,g,z,u波段的樣本數(shù)分別是6465,6547,6439,6373和6092.

圖3 SDSS J012228.72-001332.0的光變曲線(xiàn)Fig.3.Light curve of SDSS J012228.72-001332.0.

將模式維數(shù)和相似容限分別取m=1,r=0.2對(duì)選取的樣本進(jìn)行近似熵分析,獲得的結(jié)果如圖5和表1所示.可以看出,類(lèi)星體光變的近似熵值都處于0—0.6范圍內(nèi),最大值分別是0.552,0.553,0.575,0.571和0.580,最小值都為0.008,接近于0.白噪聲序列和混有白噪聲序列的混沌序列都比光變的近似熵最大值超過(guò)2倍,因此類(lèi)星體樣本光變不可能完全是白噪聲序列,這是因?yàn)轭?lèi)星體觀測(cè)總會(huì)有噪聲干擾,但其占的比例不大.近一半的近似熵值分布在0.3—0.4的范圍內(nèi),遠(yuǎn)大于周期序列,又遠(yuǎn)小于白噪聲序列,和混沌序列在初始值取為 x (0)=0.7,ω=0.75 時(shí)的近似熵值非常接近,說(shuō)明近一半的類(lèi)星體光變既不是完全周期性的也不是完全白噪聲的系統(tǒng),而是一個(gè)與混沌序列近似的、帶有非線(xiàn)性結(jié)構(gòu)的、具有混沌特性的動(dòng)力學(xué)系統(tǒng).

圖4 SDSS J012228.72-001332.0的近似熵Fig.4.The approximate entropy of SDSS J012228.72-001332.0.

近似熵值的大小和被分析的時(shí)間序列復(fù)雜程度存在成正比例的關(guān)系,時(shí)間序列越復(fù)雜,其近似熵值也就越大,時(shí)間序列越趨近于隨機(jī)性,反之,近似熵值越小,越趨于周期性.類(lèi)星體光變的近似熵值低于0.2的樣本數(shù)很少,最多的r波段才占到7.3%,最少的u波段只有3.56%,說(shuō)明周期性光變占主導(dǎo)的類(lèi)星體比例很少.近似熵值高于0.5的類(lèi)星體樣本數(shù)也少,除z波段剛超過(guò)1%,其余波段1%都不到,說(shuō)明所有類(lèi)星體光變的近似熵值都偏小,類(lèi)星體光變趨于隨機(jī)性概率較少.同一類(lèi)星體樣本不同波段,近似熵值都比較接近,說(shuō)明同一類(lèi)星體光變的復(fù)雜性基本一致.不同類(lèi)星體近似熵值不盡相同,說(shuō)明類(lèi)星體光變的非線(xiàn)性結(jié)構(gòu)是有區(qū)別的,可能光變產(chǎn)生的原因不一樣引起的.

圖5 SDSS類(lèi)星體的近似熵分布Fig.5.Approximate entropy distribution of SDSS quasars.

表1 SDSS類(lèi)星體的近似熵Table 1.Approximate entropy of SDSS quasars.

5 結(jié) 論

自從類(lèi)星體發(fā)現(xiàn)以來(lái),劇烈的光變現(xiàn)象引起了天文學(xué)家的廣泛關(guān)注.盡管現(xiàn)在觀測(cè)技術(shù)有很大進(jìn)步,但現(xiàn)有的觀測(cè)設(shè)備仍然無(wú)法分辨類(lèi)星體的中心區(qū)域,對(duì)類(lèi)星體光變的深入研究有助于了解其中心區(qū)域輻射的物理機(jī)制.常用分析類(lèi)星體光變資料的方法如結(jié)構(gòu)函數(shù)[24]、功率譜[25]、經(jīng)驗(yàn)?zāi)B(tài)分解[26]等,應(yīng)用這些方法如要獲得可信度高的結(jié)論,對(duì)觀測(cè)數(shù)據(jù)量要求較高,這就對(duì)數(shù)據(jù)分析方法提出了要求.本文采用的近似熵只需要100個(gè)以上的數(shù)據(jù)點(diǎn)就能獲得比較可靠的復(fù)雜度,完全滿(mǎn)足數(shù)據(jù)分析的需要.

本文利用近似熵方法分析了混沌序列、周期序列、白噪聲序列和不同成分組合序列的近似熵,發(fā)現(xiàn)近似熵值能有效區(qū)分不同隨機(jī)程度的時(shí)間序列,驗(yàn)證了近似熵對(duì)不同復(fù)雜性的時(shí)間序列有較好的識(shí)別能力,金寧德等[20]也用模擬信號(hào)獲得了跟我們一樣的結(jié)論.本文進(jìn)一步基于SDSS stripe 82天區(qū)的類(lèi)星體測(cè)光數(shù)據(jù)分析了類(lèi)星體光變的復(fù)雜性,從單個(gè)源的分析結(jié)果看來(lái),5個(gè)波段光變曲線(xiàn)波形非常相似,它們的近似熵值相差較小,同一類(lèi)星體不同波段光變曲線(xiàn)復(fù)雜性基本相同.對(duì)6000多個(gè)大樣本的類(lèi)星體分析發(fā)現(xiàn),所有的近似熵值都低于0.6,最大值為0.58,說(shuō)明類(lèi)星體光變中白噪聲占的比例都不高.近一半的類(lèi)星體光變和Logistic方程產(chǎn)生的混沌序列的復(fù)雜性是一樣的,說(shuō)明一部分類(lèi)星體光變以混沌成分占主導(dǎo).小部分近似熵值接近于0,這說(shuō)明周期性成分占主導(dǎo)的類(lèi)星體比例非常少.類(lèi)星體光變不可能是完全周期性的,也不可能是白噪聲,因?yàn)槿羰峭耆芷谛曰虬自肼?近似熵值是0或超過(guò)1,也不可能是完全混沌的,因?yàn)橛^測(cè)中總夾雜著噪聲的影響.因此類(lèi)星體光變可能是周期性成分,混沌成分和噪聲疊加在一起,并以3種成分中的一種占主導(dǎo).

用結(jié)構(gòu)函數(shù)、功率譜等周期分析方法分析類(lèi)星體長(zhǎng)期光變資料,發(fā)現(xiàn)光變具有周期性[1,8-13],這和我們的結(jié)論并不矛盾,這些文獻(xiàn)只提取了類(lèi)星體光變中隱含的部分周期信息,忽略了其他部分的有效信息.MacLeod等[4]用隨機(jī)游走模型來(lái)描述光變,發(fā)現(xiàn)該模型能較好地?cái)M合類(lèi)星體光變曲線(xiàn),但Kasliwal等[5]和Guo等[6]認(rèn)為隨機(jī)游走模型在描述類(lèi)星體光變中還是存在部分缺陷.Misra等[14]運(yùn)用混沌動(dòng)力學(xué)理論分析了GRS 1915+105 的光變曲線(xiàn),獲得的結(jié)果是GRS 1915+105光變也存在混沌性疊加隨機(jī)性的現(xiàn)象.唐潔[17]應(yīng)用集合經(jīng)驗(yàn)?zāi)B(tài)分解方法將類(lèi)星體光變資料分解成周期項(xiàng)、混沌項(xiàng)和趨勢(shì)項(xiàng),認(rèn)為類(lèi)星體光變是由周期成分、混沌成分和趨勢(shì)成分疊加而成.這些已有的研究成果與我們分析的結(jié)論是比較一致的.

感謝美國(guó)SDSS項(xiàng)目提供的類(lèi)星體光變數(shù)據(jù).

猜你喜歡
分析方法
隱蔽失效適航要求符合性驗(yàn)證分析
學(xué)習(xí)方法
電力系統(tǒng)不平衡分析
電子制作(2018年18期)2018-11-14 01:48:24
電力系統(tǒng)及其自動(dòng)化發(fā)展趨勢(shì)分析
用對(duì)方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢(qián)方法
捕魚(yú)
中西醫(yī)結(jié)合治療抑郁癥100例分析
在線(xiàn)教育與MOOC的比較分析
主站蜘蛛池模板: 亚洲第一区精品日韩在线播放| 在线综合亚洲欧美网站| 成人中文在线| 国产黄在线免费观看| 国产小视频免费| 亚洲第一区在线| 国产日韩欧美一区二区三区在线| 精品久久高清| 不卡网亚洲无码| 高清色本在线www| 亚洲成A人V欧美综合天堂| 99re这里只有国产中文精品国产精品 | 高清精品美女在线播放| 成人毛片免费观看| 国产黑人在线| 国产区精品高清在线观看| 波多野结衣AV无码久久一区| 日韩欧美中文| 亚洲成人一区在线| 男女男免费视频网站国产| 色妞永久免费视频| 91蝌蚪视频在线观看| 秘书高跟黑色丝袜国产91在线| 成人亚洲天堂| 国产精品极品美女自在线网站| 91小视频在线观看免费版高清| 欧美一区国产| 国产一区二区三区视频| 97综合久久| 高清色本在线www| 波多野结衣一区二区三区四区| 成人国产精品2021| 孕妇高潮太爽了在线观看免费| 国内精品伊人久久久久7777人| 伊人久久大香线蕉影院| 国产丰满大乳无码免费播放| 亚洲天堂区| 少妇精品网站| 九色视频最新网址| 国产精品久久久久久久久| www.亚洲一区二区三区| 3344在线观看无码| 国产裸舞福利在线视频合集| 亚洲国产理论片在线播放| 伊人久久福利中文字幕| 久久久久久高潮白浆| 成人在线不卡| 中文字幕一区二区人妻电影| 免费午夜无码18禁无码影院| 88av在线播放| 免费观看国产小粉嫩喷水| 天堂成人av| 亚洲欧美成人综合| 国产亚洲精品无码专| 久久精品aⅴ无码中文字幕| 日本不卡视频在线| 真人免费一级毛片一区二区| 91小视频在线观看| 中文字幕在线不卡视频| 九九视频免费在线观看| 国产无遮挡猛进猛出免费软件| 孕妇高潮太爽了在线观看免费| 亚洲色图另类| 久久国产精品夜色| 2021国产精品自产拍在线| 天天爽免费视频| 国产毛片网站| 婷婷丁香色| 国产在线自揄拍揄视频网站| 伊人久久精品亚洲午夜| 精品天海翼一区二区| 亚洲久悠悠色悠在线播放| 亚洲精品在线观看91| 成人夜夜嗨| 久久国产拍爱| 欧美黄色a| 18禁黄无遮挡免费动漫网站| 亚洲青涩在线| 欧美三级不卡在线观看视频| 亚洲乱码在线播放| 亚洲乱码视频| 国产小视频免费|