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

基于降維稀疏重構的相干信源二維DOA估計方法

2016-08-15 07:28:38王秀紅毛興鵬張乃通
系統工程與電子技術 2016年8期
關鍵詞:方法

王秀紅, 毛興鵬, 張乃通

(1. 哈爾濱工業大學電子與信息工程學院, 黑龍江 哈爾濱 150001;2. 哈爾濱工業大學(威海)信息與電氣工程學院, 山東 威海 264209;3. 哈爾濱工業大學信息感知技術協同創新中心, 黑龍江 哈爾濱 150001)

?

基于降維稀疏重構的相干信源二維DOA估計方法

王秀紅1,2, 毛興鵬1,3, 張乃通1

(1. 哈爾濱工業大學電子與信息工程學院, 黑龍江 哈爾濱 150001;2. 哈爾濱工業大學(威海)信息與電氣工程學院, 山東 威海 264209;3. 哈爾濱工業大學信息感知技術協同創新中心, 黑龍江 哈爾濱 150001)

直接將壓縮感知(compressed sensing,CS)思想應用到相干信源二維波達方向(direction of arrival,DOA)估計中會帶來高計算復雜度的問題。為了解決這一問題,提出了一種基于降維稀疏重構的二維DOA估計方法,該方法利用特殊陣列結構將二維冗余字典構建問題轉化為一維冗余字典的構建,同時提出了一種基于子字典空間譜重構的配對算法,從而在極大降低算法計算復雜度的同時,提高了配對成功概率。仿真結果表明,該方法對相干信源具有接近于克拉美羅下界(Cramér-Rao lower bound, CRLB)的估計性能,即使是在低信噪比、少快拍數和小角度間隔的情況下,仍有良好的估計性能。

波達方向; 二維波達方向估計; 相干信源; 降維稀疏重構; 冗余子字典

0 引 言

近年來高精度、高分辨率的波達方向(direction of arrival,DOA)估計方法得到了快速的發展并且成功應用于雷達、聲納和無線通信等領域[1]。這些方法絕大多數都是針對多信源的一維方向(比如方位角)估計。但是在很多應用中目標到達角卻具有二維特性(比如方位角和俯仰角),因此二維DOA估計方法得到了廣泛的關注和深入的研究[2-13]。

經典的子空間類方法已經成功地應用到了二維DOA估計中[4-5],比如二維多重信號分類(multiple signal classification,MUSIC)算法[4]和二維旋轉不變技術估計信號參數(estimate signal parameter via rotational invariance technique,ESPRIT)算法[5]等,雖然該類方法具有高分辨率、高估計精度等優點,但它們的應用前提也較為苛刻,比如要求大量的快拍數、信源之間的獨立性、噪聲為不相關的白噪聲等。很多實際應用環境都無法滿足上述條件,這使得很多基于子空間類的估計方法因性能嚴重下降而失效。因此,在各種復雜的環境和復雜電磁干擾的背景下,如何抑制強相關甚至相干信源的影響并且快速估計出目標信源的二維DOA角度,是一個具有重要的實際應用價值的課題方向。在這方面國內外的學者已經取得了一定的研究進展[6-21]。文獻[7]采用了前后向空間平滑(forward backward spatial smoothing, FBSS)[6]技術對信號去相干,而后利用二維ESPRIT 算法進行角度估計;文獻[8-9]通過對接收數據的互協方差矩陣重構實現解相干,重建信號子空間后采用ESPRIT等算法估計目標DOA;文獻[10]在三維垂直陣列上利用波達方向矩陣法實現解相干,并將二維DOA估計問題轉化為3個一維DOA估計問題進行處理;文獻[11]通過對接收數據構建成Toeplitz結構的矩陣實現了信源的解相干;文獻[12]通過將接收數據的四階累積量構建成Toeplitz矩陣結構實現信源的解相干;文獻[13]提出利用斜投影的方法將相干信源和非相干信源分離,從而實現所有信源的二維角度估計。

