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

基于IRS輔助多用戶通信系統的信道容量優化

2022-11-19 06:53:38劉金枝梅志強梁家敏
系統工程與電子技術 2022年12期
關鍵詞:優化用戶

王 丹, 劉金枝, 梅志強, 梁家敏

(重慶郵電大學通信與信息工程學院, 重慶 400065)

0 引 言

隨著5G無線通信網絡的商業化,人們渴望獲得更快的數據傳輸速度。大規模多輸入多輸出(multiple input multiple output, MIMO)、毫米波通信等5G關鍵技術為實現此功能做出了有效貢獻。然而,額外的高硬件成本,巨大的功耗/能源消耗以及5G基站的選址是其在實踐中實施的主要障礙[1-2]。因此,為了在5G和無線網絡之外實現綠色和可持續發展,尋找頻譜和高能效技術的研究對于可持續容量增長仍然至關重要。為了解決上述挑戰,在5G和下一代移動通信系統中,智能反射面(intelligent reflecting surface,IRS)被認為是一種有前途的、綠色、高成本效益的技術。IRS是由電磁材料組成的人造表面[3]。具體來說,IRS可以通過在表面上集成大量低成本的無源反射元件來智能地調整無線傳播環境[4]。此外,IRS外形小巧,重量輕,可以輕松地從墻壁,天花板,廣告面板甚至衣服上安裝/拆卸[5-7],為實際實施提供了很高的靈活性。并且,與傳統的多用戶通信系統相比,多天線基站若使用線性預編碼器去同時服務不同的用戶,這比單天線的基站極大地提高了頻譜效率,但是當不同的用戶靠得比較近時,會產生較嚴重的共信道干擾,那么傳統的線性預編碼器對消除用戶間的干擾無效[8]。但是,通過有效地布置IRS可以解決上述問題,因為IRS通過優化其反射系數來提供附加的控制信號路徑,可以有效地減少多用戶之間不期望的信道干擾[9]。因此,IRS具有良好的干擾相消能力。基于上述優點,IRS在未來的無線網絡中具有廣闊的前景。

針對先前有關IRS輔助無線通信系統中IRS反射優化的工作。在文獻[7]中,針對IRS輔助的多用戶多輸入單輸出(multi input single output, MISO)通信系統,通過利用半正定松弛(semi-definite relaxation, SDR)方法和交替優化技術來最大化每個用戶的接收功率。但是,SDR方法不僅具有高復雜度,并且只能獲得近似解。此外,該系統中IRS上的每個反射元素都具有連續的相移,由于硬件的限制,這會使制造成本高昂,甚至無法實際應用。在文獻[10]中,作者在設計資源分配時采用了梯度下降搜索和順序分數規劃(sequential fraction programming, SFP),以最大程度地同時提高能量或頻譜效率。然而,IRS引起的單位模量約束是通過SFP方法中的一系列近似子問題解決的,這可能會導致一些誤差并導致性能損失。在文獻[11]中,通過利用連續凸逼近來增強物理層的安全性,研究了IRS輔助通信中用于安全通信的資源分配設計。在文獻[12]中,作者提出基于梯度下降算法的復圓流形(complex circle manifold, CCM)法來解決單位模量約束。雖然可以獲得相移的全局解,但是CCM法是黎曼一階算法,越接近目標值,步長會越小時,進而收斂速度會變慢,求解會需要多次迭代。綜上所述,在IRS輔助通信系統中,目前沒有一種方法可以在具有較低復雜度和快速收斂速度時保持優異的性能。因此,提供有效的IRS反射優化算法是對于提高IRS輔助通信系統的性能是很有必要的。

因此,在基站(base station, BS)功率限制和IRS單位模量約束下,本文通過共同優化BS預編碼矩陣和IRS相移向量來實現所有用戶信道容量的最大化。但是由于優化變量高度耦合,以及單位模量約束是高度非凸的,容量最大化問題較為棘手。本文的貢獻總結如下:

