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

譜分析與啟發式遺傳算法相結合的多尺度社區檢測方法

2015-06-14 07:38:36郭玉泉李雄飛
吉林大學學報(工學版) 2015年5期
關鍵詞:結構檢測

郭玉泉,李雄飛,劉 昕

(1.吉林大學 計算機科學與技術學院,長春130012;2.吉林大學 網絡中心,長春130022)

0 引 言

目前,對于網絡社區還沒有統一定義,通常認為社區是一個內部連接緊密的節點集合,而該節點集合與網絡的其他部分連接稀疏。在現實網絡中,網絡社區表現出多尺度結構特征,即在不同尺度下有不同的社區劃分結果。常用網絡社區檢測方法可以分為三類:第一類是基于優化的方法,典型算法有GN[1]、FN[2]、LPAm[3]和LGA[4]方法。優化方法多以Newman等[5]提出的Q 函數為優化目標,而最大化Q 函數是NP 完全問題,所以這些優化算法都是近似的優化算法。第二類是層次聚類方法,典型算法有EAGLE[6]和CPM[7]。層次聚類方法中樹狀圖生成與分割的時間復雜度都很高,因此層次聚類方法效率較低。第三類是網絡動力學方法,如Jin等[8]提出的基于馬爾可夫方法的社區檢測方法。以上方法多數都是單一尺度的網絡社區結構檢測,不能檢測多尺度網絡社區結構。網絡社區結構的多尺度特征是通過網絡傳播特性與譜特征來體現的,目前檢測都是采用基于網絡動力學與譜方法的檢測算法。Alex等[9]利用同步過程與拉普拉斯矩陣譜檢測多尺度社區;Delvenne等[10]通過在時間尺度上社區穩定性檢測多尺度社區;Cheng等[11]通過復雜網中傳播過程的穩定狀態檢測多尺度社區。

為評估社區檢測結果的優劣,2004 年Newman[5]提出了網絡社區的度量標準,稱為模塊 化 函 數。隨 著 研 究 的 深 入,Fortunato[12]、Arenas[13]等發現模塊化函數存在分辨率限制(Resolution limit),無法發現特定尺寸的社區結構。2010年,Cheng 等[11]提出的傳導率函數(C函數)作為網絡社區評價標準,傳導率函數從網絡動力學角度評價網絡社區劃分結果,克服了模塊化函數的分辨率限制問題。

本文提出的HGASA 算法(Heuristic genetic algorithm with spectral analysis,HGASA)通過譜分析獲得復雜網絡多尺度結構信息,將結構信息應用于遺傳算法優化過程。算法從網絡動力學角度分析了個體變異過程,提出基于網絡動力學的變異操作啟發函數,用于指導變異操作,同時證明該啟發函數與優化目標函數(C函數)存在單調關系。

1 譜分析

矩陣的譜特性已經被廣大學者深入研究,將其應用到復雜網絡社區檢測中時,可以根據復雜網絡矩陣的譜特性揭示出復雜網絡的社區結構[14]。復雜網絡通常使用圖G=(V,E)來描述,其中V 表示網絡的節點集合,E 表示網絡中邊的集合。鄰接矩陣是圖的主要表示方法,當使用鄰接矩陣A 描述圖時,Aij定義如下:

鄰接矩陣有多種擴展形式,如標準拉普拉斯矩陣、規范化拉普拉斯矩陣、模塊化矩陣、關聯矩陣。本文只討論無權重無向網絡,采用規范化拉普拉斯矩陣[14]用于社區結構檢測。規范化拉普拉斯矩陣定義為N =I-T,其中I為單位陣,T=D-1A,D 為對角陣,Dii為節點i的度,A 為網絡鄰接矩陣。

定義1 假設λ1,…,λn是規范化拉普拉斯矩陣N 的n 個特征值,將特征值按遞增排序并去除重復的特征值得到一個遞增序列ω1,…,ωm(m <n),稱此序列為矩陣N 的譜,記為S(N)={ω1,…,ωm}。

定義2 S(N)為規范化拉普拉斯矩陣N 的譜,設ei=ωi+1-ωi,則ei稱為矩陣N 的第i個譜隙。

設定一個閾值ep,當譜隙大于此閾值(即ei>ep)時,認為ei是一個有效譜隙,實際應用中ep取值0.5[15]。將矩陣N 的有效譜隙按遞減順序排列,每個譜隙對應一個社區尺度,譜隙的下標值與網絡中社區的數量一致,由此得到了復雜網絡的多尺度社區結構。多尺度社區結構反映了網絡在傳播過程中的動力學特征[11],矩陣的譜分析方法可以有效地揭示出網絡的社區結構。