近年來發展起來的壓縮感知(compressed sensing,CS)技術也已經成功應用到了DOA估計中[14-19]。該類方法由于不需要利用接收數據的統計特性,因而無需任何預處理即可對相干信源實現DOA估計;同時由于不需要進行特征分解和空間平滑技術,并且可單快拍處理,所以具有低計算復雜度的優點;因此,該類方法具有重要的實用價值和廣泛的應用前景。目前基于CS的DOA估計方法主要應用于均勻線陣[14-15],也即一維DOA估計中。而針對面陣的二維CS DOA估計目前文獻較少,其主要原因之一就是CS方法中需要構建的二維角度冗余字典長度過長,使得算法的復雜度太高,難以實現,從而無法直接將其應用于二維DOA估計中。為了解決這個問題,文獻[16]利用L型陣列的特殊結構,提出了一種空間角度稀疏表示的方法,將二維DOA角度解耦并轉變為兩個一維DOA估計的問題,但是該方法的配對算法在信源功率相近甚至相等時將會失效。文獻[17]對接收數據的互相關矩陣進行稀疏表示來實現二維角度估計,但是無法實現對相干信源的估計。

為了解決相干信源二維CS DOA估計的低復雜度問題,本文提出了一種基于降維稀疏重構(reduction dimension sparse reconstruction,RDSR)的DOA估計方法。該方法利用十字型陣列結構特點,將二維DOA角度冗余字典的構建轉化為兩個一維的角度冗余字典的構建,從而大大降低算法的計算復雜度,并且提出了基于冗余子字典的二維空間譜重構的配對方案,從而克服了文獻[16]配對算法的缺陷。

因此,本文方法的優勢主要體現在以下幾個方面:①無需空間平滑等預處理即可實現相干信源的估計;②通過降維處理將二維問題轉化為兩個一維問題并行處理,大大降低算法復雜度和運行時間;③該方法對快拍數目要求不高,即使單快拍數據也可以完成二維DOA估計。

1 數學模型

1.1系統模型

考慮有K個遠場信源入射到一個十字型陣列,如圖1所示,該陣列位于x-z平面,由兩個分別位于x軸和z軸上的均勻線陣垂直交叉構成,其中位于x軸和z軸子陣列的陣元數分別為M和N,陣元間距均勻且分別為dx和dz。假設第k個信源入射到天線陣列的方位角為θk,俯仰角為φk,其中方位角定義為入射信號在xoy平面上的投影與x軸的夾角,俯仰角定義為入射信號與其在xoy平面投影之間的夾角。

圖1 陣列系統模型

該十字陣列(共M+N-1個陣元)的接收信號可以表示為

(1)

式中,S(t)=[s1(t),s2(t),…,sK(t)]T為K個入射信號的矢量;N(t)=[n1(t),n2(t),…,n(M+N-1)(t)]T為陣列接收噪聲矢量;A為(M+N-1)×K維陣列流形矩陣,可表示為

(2)

式中,a(θk,φk)為導向矢量,其各元素表示為

(3)

式中,xi和zi為第i個陣元在xoz平面的坐標位置。

1.2基于CS的二維DOA估計

首先將式(1)的接收信號模型轉化為稀疏化模型。利用信源在空域上具有稀疏性的特點,信源矢量S(t)可按照空域角度進行稀疏化表示為稀疏的(Lθ×Lφ)×1維矢量P(t),其中Lθ為待搜索的方位角數目,Lφ為待搜索的俯仰角數目。S(t)稀疏化表示后則需要將陣列流形矩陣A進行列擴展,生成角度冗余字典Ψ,其由所有可能的Lθ×Lφ個二維DOA角度的導向矢量組成。

(4)

因此,二維DOA估計中接收信號的稀疏化模型可表示為

(5)

式中,R=[R(t1),R(t2),…,R(tT)]∈C(M+N-1)×T和P=[P(t1),P(t2),…,P(tT)]∈CL×T分別為陣列接收信號與入射信號的多快拍數據矩陣(其中LLθ×Lφ,T為快拍數)。這樣,式(1)的DOA估計問題則轉化為對式(5)的稀疏解求解問題,可以通過無約束最小化問題來求解。