(1) 針對非凸非確定性多項式(non-deterministic polynomial, NP)難優化問題,通過利用基于信道容量與加權最小均方誤差(weighted minimum mean square error, WMMSE)之間的等價關系,將原問題轉換成一個易解決的問題形式。再利用交替優化算法對預編碼矩陣和相移向量進行交替優化。

(2) 當固定相移向量時,預編碼矩陣優化問題變成一個凸優化問題,由約束條件將其看作一個二階錐規劃問題(second order cone problem, SOCP),然后使用標準優化包[13]求解最優預編碼矩陣。當固定預編碼矩陣時,經過復雜的矩陣變換將相移優化問題轉換為受單位模量約束的非凸二次規劃(quadratic programming, QP)問題。然后,根據黎曼流形優化技術[14],單位模量約束條件可以構成一個復圓流形并嵌入搜索空間后,使問題變成一個無約束問題。最后,然后,利用黎曼信賴域(Riemannian trust-region, RTR)算法獲得局部最優解。

(3) 數值結果表明,相比于基準算法,RTR算法在保持較低復雜度和快速收斂速度的同時提高了性能。

1 系統模型

如圖1所示,本文研究了一個IRS輔助多用戶MIMO通信系統。

圖1 IRS輔助多用戶MIMO通信系統模型

該系統中,由N個反射元件組成的IRS協助帶有M1根天線的BS和帶有M2根天線的K個用戶之間的下行鏈路通信。由于路徑損耗的影響,本文僅考慮由IRS首次反射的路徑,不考慮二次及之后的反射路徑。此外,假設所有信道都是準靜態平坦衰落信道,并且其信道狀態信息(channel state information, CSI)在BS處是已知的。因此,用戶k處的基帶接收信號yk可以寫成:

(1)

此外,本文考慮一個實用的離散相移的IRS模型,為了利于實現,設置振幅系數α=1,θn僅可取值為Lps=2b個離散值,其中b為位數[15]。那么,θn在間隔[0,2π)將被均勻量化,滿足集合θn∈F={0,Δθ,…,(Lps-1)Δθ},其中Δθ=2π/Lps。那么,對于用戶k而言,信道容量表示為

(2)

(3)

式中:Pmax>0是BS的最大發射功率閾值。然后,針對IRS處的反射優化約束,已知IRS反射矩陣φ中的φn=ejθn,注意到離散相移θn通常很難解決。因此,取出φ的對角元素定義為向量ζ=[ζ1,ζ2,…,ζN]H,其中ζn=ejθn,?n∈N。那么,離散相移約束θn∈F可以等價于單位模量約束[12]:

|ζn|=1,?n∈N

(4)

因此,信道容量最大化問題可以寫成:

(5)

式中:ωk表示用戶k權重因子,限制條件C1表示所有用戶的發射功率之和不能超過其最大發射功率Pmax,C2表示IRS處的單位模量約束。

2 問題變形

從式(5)中可以觀察到,此時的預編碼矩陣和相移向量是耦合的,并且約束條件C2是高度非凸的,故該問題是一個非凸NP難問題。因此,本文首先利用最優解碼矩陣的均方誤差(mean square error, MSE)和信道容量表達式之間的關系[16]對式(5)進行變形。為利于實現,本文考慮使用線性解碼矩陣,那么估計信號向量表示為

(6)

(7)

rk(W,D,G,ζ)=log2|Wk|-tr(WkEk)+q

(8)

因此,通過利用信道容量表達式和MSE矩陣之間的關系[16-17],可以建立式(5)與上述函數之間的聯系。那么,式(5)重新表示為

s.t. C1, C2

(9)

與式(5)相比,盡管式(9)引入了兩個優化變量,但更易于解決。因為對于式(9),當給定G或ζ時,如果輔助矩陣W和解碼矩陣D也固定時,則rk(W,D,G,ζ)對于每組優化矩陣都是一個凹函數,可以利用交替優化算法對預編碼矩陣G或相移向量ζ進行求解。

3 預編碼矩陣優化

