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

基于二次虛擬擴展的高分辨率波達方向估計方法

2017-07-05 14:27:18李文興毛曉軍
哈爾濱工業大學學報 2017年5期
關鍵詞:方法

趙 宇,李文興,毛曉軍

(哈爾濱工程大學 信息與通信工程學院, 哈爾濱 150001)

基于二次虛擬擴展的高分辨率波達方向估計方法

趙 宇,李文興,毛曉軍

(哈爾濱工程大學 信息與通信工程學院, 哈爾濱 150001)

為進一步提高天線陣波達方向估計的分辨率,在四階量多重信號分類方法的基礎上,提出一種高分辨率的波達方向估計方法. 利用陣列接收數據的四階矩量進行虛擬陣列擴展,再利用陣列接收數據的共軛進行虛擬陣列擴展,實現二次虛擬擴展;將擴展后的陣列導向矢量和協方差矩陣用于波達方向估計,與原陣列導向矢量和協方差矩陣相比,相當于構造了更多的虛擬陣元,并擴展了陣列的孔徑. 仿真結果表明:與四階量多重信號分類等波達方向估計方法相比,所提出的方法在波達方向估計中成功概率更高,均方誤差更低,具有更高的分辨率. 所提方法通過二次虛擬擴展,構造了更多的虛擬陣元,有效地提高了天線陣波達方向估計的分辨率.

波達方向估計;多重信號分類;四階量法;虛擬陣列;分辨率

波達方向(direction-of-arrival, DOA)估計在雷達、聲吶、無線通信等眾多應用領域有著重要作用[1-2]. Schmidt等[3]提出的多重信號分類(multiple signal classification, MUSIC)是最具有代表性的DOA估計方法,是一種基于信源協方差矩陣特征分解的子空間類方法,具有較高的分辨率.

隨著科技的快速發展,對DOA估計的分辨率也提出更高的要求. 虛擬天線陣是一種研究構造虛擬陣元的方法和將實際天線陣轉換為虛擬天線陣的先進技術,主要包括內插變換法[4-5]和四階量法[6-8]. 將虛擬天線陣技術與MUSIC方法結合,可以增加天線陣的自由度,有效提高天線陣DOA估計的分辨率. 內插變換法通過在實際陣元間插入虛擬陣元,提高了MUSIC方法的分辨率[5],但是內插變換法不利于陣列孔徑擴展,對分辨率的提高有限. 基于四階量的MUSIC(fourth-order MUSIC, FOMUSIC)方法能夠有效地擴展陣列孔徑,并提高陣列的分辨率[6],是一種高分辨率的DOA估計方法.

FOMUSIC方法能夠擴展陣列孔徑,但是運算量較大[9]. Akkar等[10]利用正交分解代替特征分解; Lie等[11]用二階和四階協方差矩陣進行混合估計;Wang等[12]則去除了四階量矩陣中的部分冗余信息,都減小了FOMUSIC方法的運算量. Liao等[13]利用基于行列式的空間譜,實現了未知互耦情況下的DOA估計;宋海巖等[14]用一組四階量矩陣進行聯合對角化,能處理相干信號,且與單個四階量矩陣相比獲得了更高的分辨率;Li[15]等通過增強信號子空間,減小了內插變換法的變換誤差,與常規內插變換法相比分辨率更高,且運算量較小;Shan等[16]提出了共軛擴展多重信號分類(conjugate augmented MUSIC, CAM)方法,通過擴展陣列導向矢量以及協方差矩陣所含信息,大幅提高了分辨率,是一種分辨率極高的虛擬天線陣DOA估計方法.

為了進一步提高天線陣DOA估計的分辨率,本文提出一種基于二次虛擬擴展的多重信號分類(twice virtual expansion MUSIC, TVEM)方法. 該方法先求出陣列接收數據的四階矩量,以及相應擴展的陣列導向矢量,再求出陣列接收數據及其共軛,以及相應擴展的陣列導向矢量,實現二次虛擬擴展,構造了更多的虛擬陣元,進一步提高了天線陣DOA估計的分辨率.