算法1 多尺度社區結構的算法

Input:網絡鄰接矩陣A,譜隙閾值ep,輸出列表中元素的數量m

Output:社區數量值列表L

1 T=D-1A

2 N=I-T

3 計算N 的特征值λ1…λn;

4 計算N 的譜隙e1…en-1;

5 對大于閾值ep的譜隙按遞減進行排序,得到譜隙序列EP;

6 取EP中前m 個元素的下標值放入到L中;

7 Return

2 HGASA 算法

2.1 編解碼

HGASA 算法采用字符串編碼方式[16],每個字符串代表一個個體,每個字符位置代表節點,每個字符代表節點所屬的社區。Ii=(Li1,Li2,…,Lin),Ii表示第i個個體,Lin表示第i個個體中第n個基因表達,這里表示節點n的社區標識符。由上面的定義知,如果Lip=Liq,則表明Ii所對應的社區結構中,節點p 與節點q 在同一個社區中。

2.2 群體初始化

遺傳算法在群體初始化時通常采用隨機生成個體的方式,這樣可以使初始群體具有多樣性。HGASA 算法中采用標簽傳播[17]方法生成個體,初始時個體每個節點分配一個社區標識符,然后經過數次異步標簽傳播產成個體。標簽傳播方法的隨機性保障生成個體的多樣性,同時具有很高的效率。但是標簽傳播方法生成個體的社區結構具有不確定性,使個體需要較長的進化時間才能達到目標狀態。標簽傳播算法生成的個體雖然具有一定的社區結構,但與社區檢測的目標有一定的差距。HGASA 算法利用矩陣譜分析的結果作為標簽傳播的約束條件,在標簽傳播過程中控制個體中標簽的數量與譜分析的結果一致,這樣個體生成時便具有了基本的社區結構,從而保障后續優化過程的效率和精度。

算法2 群體初始化算法

Input:個體包含社區的數量k,群體包含個體的數量m

Output:群體p

1 for i=1to m

2 生成個體I;

3 While(個體I社區標識數量>k)

4 隨機選擇個體I的一個基因位置,進行標簽傳播操作;

5 將I加入到群體p中;

6 Return p

2.3 交叉操作

采用遺傳算法解決社區檢測問題時,交叉操作中有一定的特殊性,進行交叉操作時要將一組關聯密切的基因交叉給下一代而不是個別基因,因此,不能采用傳統的單點、多點交叉方法,而是采用單路交叉、多路交叉。本文對Tasgin[16]提出的單路交叉算法進行了改進,稱為合并式單路交叉(One-way incorporating crossing)。改進的目的是使交叉后新個體保留兩個父母個體交叉位置的最大社區結構特征。合并式單路交叉操作過程定義如下:設A、B是任意兩個體,A 作為源個體,B作為目的個體,CB(n)為個體B 中節點n 所屬社區的節點集合,CA(n)為個體A 中與節點n 在同一社區節點的集合,L 為初始狀態(每個節點一個社區)社區標識符的集合,L(B)為個體B 中社區標識符的集合,LB(n)為個體B 中節點n 的社區標識符。進行交叉操作時,首先在個體A 中選擇一個節點作為交叉位置,設節點n 被選為交叉位置;然后,在個體B 中所有包含于CB(n)∪CA(n)中的節點進行交叉操作,對每個節點取一個不同于個體B 中現有標簽進行賦值,即LB(m)←L?m∈CB(n)∪CA(n)且L{L-LB}。

交叉操作如圖1 所示。在交叉后個體B 中標簽6所代表的社區結構繼承了個體A 與個體B的結構特征。

圖1 交叉操作Fig.1 Cross operation

合并單路交叉能保留上一代個體的社區結構特征,在新個體中社區結構特征更加突出。

2.4 變異操作

對于復雜網絡社區檢測問題,遺傳算法通常采用隨機方式變異,變異操作缺乏必要的啟發信息,對于變異操作的進化方向缺乏有效控制,從而使個體進化速度緩慢,優化效果不佳?;谏鲜鲈?,遺傳算法難以應用于大規模網絡的社區檢測。本節針對復雜網絡社區檢測過程,構建了局部動力學啟發函數,以啟發變異操作進化方向,并且證明啟發函數和HGASA 算法目標函數單調,從而提高個體進化速度和檢測精度。