(6)

為了描述方便,定義CS框架下的空間譜為

(7)

這樣,估計出的二維DOA角度可表示為

(8)

綜上可以發現,二維CSDOA估計的關鍵和主要的計算量都在于式(6)的求解,而該式求解的復雜度與冗余字典Ψ的維數Lθ×Lφ密切相關,如采用內點法求解,式(6)的求解復雜度約與O(Lθ×Lφ)3[15]成正比。當需要搜索的方位角和俯仰角的數目較多時,冗余字典Ψ的維數Lθ×Lφ很大,此時式(6)的求解復雜度將會是巨大的,并且在實際系統中難以實現,因此限制了該類方法的應用。

2 基于降維稀疏重構二維DOA估計

為了降低上述二維DOA估計算法的復雜度,可以考慮縮減冗余字典Ψ的維數。為此,本文通過將二維DOA角度解耦,把二維DOA估計問題轉化為兩個一維DOA估計問題,從而大大減小冗余字典Ψ的維數。

2.1降維處理

考慮將圖1所示的十字陣看做是兩個獨立的均勻線陣,即將z軸的子陣列視為一個獨立的均勻線陣,可以使用該子陣列估計出入射信號與該子陣列之間的夾角,該夾角為π/2-φ,因此可以利用z軸子陣列估計出俯仰角φ∈[-π/2,π/2];同理,可以使用x軸的子陣列估計出入射信號與該子陣列之間的夾角,即入射信號與x軸的夾角α(為了表述方便,將夾角α稱為合成角),其與方位角θ和俯仰角φ之間有如下關系

(9)

因此,只要能夠利用子陣列求出俯仰角φ和合成角α,就可根據式(9)求出方位角θ∈[0,π]。

將z軸和x軸子陣列接收的信號分別表示為

(10)

類似于第1.2節中的稀疏化過程,由式(10)得到其稀疏表示模型

(11)

式中,ΦZ和ΦX分別為z軸和x軸子陣列的冗余字典,由AZ和AX擴展而成,即

(12)

(13)

可將式(11)的單快拍模型推廣到多快拍模型

(14)

(15)

式(15)的稀疏解可表示為

(16)

需要注意的是,式中F范數部分代表重構誤差, l1范數部分代表稀疏性,即DOA估計誤差,而參數μ1和μ2則的選擇方法同式(6),它們用來控制重構誤差與稀疏性之間的折中關系,即調節兩者在整體誤差中所占的比重。

利用求解出的稀疏解分別生成一維空間譜

(17)

找出K個最大的譜峰對應的角度,即為各信源的俯仰角和合成角

(18)

2.2角度配對存在的問題

(19)

同時,我們還發現在求解式(16)的問題中,參數μ1和μ2的選擇對估計誤差和配對成功概率影響較大。對于文獻[16]中基于幅度排序的配對方案,為了保證盡可能大的配對概率,因此在選擇參數μ時要求幅度估計誤差項所占比重較大(即要使幅度估計誤差最小),這樣位置估計誤差項比重就會較小,從而導致DOA估計誤差較大。所以這種配對算法會出現在信噪比較高的情況下DOA估計誤差較大的問題。

2.3基于子字典的角度配對方案

為進一步提高角度配對的成功概率和DOA估計性能,提出一種基于冗余子字典空間譜重構的角度配對方案。

首先需要構建冗余子字典(為了方便區分,將上節提到的冗余字典稱為完全冗余字典,簡稱全字典,而此處構建的冗余子字典簡稱為子字典),其具體構建方法如下。

假設系統中有K個入射信源,從第2.1節的方法中估計出它們的俯仰角為{φ1,…,φK}和合成角為{α1,…,αK},那么按照排列組合方案將形成I=K2個組合方案,將其表示為

(20)

(21)

(22)