1 天線陣模型及MUSIC方法

(1)

式中:M為信源數,(·)T表示轉置,A=[a1(θ),a2(θ), …,aM(θ)]為陣列導向矢量矩陣,S(t)為信號復包絡,N(t)為天線陣的噪聲.

陣列接收數據的協方差矩陣為

協方差矩陣Rxx是Hermitian矩陣,特征分解后為

Rxx=UΣUH,

式中Σ=diag[λ1,…,λN]是Rxx的特征值按降序排列組成的對角陣,U為特征向量矩陣.

將Rxx的特征向量根據特征值的大小分為信號和噪聲兩個部分:

式中US為大特征值所對應的信號子空間,UN為小特征值所對應的噪聲子空間.

根據協方差矩陣信號子空間與噪聲子空間的正交性,可以由噪聲子空間估計出信號的波達方向,在實際計算中,Rxx通常用有限次快拍數K的采樣協方差矩陣代替:

MUSIC算法的空間譜為

式中a(θ)為陣列導向矢量. 再通過對PMUSIC譜峰搜索即可估計出信號的波達方向.MUSIC方法是一種分辨率較高的DOA估計方法,但是當信源數超過實際陣元數時,MUSIC方法將會失效,此外,MUSIC方法的分辨率也受到實際陣列孔徑的限制.

2 FOMUSIC方法

用四階量替代二階量進行DOA估計,可以保留更豐富的波達方向信息,獲得更好的DOA估計性能[7]. 設x(t)為零均值平穩實隨機過程,則其四階累積量矩陣為

(2)

式中(·)*表示共軛. 式(2)可以表示為[7]C4x=E{[X(t)?X*(t)][X(t)?X*(t)]H}- E{X(t)?X*(t)}E{[X(t)?X*(t)]H}-E{X(t)XH(t)}?E{[X(t)XH(t)]*}.

(3)

式中?表示克羅內克積,第一項為四階矩量協方差矩陣:

D4x=E{[X(t)?X*(t)][X(t)?X*(t)]H}.

式中:X(t)為N×K維矩陣,X(t)?X*(t)為N2×K維矩陣,D4x為N2×N2維矩陣. 由式(3)可以得到陣列接收數據的四階協方差矩陣.

根據克羅內克積的性質[17],將式(1)代入X(t)?X(t)*:

X(t)?X*(t)=[AS(t)+N(t)]?

[AS(t)+N(t)]*=[AS(t)]?

[AS(t)]*+N(t)?N*(t)=

(A?A*)[S(t)?S*(t)]+

N(t)?N*(t).

(4)

再將式(4)代入D4x[7]:

D4x=E{[X(t)?X*(t)][X(t)?X*(t)]H}=

(A?A*)E{[S(t)?S*(t)]·

[S(t)?S*(t)]H}(A?A*)H+

E{[N(t)?N*(t)][N(t)?N*(t)]H}=

(A?A*)D4S(A?A*)H+D4N.

(5)

式中:A?A*為擴展后的陣列導向矢量矩陣,D4S和D4N分別為信號和噪聲的四階矩量矩陣.

與推導四階矩量矩陣的式(5)類似,可由式(3)推導出四階累積量矩陣為

C4x=(A?A*)C4S(A?A*)H+C4N,

(6)

式中C4S和C4N分別為信號和噪聲的四階累積量矩陣. 由式(5)和式(6)可知,與二階協方差矩陣相比,四階協方差矩陣的維度得到了擴展,所含波達方向信息也更豐富,實際陣列通過克羅內克積運算被轉換為虛擬陣列,陣列孔徑得到了擴展.

由式(3)得到C4x,再對其特征分解得到噪聲子空間U4N,根據噪聲子空間與信號子空間的正交性,可知,FOMUSIC方法的空間譜為

FOMUSIC方法將實際陣列轉換為虛擬陣列,擴展了陣列孔徑,能夠在信源數超過實際陣元數的情況下進行DOA估計,與MUSIC方法相比,分辨率也更高.

3 TVEM方法