利用交替優化算法,首先考慮解碼矩陣D與輔助矩陣W,因為這兩個優化變量僅與函數rk(W,D,G,ζ)有關。當求解碼矩陣Dk的最優解時,通過固定W,G,ζ,并將函數rk(W,D,G,ζ)關于Dk求一階導數后令其等于零。最優解碼矩陣Dk為

(10)

同理,最優輔助矩陣Wk為

(11)

然后,在對發射預編碼矩陣G的優化時。對于任何給定的W,D,ζ,將Ek代入式(9)后,通過忽略常數項,可以將發射預編碼矩陣優化問題表示為

(12)

s.t. C1

由式(12)易知,該問題是一個凸優化問題,并且約束條件C1是一個二階錐約束條件。因此,該問題可看作是一個SOCP問題,可以利用現有的標準優化工具包例如CVX[13],直接求得發射預編碼矩陣的最優解。

4 相移向量優化

同理,在固定矩陣W,D,G時,將Ek代入式(9)并忽略常數項,關于相移向量的優化問題可以表示為

s.t. C2

(13)

s.t. C2

(14)

式中:

ζ是由φ的對角元素的收集向量,那么根據文獻[18]中的矩陣等式,有以下等式:

tr(φHBφC)=ζH(B⊙C)ζ

(15)

同理,令v=[[V]1,1,[V]2,2,…,[V]N,N]H為矩陣V對角元素的收集向量,有tr(φHVH)=vTζ,tr(Vφ)=ζHv*。那么,式(15)中的問題可以變成以下等效問題:

s. t. C2

(16)

式中:矩陣U=B⊙C。顯然,約束條件C2中的單位模量約束是解決式(16)中問題的主要障礙。針對該非凸約束條件,本文決定利用黎曼流形優化工具將該約束條件嵌入搜索空間[19],求解基于約束空間的無約束優化問題。

首先,該問題的單位模量約束可看作是N個復數圓的乘積。那么,定義一個如圖2所示的復圓乘積流形M:

圖2 黎曼流形幾何釋義

MSN={ζ∈CN:|ζn|=1,?n∈N}

(17)

式中:S{ζ∈C:ζ*ζ=Re{ζ}2+lm{ζ}2=1}表示一個復數圓。在算法闡述之前,首先提供關于流形優化的背景知識[14]。

根據圖2,設當前迭代點為ζ(i)。可以看出對于乘積流形M,切空間Tζ(i)M對算法搜索方向ηi∈Tζ(i)M的確定意義重大。對于任何點ζ(i)∈M,切空間由所有經過該點的切向量組成。那么,在點ζ(i)∈M上的切空間Tζ(i)M如下所示:

(18)

式中:z是在點ζ(i)上的切向量;0N是N維零向量。

在歐氏空間中,歐氏梯度表示某一函數在該點處變化最快的方向。那么在點ζ(i)∈M,式(16)目標函數f(ζ(i))的歐氏梯度Gradζ(i)f具體如下所示:

Gradζ(i)f=2(Uζ(i)+v)

(19)

與歐式梯度相似,黎曼梯度是目標函數增長最快的切線方向。并且目標函數的黎曼梯度gradf(ζ(i))被定義為滿足以下條件的切空間中的獨特元素[14]:

Df(ζ(i))[ζ(i)]=〈gradf(ζ(i)),ζ(i)〉

(20)

式中:Df(ζ(i))為在點ζ(i)時f(ζ(i))的方向導數。如圖2所示,在切空間Tζ(i)M上,黎曼梯度是歐氏梯度的正交投影,如下所示:

(21)

如圖2所示,在流形上,沿著點ζ(i)∈M的切線向量方向移動來找到在流形上的目標點定義為收縮Rζ(i)(ηi):Tζ(i)M→M,是從切線空間Tζ(i)M到流形M的映射[14]。黎曼指數收縮是在黎曼流形優化時最常用的收縮,但其計算成本可能很高。因此,本文決定使用二階收縮替代指數收縮,因為其不僅具有較低的計算成本,并且是指數映射的近似值[20]。因此,點ζ(i)∈M處的收縮Rζ(i)(ηi)定義為