HGASA 算法采用傳導率函數C(P)作為目標函數,其定義如下:

式中:P 為復雜網絡的一個社區檢測結果,P ={V1,…,Vk},其中k為社區個數,Vi表示第i個社區中的節點集合。

定義3 社區Vi的內部度in_vol(Vi)定義為

定義4 社區Vi的度vol(Vi)定義為vol(Vi)

從傳導率定義(式(2))可以看出,傳導率是每個社區分離概率的平均值,反映了復雜網絡中社區間傳播的能力,體現了復雜網絡的一種動力學特征。傳導率越小,社區結構劃分越合理。

命 題1 n 是 一 個 節 點,C 是 一 個 社 區,n ?C,n與C 中節點存在一條邊,當節點n社區標識符由l變化為l′時,社區C的內部度in_vol(C)不變化。

證明 如圖2所示,在圖2(a)中,n∈Vi,節點n與社區Vk中某一節點間存在一條邊。當節點n由社區Vi劃分到Vj后,即Vi-{n},Vj+{n},如圖2(b)所示,社區Vk的內部度與外部度都沒有變化,所以命題1成立,證畢。

圖2 節點社區變化示意圖Fig.2 Community of node

命題1闡明了節點社區標識變化和其相鄰社區內部度的關系,在變異操作中性質1表現為將一個節點由包含它的當前社區分離,然后融入到與節點連接更緊密的相鄰社區。從網絡動力學角度分析此過程,可以認為是“社區標識符”從鄰居社區傳遞到當前節點過程,即節點i的社區標識符由L(V(i))變為L(V(j)),其中L(V(i))表示包含節點i的社區,L(V(i))表示節點i的社區標識符。由于社區標識符傳播導致復雜網絡社區劃分的傳導率發生變化。社區標識符由一個社區傳播到相鄰社區邊緣節點時連接兩個社區間邊的數量發生變化,因此引入函數h(Ci,Cj)(簡稱h 函數)評估兩個社區之間社區標識符的傳播能力,其定義如下:

式中:Ci、Cj代表兩個相鄰的社區。

h(Ci,Cj)函數代表兩個相鄰社區的凝聚概率的平均值。當社區標識傳播達到穩定狀態時,社區標識符所表達的社區結構是網絡的社區結構。

命題2 當社區數量不變時,傳導率函數C(P)與h(Vi,Vj)單調遞減,P 為一個分區,Vi、Vj為相鄰的社區,Vi∈P,Vj∈P,P ={V1,…,Vk}。

證明 設P 為一個社區劃分,P ={V1,…,Vk},節點n∈Vi。社區標識符通過網絡傳播,節點n當前標識符為L(Vi),當社區標識符L(Vj)傳播至節點n時,引發社區劃分由P變為P′,P′={V′1,…,V′k}。社區標識符L(Vj)傳播到節點n后引起的函數h的增量如下:

由n 的社區標識變化引起的傳導率函數C(P)的增量如下:

根據命題1 除社區Vi與Vj外其他社區有in_vol(V′i)=in_vol(Vi),于是有:

綜上所述,C(P)與h(Ci,Cj)單調遞減,命題2成立。需要特別說明的是,命題2 中在社區標識符傳遞的過程中不討論社區標識符減少的情況,因為這種情況不滿足優化目標的要求。

從命題2 可知:傳導率函數C(P)與函數h(Vi,Vj)單調遞減,因此使用h(Vi,Vj)作為社區標識符傳播的啟發信息,當社區標識符向著h(Vi,Vj)增大的方向傳播時,C(P)將變小,此時得到的社區結構更優良。

當社區標識符在網絡上傳播時,它到達的下一個節點與它在同一社區的概率最大?;诖吮疚奶岢鋈缦碌淖儺惙椒ǎ簩τ趥€體I的第i個基因(即節點i被選擇進行變異),節點i的鄰接社區標識符集合為Ln(i),只需在Ln(i)中選擇一個社區標識符,不需要考慮節點i的所有鄰居節點的社區標識符,而是以h(Vi,Vj)作為啟發函數,選擇具有最大h 值的標簽傳播到當前節點。Li←arg maxL{h(Vi,VL)L ∈Ln(i)}。Li為節點i的社區標識符。

算法3 基于h函數的局部啟發式變異算法(Local heuristic mutation algorithm,LHMA)

Input:待變異個體Ind,變異基因位置Pos

Output:變異后個體Ind