為了進一步提高天線陣DOA估計的分辨率,本文在FOMUSIC方法的基礎上,提出了TVEM方法,該方法通過二次虛擬擴展,構造了更多的虛擬陣元,能夠獲得比FOMUSIC方法更高的分辨率.

先求出陣列接收數據的四階矩量,進行第一次虛擬擴展,陣列接收數據X(t)的四階矩量為

Xu1(t)=X(t)?X*(t).

(7)

根據克羅內克積的性質[17],將式(1)代入式(7)得

Xu1(t)=X(t)?X*(t)=[AS(t)+N(t)]?

[AS(t)+N(t)]*=[AS(t)]?[AS(t)]*+

N(t)?N*(t)=(A?A*)[S(t)?S*(t)]+

加入死亡話題的親子談話這個中介變量后,采用線性回歸分析來考察父母來生信念和死亡話題的親子談話與兒童來生信念的關系。依次將父母的來生信念(BA量表得分)、死亡話題親子談話(PDCA量表得分)納入回歸方程,進行回歸。采用Sobel檢驗方法對心理理論的中介效應進行檢驗,Z=-2.36,p<0.05。當將親子談話納入回歸模型時,父母來生信念的預測作用不再顯著,表明親子談話起到顯著、完全的中介作用(具體指標見表7)。

N(t)?N*(t),

式中Xu1為N2×K維矩陣,A?A*為擴展后的陣列導向矢量矩陣. 擴展后的陣列導向矢量為b(θ)=a(θ)?a*(θ)=[1,e,…,e]T?[1,e,…,e]H=

[1,e,…, e,e,…,e]T.

(8)

由式(8)可知,b(θ)中有大量冗余,只有2N-1個非重復元素,即第一次虛擬陣列擴展的陣元位置,如圖1所示.

圖1 第一次虛擬擴展的虛擬陣元位置

Fig.1 Position diagram of virtual array elements of the first extension

Fig.2 Position diagram of virtual array elements of the second extension

求出二次虛擬擴展后的擴展協方差矩陣:

(9)

對Ru進行特征分解:

式中UuN為噪聲子空間.

二次虛擬擴展后的虛擬陣列導向矢量為

(10)

由式(10)和圖1、圖2可知,通過二次虛擬擴展,實際陣列被轉換為虛擬陣列,陣列的導向矢量得到了進一步的擴展,構造了更多的虛擬陣元.TVEM方法的空間譜為

TVEM方法先求出陣列接收數據的四階矩量,以及相應擴展的陣列導向矢量,再求出陣列接收數據及其共軛,以及相應擴展的陣列導向矢量,通過二次虛擬擴展,將實際陣列轉換為虛擬陣列.TVEM方法分辨率提高的原因可以從數學和物理意義兩方面進行分析: 從數學上看,如式(6)、(9)所示,TVEM方法中協方差矩陣得到擴展,其中包含了更多的波達方向信息和冗余,從而提高了分辨率;從物理意義上看,如圖1、2所示,TVEM方法構造了更多的虛擬陣元,擴展了陣列孔徑,因此獲得了比FOMUSIC方法更高的分辨率.

4 運算量分析

將常規MUSIC方法、FOMUSIC方法、CAM方法[16]以及本文提出的TVEM方法的計算復雜度進行比較,如表1所示.

表1 4種方法計算復雜度對比

由表1可知,在未去冗余的情況下,FOMUSIC法、CAM法和TVEM法這3種基于四階量的虛擬天線陣DOA估計方法的計算復雜度遠高于MUSIC法,這是因為四階量矩陣中包含大量冗余,其中CAM法計算復雜度比TVEM法高,TVEM法計算復雜度比FOMUSIC法高. 在天線陣為均勻直線陣的情況下,能夠去除四階量矩陣中的冗余[16],此時TVEM法的計算復雜度與CAM方法十分接近,而TVEM法的計算復雜度明顯高于FOMUSIC法,這是因為TVEM法通過二次虛擬擴展,構造了比FOMUSIC法更多的虛擬陣元,從而獲得了更高的分辨率,其代價是運算量增大.