(22)

黎曼黑塞是梯度向量場相對于Levi-Civita連接的協變導數,當在點ζ(i)處時,是從切空間Tζ(i)M到自身的線性映射Hessf(ζ(i))[14],如下所示:

(23)

(24)

式中:DGradf(ζ(i))[ηi]是Gradf(ζi)在ζ(i)上沿搜索方向ηi∈Tζ(i)M的方向導數,根據式(20)計算如下:

DGradf(ζ(i))[ηi]=〈2U,ηi〉

(25)

在文獻[12]提出的CCM方法,是一種一階黎曼流形優化方法,主要利用的是黎曼最速下降方法。該方法得到的是全局解,但不能保證全局最優解,并且越接近目標值,步長越小時,收斂速度會變慢,求解會需要多次迭代。因此,本文提出一種二階黎曼流形優化方法——RTR算法。該算法利用二階幾何形狀能夠避開鞍點并獲得更準確的結果,同時具有低復雜度和快速收斂速度。

首先,需要構造一個信賴域子問題。在第i次迭代時利用二次逼近,在流形M上將目標函數f(ζ(i))的二階泰勒展開定義一個二次模型是足夠的[14]。那么,在點ζ(i)∈M的由二次模型構成的信賴域子問題,如下所示:

(26)

式中:ηi∈Tζ(i)M是當前迭代的搜索方向;Δi是當前迭代的信賴域半徑。對于式(26)中的信賴域子問題,截斷共軛梯度(truncated conjugate-gradient, TCG)算法特別適用于解決該問題[14]。因為盡管隨著IRS反射元件數量N增加時,問題規模也逐漸變大,TCG算法依舊可以快速地處理每個信賴域子問題。

然后,是否調整信賴域半徑Δi和更新下一迭代變量ζ(i+1),由以下評價函數確定:

(27)

根據式(27)可以看出,該評價函數[14]定義是在第i次迭代時目標函數的實際下降量f(ζ(i))-f(Rζ(i)(ηi))與二次模型函數的預測下降量mζ(i)(0)-mζ(i)(ηi)的比值。根據計算的ρi,如果ρi<1/4,那么需要減小信賴域半徑Δi和重新計算搜索方向ηi;如果ρi→1時,這表明二次模型和目標函數在信任區域中具有良好的近似性,可以在下一次迭代時擴展信賴域半徑。RTR算法的關鍵步驟如算法1所示。

算法 1 RTR算法步驟 1 初始化:最大信賴域半徑Δ>0,初始信賴域半徑Δ0∈(0,Δ), 步長接受閾值ρ'∈0,14 ,初始迭代點ζ(0)∈M和算法迭代停止門限值ε。步驟 2 for i=0,1,… do

步驟 3 求解式(26)中的信賴域子問題獲得搜索方向ηi步驟 4 根據式(27)計算比值ρi步驟 5 調整信賴域半徑Δi+1=min(2Δi,Δ), ρζi>34; ηi=Δi14Δi, ρζi<14Δi, 其他步驟 6 計算下一迭代變量ζ(i+1)if ρi>ρ', then ζ(i+1)=Rζ(i)(ηi)elseζ(i+1)=ζ(i)步驟 7 i←i+1步驟 8 untilgradf(ζ(i))2≤ε步驟 9 end for步驟 10 輸出:當前迭代點ζ(i)

參考文獻[14]中有關RTR算法的收斂性分析,表明RTR算法產生的序列{ζ(i)}將收斂至目標函數f(ζ(i))的固定點集合,即利用該算法將最終收斂于目標函數的局部最優值。

5 交替優化算法與復雜度分析

5.1 交替優化算法

基于上述對預編碼矩陣和相移向量的優化分析,解決式(9)的關鍵步驟總結在算法2中。

