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

        ?

        GNSS多徑信道模擬的聚類稀疏擬合方案*

        2023-09-28 07:21:24唐小妹
        國防科技大學學報 2023年5期
        關(guān)鍵詞:方法模型

        周 順,歐 鋼,唐小妹

        (1. 國防科技大學 電子科學學院, 湖南 長沙 410073; 2. 國防科技大學 第六十三研究所, 江蘇 南京 210007)

        在非理想環(huán)境下,多徑信號仍然是全球?qū)Ш叫l(wèi)星系統(tǒng)(global navigation satellite system, GNSS)應(yīng)用中最主要的精度降級源之一[1]。因此,如何在多徑環(huán)境中提高定位性能成為一項具有挑戰(zhàn)性的課題。研究多徑環(huán)境對GNSS性能的影響,可以利用信道模擬器模擬各種多徑信道和信號組合。同時,無線系統(tǒng)信號級仿真需要模擬電磁波及傳播環(huán)境之間的相互作用,計算負擔的平衡是無線仿真系統(tǒng)最大的挑戰(zhàn)。

        文獻[2]介紹了全球最大的基于射頻信道模擬器的頻譜協(xié)作挑戰(zhàn)賽測試平臺,該平臺采用稀疏有限長沖擊響應(yīng)(finite impulse response, FIR)濾波器近似地模擬多徑信道,目的在于降低大規(guī)模信道模擬的硬件資源消耗。但是,其公開的文獻內(nèi)容并未涉及設(shè)計稀疏濾波器抽頭的方法。同樣,在GNSS星座信道模擬器研究中,衛(wèi)星導航信道沖擊響應(yīng)(channel impulse response, CIR)多徑參數(shù)數(shù)量眾多,難以實時仿真計算。因此,需要研究減少CIR參數(shù)的優(yōu)化方法以降低計算負擔,同時盡可能地保持信道模型對導航接收機的影響[3-4]。

        本文針對衛(wèi)星導航系統(tǒng)測試與評估中大規(guī)模信道模擬實時計算復雜度高的問題,提出一種GNSS多徑信道模擬的聚類稀疏擬合方案。從整體結(jié)構(gòu)角度,提出GNSS星座信道的稀疏抽頭延遲線(tapped-delay-line, TDL)結(jié)構(gòu)模擬框架;在方法上,提出了基于K中心聚類CIR參數(shù)萃取的信道稀疏擬合方法,提取等效精簡CIR參數(shù),再利用稀疏TDL結(jié)構(gòu)實現(xiàn)等效精簡CIR參數(shù)模擬,使信道模擬的復雜度不依賴于傳播路徑數(shù)目。該方案采用TDL結(jié)構(gòu)FIR濾波器稀疏擬合原始GNSS多徑信道模型,簡化了GNSS信道模擬的復雜度,同時保持了原始模型的大部分精度。這對于GNSS系統(tǒng)設(shè)計優(yōu)化,提供能夠?qū)崟r計算的GNSS算法測試評估環(huán)境具有重要意義。

        1 問題描述與方案框架

        GNSS信道模擬系統(tǒng)的核心部分是采用 TDL結(jié)構(gòu)的 FIR濾波器[5-6],每個濾波器的抽頭代表一個多徑信號拷貝,通常使用系統(tǒng)采樣速率的倒數(shù)Ts作為抽頭延遲間隔,以使信道模型仿真與接收機信號處理能夠以相同的更新速率連接。TDL結(jié)構(gòu)是實現(xiàn)多徑信道模型的關(guān)鍵硬件結(jié)構(gòu),它占用了大量實時硬件資源,因此需要優(yōu)化它的實現(xiàn)方式。

        圖1為密集的N抽頭FIR濾波器結(jié)構(gòu)圖,其具有密集復系數(shù)抽頭,抽頭之間是一個單位延遲單元,每個抽頭對應(yīng)一個復數(shù)乘累加。當需要模擬的多徑數(shù)目較多時,所耗硬件資源非常龐大,特別是對于GNSS星座的信道模擬,將很快耗盡有限的硬件乘法器和加法器資源。文獻[2]提出用4個抽頭FIR濾波器代替具有密集單位延遲抽頭的FIR濾波器結(jié)構(gòu)。由于在FIR濾波器中對乘法器、加法器進行裁剪,大大減少了硬件資源的消耗。然而,該文獻并未公開剪裁方法及效果。

        圖1 密集的N抽頭FIR濾波器結(jié)構(gòu)Fig.1 Structure of dense N taps FIR filter

        如圖2所示,本文提出GNSS星座信道的稀疏TDL結(jié)構(gòu)模擬框架,其中使用稀疏FIR濾波器作為GNSS星座信道模擬的核心模擬單元,將一組稀疏FIR濾波器和加法器組合在一起。進一步地,除了采用以上稀疏FIR濾波器結(jié)構(gòu)以減少硬件乘法器的使用,還提出稀疏擬合CIR參數(shù)的方法,方法在第2節(jié)詳細描述。

        圖2 GNSS星座信道的稀疏TDL結(jié)構(gòu)模擬框架Fig.2 Framework of sparse TDL structure emulation for GNSS satellite constellation channel

        2 稀疏擬合方法

        本文進一步提出了基于K中心聚類CIR參數(shù)萃取的稀疏擬合方法,得到了用于稀疏擬合CIR參數(shù)的等效精簡參數(shù)。稀疏擬合方法框圖如圖3所示,包括四個主要部分:信道重采樣、K中心多徑聚類分析、最鄰近多徑分量聚合和稀疏FIR濾波。其中K中心多徑聚類分析與最鄰近多徑分量聚合又合稱為K中心聚類CIR參數(shù)萃取。信道模擬的過程中,采用GNSS信號發(fā)生器產(chǎn)生標準測試激勵信號,GNSS接收機作為被測試對象。圖中所示的GNSS參考信道模型是用來根據(jù)仿真條件計算產(chǎn)生CIR參數(shù)的。

        圖3 稀疏擬合方法框圖Fig.3 Diagram of the sparse fitting method

        該方法采用信道重采樣、K中心多徑聚類分析、最鄰近多徑分量聚合等步驟萃取等效精簡CIR參數(shù),作為稀疏FIR濾波器的抽頭系數(shù)。它的技術(shù)效果是擬合原始GNSS信道模型中的多徑效應(yīng),保持GNSS接收機引入的多徑誤差不變,從而通過稀疏FIR濾波器模擬導航信號的多徑效應(yīng)。

        2.1 信道重采樣

        通常,由GNSS參考信道模型生成CIR參數(shù)作為TDL結(jié)構(gòu)FIR濾波器的抽頭系數(shù),存在抽頭之間的差分延遲不是仿真采樣時間倍數(shù)的問題[7]。為了克服這個問題,通過信道重采樣將連續(xù)時間時延域上CIR參數(shù)轉(zhuǎn)換為均勻采樣的離散時間形式[4]。離散多徑信道沖擊響應(yīng)函數(shù)hd(t,τ)可表示為

        (1)

        其中,N為離散多徑的數(shù)量,ai(t)、φi(t)和τi(t)分別為第i路多徑的幅度、相位和延遲。FIR濾波器的抽頭系數(shù)取值為復加權(quán)系數(shù)

        Ai(t)=ai(t)exp(-jφi(t))

        (2)

        延遲τi(t)與濾波器抽頭延遲間隔無法對應(yīng)。因此,需要將hd(t,τ)與sinc函數(shù)卷積,即

        hc(t,τ)=hd(t,τ)*sinc(πBτ)

        (3)

        通過這個信道重采樣步驟,hd(t,τ)變換為以B為采樣率的CIR函數(shù)hc(t,τ)。提取hc(t,τ)的CIR參數(shù),包括復加權(quán)系數(shù)和對應(yīng)的時延作為重采樣后的CIR參數(shù)輸出到后續(xù)處理步驟。

        2.2 基于K中心聚類的CIR參數(shù)萃取

        無監(jiān)督學習通過對無標記的數(shù)據(jù)樣本進行學習,揭示數(shù)據(jù)樣本的內(nèi)在特征,聚類分析被認為是一種最重要的無監(jiān)督學習任務(wù)[8]。一種常用的聚類算法是K均值(K-means)算法,在文獻[9]中,作者已經(jīng)提出了采用K均值聚類算法簡化CIR參數(shù)。作為一種改進,本文選擇K中心聚類算法,它是一種針對K均值聚類算法的改進,作為聚類分析的候選算法。主要是基于以下考慮:

        1)K均值算法對異常值非常敏感,具有極大值的對象可能會嚴重扭曲數(shù)據(jù)分布。

        2)K中心聚類算法采用數(shù)據(jù)集的實際樣本作為簇中心,而不是像K均值算法一般使用簇中對象的均值[10],從而降低了算法對噪聲和孤立點的敏感性。

        3)在正態(tài)分布和均勻分布的情況下,K中心聚類算法所花費的平均時間小于K均值算法所耗費的時間[11]。

        本文在對CIR參數(shù)實施K中心聚類分析的基礎(chǔ)上,進行最鄰近多徑分量聚合,使衛(wèi)星導航多徑信道等效精簡。下面對K中心聚類分析算法的具體細節(jié)進行介紹。

        首先將重采樣后的CIR參數(shù)視為一個數(shù)據(jù)集X=[x1,x2,…,xN],其中有N個CIR參數(shù)樣本,每個樣本由兩個屬性的特征向量組成。兩個屬性分別為復加權(quán)系數(shù)g和多徑分量時延τ,可表示為2×N的矩陣形式

        (4)

        其中:g1,g2,…,gN為抽頭增益;τ1,τ2,…,τN為時延。K中心聚類是通過選擇數(shù)據(jù)集X的M個對象作為中心(medoids),并將X剩余的未被選擇成員分配給其最近的中心來生成的。具體來說,若采用K中心聚類方法將數(shù)據(jù)集X劃分為一組簇C={C1,C2,…,CM},C滿足以下三條屬性[12-13]:

        1)Ci≠?,i=1,…,M;

        3)Ci∩Cj=?,i≠j,1≤i≤M,1≤j≤M。

        此外,簇應(yīng)當反映數(shù)據(jù)的結(jié)構(gòu),即同一簇中的對象彼此相似,不同簇中的對象彼此不同。為了解決聚類問題,需要一種定量的方法來區(qū)分相似和不同的對象。一個常用的定量準則是簇內(nèi)變分和S來表示

        (5)

        其中,mj?X代表第j簇的中心點,d(xi,mj)是第i個對象xi和第j個簇中心mj之間的相異性度量。本研究中,采用歐幾里得距離作為相異性度量,對象i和對象j之間的歐幾里得距離由式(6)給出。

        (6)

        通過K中心聚類算法對數(shù)據(jù)集分簇,劃分出子集,子集的中心稱為中心點,聚類分析的作用是提取各個簇中心對應(yīng)的多徑時延。該算法使用簇中的真實樣本,而不是如K均值算法那樣使用均值作為簇中心,從而降低了算法對噪聲和異常值的敏感性。

        值得注意的是,雖然文獻[4]已經(jīng)提出了采用K均值聚類算法進行多徑聚合,并且認為K均值聚類算法具有最優(yōu)的性能,但忽略了K均值聚類得到的聚類中心不在采樣時刻,因而無法直接應(yīng)用于TDL結(jié)構(gòu)信道模擬。因此,為了能夠在信道模擬中使用K均值聚類,本文首先處理了聚類中心集合,然后將聚類中心移動到最近的實際采樣時刻,從而與TDL結(jié)構(gòu)的信道模擬機制相匹配。采用本文提出的K中心聚類算法是以實際樣本作為簇中心,因此無須附加額外的搬移步驟。

        需要說明的是,K中心聚類可能在一些場景下與真實信道的聚類結(jié)構(gòu)不匹配,而且聚類算法中的隨機處理等原因也會給算法的穩(wěn)定收斂帶來影響,導致不能很好反映多徑信道的原始形態(tài)。為了減少這些問題的發(fā)生,可以通過集成多個聚類學習器等方法[9]來嘗試解決。

        2.3 稀疏FIR濾波

        信道模擬器本質(zhì)為復數(shù)FIR濾波器,采用抽頭延遲線結(jié)構(gòu),由復系數(shù)表示濾波器的單個抽頭。本文采用稀疏抽頭數(shù)字FIR濾波器模擬GNSS多徑信道,對送入該濾波器的輸入信號進行稀疏抽頭數(shù)字FIR濾波,由于未采用完全填充的FIR濾波器,大多數(shù)抽頭上的乘法器和加法器都不需要配置,因此硬件資源成本大大降低。

        圖4為稀疏抽頭數(shù)字FIR濾波器的結(jié)構(gòu)示意圖。稀疏抽頭數(shù)字FIR濾波器包括一組橫向排列的移位寄存器延遲單元、一組加權(quán)系數(shù)單元和一組對應(yīng)的加法器。輸入信號x(i)在每個延遲單元中延時,延遲單元具有相同單位延遲時間。與均勻間隔數(shù)字FIR濾波相區(qū)別的是,各個加權(quán)系數(shù)之間的延時不再是均勻間隔的單位延遲,各延時的信號在加權(quán)系數(shù)單元中用所對應(yīng)的加權(quán)系數(shù)進行加權(quán)。也可以說,加權(quán)系數(shù)是多徑信道的稀疏擬合值[14],本文中就叫作稀疏抽頭FIR濾波器的抽頭系數(shù)。通常加權(quán)系數(shù)是短時平穩(wěn)的,與信號的變化速度相比,加權(quán)系數(shù)的變化是比較緩慢的。延時并加權(quán)之后的信號在一組加法器中逐級相加,這樣的好處是可把加法計算時間分攤到分散的時延內(nèi),最后得到輸出信號y(i)。

        圖4 稀疏抽頭FIR濾波器示意圖Fig.4 Sparse taps FIR filter

        將圖4和圖1所示濾波器結(jié)構(gòu)特征進行比較,可以看出這里描述的信道模擬濾波器是一種稀疏抽頭FIR濾波器。與參考信道模擬濾波器相比,其模擬GNSS多徑信道實時計算復雜度大大簡化。

        3 性能評估

        3.1 評估準則

        由于多徑導致信號自相關(guān)函數(shù)失真,延遲鎖相環(huán)S曲線的過零點被偏移,并且該偏移與接收機偽距誤差直接相關(guān)。因此,本文所用的評估準則著重于鑒別器開環(huán)誤差評估方法,將針對每個CIR樣本獨立考慮鑒別器誤差。采用前后功率差(early minus late power,EMLP)鑒別器。EMLP鑒別器表達式為

        (7)

        式中,符號RS,H*S表示經(jīng)過信道加權(quán)的輸入信號與本地復現(xiàn)信號的復相關(guān)函數(shù),CS是相關(guān)器間隔,D(τ)構(gòu)成了S曲線。如果只有直射信號到達接收機,則方程D(τ)=0的解為τ=0。當多徑信號疊加到直射信號上,使接收機自相關(guān)函數(shù)失真時,同樣意味著S曲線失真,方程D(τ)=0的解會偏離0。鑒別器誤差Ed由式(8)給出

        Ed=c×τ′

        (8)

        式中,c為光速,τ′為方程D(τ)=0的解。

        在仿真結(jié)果中,利用鑒別器開環(huán)誤差來評估本文所提算法的性能。選擇這種方法是因為:鑒別器誤差量化了多徑信道對GNSS偽距誤差的影響;此外,開環(huán)特性允許獨立于任何環(huán)路效應(yīng)來研究多路徑的影響。由簡化后的CIR快照引起的鑒別器誤差相對于原始CIR快照的偏差(以下簡稱為“鑒別器誤差偏差”),可以量化使用本文方法模擬多徑信道時對于衛(wèi)星導航接收機的影響。

        3.2 仿真設(shè)置

        下面介紹評估多徑信道稀疏擬合效果的仿真設(shè)置。仿真評估框圖如圖5所示,信道模型產(chǎn)生一組CIR參數(shù)數(shù)據(jù)集,對應(yīng)于仿真時間長度,在仿真運行過程中讀取數(shù)據(jù)集,每個仿真快拍的CIR參數(shù)需要經(jīng)過信道重采樣,然后用K中心聚類CIR參數(shù)萃取進行精簡。GNSS信號發(fā)生單元產(chǎn)生基帶仿真信號和CIR參數(shù)進行卷積運算,然后輸出給多徑誤差評估單元。

        圖5 仿真評估框圖Fig.5 Diagram of simulation evaluation

        本文所用的信道模型是基于最新版本ITU-R P.681-11建議書[15]中的混合傳播條件下物理統(tǒng)計寬帶模型,該模型從2002年開始開展信道測量活動,到2009年被納入ITU建議書中,旨在逼真和精確建模GNSS應(yīng)用中的多徑傳播效應(yīng)。圖6展示了物理統(tǒng)計寬帶模型概念。文獻[16]中提供了有關(guān)該信道模型的詳細信息。

        圖6 物理統(tǒng)計寬帶模型概念框圖Fig.6 Structure of the physical-statistical wideband model

        設(shè)置信道模型為城市中的車輛運動場景,所定義的場景針對1 575.42 MHz頻點生成。設(shè)仿真采樣率為100 MHz,信道的最大延遲約為250 ns,信道參數(shù)每10 ms更新一次,人工虛擬環(huán)境中的車輛運動速度設(shè)為10 m/s。

        3.3 仿真結(jié)果

        3.3.1 單次CIR快照結(jié)果

        由仿真設(shè)置中所介紹的信道模型產(chǎn)生的單次CIR快照如圖7所示,圖中的紅色圓圈標示出了原始信道的CIR,圖中“?”標記的是信道重采樣后的CIR,在最大時延處得到24個TDL等效抽頭,對應(yīng)24個抽頭的復系數(shù)FIR濾波器。單次CIR快照涉及的仿真將在此基礎(chǔ)上展開優(yōu)化。

        圖7 原始信道模型CIR與信道重采樣CIRFig.7 CIR snapshot of original channel model and re-sampled channel

        從圖7中可以看到,重采樣得到的CIR中有24個抽頭,每個抽頭的幅度(相對于視距(line of sight, LOS)分量的比值)在圖中顯示,如若將CIR的抽頭數(shù)量從24個降低到9個,可以減少62.5%的抽頭數(shù)目。如果采用前述的稀疏FIR濾波器來模擬信道,可以節(jié)約大量的硬件乘法器和加法器資源。

        使用基于K均值聚類的方法與本文所提出的基于K中心聚類算法的方法,分別來簡化這24個抽頭的CIR。如圖8所示,分別將K均值聚類和K中心聚類兩種方法得到的CIR與信道重采

        (a) K均值聚類(a) K-means clustering

        (b) K中心聚類(b) K-medoids clustering圖8 兩種聚類方法得到的CIR與信道重采樣后CIR的比較Fig.8 Comparison of CIR given by two clustering methods and re-sampled channel

        樣后的CIR進行比較,發(fā)現(xiàn)兩種算法的精簡效果基本一致。另外,如圖9所示,采用L1 C/A碼,將兩種方法得到的信道自相關(guān)函數(shù)(autocorrelation function,ACF)與原始信道ACF進行比較,圖中橫坐標的時延已歸約到L1 C/A碼的碼片長度。由圖可以看出,兩種方法簡化得到的ACF同原始信道ACF的區(qū)別很小。

        (a) K均值聚類(a) K-means clustering

        (b) K中心聚類(b) K-medoids clustering圖9 兩種聚類方法得到的ACF與原始信道ACF比較Fig.9 Comparison of ACF given by two clustering methods and original channel

        3.3.2 小型CIR快照數(shù)據(jù)集結(jié)果

        由仿真設(shè)置中所介紹的信道模型產(chǎn)生1 500組CIR快照,其功率延遲剖面概率密度如圖10所示,將這1 500組數(shù)據(jù)稱為小型CIR快照數(shù)據(jù)集。以下仿真首先分析小型CIR快照數(shù)據(jù)集的不同期望抽頭下的均方根性能,然后借助箱線圖工具分析結(jié)果的統(tǒng)計性能。

        圖10 小型CIR快照數(shù)據(jù)集的功率延遲剖面概率密度圖Fig.10 Probability density for power delay profile of the small CIR snapshots dataset

        首先,比較了在設(shè)置4至20個期望抽頭數(shù)的目標下的鑒別器誤差偏差。取鑒別器誤差偏差的均方根值,得到如圖11所示的期望抽頭數(shù)與鑒別器誤差偏差均方根值的關(guān)系曲線。觀察到采用K中心聚類對應(yīng)的曲線,設(shè)置抽頭數(shù)減少到4個時,偏差均方根值約為0.77 m;期望目標變?yōu)?0個抽頭時,偏差均方根值約為0.13 m。由此可見,隨著抽頭數(shù)量的增加,稀疏擬合信道逼近原始信道的程度也在不斷提升。同時,比較圖11中兩條曲線,顯示本文提出的方法性能逼近傳統(tǒng)的K均值聚類方法,而且在抽頭數(shù)大于9個時,所提方法的均方根值逐漸小于K均值聚類方法。注意到,在抽頭數(shù)為8個時,鑒別器誤差偏差的均方根值小于0.5 m,此時已經(jīng)比大部分GNSS用碼偽距定位的精度小了一個數(shù)量級,因此建議將稀疏擬合的抽頭數(shù)目標定為8個。

        圖11 期望抽頭數(shù)與鑒別器誤差偏差均方根值的關(guān)系Fig.11 Relation of desired taps and root mean square values of deviation from discriminator errors

        然后,如圖12所示,給出了具有不同期望抽頭數(shù)的鑒別器誤差偏差箱線圖,比較在設(shè)置4~20個期望抽頭數(shù)的目標下兩種方法的統(tǒng)計性能。圖12中箱子的上、下四分位數(shù)和中值處有一條線段,箱子末端延伸出去的線稱為須,須線會延伸到不是離群值的最遠端數(shù)據(jù)點,離群值以紅色圓點符號單獨繪制。結(jié)果顯示,不論是K均值聚類還是K中心聚類,隨著抽頭數(shù)量的增加,稀疏擬合的效果都越來越好。將圖12 (a)、(b)的結(jié)果比較,可見本文提出的K中心聚類方法統(tǒng)計意義上的性能在抽頭數(shù)大于9之后開始優(yōu)于傳統(tǒng)的K均值聚類方法。

        (a) K均值聚類(a) K-means clustering

        (b) K中心聚類(b) K-medoids clustering圖12 具有不同期望抽頭數(shù)的鑒別器誤差偏差箱線圖Fig.12 Boxplot for deviation of discriminator errors with different numbers of desired taps

        值得一提的是,觀察箱線圖的離群值,發(fā)現(xiàn)有部分數(shù)據(jù)點偏差比較大。分析其原因:在時變運動場景下可能會存在聚類結(jié)構(gòu)偏差,不能很好地匹配信道CIR快照數(shù)據(jù)形態(tài),單一種類的聚類算法可能難以穩(wěn)定有效地聚類分析所有的CIR數(shù)據(jù)形態(tài)。這也是后續(xù)值得深入研究的問題。

        4 結(jié)論

        衛(wèi)星導航信道模型的電磁場形式傳播模型通常是極其煩瑣和難以模擬的。信道的沖擊響應(yīng)和傳遞函數(shù)比較清晰,但是還需要簡化以便工程應(yīng)用。如何把復雜信道模型簡化成無失真?zhèn)尉嗾`差、計算代價低的信道模型是信道建模與仿真的重要組成部分。本文采用稀疏抽頭濾波器結(jié)構(gòu)對模擬信道進行仿真,使仿真復雜度不受傳播路徑數(shù)目的影響,克服了復雜度較高的計算困難;提出了GNSS多徑信道模擬的聚類稀疏擬合方案,其中包括GNSS星座信道的稀疏TDL結(jié)構(gòu)模擬框架,采用較少抽頭數(shù)目的TDL結(jié)構(gòu)FIR濾波器稀疏擬合原始GNSS多徑信道模型,簡化了GNSS信道仿真的復雜度,同時保持了大部分模型的精確度;提出了基于K中心聚類CIR參數(shù)萃取的信道稀疏擬合方法,得到了等效精簡參數(shù)。上述研究結(jié)果為大規(guī)模低復雜度實時GNSS信道模擬提供了技術(shù)基礎(chǔ)。

        猜你喜歡
        方法模型
        一半模型
        重要模型『一線三等角』
        重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
        學習方法
        可能是方法不對
        3D打印中的模型分割與打包
        用對方法才能瘦
        Coco薇(2016年2期)2016-03-22 02:42:52
        FLUKA幾何模型到CAD幾何模型轉(zhuǎn)換方法初步研究
        四大方法 教你不再“坐以待病”!
        Coco薇(2015年1期)2015-08-13 02:47:34
        賺錢方法
        久久天天躁夜夜躁狠狠85麻豆| 亚洲AV色欲色欲WWW| 三级黄片一区二区三区| 国内精品少妇高潮视频| 久久综合九色综合97欧美| 亚洲精品国产美女久久久| 精品黑人一区二区三区| 美女免费视频观看网址| 国产麻豆精品一区二区三区v视界| 欧美国产成人精品一区二区三区| 久久99亚洲网美利坚合众国| 亚洲精品在线一区二区| 丰满熟女高潮毛茸茸欧洲视频| 破了亲妺妺的处免费视频国产| 2021最新久久久视精品爱| 亚洲日本中文字幕乱码在线| 国产精品a免费一区久久电影| 久久精品国产亚洲av电影| 日韩美女av二区三区四区| 免费一区二区在线观看视频在线| 亚洲乱码日产精品一二三| 一级免费毛片| 青青草手机成人自拍视频| 国产av自拍视频在线观看| 97成人碰碰久久人人超级碰oo | 亚洲无码视频一区:| 国产天堂av在线播放资源| 鲁丝片一区二区三区免费| 亚洲人成7777影视在线观看| 日美韩精品一区二区三区| 第一次处破女18分钟高清| 久久午夜伦鲁片免费无码| 亚洲AV成人综合五月天在线观看| 日本一区二区免费在线看| 国产99久久久国产精品免费看 | 在线国产丝袜自拍观看| 97碰碰碰人妻无码视频| 国产精品片211在线观看| 男女啪啪免费视频网址| 久久影院午夜理论片无码| 中日av乱码一区二区三区乱码|