5 仿真分析

考慮遠場窄帶情況下的來波信號,信號之間互不相關,天線陣為四元均勻直線陣,陣元間距為半波長,快拍數為200. 下面通過三組仿真實驗,比較TVEM方法、FOMUSIC方法和CAM方法[16]的分辨率,三組仿真結果均通過200次重復獨立實驗得到. 分辨率的評價指標采用均方誤差和成功概率. 成功概率的計算方法為:獨立重復實驗中,波達方向估計的誤差低于門限值的次數與總次數之比,仿真中門限值設為0.2°.

5.1 在信源數少于實際陣元數的情況下,比較3種方法的分辨率

仿真1 信源數為3,來波方向為-45°、0°和40°,其中天線陣法線方向設為0°,左側為-90°到0,右側為0°到90°.

圖3為3種方法的均方誤差(root mean square error, RMSE)隨輸入信噪比RSN的變化曲線,其中RSN表示天線陣輸入的期望信號功率與噪聲功率之比;圖4為3種方法的成功概率隨輸入信噪比的變化曲線.

圖3 輸入信噪比變化時3種方法的均方誤差曲線

圖4 輸入信噪比變化時3種方法的成功概率曲線

Fig.4 Probability of target resolution versusRSNofthethreemethods

從圖3可以看出:當RSN低于3dB時,TVEM法的均方誤差高于CAM法,當RSN高于3dB時,TVEM法的均方誤差低于CAM法;總體上,TVEM法的均方誤差略低于CAM方法,且二者的均方誤差都明顯低于FOMUSIC法. 從圖4可以看出:當RSN低于3dB時,TVEM法的成功概率低于CAM法,當RSN高于3dB時,TVEM法的成功概率高于CAM法;總體上,TVEM法的成功概率略高于CAM法,且二者的成功概率都明顯高于FOMUSIC法. 由此可知,在信源數少于實際陣元數的情況下,TVEM法的分辨率略高于CAM法,明顯高于FOMUSIC法.

5.2 在波達角度較為接近的情況下,比較3種方法的分辨率

仿真2 信源數為2,來波方向為-3°和3°.

圖5為3種方法的均方誤差隨輸入信噪比的變化曲線. 圖6為3種方法的成功概率隨輸入信噪比的變化曲線.

從圖5可以看出,在波達角度較為接近的情況下,隨著RSN不斷提高,3種DOA估計方法的均方誤差都不斷下降,3種方法均方誤差曲線的變化過程都可以大致分為3個階段. 第一階段,RSN相對較低時,均方誤差較大且下降較為緩慢;第二階段,隨著RSN的提高,均方誤差快速下降;第三階段,RSN提高到一定值后,均方誤差相對較小且平穩下降. 這是因為當RSN較低時,DOA估計方法無法分辨波達角度較為接近的信源,因而均方誤差很大且下降較為緩慢;隨著RSN不斷提高,DOA估計方法可以在一定程度上分辨波達角度較為接近的信源,此時均方誤差快速下降;當RSN達到一定值后,DOA估計方法可以相對有效地分辨波達角度較為接近的信源,則均方誤差相對較小且平穩下降. 對比圖5中3種方法的均方誤差曲線可以知道,TVEM法的均方誤差曲線在RSN超過9dB后進入第三階段,CAM法的均方誤差曲線在RSN超過11dB后進入第三階段,FOMUSIC法的均方誤差曲線在RSN超過17dB后進入第三階段. 此外,在第二、三階段,TVEM法的均方誤差略低于CAM法,明顯低于FOMUSIC法. 由此可知,TVEM法的均方誤差略低于CAM法,明顯低于FOMUSIC法.

圖6 輸入信噪比變化時3種方法的成功概率曲線

Fig.6 Probability of target resolution versusRSNofthethreemethods

從圖6可以看出:當RSN低于5dB時,TVEM法的成功概率略低于CAM法,當RSN高于5dB時,TVEM法的成功概率高于CAM法;總體上,TVEM法的成功概率略高于CAM法,且二者的成功概率都明顯高于FOMUSIC法. 由圖5和圖6可知,在波達角度較為接近的情況下,TVEM法的分辨率略高于CAM法, 明顯高于FO-MUSIC法.