由集合Ω內所有角度生成的導向矢量構成矩陣

(23)

(24)

由于子字典是全字典的一部分,仍然滿足式(5)的DOA稀疏模型,故有

(25)

因此,可重構出壓縮感知框架下的二維空間譜為

(26)

(27)

需要注意的是,由于在每個子字典下至多只有一個信源存在,因此對于子字典半徑Dθ,Dφ和間隔Δθ,Δφ的選取只需保證空間譜的稀疏性即可;同時,由于至多只有一個信源存在,所以對式(25)求解稀疏解算法的精度和分辨率要求也不高,因此除了式(6)的方法之外,還可以采用低復雜度的貪婪算法,比如正交匹配追蹤(orthogonalmatchingpursuit,OMP)[14]算法等。

圖2 基于子字典的配對方案舉例(兩個相干信源)

以上所述的基于子字典的配對方案,基本不受2.1節重構算法誤差的影響,并且對信噪比要求不高,對于信源功率相等或相近的情況仍適用,從而可以大大提高角度配對成功概率。

綜上,本文提出的RDSR算法的流程歸納如下:

步驟 1根據式(12)和式(13)構建子陣列的一維冗余字典ΦZ和ΦX;

3 數值仿真

為了驗證本文方法的有效性,本文將從下面幾個方面來仿真分析算法的性能。仿真中采用圖1所示的十字型陣列,其中每個子陣列的陣元數M=N=11,陣元間隔為半波長。仿真中對信源的估計精度采用均方根誤差(rootmeansquareerror,RMSE)來表征,其具體定義為

3.1信源幅度相同時DOA估計結果

圖3為100次獨立實驗的3個相干信源二維DOA估計結果,其中每個信源的幅度相同,SNR為0dB,快拍數為100。從圖3可以看出,當信源的幅度相同時,文獻[16]的方法會出現錯誤配對的情況,而導致算法失效,而本文方法不存在這種情況,可以較準確的估計出信源的二維DOA角度。

圖3 二維DOA估計結果對比

3.2信噪比對算法性能影響

圖4給出了2個相干信源在幅度比為1∶1情況下各種算法的估計誤差和配對成功概率比較。仿真中兩個信源角度分別為[θ1,φ1]=[63°,15°],[θ2,φ2]=[107°,42°],快拍數為100。從圖4可以看出,當信源幅度相同時,文獻[16]的方法由于配對失敗而導致估計性能很差,即這種情況下該方法已失效,而文獻[7]中基于FBSS的方法和本文方法都能較準確地估計出信源角度,但在低信噪比下,本文方法的估計性能則優于FBSS方法,且性能更接近CRB。同時,相對于FBSS方法只適用于由均勻線陣構成的陣列這一缺點,本文算法則也可以適用于由非均勻線陣構成的陣列,但由于篇幅有限本文不做深入討論。

圖4 相干信源DOA估計性能隨信噪比變化曲線

3.3信源角度間隔對算法性能影響

圖5給出了兩個相干信源在幅度比為0.8∶1情況下二維DOA估計誤差與信源角度間隔的變化關系。仿真中角度間隔Δ定義如下:如果設第一個信源的DOA角度為(θ1,φ1),那么與其角度間隔Δ的第二個信源的DOA角度(θ1+Δ,φ1+Δ)。仿真中SNR為10 dB,快拍數為100。從仿真結果來看,在3種不同方法中,本文方法的估計性能最好,且接近于CRB界;文獻[16]方法在角度間隔小于10°以內時估計性能優于FBSS方法,而在角度間隔大于10°以上時不如FBSS。

圖5 相干信源DOA估計誤差隨信源角度間隔變化曲線

3.4快拍數對算法性能影響

圖6給出了SNR分別為-10 dB、0 dB和10 dB情況下二維DOA估計誤差與快拍數的變化關系曲線。仿真中兩個相干信源的幅度比為1∶1,信源角度同上。從仿真結果來看,在不同的信噪比下和不同的快拍數目條件下,本文算法都具有接近于CRB的估計性能,而FBSS方法則在低信噪比情況下估計性能很差,即使增加快拍數,性能也得不到明顯改善。同時值得注意的是,本文算法在單快拍的情況下也可工作,而FBSS算法不能。