算法 2 交替優化算法步驟 1 初始化:當前迭代次數i=0,最大迭代次數imax,初始預編碼矩陣G(0)和初始相移向量ζ(0),誤差容忍值ε,計算式(5)中的目標函數值R(G(0),ζ(0))步驟 2 給定G(i)和ζ(i)時,計算式(10)中的最優解碼矩陣D(i)步驟 3 給定G(i),ζ(i),D(i)時,計算式(11)中的最優輔助矩陣W(i)步驟 4 給定W(i),ζ(i)和D(i)時,利用CVX工具求解獲得最優預編碼矩陣G(i+1)步驟 5 給定G(i),W(i)和D(i)時,利用RTR算法求解相移向量ζ(i+1)步驟 6 if|R(G(i+1),ζ(i+1))-R(G(i),ζ(i))|R(G(i),ζ(i))<ε or i>imaxBreak;else i=i+1;返回步驟2;end if

5.2 復雜度分析

(28)

此外,與本文提出的RTR算法相比較的CCM算法,其計算復雜度為O(TCCMN2+N3)。其中TCCM代表CCM算法的總迭代次數,下一節的數值結果將比較其收斂速度。

6 仿真與數值結果分析

在本節中,提供數值結果來證明本文的理論分析以及所提出算法的出色性能。參照文獻[7],本文考慮一個IRS輔助的下行鏈路通信系統,該系統由具有M1=4根天線的BS和M2=2根天線K=4個用戶組成,數據流q=2,IRS反射元件數量N=20。假定所有涉及的信道都服從獨立瑞利分布,呈現平坦衰落,其中以dB為單位的路徑損耗模型[12]表示為

(29)

式中:PL0為參考距離d0=1 m時的路徑損耗,設定為-30 dB;α為路徑損耗指數,均設定為3。除稍后指定的某些特定參數外,其他參數如表1所示。

表1 仿真參數

根據表1,已知BS和IRS的固定位置,那么假設用戶隨機地分布在以(xu,0)為中心,半徑為10 m的圓內,如圖3所示。為了便于實施,假設xu=90 m,即用戶位于IRS附近。

圖3 IRS輔助MIMO多用戶情況 (俯視圖)

然后,本文將以下方案作為基準方案與RTR算法進行比較。

(1) CCM算法:交替優化算法的步驟5用文獻[12]中的CCM算法代替求解。

(2) 隨機相移:假設每個反射元素的相移是從集合θn∈F ={0,Δθ,…,(Lps-1)Δθ}獨立隨機生成的。那么在交替優化算法中,可以直接跳過步驟5。

(3) Without-IRS:通過設置不使用IRS的系統。

在圖4中,比較了RTR算法與CCM算法進行的收斂性能。從圖4中可以看出,RTR算法的收斂速度比CCM算法的收斂速度要快得多,并且RTR算法收斂后達到信道容量值比CCM算法的值要大。并且,在表2中比較了RTR算法與CCM算法在圖4中達到收斂時對應迭代次數的運行時間。結合第5節中分析了RTR算法與CCM算法的復雜度,已知兩者復雜度都較低。可以得出結論,RTR算法具有更快的收斂速度、更短的運行時間和更高的信道容量值,這說明了RTR算法比CCM算法具有更好的優化性能。

圖4 RTR算法與CCM算法的收斂性比較

表2 RTR算法與CCM算法的仿真時間

在圖5中,可以看到所有方案的信道容量值與IRS上反射元件數量N的關系。可以觀察到,隨著反射元件數量的增加,相比較于其他3種基準方案,本文提出的算法的性能增益越來越明顯。例如,相比較于Without-IRS情況,當N=20時RTR算法的性能增益僅為14.2%,但當N=70時,性能增益提升至43.8%。這主要是因為隨著N變大,可以增強在IRS處接收信號功率,從而獲得更高的陣列增益;并且通過適當地設計相移,用戶接收到的反射信號功率將隨著N的增加而增加。因此,所提出的IRS輔助系統不僅可以利用陣列增益,還可以利用IRS處的反射波束成形增益。此外,隨著N值的增加,可以注意到所有With-IRS的方案都比Without-IRS方案具有更高的信道容量值,并且由于IRS是無源的,不需要額外的射頻鏈和額外的功耗。因此,可以得出結論,在無線通信系統中部署IRS的有效和可行的。