5.3 在信源數超過實際陣元數的情況下,比較3種方法的分辨率

仿真3 信源數為5,來波方向為-50°、-25°、0°、15°和45°.

圖7為3種方法的均方誤差隨輸入信噪比的變化曲線,圖8為3種方法的成功概率隨輸入信噪比的變化曲線.

圖7 輸入信噪比變化時3種方法的均方誤差曲線

圖8 輸入信噪比變化時3種方法的成功概率曲線

Fig.8 Probability of target resolution versusRSNofthethreemethods

從圖7可以看出,TVEM法的均方誤差低于CAM法,CAM法的均方誤差低于FOMUSIC法. 從圖8可以看出,TVEM法的成功概率高于CAM法,CAM法的成功概率高于FOMUSIC法. 由此可知,在信源數超過實際陣元數的情況下,TVEM法的分辨率略高于CAM法,明顯高于FOMUSIC法.

由上述仿真結果可知,TVEM方法的分辨率略高于CAM方法,且TVEM方法的分辨率明顯高于FOMUSIC方法. 因此,TVEM方法是進一步提高天線陣DOA估計分辨率的一種有效方法.

6 結 論

本文提出的基于二次虛擬擴展的高分辨率波達方向估計方法,通過陣列接收數據的四階矩量和共軛,擴展了陣列的導向矢量和協方差矩陣,構造出更多的虛擬陣元,并擴展了陣列的孔徑. 將實際天線陣轉換為虛擬天線陣,實現了二次虛擬擴展,進一步提高了天線陣DOA估計的分辨率. 仿真結果表明,本文所提出的方法與四階量多重信號分類等DOA估計方法相比,具有更低的均方誤差和更高的成功概率,是一種提高天線陣DOA估計分辨率的有效方法.

[1] WANG Yi, YANG Minglei, CHEN Baixiao, et al. Improved DOA estimation based on real-valued array covariance using sparse Bayesian learning [J]. Signal Processing, 2016, 129(C): 183-189. DOI: 10.1016/j.sigpro.2016.06.002.

[2] 閆鋒剛, 張薇, 金銘. 求根MUSIC初值設置和更新算法 [J]. 哈爾濱工業大學學報, 2015, 47(3): 88-92. DOI: 10.11918/j.issn.0367-6234.2015.03.015.

YAN Fenggang, ZHANG Wei, JIN Ming. A new method for setting and updating theinitiation of root-MUSIC [J]. Journal of Harbin Institute of Technology, 2015, 47(3): 88-92. DOI: 10.11918/j.issn.0367-6234.2015.03.015.

[3] SCHMIDT R O. Multiple emitter location and signal parameter estimation [J]. IEEE Transactions on Antennas and Propagation, 1986, 34(3): 276-280. DOI: 10.1109/TAP.1986.1143830.

[4] BELLONI F, RICHTER A, KOIVUNEN V. DOA estimation via manifold separation for arbitrary array structures [J]. IEEE Transactions on Signal Processing, 2007, 55(10): 4800-4810. DOI: dx.doi.org/10.1109/TSP.2007.896115

[5] WANG Yongliang, CHEN Hui, WAN Shanhu. An effective DOA method via virtual array transformation [J]. Science in China Series E, 2001, 44(1): 75-82. DOI: 10.1007/BF02916727.

[6] LI Xiukun, LI Tingting, GU Xinyu, et al. Array gain of fourth-order cumulants beamforming under typical probability denisity background [J]. Chinese Journal of Acoustics, 2015, 34(1): 15-26. DOI: 10.15949/j.cnki.0217-9776.2015.01.002.

[7] 丁齊, 魏平, 肖先賜. 基于四階累積量的DOA估計方法及其分析 [J]. 電子學報, 1999, 27(3): 25-28. DOI: 10.3321/j.issn:0372-2112.1999.03.008.