圖6 相干信源DOA估計誤差隨快拍數變化曲線

4 復雜度分析

在采用CS理論進行二維DOA估計時,如果直接進行二維角度搜索,那么DOA估計算法的復雜度將與角度冗余字典的規模大小密切相關。同時,估計中采用不同的稀疏重構算法復雜度也不相同。結合本文二維DOA估計背景,將幾種常見的稀疏分解算法(OMP算法[14]和L1-SVD算法[15])的復雜度如表1所示。

表1 幾種常用算法與本文算法的復雜度對比

注:K為信源數;M為陣元總數;T為快拍數;Lθ,Lφ和Lsd分別為方位角、俯仰角和子字典的搜索角度數。

表2給出了不同DOA估計算法的運行時間和估計誤差的對比,其中“運行時間”是由計算機通過Matlab軟件進行100次統計算法的運行時間并平均后得到的。假設2個相干信源,信源角度同第3.2節,子陣元數為M=N=11,SNR=10dB,快拍數為100,方位角和俯仰角的冗余字典大小都為400,角度間隔均為0.2°。從表2中可以看出,降維處理可以極大地降低算法的復雜度,并且幾乎不會帶來估計精度的損失,同時,本文方法與文獻[16]的方法比較,算法復雜度相差不多,但是本文方法的估計精度更高。

表2 不同DOA估計算法的運行時間和估計誤差對比

5 結 論

本文提出了一種RDSR的DOA估計方法,以實現低復雜度的相干信源二維DOA估計問題,同時該方法也成功解決了當信源幅度或功率相同或相近時配對算法失效的問題。仿真結果表明,本文方法在低信噪比、小角度間距以及少快拍數情況下都具有良好的估計性能,并且估計誤差可以接近于克拉美羅界。另外,與直接應用CS進行二維DOA搜索的估計方法相比,本文方法具有更低的計算復雜度。

[1] Krim H, Viberg M. Two decades of array signal processing research: The parametric approach[J].IEEESignalProcessingMagazine, 1996, 13(4): 67-94.

[2] Nie X, Li L P. A computationally efficient subspace algorithm for 2-D DOA estimation with L-Shaped array[J].IEEESignalProcessingLetters, 2014, 21(8): 971-974.

[3] Gu J F, Zhu W P, Swamy M N S. Joint 2-D DOA estimation via sparse L-shaped array[J].IEEETrans.onSignalProcessing, 2015, 63(5): 1171-1182.

[4] Chan A Y J, Litva J. MUSIC and maximum techniques on two-dimensional DOA estimation with uniform circular array[J].IEEProceedings-Radar,SonarandNavigation, 1995, 142(3): 105-114.

[5] Kedia V S, Chandna B. A new algorithm for 2-D DOA estimation[J].SignalProcessing, 1997, 60(3): 325-332.

[6] Pillai S U, Won B H K. Forward/backward spatial smoothing techniques for coherent signals identification[J].IEEETrans.onAcousticsSpeechandSignalProcessing, 1989, 37(1): 8-15.

[7] Cheng Q, Zhao Y, Yang J. A Novel 2-D DOA estimation for coherent signals based on L-shaped array[C]∥Proc.oftheInternationalConferenceonWirelessCommunicationandSignalProcessing, 2012, 7363(1): 1-5.

[8] Palanisamy P, Kalyanasundaram N, Swetha P M. Two-dimensional DOA estimation of coherent signals using acoustic vector sensor array[J].SignalProcessing, 2012, 92(7): 19-28.

[9] Liu Z T, Ruan X Y, He J. Efficient 2-D DOA estimation for coherent sources with a sparse acoustic vector-sensor array[J].MultidimensionalSystemsandSignalProcessing, 2013, 24(1): 105-120.