圖5 信道容量與IRS反射原件數量關系

在圖6中,比較了所有方案實現的信道容量值與基站到用戶分布圓心點(xu,0)之間的水平距離xu的關系。可以觀察到,與其他3個基準方案相比,隨著xu的增加,即當用戶遠離BS,越來越靠近IRS時,RTR算法的性能增益也越明顯。這是因為用戶會從IRS接收到強烈的反射信號,進而可以減輕小區內用戶間干擾。

圖6 信道容量與用戶位置的關系

假定圖3中IRS坐標為(xIRS,0)。圖7中,比較了IRS的位置從xIRS=20 m到xIRS=100 m對信道容量的影響。可以觀察到,當2050時隨著xIRS增大而增加。這是因為當IRS位于BS與用戶的中間點時,組合的信道增益也將達到其最小值。并且可以再次觀察到,RTR算法與其他基準方案相比,大大提高了信道容量性能。例如,相比較于Without-IRS情況,當xIRS=20 m時RTR算法的性能增益僅為2.1 bps/Hz;但當xIRS=100 m時RTR算法的性能增益為4.8 bps/Hz,相比于xIRS=20 m時的性能增益提升了128.6%。這歸因于IRS與用戶之間的良好鏈路。同時應仔細設計IRS相移,因為當xIRS=50 m時,注意到RandPhase的實際性能比Without-IRS時的性能還差。

圖7 信道容量與IRS位置的關系

在圖8中繪制了信道容量性能與最大發射功率的關系。從圖8可以看出,當最大發射功率太小而無法實際實現時,即如果從BS到用戶的發射功率不足以滿足用戶的速率需求,那么所有方案的性能都對應于非常低的值。然后,隨著最大發射功率值的增加,曲線的斜率顯然會增加,并且RTR算法的性能依舊是最優異的。因此,可以結合圖5和圖8的分析得出結論,IRS輔助通信系統應該在IRS的反射元件數量和BS的發射功率之間做出最佳折衷。

圖8 信道容量與最大發射功率的關系

如式(5)中所述,權重ωk用于控制用戶之間的公平性。假設4個用戶的位置分別在圖3中坐標為(30,0),(50,0),(70,0)和(90,0)的位置。其中,第1個用戶比第2個用戶離BS更近,第4個用戶比第3個用戶離IRS更近。測試了兩組權重值:①ωk=0.25,?k;②ω1=0.15,ω2=0.35,ω3=0.3,ω4=0.2。從圖9中可以發現,在權重相同時,用戶1的數據速率值最高,這是因為用戶1與基站離得最近;用戶4的數據速率值比用戶1的略低一點,但比其他兩個用戶都要高,這是因為盡管用戶4離BS較遠但是離IRS最近,進而使該用戶信道增益增加;在權重不等時,為了保證用戶之間的公平性,通過給信道增益低的用戶分配更高的權重值,使得用戶的數據速率分布更平衡。但是,還可從圖9中發現,信道容量增益存在損耗,權重相等時,該系統的信道容量約為28.2 bps/Hz;在權重不等時,該系統信道容量約為23.1 bps/Hz。這說明存在最大化信道容量和用戶之間的公平性的矛盾。因此,若要使得信道容量與用戶之間的公平性達到平衡,在分配權重時盡量在滿足各個用戶的最低質量要求時,保證達到信道容量最大化。

圖9 權重影響

7 結 論