1 for k=1to n //n為節點數量

2 List =NeiLabel(pos,Ind)//求節點相鄰社區標識符集合

3 m=length(List) //計算隊列長度

4 for i =1to m

5 t=h(Vpos,VList[i]) //計算h函數

6 if(t>max_h)

7 max_h=t

8 L=Label(VList[i])

9 將Ind中pos位置基因變為L

10 Return Ind

性質1 LHMA 的時間復雜度為O(Cn)。

證明 假設網絡的社區數量為k,節點的平均外度為dout,網絡節點數量為n,對于個體I發生變異節點的概率為β。

設計算h函數的平均時間為t,每個節點直接相連外部社區數量為dout,則每個節點計算所有h函數的時間的最大值為O(doutt),總的時間復雜度為O(βdouttn)。將tdoutβ看作常數C,則LHMA的時間復雜度為O(Cn)。

2.5 HGASA 框 架

HGASA 的流程圖如圖3所示。算法首先對復雜網絡的拉普拉斯矩陣進行譜分析,得到不同尺度下復雜網絡的社區數量;然后根據社區數量進行群體初始化,得到具有特定尺度社區結構的個體;最后使用局部啟發式變異算法與合并式單路交叉算法對群體進行優化,從而得到復雜網絡的社區結構。HGASA 算法選擇在所有個體進行完交叉或變異操作后進行,將交叉和變異產生的新個體加入群體中,然后在新群體中選擇傳導率最小的個體作為下一代種群。

圖3 HGASA算法流程圖Fig.3 Flow diagram of HGASA

算法4 HGASA

Input:鄰接矩陣A,變異概率β,進化代數L,種群數量I

Output:為社區檢測結果C

1 spectrumanalysis(A) //對矩陣進行譜分析

2 Initialization() //群體初始化

3 While(進化代數<L)

4 for i=1to I

5 If rand()>β

6 LHMA();//局部啟發式變異算法

7 Else

8 OWICA();//合并式單路交叉算法

9 Select();//選擇操作

10 C←社區檢測結果

11 Return C

性質2 HGASA 算法的時間復雜度為O(Mn)。

證明 步驟1的時間復雜性為O(kn)。步驟2對于每個個體經2~3次的標簽傳播可以得到滿足要求的個體,步驟2時間復雜性為O(2In),步驟5 至 步 驟8 時 間 復 雜 度 為O(LβCn)。所 以HGASA 算 法 時 間 復 雜 度 為O(kn +2In +LβCn),即O((k+2I+LβC)n),設k+2I+LC =M,所以算法時間復雜度為O(Mn)。

3 實 驗

采用人工網絡和現實網絡對HGASA 的性能進行測試。算法采用VC++6.0、Matlab7.0在Windows XP 下實現。VC++6.0實現遺傳算法部分,Matlab7.0實現矩陣分析。

HGASAd的主要參數有5 個,群體規模I、進化代數G、個體變異概率α、個體基因發生變異的概率β和隙閾值ep。前四個參數是遺傳算法常用參數。第五個參數是用于控制譜隙的有效值,一般取0.5。實驗采用召回率和精度來比較不同算法的分區結果。

3.1 人工網絡

目前人工生成網絡方法有多種,本文采用Lancichinetti[18]提出的方法進行人工網絡生成。這種方法可以根據節點的度、社區尺寸等多種方式生成社區,使生成的網絡能夠更深入地檢查社區檢測方法的性能?;旌下蔒ix(0≤Mix≤1)是Lancichinetti方法控制網絡生成的重要參數,它反映了社區結構的清晰程度,Mix越小網絡社區結構越清晰,反之社區結構模糊。選擇GCE[19]、CPM[7]、COPRA[17]、EAGALE[6]作 為 對 比 算 法。從圖4可以看出:0.1≤Mix≤0.2時,網絡社區結構明顯,對比算法和HGASA 都有較高的召回率和準確率,隨著Mix 增大網絡社區結構變得模糊,對比算法的準確率和召回率顯著降低,但HGASA 仍可較好地發現社區結構,說明HGASA 在檢測模糊社區結構方面性能明顯優于對比算法。

圖4 HGASA算法的召回率和準確率Fig.4 Recall rate and precision of HGASA

3.2 現實網絡