[10] Wang L, Li G L. A new method for estimating 2-D DOA in coherent sources environment using single snapshot[J].TransactionsofBeijingInstituteofTechnology, 2015, 35(5): 512-518. (王凌, 李國林. 利用單次快拍實現相干信源二維測向的新算法[J].北京理工大學學報, 2015, 35(5): 512-518.)

[11] Chen H, Zhang X. Two-dimensional DOA estimation of coherent sources for acoustic vector-sensor array using a single snapshot[J].WirelessPersonalCommunications, 2013, 72(1): 1-13.

[12] Chen H, Hou C P, Wang Q, et al. Cumulants-based toeplitz matrices reconstruction method for 2-D coherent DOA estimation[J].IEEESensorsJournal, 2014, 14(8): 2824-2832.

[13] Tao H, Xin J M, Wang J S, et al. Two-dimensional direction estimation for a mixture of no coherent and coherent signals[J].IEEETrans.onSignalProcessing, 2015, 63(2): 318-333.

[14] Tropp J A, Gilbert A C. Signal recovery from random measurements via orthogonal matching pursuit[J].IEEETrans.onInformationTheory, 2007, 53(12): 4655-4666

[15] Malioutov D, Cetin M, Willsky A. S. A sparse signal reconstruction perspective for source location with sensor arrays[J].IEEETrans.onSignalProcessing, 2005, 53(8): 3010-3022.

[16] Li P F, Zhang M, Zhong Z F. Two dimensional DOA estimation based on sparse representation of space angle[J].JournalofElectronics&InformationTechnology, 2011, 33(10): 2402-2406. (李鵬飛,張旻, 鐘子發. 基于空間角稀疏表示的二維DOA估計[J].電子與信息學報,2011,33(10):2402-2406.)

[17] Luo X Y, Fei X C, Wei P, et al. Sparse representation-based method for two-dimensional direction-of-arrival estimation with L-shaped array[C]∥Proc.oftheIEEEInternationalConfe-renceonCommunicationProblem-solving, 2014: 277-280.

[18] Northardt E T, Bilik I, Abramovich Y I. Spatial compressive sensing for direction-of-arrival estimation with bias mitigation via expected likelihood[J].IEEETrans.onSignalProcessing, 2013, 61(5): 1183-106.

[19] Massa A, Rocca P, Oliveri G. Compressive sensing in electromagnetics-a review[J].IEEEAntennasandPropagationMagazine, 2015, 57(1): 224-238.

[20] Hansen P C, O’ Leary D P. The use of the L-curve in the regularization of discrete ill-posed problems[J].SiamJournalonScientificComputing, 1993, 14(6): 1487-1503.

[21] Heng Y, Lu S, Mhamdi A, Pereverzev S V. Model functions in the modified L-curve method case study: the heat flux reconstruction in pool boiling[J].InverseProblems, 2010, 26(5): 55006-55018.

Two-dimensional DOA estimation for coherent sources based on reduction dimension sparse reconstruction

WANG Xiu-hong1,2, MAO Xing-peng1,3, ZHANG Nai-tong1

(1. School of Electronics and Information Engineering, Harbin Institute of Technology, Harbin 150001, China; 2. School of Information and Electrical Engineering, Harbin Institute of Technology(Weihai), Weihai 264209, China; 3. Collaborative Innovation Center of Information Sensing and Understanding at Harbin Institute of Technology, Harbin 150001, China)

The problem of high computational complexity will be caused if compressed sensing (CS) is directly applied to two-dimensional (2-D) direction of arrival (DOA) Estimation of coherent sources. To solve this problem, a 2-D DOA estimation method based on reduction dimension sparse reconstruction (RDSR) is proposed. The proposed method converts the construction of a 2-D redundancy dictionary into that of a 1-D dictionary by using the special array structure. In addition, a pair-matching scheme is proposed based on spatial spectrum reconstruction of the sub-dictionary. Therefore, the proposed method not only reduces the computational complexity but also improves the pairing probability of success. Simulation results show that the estimated performance of the method is close to the Cramér-Rao lower bound (CRLB), even in the case of low signal-to-noise ratio (SNR), small number of snapshots and small angle interval, the estimation performance is still good.