本文針對IRS輔助多用戶通信系統中的信道容量最大化問題,首先由于優化變量是高度耦合,所產生的優化問題是一個非凸NP難問題。因此,通過利用基于信道容量與最優解碼矩陣的均方誤差之間的等價關系,將原問題轉換成一個易解決的形式。再利用交替優化算法對預編碼器和相移器進行交替優化。為求解在單位模量約束下的相移向量,本文提出使用RTR算法,將單位模量約束條件構成一個CCM并嵌入搜索空間后,使問題變成一個無約束問題再獲得局部最優解。數值結果表明,相比于基準算法,RTR算法在保持較低復雜度和快速收斂速度的同時提高了性能。

猜你喜歡
優化用戶
超限高層建筑結構設計與優化思考
房地產導刊(2022年5期)2022-06-01 06:20:14
民用建筑防煙排煙設計優化探討
關于優化消防安全告知承諾的一些思考
一道優化題的幾何解法
由“形”啟“數”優化運算——以2021年解析幾何高考題為例
關注用戶
商用汽車(2016年11期)2016-12-19 01:20:16
關注用戶
商用汽車(2016年6期)2016-06-29 09:18:54
關注用戶
商用汽車(2016年4期)2016-05-09 01:23:12
基于低碳物流的公路運輸優化
現代企業(2015年2期)2015-02-28 18:45:09
Camera360:拍出5億用戶
創業家(2015年10期)2015-02-27 07:55:08
主站蜘蛛池模板: 中美日韩在线网免费毛片视频| 麻豆精品久久久久久久99蜜桃| 九九九精品视频| 国产成人精品第一区二区| 国产精品成人免费综合| 国产成人精品在线1区| 国语少妇高潮| 亚洲欧洲日本在线| 亚洲国产中文精品va在线播放| 久久久精品无码一区二区三区| 欧美成人综合视频| 看看一级毛片| 色综合a怡红院怡红院首页| 婷婷色一区二区三区| 久久6免费视频| 亚洲人成网站在线播放2019| 亚洲精品爱草草视频在线| 久久人人97超碰人人澡爱香蕉 | 免费全部高H视频无码无遮掩| 高清无码手机在线观看| 亚洲综合中文字幕国产精品欧美| 久久精品欧美一区二区| 国产麻豆精品久久一二三| 欧美性猛交xxxx乱大交极品| 精品国产成人av免费| 色妺妺在线视频喷水| 久久精品66| 欧美黄色网站在线看| 婷婷色狠狠干| 2022国产91精品久久久久久| 91国内外精品自在线播放| 国产自产视频一区二区三区| 国产97视频在线观看| 日韩在线网址| 中国丰满人妻无码束缚啪啪| 国产欧美在线观看精品一区污| 国产在线91在线电影| 亚洲精品无码AV电影在线播放| 国产一级做美女做受视频| 日本欧美视频在线观看| 米奇精品一区二区三区| 女人18毛片水真多国产| 蜜桃臀无码内射一区二区三区| 怡春院欧美一区二区三区免费| 操美女免费网站| 91视频国产高清| 国产麻豆福利av在线播放| 国产一区二区三区在线无码| 亚洲精品人成网线在线 | 欧美一级在线播放| 亚洲国产综合精品一区| 在线观看亚洲精品福利片| 人妻免费无码不卡视频| 国产欧美网站| 日韩国产 在线| 国内毛片视频| 亚洲美女视频一区| 欧美日韩中文国产| 毛片免费观看视频| 中文字幕无线码一区| av在线无码浏览| 在线观看网站国产| 国产精品成人啪精品视频| 亚洲一区色| 在线毛片免费| av午夜福利一片免费看| 国产亚洲视频中文字幕视频| 亚洲欧美一区二区三区蜜芽| 日本一区二区不卡视频| 国产精品网址在线观看你懂的| 免费在线观看av| 亚洲国产中文欧美在线人成大黄瓜 | 激情国产精品一区| 亚洲a级毛片| 亚洲一区二区成人| 亚洲精品成人福利在线电影| 国产午夜看片| 国产精品lululu在线观看| 美女一区二区在线观看| 伊人久久精品无码麻豆精品| 好久久免费视频高清| 国产激情无码一区二区免费|