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

        ?

        基于遺傳-最大似然方法高速旋轉(zhuǎn)彈丸阻力系數(shù)辨識

        2016-12-14 01:25:15易文俊常思江史繼剛劉世平
        彈道學(xué)報(bào) 2016年4期
        關(guān)鍵詞:方法

        管 軍,易文俊,劉 海,常思江,史繼剛,劉世平

        (1.南京理工大學(xué) 瞬態(tài)物理國家重點(diǎn)實(shí)驗(yàn)室,南京 210094;2.重慶房地產(chǎn)職業(yè)學(xué)院 成本控制系,重慶 401331;3.南京理工大學(xué) 能源與動(dòng)力工程學(xué)院,南京 210094)

        ?

        基于遺傳-最大似然方法高速旋轉(zhuǎn)彈丸阻力系數(shù)辨識

        管 軍1,易文俊1,劉 海2,常思江3,史繼剛1,劉世平3

        (1.南京理工大學(xué) 瞬態(tài)物理國家重點(diǎn)實(shí)驗(yàn)室,南京 210094;2.重慶房地產(chǎn)職業(yè)學(xué)院 成本控制系,重慶 401331;3.南京理工大學(xué) 能源與動(dòng)力工程學(xué)院,南京 210094)

        為了進(jìn)一步提高高速旋轉(zhuǎn)彈丸氣動(dòng)參數(shù)辨識技術(shù)的精度,為射表編制、彈箭飛行控制技術(shù)等提供更加可靠的基礎(chǔ)數(shù)據(jù),對參數(shù)的可辨識性問題進(jìn)行了理論上的定量分析,利用遺傳-最大似然方法,通過彈丸空中自由飛行的速度數(shù)據(jù)對零升阻力系數(shù)進(jìn)行了辨識。通過仿真實(shí)驗(yàn)對算法的精確性和可靠性進(jìn)行了充分的驗(yàn)證。用該算法對實(shí)際飛行數(shù)據(jù)進(jìn)行了處理,得到了較為精確的彈丸零阻系數(shù),并應(yīng)用于工程實(shí)際問題,分析了仿真實(shí)驗(yàn)和實(shí)際實(shí)驗(yàn)的數(shù)據(jù)。結(jié)果表明,該算法具有較高的數(shù)據(jù)辨識精度并能有效地解決實(shí)際工程問題。

        旋轉(zhuǎn)彈丸;氣動(dòng)參數(shù)辨識;遺傳算法;最大似然估計(jì)

        空氣動(dòng)力是影響彈丸運(yùn)動(dòng)狀態(tài)的重要因素,氣動(dòng)力和力矩系數(shù)的精確與否又直接決定了空氣動(dòng)力的精確度,因此氣動(dòng)系數(shù)的辨識是彈箭飛行控制等領(lǐng)域重要的研究課題。目前,確定彈箭氣動(dòng)系數(shù)的方法主要有:理論計(jì)算法、風(fēng)洞實(shí)驗(yàn)法和射擊實(shí)驗(yàn)法。理論計(jì)算法因?yàn)椴荒軌蛲耆紤]彈箭飛行的實(shí)際情況,或多或少地進(jìn)行了一些近似或簡化,所以計(jì)算結(jié)果和實(shí)際結(jié)果存在偏差。風(fēng)洞實(shí)驗(yàn)法相比理論計(jì)算法精度有大幅的提高,但也存在一些問題,比如吹風(fēng)實(shí)驗(yàn)中難以模擬彈丸的各種運(yùn)動(dòng),這使動(dòng)態(tài)條件下的氣動(dòng)力測定較為困難;同樣由于模型支撐桿件的干擾和洞壁反射,測量結(jié)果也會存在一些誤差。射擊實(shí)驗(yàn)法是針對彈箭的真實(shí)自由飛行數(shù)據(jù)進(jìn)行相關(guān)的參數(shù)辨識,通過其辨識的結(jié)果可以對理論計(jì)算數(shù)據(jù)或者風(fēng)洞數(shù)據(jù)進(jìn)行適當(dāng)?shù)男拚蛘咧苯佑迷诠こ虒?shí)際問題中,進(jìn)一步提高了氣動(dòng)參數(shù)的精確度。對于射擊實(shí)驗(yàn)法,國外大部分是通過專門的室內(nèi)實(shí)驗(yàn)靶道來獲取彈丸自由飛行的相關(guān)數(shù)據(jù),并進(jìn)行氣動(dòng)參數(shù)的辨識,例如姿態(tài)、速度、位置等信息。據(jù)公開報(bào)道的資料顯示,目前國內(nèi)還沒有大口徑彈丸的室內(nèi)靶道。

        基于最大似然準(zhǔn)則的氣動(dòng)參數(shù)辨識技術(shù)是目前主要的參數(shù)辨識方法。汪清[1]等人利用最大似然估計(jì)對高速自旋飛行器進(jìn)行了相關(guān)參數(shù)的辨識;崔平遠(yuǎn)[2]等人用最大似然方法對有控飛行器進(jìn)行了參數(shù)辨識;程振軒[3]、祁載康等人利用雙加速度計(jì)的方法對彈體氣動(dòng)參數(shù)進(jìn)行了辨識。另外,國外的相關(guān)學(xué)者[4-7]分別根據(jù)飛行器的自由飛行數(shù)據(jù)對相關(guān)參數(shù)進(jìn)行了辨識;錢煒琪[8]、王雅平[9]、杜昌平[10]等人分別用遺傳算法對相關(guān)飛行器進(jìn)行了參數(shù)辨識。遺傳算法是求解復(fù)雜系統(tǒng)優(yōu)化問題的一種有效方法,具有較強(qiáng)的魯棒性和全局尋優(yōu)能力,將其與最大似然方法組合使用,能夠有效提高收斂速度,更快速地找到最優(yōu)解。

        本文通過將遺傳算法與最大似然方法相結(jié)合的方法來解決高速旋轉(zhuǎn)彈丸氣動(dòng)參數(shù)中零升阻力系數(shù)的辨識問題。

        1 遺傳-最大似然法的基本原理

        本文利用全彈道坐標(biāo)跟蹤雷達(dá)對彈丸的速度和位置信息進(jìn)行了測量,所用的彈道跟蹤雷達(dá)如圖1所示。

        最大似然方法是目前在參數(shù)辨識領(lǐng)域應(yīng)用最為廣泛的方法之一,本文將遺傳算法與最大似然方法相結(jié)合來獲取待辨識參數(shù)的最優(yōu)解[11]。

        圖2給出了遺傳-最大似然算法的基本原理。

        圖2 遺傳-最大似然方法原理圖

        下面將針對彈丸的自由飛行試驗(yàn),分別給出氣動(dòng)參數(shù)辨識所需要的數(shù)學(xué)模型、辨識準(zhǔn)則以及尋優(yōu)方法。

        1.1 高速旋轉(zhuǎn)彈丸運(yùn)動(dòng)模型

        為了更加準(zhǔn)確地描述彈丸在空中的運(yùn)動(dòng)狀態(tài),并考慮到高速旋轉(zhuǎn)彈丸的縱橫向運(yùn)動(dòng)相互耦合,采用6自由度彈道方程[12]作為參數(shù)辨識的基本模型,如式(1)~式(12)所示:

        mgsinθacosψ2]

        (1)

        (2)

        (3)

        (4)

        (5)

        ωηωζtanφ2

        (6)

        (7)

        (8)

        (9)

        (10)

        (11)

        (12)

        1.2 參數(shù)辨識準(zhǔn)則

        本文采用最大似然準(zhǔn)則作為目標(biāo)函數(shù),其形式可表示為

        (13)

        1.3 基于遺傳算法的尋優(yōu)方法

        由旋轉(zhuǎn)彈丸的運(yùn)動(dòng)方程可知,當(dāng)彈丸參數(shù)、發(fā)射諸元、氣象諸元以及氣動(dòng)參數(shù)等確定以后,彈丸的運(yùn)動(dòng)狀態(tài)是確定的。反之,若已知彈丸的運(yùn)動(dòng)狀態(tài),通過優(yōu)化的方法,總可以找到一個(gè)合適的阻力系數(shù),使得計(jì)算值以最佳效果逼近實(shí)測數(shù)據(jù)。

        1.3.1 目標(biāo)函數(shù)

        本文以最大似然準(zhǔn)則式(13)作為遺傳算法的目標(biāo)函數(shù),即通過遺傳算法搜索零升阻力系數(shù)cx0的最優(yōu)解,使得目標(biāo)函數(shù)最小。

        1.3.2 編碼及初始化種群

        由于二進(jìn)制編碼具有簡便、易操作等優(yōu)點(diǎn),所以二進(jìn)制編碼在遺傳算法中得到了廣泛的應(yīng)用。本文同樣適用二進(jìn)制編碼的策略對種群中的個(gè)體進(jìn)行編碼。根據(jù)文獻(xiàn)[13],種群的大小通常選為20~200,結(jié)合本課題的實(shí)際情況,選擇種群大小為50。

        1.3.3 編碼及個(gè)體適應(yīng)度函數(shù)

        通過前期的理論計(jì)算、風(fēng)洞吹風(fēng)試驗(yàn)等數(shù)據(jù),所用彈丸的零升阻力系數(shù)大概在0.1~0.5之間,文中選擇每個(gè)個(gè)體染色體的長度為10位,則其最小分辨率可表示為rmin=0.4/(210-1)。個(gè)體的解碼方式為:先將個(gè)體所對應(yīng)的二進(jìn)制數(shù)據(jù)d2轉(zhuǎn)換成十進(jìn)制數(shù)據(jù)d10,則該個(gè)體所對應(yīng)的零升阻力系數(shù)為cx0=d10rmin+0.1。為了達(dá)到尋優(yōu)的目標(biāo),適應(yīng)度函數(shù)一般是通過目標(biāo)函數(shù)的某種轉(zhuǎn)換而來,本課題中最優(yōu)的結(jié)果是使目標(biāo)函數(shù)往最小的方向發(fā)展,又因?yàn)楸菊n題中的目標(biāo)函數(shù)J有可能會出現(xiàn)負(fù)數(shù)的情況,綜合考慮以后選擇適應(yīng)度函數(shù)為f=1 000-J。每代的平均適應(yīng)度用fa表示。

        1.3.4 遺傳算子的確定

        1)選擇算子。

        在選擇算子的操作過程中,選用輪盤賭的策略對個(gè)體進(jìn)行選擇。

        2)交叉算子。

        交叉算子是以概率Pc進(jìn)行交叉運(yùn)算,文中采用相鄰個(gè)體之間單點(diǎn)交叉的方式進(jìn)行運(yùn)算。

        3)變異算子。

        變異是增加種群多樣性的搜索算子。變異概率Pm的選取對遺傳算法的性能具有較大的影響,概率過小有可能會陷入局部收斂,得不到全局最優(yōu)解,概率過大有可能會使迭代周期變長,一般取Pm=0.001~0.1,本文取Pm=0.01。

        4)終止條件。

        當(dāng)整個(gè)種群收斂,即各個(gè)體的適應(yīng)度相等,或者達(dá)到進(jìn)化代數(shù)的上限,停止運(yùn)算。本文以運(yùn)行到進(jìn)化代數(shù)的上限作為終止條件,進(jìn)化代數(shù)NG選為100。

        2 可辨識性問題分析

        在參數(shù)辨識的研究領(lǐng)域,參數(shù)的可辨識性問題是工程可行性分析的重要依據(jù)。若利用已有的測量數(shù)據(jù)很難辨識出或者根本不能夠辨識出待辨識參數(shù),那么該測量數(shù)據(jù)對該待辨識參數(shù)將不具有可辨識性。敏感系數(shù)可以衡量測量數(shù)據(jù)和待辨識參數(shù)之間相關(guān)性的大小,相關(guān)性越大,則敏感系數(shù)越大,反之越小,即通過判斷敏感系數(shù)的大小來判斷參數(shù)的可辨識性問題。在本文中,利用雷達(dá)測量彈丸的飛行速度-時(shí)間數(shù)據(jù)來辨識彈丸的零升阻力系數(shù),其敏感系數(shù)的推導(dǎo)過程如下。

        式(1)-式(12)中各微分方程中左側(cè)的物理量對cx0的敏感系數(shù)分別為

        以速度v對零升阻力系數(shù)cx0的敏感系數(shù)P11的推導(dǎo)來說明推導(dǎo)過程,由于:

        (14)

        所以只要推出了所有的敏感系數(shù),然后和六自由度彈道方程聯(lián)合求解,即可得到各敏感系數(shù)隨時(shí)間的數(shù)值解。式(14)中:

        cosδ2sin′δ1|cx0)]

        A17=-mg(cosθacosψ2P21-sinθasinψ2P31)

        圖3為敏感系數(shù)P11、P121隨時(shí)間變化的數(shù)值解。

        從圖3中可以看出,不同的測量數(shù)據(jù)對待辨識參數(shù)的敏感系數(shù)具有很大的差別。圖3(a)為速度v對cx0的敏感系數(shù)P11,圖3(b)為側(cè)偏z對cx0的敏感系數(shù)P121,由于側(cè)偏z對cx0的敏感系數(shù)P121的數(shù)值過小,即兩者之間的相關(guān)性過小,所以基本不能夠用側(cè)偏z來對cx0進(jìn)行辨識。反之,速度v對cx0的敏感系數(shù)P11的數(shù)值很大,v和cx0之間具有很強(qiáng)的相關(guān)性,即可以通過v對cx0進(jìn)行辨識,后續(xù)的實(shí)驗(yàn)也證明了該結(jié)論的正確性。

        圖3 敏感系數(shù)P11和P121

        3 實(shí)驗(yàn)驗(yàn)證

        3.1 仿真校驗(yàn)

        本節(jié)以某型彈丸為例進(jìn)行了算法的驗(yàn)證實(shí)驗(yàn),在已知?dú)鈩?dòng)力的情況下計(jì)算一組彈道數(shù)據(jù),利用這組彈道數(shù)據(jù)中的速度-時(shí)間數(shù)據(jù),根據(jù)所設(shè)計(jì)的算法檢驗(yàn)?zāi)芊竦玫綔?zhǔn)確的零升阻力系數(shù),算例中cx0的真實(shí)值為0.4。圖4為遺傳算法的搜索過程,顯示了適應(yīng)度函數(shù)隨代數(shù)的變化過程,隨著代數(shù)的不斷增加,各代平均適應(yīng)度也越來越大,表明個(gè)體在向好的方向發(fā)展。表1為不同射角β和不同誤差條件下的測試結(jié)果。表中,cx0為真實(shí)值;情況1:仿真值作為測量值;情況2:仿真值加上隨機(jī)白噪聲作為測量值,信噪比為100∶5。

        圖4 遺傳算法各代平均適應(yīng)度的迭代收斂過程

        β/(°)cx0情況1cx0,1Δcx0/%情況2cx0,2Δcx0/%50.40000.40621.540.3979-0.53250.40000.39990.010.4012-0.40450.40000.40050.140.3993-0.18550.40000.3989-0.270.3935-1.64650.40000.40140.360.40581.45

        上述測試結(jié)果表明,利用速度數(shù)據(jù),基于遺傳算法的最大似然估計(jì)算法可以有效地對零阻系數(shù)進(jìn)行準(zhǔn)確的辨識。

        3.2 實(shí)際數(shù)據(jù)處理

        在實(shí)際數(shù)據(jù)的處理過程中,零升阻力系數(shù)cx0是馬赫數(shù)的函數(shù),全彈道馬赫數(shù)在不斷變化,所以阻力系數(shù)也在不斷變化。根據(jù)文獻(xiàn)[14],利用小區(qū)間常數(shù)法,認(rèn)為在小區(qū)間內(nèi)零升阻力系數(shù)是不變的,所以將全彈道的測量數(shù)據(jù)分成若干個(gè)數(shù)據(jù)段,分別對每段進(jìn)行辨識,然后再插值得到全彈道的阻力系數(shù)。在本課題的處理中每段的數(shù)據(jù)點(diǎn)數(shù)選為50個(gè)。

        3.2.1 對一段數(shù)據(jù)的辨識

        表2 實(shí)際數(shù)據(jù)處理結(jié)果

        3.2.2 全彈道零升阻力系數(shù)的辨識

        全彈道的實(shí)際數(shù)據(jù)處理是分別對每個(gè)小區(qū)間的數(shù)據(jù)進(jìn)行處理,得到結(jié)果以后再通過插值的方法得到隨馬赫數(shù)變化的零升阻力數(shù)據(jù),如圖5所示。

        圖5 實(shí)際數(shù)據(jù)辨識的零阻隨馬赫數(shù)變化的曲線

        用該零阻系數(shù)代入彈道方程,計(jì)算得到的速度數(shù)據(jù)和實(shí)測雷達(dá)數(shù)據(jù)的對比如圖6所示,通過數(shù)據(jù)的分析,全彈道的擬合概率誤差在1.58m/s左右,有效地驗(yàn)證了算法和氣動(dòng)辨識結(jié)果的正確性。

        圖6 速度計(jì)算值和測量值的曲線

        4 結(jié)束語

        通過數(shù)值仿真和實(shí)際數(shù)據(jù)處理的結(jié)果可以看出,通過基于遺傳算法的最大似然方法能夠有效地辨識出高速旋轉(zhuǎn)彈丸的零升阻力系數(shù)。阻力系數(shù)的精確辨識為有控彈箭的自適應(yīng)控制提供了可靠的在線辨識技術(shù),為無控彈箭的射表編制等技術(shù)提供了必要的基礎(chǔ)數(shù)據(jù)。通過本文的研究,對推動(dòng)彈箭氣動(dòng)參數(shù)辨識技術(shù)的發(fā)展具有一定的意義。

        [1] 汪清,萬宗國,錢煒琪,等.高速自旋飛行器氣動(dòng)參數(shù)辨識[J].空氣動(dòng)力學(xué)學(xué)報(bào),2004,22(1):1-6. WANG Qing,WAN Zong-guo,QIAN Wei-qi,et al.Aerodynamic parameter identification of rapid self-rotating flight vehicle[J].Acta Aerodynamic Sinica,2004,22(1):1-6.(in Chinese)

        [2]崔平遠(yuǎn),楊滌.極大似然法及在有控飛行器動(dòng)參數(shù)辨識中的應(yīng)用[J].航空學(xué)報(bào),1991,12(11):A644-A649. CUI Ping-yuan,YANG Di.Maximum likelihood algorithm and its application to parameters identification[J].Acta Aeronautica Et Astronautica Sinica,1991,12(11):A644-A649.(in Chinese)

        [3]程振軒,祁載康,林福德,等.基于雙加速度計(jì)的彈體氣動(dòng)參數(shù)辨識[J].北京理工大學(xué)學(xué)報(bào),2009,29(4):283-285. CHEN Zhen-xuan,QI Zai-kang,LIN Fu-de,et al.Aerodynamic parameter identification based on dual-accelerometer[J].Transactions of Beijing Institute of Technology,2009,29(4):283-285.(in Chinese)

        [4]GIRISH C,RAVINDRA J.Aerodynamic parameter estimation from flight data applying extended and unscented Kalman filter[C]//2006 AIAA Atmospheric Flight Mechanics Conference and Exhibit.Keystone,Colorado:AIAA,2006:106-117.

        [5]ARDA A.Aerodynamic parameter estimation of a missile without wind angle measurements[C]//2014 AIAA Atmospheric Flight Mechanics Conference.Atlanda,GA:AIAA,2014:1-10.

        [6]LAWRENCE E H,MAYURESH P,CHRISTOPHER J R.Aerodynamic parameter identification and uncertainty quantification for small unmanned aircraft[C]//2015 AIAA Guidance,Navigation and Control Conference.Kissimmee,Florida:AIAA,2015:1 538-1 550.

        [7]BRADLEY T B.Aerodynamic parameter identification for symmetric projectiles:comparing gradient based and evolutionary algorithms[C]//2014 AIAA Atmospheric Flight Mechanics Conference.Portland,Oregon:AIAA,2011:235-248.

        [8]錢煒祺,汪清,何開鋒,等.混合遺傳算法在氣動(dòng)力參數(shù)辨識中的應(yīng)用[J].飛行力學(xué),2004,22(1):33-36. QIAN Wei-qi,WANG Qing,HE Kai-feng,et al.Application of hybrid genetic algorithms for aerodynamic parameter estimation[J].Flight Dynamics,2004,22(1):33-36.(in Chinese)

        [9]王雅平,齊曉慧,張昊,等.基于改進(jìn)遺傳算法的動(dòng)力傘氣動(dòng)參數(shù)辨識研究[J].計(jì)算機(jī)仿真,2015,32(6):31-34. WANG Ya-ping,QI Xiao-hui,ZHANG Hao,et al.Identification for aerodynamic coefficients of powered paraglider[J].Computer Simulation,2015,32(6):31-34.(in Chinese)

        [10]杜昌平,周德云,宋筆鋒,等,基于遺傳算法的彈道參數(shù)辨識方法研究[J].西北工業(yè)大學(xué)學(xué)報(bào),2008,26(3):373-376. DU Chang-ping,ZHOU De-yun,SUN Bi-feng,et al.Weapon parameter identification algorithm based on genetic algorithm[J].Journal of North Western Poly Technical University,2008,26(3):373-376.(in Chinese)

        [11]張?zhí)戽?錢煒琪,何開峰,等.基于最大似然法的風(fēng)洞自由飛試驗(yàn)氣動(dòng)力參數(shù)辨識技術(shù)研究[J].實(shí)驗(yàn)流體力學(xué),2015,29(5):8-14. ZHANG Tian-jiao,QIAN Wei-qi,HE Kai-feng,et al.Research on aerodynamic parameter identification technology in wind tunnel free-flight test based on maximum likelihood estimation[J].Journal of Experiments in Fluid Mechanics,2015,29(5):8-14.(in Chinese)

        [12]韓子鵬.彈箭外彈道學(xué)[M].北京:北京理工大學(xué)出版社,2014. HAN Zi-peng.Exterior ballistics of projectile and rocket[M].Beijing:Beijing Institute of Technology Press,2014.(in Chinese)

        [13]王小平.遺傳算法理論、應(yīng)用與軟件實(shí)現(xiàn)[M].西安:西安交通大學(xué)出版社,2002. WANG Xiao-ping.Theory,application and software realization of genetic algorithm[M].Xi’an:Xi’an Jiaotong University Press,2002.(in Chinese)

        [14]閆章更,祁載康.射表技術(shù)[M].北京:國防工業(yè)出版社,2000. YAN Zhang-geng,QI Zai-kang.Firing table technology[M].Beijing:National Defense Industry Press,2000.(in Chinese)

        Drag Coefficient Identification of Spinning Projectile Using GA-MSE

        GUAN Jun1,YI Wen-jun1,LIU Hai2,CHANG Si-jiang3,SHI Ji-gang1,LIU Shi-ping3

        (1.National Key Laboratory of Transient Physics,Nanjing University of Science and Technology,Nanjing 210094,China;2.Department of Cost Control,Chongqing Real Estate College,Chongqing 401331,China;3.School of Power and Engineering,Nanjing University of Science and Technology,Nanjing 210094,China)

        In order to improve the accuracy of high-speed spinning projectile and to supply more precise base data for compiling firing tables and designing control systems for guided projectile,the analysis of identifiability for aerodynamic parameters in theory was carried out.Based on Genetic Algorithm-Maximum Likelihood Estimation(GA-MSE),the drag coefficients of spinning projectile were identified by using radar data of velocity.The accuracy and reliability of GA-MSE were validated by using simulated data.The actual data was processed by using GA-MSE,and a satisfied result was obtained.The result of simulation experiment data and actual experiment data were analyzed.The result shows that GA-MSE has high precision for parameter identification.

        spinning projectile;aerodynamic-parameter identification;genetic algorithm;maximum likelihood estimation

        2016-07-06

        國家自然科學(xué)基金項(xiàng)目(11472136,11402117)

        管軍(1987- ),男,博士研究生,研究方向?yàn)閺椉w行控制。E-mail:guanjun8710@163.com。

        易文俊(1970- ),男,教授,博士生導(dǎo)師,研究方向?yàn)閺椉w行控制。E-mail:yiwenjun0@163.com。

        V212

        A

        1004-499X(2016)04-0001-06

        猜你喜歡
        方法
        中醫(yī)特有的急救方法
        中老年保健(2021年9期)2021-08-24 03:52:04
        高中數(shù)學(xué)教學(xué)改革的方法
        化學(xué)反應(yīng)多變幻 “虛擬”方法幫大忙
        變快的方法
        兒童繪本(2020年5期)2020-04-07 17:46:30
        學(xué)習(xí)方法
        可能是方法不對
        用對方法才能瘦
        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
        賺錢方法
        91人妻人人做人人爽九色| 久久精品亚洲乱码伦伦中文| 图图国产亚洲综合网站| 日本一区不卡高清在线观看| 亚洲国产系列一区二区| 人妻av无码一区二区三区| 娇妻玩4p被三个男人伺候电影| 狼色在线精品影视免费播放| 精品人妻在线一区二区三区在线| 在线观看av片永久免费| 国产suv精品一区二人妻| 男女好痛好深好爽视频一区| 在线观看av不卡 一区二区三区| 狠狠躁天天躁无码中文字幕图| 最新高清无码专区| 五月天综合社区| 国产大屁股熟女流白浆一区二区| 国产一区二区精品久久岳| 婷婷丁香五月中文字幕| 国产欧美激情一区二区三区| av影片手机在线观看免费网址| 日本大骚b视频在线| 日本一区二区不卡视频 | 欧美aaaaaa级午夜福利视频| 九色91精品国产网站| 色偷偷亚洲精品一区二区| 日本真人边吃奶边做爽电影| 国产嫖妓一区二区三区无码| 亚洲欧美香港在线观看三级片| av在线高清观看亚洲| 国产精品成人aaaaa网站 | 亚洲综合久久精品无码色欲| 国产91 对白在线播放九色| 亚洲国产精品av麻豆一区| 三个男吃我奶头一边一个视频| 18禁高潮出水呻吟娇喘蜜芽| 无码啪啪熟妇人妻区| 久久av粉嫩一区二区| 亚瑟国产精品久久| 亚洲VA中文字幕无码毛片春药 | 婷婷丁香91|