HGASA 對現實世界網絡的檢測選用空手道俱樂部網絡和海豚網絡來進行??帐值谰銟凡烤W絡(Karate)[20]展示了美國一所大學空手道俱樂部34名成員間的社會關系,每個節點代表一名成員,兩個節點間有一條邊則意味著相應的兩個成員之間是交往頻繁的朋友關系。用HGASA 對空手道俱樂部網絡進行檢測的結果如圖5、圖6所示。圖5是空手道俱樂部網絡的譜隙,e2、e4是兩個譜隙,對應的社區數量為2 和4。圖6 為HGASA 對空手道俱樂部網絡的檢測結果,此社區結構與現實觀察的檢測結果一致。通過實驗可以看到HGASA可以有效地檢測出Karate的多尺度社區結構。

圖5 Karate網絡譜隙Fig.5 Spectral of Karate

圖6 Karate網絡的多尺度社區結構Fig.6 Multiple-scale community structure of Karate

海豚網絡[21]是Lusseau通過對新西蘭神奇灣中62只不同性別海豚觀測而構建的動物社會網。網絡中每個節點代表一只海豚,如果兩個海豚聯系密切,那么海豚對應的頂點之間就會有一條邊連接。這些海豚被天然地分為雄性海豚組和雌性海豚組兩個社區。由圖7可以看到海豚網絡的e2、e5對應兩個社區結構。圖8(a)為海豚網絡對應e2的社區結構,海豚網絡被分成兩個大的社區,圖8(b)為海豚網絡對應e5的社區結構,可以看到在圖8(a)中右側社區又被細分為4個社區。

圖7 海豚網絡譜隙Fig.7 Spectra of dolphin

通過在人工網絡和現實網絡上對HGASA的測試可以看出,HGASA 具有較強的社區檢測能力,在社區結構不明顯時仍有較好的檢測結果,并且可以有效地檢測出多尺度網絡社區結構。

圖8 海豚網絡多尺度社區結構Fig.8 Multiple-scale community structure of dolphin

4 結束語

利用矩陣譜分析方法揭示出復雜網絡的多尺度特征,并且結合局部網絡動力學啟發函數提出了遺傳算法HGASA。計算機生成網絡和現實網絡的測試結果表明,HGASA 可有效地檢測多尺度網絡社區結構。在后續的研究中,作者將進一步研究多尺度社區結構與網絡功能演化的關系。

[1]Girvan M,Newman M E J.Community structure in social and biological networks[J].Proceedings of the National Academy of Sciences of the United States of America,2002,99(12):7821-7826.

[2]Newman M E J.Fast algorithm for detecting community structure in networks[J].Physical Review E,2004,69(6):066133.

[3]Barber M J,Clark J W.Detecting network communities by propagating labels under constraints[J].Physical Review E,2009,80(2):026129.

[4]Liu D Y,Jin D,Baquero C,et al.Genetic algorithm with a local search strategy for discovering communities in complex networks[J].International Journal of Computational Intelligence Systems,2013,6(2):354-369.

[5]Newman M E J.Detecting community structure in networks[J].European Physical Journal B,2004,38(2):321-330.

[6]Shen H W,Cheng X Q,Cai K,et al.Detect overlapping and hierarchical community structure in networks[J].Physica A:Statistical Mechanics and Its Applications,2009,388(8):1706-1712.

[7]Palla G,Derenyi I,Farkas I,et al.Uncovering the overlapping community structure of complex networks in nature and society[J].Nature,2005,435(7043):814-818.

[8]Jin D,Yang B,Baquero C,et al.A Markov random walk under constraint for discovering overlapping communities in complex networks[J].Journal of Statistical Mechanics-Theory and Experiment,2011(5):P05031.

[9]Alex Arenas,Albert Diaz-Guilera,Conrad J.Synchronization reveals topological scales in complex networks[J].Physical Review Letters,2006,96(11):114102.

[10]Delvenne J C,Yaliraki S N,Barahona M.Stability of graph communities across time scales[J].Proceedings of the National Academy of Sciences,2010,107(29):12755-12760.

[11]Cheng X Q,Shen H W.Uncovering the community structure associated with the diffusion dynamics on networks[J].Journal of Statistical Mechanics-Theory and Experiment,2010(4):P04024,

[12]Fortunato S,Barthelemy M.Resolution limit in community detection[J].Proceedings of the National Academy of Sciences,2007,104(1):36-41.

[13]Arenas A,Fernandez A,Gomez S.Analysis of the structure of complex networks at different resolution levels[J].New Journal of Physics,2008,10(5):053039.

[14]Shen H W,Cheng X Q.Spectral methods for the detection of network community structure:a comparative analysis[J].Journal of Statistical Mechanics-Theory and Experiment,2010(10):P10020,