DING Qi, WEI Ping, XIAO Xianci. Estimation and analysis of DOA based on fourth-order cumulant [J]. ACTA Electronica SINICA, 1999, 27(3): 25-28. DOI: 10.3321/j.issn:0372-2112.1999.03.008.

[8] ZHANG Xin, LIU Xiaoming, YU Haixia, et al. Improved MUSIC algorithm for DOA estimation of coherent signals via toeplitz and fourth-order-cumulants [J]. International Journal of Control and Automation, 2015, 8(10): 261-272. DOI: 10.14257/ijca.2015.8.10.25.

[9] 陳建, 王樹勛. 基于高階累積量虛擬陣列擴展的DOA估計 [J]. 電子與信息學報, 2007, 29(5): 1041-1044.DOI: 10.3969/j.issn.1671-5896.2006.04.001

CHEN Jian, WANG Shuxun. DOA estimation of virtual array extension based on fourth-order cumulant [J]. Journal of Electronics and Information Technology, 2007, 29(5): 1041-1044. DOI: 10.3969/j.issn.1671-5896.2006.04.001.

[10]AKKAR S, HARABI F, GHARSALLAH A. MUSIC like algorithms for fast direction of arrival estimation [C]// 1stInternational Conference on Advanced Technologies for Signal and Image Processing. Sousse: IEEE, 2014: 550-554. DOI: 10.1109/ATSIP.2014.6834675.

[11]LIE S, LEYMAN A R, CHEW Y H. Fourth-order and weighted mixed order direction-of-arrival estimators [J]. IEEE Signal Processing Letters, 2006, 13(11): 691-694. DOI: 10.1109/LSP.2006.879456.

[12]WANG Qing, CHEN Hua, ZHAO Guohuang,et al. An improved direction finding algorithm based on toeplitz approximation [J]. Sensors, 2013, 13(1): 746-757. DOI: 10.3390/s130100746.

[13]LIAO Bin, CHAN S C. A cumulant-based method for direction finding in uniform linear arrays with mutual coupling [J]. IEEE Antennas and Wireless Propagation Letters, 2014, 13: 1717-1720. DOI: 10.1109/LAWP.2014.2352939.

[14] 宋海巖, 樸勝春. 基于高階累積量矩陣組正交聯合對角化的高分辨率方位估計方法 [J]. 電子與信息學報, 2010, 32(4): 967-972. DOI: 10.3724/SP.J.1146.2008.01176.

SONG Haiyan, PIAO Shengchun. DOA estimation method based on orthogonal joint diagonalization of high-order cumulant [J]. Journal of Electronics and Information Technology, 2010, 32(4): 967-972. DOI: 10.3724/SP.J.1146.2008.01176.

[15]LI Wenxing, MAO Xiaojun, YU Wenhua, et al. An effective technique for enhancing direction finding performance of virtual arrays [J]. International Journal of Antenna and Propagation, 2014(3):1-7. DOI: 10.1155/2014/728463.

[16]SHAN Zhilong, YUM T S P. A conjugate augmented approach to direction-of-arrival estimation [J]. IEEE Transactions on Signal Processing, 2005, 53(11): 4104-4109. DOI: 10.1109/TSP.2005.857012.

[17]BREWER J W. Kronecker products and matrix calculus in system theory [J]. IEEE Transaction on Circuits and Systems, 1978, 25(9): 772-781. DOI: 10.1109/TCS.1978.1084534.

(編輯 王小唯, 苗秀芝)

High resolution approach of direction-of-arrival estimation based on twice virtual extension

ZHAO Yu, LI Wenxing, MAO Xiaojun

(College of Information and Communications Engineering, Harbin Engineering University, Harbin 150001, China)

In order to further improve the resolution of direction-of-arrival (DOA) estimation of antenna array, a high resolution approach of DOA estimation is proposed based on the fourth-order multiple signal classification (FO-MUSIC) approach. First, the antenna array is extended through the fourth-order moment of the