direction of arrival (DOA); two-dimensional (2-D) DOA estimation; coherent sources; reduction dimension sparse reconstruction (RDSR); redundant sub-dictionary

2015-11-15;

2016-02-28;網絡優先出版日期:2016-06-07。

國家自然科學基金(61171180, 61371100)資助課題

TN 953.3

A

10.3969/j.issn.1001-506X.2016.08.01

王秀紅(1978-),女,講師,博士研究生,主要研究方向為陣列信號處理、波達角度估計。

E-mail:xiuhongwang@hit.edu.cn

毛興鵬(1972-),男,教授,博士研究生導師,主要研究方向為電子偵察與電子對抗、弱信號檢測、雷達信號處理。

E-mail:mxp@hit.edu.cn

張乃通(1934-),男,中國工程院院士,教授,主要研究方向為信號處理、專用移動通信系統、衛星通信、深空通信。

E-mail:ntzhang@hit.edu.cn

網絡優先出版地址:http://www.cnki.net/kcms/detail/11.2422.TN.20160607.1605.020.html

猜你喜歡
方法
中醫特有的急救方法
中老年保健(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
賺錢方法
捕魚
主站蜘蛛池模板: 伊大人香蕉久久网欧美| 欧美成人综合在线| 久久精品国产精品青草app| 国产精品无码在线看| 国产成人你懂的在线观看| 亚洲色大成网站www国产| 波多野结衣在线se| 国产成年无码AⅤ片在线| 91成人在线免费观看| 欧美国产在线看| 欧美另类精品一区二区三区| 一级毛片免费观看不卡视频| 久久这里只有精品2| 久久免费成人| 国产精品护士| 久久中文电影| 国产乱人伦偷精品视频AAA| 原味小视频在线www国产| 高h视频在线| 1024国产在线| 久久动漫精品| 激情综合网激情综合| 九色国产在线| 久久精品国产在热久久2019| 午夜三级在线| 一区二区三区毛片无码| 一级毛片在线播放免费观看| 婷婷五月在线| 在线免费亚洲无码视频| 欧美亚洲日韩不卡在线在线观看| 狠狠色丁婷婷综合久久| 亚洲免费福利视频| 精品欧美一区二区三区久久久| 日韩毛片免费观看| 九九线精品视频在线观看| 国产精品手机在线播放| 免费激情网址| 72种姿势欧美久久久久大黄蕉| 国产黄在线免费观看| 免费a在线观看播放| 婷婷午夜天| 99久久免费精品特色大片| 91精品视频播放| 日本免费福利视频| 九色视频线上播放| 亚洲欧美不卡视频| 亚洲天堂网在线视频| 亚洲国产精品日韩专区AV| 欧美精品在线看| 亚洲成网站| 欧洲一区二区三区无码| 欧美国产日韩另类| 中文字幕在线播放不卡| 欧美性猛交一区二区三区| 国产激情第一页| 亚洲成人在线网| 黄色网页在线观看| 极品国产在线| 国产色婷婷视频在线观看| 精品伊人久久久大香线蕉欧美| 亚洲青涩在线| 香蕉视频国产精品人| 欧美午夜在线观看| 91探花国产综合在线精品| 成人午夜久久| 国产成人在线无码免费视频| 欧美精品成人一区二区视频一| 天天干天天色综合网| 国产第八页| 国产尤物视频在线| 国产午夜精品一区二区三| 伊人无码视屏| 毛片最新网址| 色综合色国产热无码一| 真人免费一级毛片一区二区| 国产青榴视频在线观看网站| 黄色三级毛片网站| 在线免费a视频| 亚洲无码高清一区二区| AV无码国产在线看岛国岛| lhav亚洲精品| 97在线免费|