[15]Shen H W,Cheng X Q,Wang Y Z,et al.A dimensionality reduction framework for detection of multiscale structure in heterogeneous networks[J].Journal of Computer Science and Technology,2012,27(2):341-357.

[16]Tasgin M,Herdagdelen A,Bingol H.Community detection in complex networks using genetic algorithms[EB/OL].[2013-07-08].http://arxiv.org/abs/0711.0491.

[17]Gregory S.Finding overlapping communities in networks by label propagation[J].New Journal of Physics,2010,12(10):103018

[18]Lancichinetti A,Fortunato S,Radicchi F.Benchmark graphs for testing community detection algorithms[J].Physical Review E,2008,78(4):046110.

[19]Lee C,Reid F,Mcdaid A,et al.Detecting highly overlapping community structure by greedy clique expansion[EB/OL].[2013-09-22].http://arxiv.org/abs/1002.1827.

[20]Zachary W.An information flow modelfor conflict and fission in small groups1[J].Journal of Anthropological Research,1977,33(4):452-473.

[21]Lusseau D.The emergent properties of a dolphin social network[J].Proceedings of the Royal Society of London Series B:Biological Sciences,2003,270(Suppl 2):186-188.

猜你喜歡
結構檢測
“不等式”檢測題
“一元一次不等式”檢測題
“一元一次不等式組”檢測題
《形而上學》△卷的結構和位置
哲學評論(2021年2期)2021-08-22 01:53:34
“幾何圖形”檢測題
“角”檢測題
論結構
中華詩詞(2019年7期)2019-11-25 01:43:04
新型平衡塊結構的應用
模具制造(2019年3期)2019-06-06 02:10:54
論《日出》的結構
小波變換在PCB缺陷檢測中的應用
主站蜘蛛池模板: 无码在线激情片| 嫩草在线视频| 一本二本三本不卡无码| 亚洲午夜福利精品无码不卡| 亚洲欧美精品日韩欧美| 日韩A∨精品日韩精品无码| 国产精品成人啪精品视频| 国产女同自拍视频| 无码日韩人妻精品久久蜜桃| 91精品国产无线乱码在线| 国产在线观看第二页| 欧美性久久久久| 亚洲天堂精品视频| 亚洲经典在线中文字幕| 国产精品尤物在线| 中文字幕日韩久久综合影院| 成人中文字幕在线| 国产精品v欧美| 国产乱人免费视频| 免费a在线观看播放| 日韩av在线直播| 免费啪啪网址| 亚洲中文字幕在线一区播放| 天堂成人av| 在线精品视频成人网| 亚洲成在线观看 | 色香蕉影院| 91色在线观看| 国产成人高清在线精品| 18禁黄无遮挡网站| 婷婷伊人五月| 欧美性色综合网| 99999久久久久久亚洲| 国产麻豆精品在线观看| 婷婷五月在线| 国产一在线观看| 国产青榴视频| 国产成人精品一区二区秒拍1o| 色综合手机在线| 国产精品手机在线观看你懂的| 日韩欧美在线观看| 激情五月婷婷综合网| 日韩欧美中文| 中国国产A一级毛片| 黄色网站在线观看无码| 99久久精品久久久久久婷婷| 国产欧美日韩另类精彩视频| 久久黄色免费电影| 国产欧美日韩资源在线观看 | 亚洲国产日韩一区| 综合五月天网| 欧美激情视频一区| 一区二区三区成人| 永久在线精品免费视频观看| 成人国产一区二区三区| 国产视频 第一页| 国产精品分类视频分类一区| a毛片免费在线观看| 嫩草影院在线观看精品视频| 久久国产精品麻豆系列| 久久黄色毛片| 2024av在线无码中文最新| 国产小视频在线高清播放| 国产尤物jk自慰制服喷水| 亚洲人视频在线观看| 露脸国产精品自产在线播| 91久久性奴调教国产免费| 最新国产你懂的在线网址| 色综合天天娱乐综合网| 亚洲日韩AV无码精品| 国产乱子精品一区二区在线观看| 亚洲精品男人天堂| 丁香五月婷婷激情基地| 国产一级在线观看www色| 日韩高清中文字幕| 亚洲天堂.com| 亚洲丝袜中文字幕| 国产99视频精品免费观看9e| 日韩免费毛片| 极品国产在线| 国产精品久久久久久久久| 久久精品视频一|