data. Then, the antenna array is extended through the conjugate value of the received data. When the expanded steering vector and the expanded covariance matrix are used for DOA estimation instead of the original steering vector and the original covariance matrix, virtual array elements are formed and the aperture of array is extended. It is shown by simulation results that, compared to the FO-MUSIC approach, the proposed method has a higher probability of target resolution, lower root mean square error and higher resolution. Hence, more virtual array elements are formed by the proposed approach based on twice virtual extension, and the resolution of DOA estimation of antenna array is effectively improved.

direction-of-arrival estimation; multiple signal classification; fourth-order statistics; virtual array; resolution

10.11918/j.issn.0367-6234.201612032

2016-12-08

國防預研項目(4010403020102, 4010103020103)

趙 宇(1987—),男,博士研究生; 李文興(1960—),男,教授,博士生導師

李文興,liwenxing@hrbeu.edu.cn

TN911.7

A

0367-6234(2017)05-0122-06

猜你喜歡
方法
中醫特有的急救方法
中老年保健(2021年9期)2021-08-24 03:52:04
高中數學教學改革的方法
河北畫報(2021年2期)2021-05-25 02:07:46
化學反應多變幻 “虛擬”方法幫大忙
變快的方法
兒童繪本(2020年5期)2020-04-07 17:46:30
學習方法
用對方法才能瘦
Coco薇(2016年2期)2016-03-22 02:42:52
最有效的簡單方法
山東青年(2016年1期)2016-02-28 14:25:23
四大方法 教你不再“坐以待病”!
Coco薇(2015年1期)2015-08-13 02:47:34
賺錢方法
捕魚
主站蜘蛛池模板: 丝袜国产一区| 麻豆精品在线| 国产成人亚洲综合A∨在线播放| 欧美在线精品怡红院| 亚洲天堂精品在线| 在线毛片免费| 成人国产精品网站在线看| 精品国产欧美精品v| 国产一级在线播放| 欧美综合成人| 99爱视频精品免视看| 欧美不卡在线视频| 狠狠干综合| 精品国产一区二区三区在线观看| 黄色福利在线| 综合天天色| V一区无码内射国产| 国产95在线 | 99视频有精品视频免费观看| 第九色区aⅴ天堂久久香| 嫩草影院在线观看精品视频| 亚洲视频无码| 伊伊人成亚洲综合人网7777| 亚洲AV成人一区二区三区AV| 日本成人在线不卡视频| av天堂最新版在线| 欧美区国产区| 精品视频在线一区| AV不卡在线永久免费观看| 国产噜噜噜| 亚洲精品第一在线观看视频| 欧美亚洲国产精品第一页| 在线无码av一区二区三区| 一区二区理伦视频| 狠狠做深爱婷婷综合一区| 中文成人无码国产亚洲| 久久毛片基地| 久久综合五月婷婷| 在线观看无码a∨| 色有码无码视频| 国产精品第5页| 特级精品毛片免费观看| 国产在线视频导航| 国产亚洲一区二区三区在线| 欧美日韩国产成人高清视频| 成人国产精品网站在线看| 色婷婷视频在线| 久久一日本道色综合久久| 麻豆精品视频在线原创| 91无码视频在线观看| 亚洲欧美一级一级a| 精品久久久久无码| 精品在线免费播放| AV无码无在线观看免费| 五月激激激综合网色播免费| 日韩欧美国产三级| 一级看片免费视频| 国产精品午夜福利麻豆| 久久国产乱子伦视频无卡顿| 97视频免费看| 国产美女视频黄a视频全免费网站| 伊人天堂网| 日韩二区三区| 五月婷婷精品| 福利国产在线| 亚洲综合网在线观看| 国产爽爽视频| 99免费视频观看| 国产激情影院| 日本午夜在线视频| 在线观看亚洲人成网站| 亚洲福利一区二区三区| 久久综合九九亚洲一区| 亚洲一区精品视频在线| 国产在线观看99| 在线无码av一区二区三区| 日韩a级片视频| 欧美日韩国产一级| 国产内射在线观看| 婷婷99视频精品全部在线观看| 69综合网| 狠狠亚洲五月天|