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

空間非合作目標慣性參數在軌辨識

2016-11-09 09:39:17孔祥龍李文龍李文峰趙毅
中國空間科學技術 2016年4期
關鍵詞:機械測量

孔祥龍,李文龍,李文峰,趙毅

上海衛星工程研究所,上海 201109

?

空間非合作目標慣性參數在軌辨識

孔祥龍,李文龍*,李文峰,趙毅

上海衛星工程研究所,上海 201109

在空間冗余機械臂抓捕非合作目標的慣性參數辨識問題中,已有方法大都基于系統動量已知的假設,且辨識過程未考慮基座姿態穩定。針對目標動量未知的問題,設計了具有增量形式的慣性參數分步辨識算法。首先基于線動量方程得到關于質量和質心位置的第一組估計方程,采用增量形式消除未知線動量更新估計方程。辨識結果收斂后根據估計參數計算線動量估計值,代入以轉動慣量為未知參數的第二組估計方程中,利用其增量表達式完成對轉動慣量的估計。辨識過程中的激勵由自適應零反作用控制輸入提供,算法在保證基座姿態不受干擾的同時還能對慣性參數精確辨識。仿真結果表明,在30 s以內算法已收斂,誤差收斂到零的同時,基座姿態角速率控制精度在10-3以下,說明算法收斂快,精度高,同時還能實現基座姿態穩定。

非合作目標;姿態穩定;冗余自由度;慣性參數辨識;反作用零空間;自適應控制

精確的慣性參數對于航天器控制系統和狀態估計系統來說都是非常重要的。以空間機械臂捕獲目標為例,當捕獲完成后,航天器控制系統必須盡快使組合體穩定下來以執行接下來的任務。如果無法獲取組合體航天器的慣性參數信息,控制系統無法精確控制組合體達到期望狀態,從而無法完成給定的任務。此外,由于服務航天器必須滿足天線對地定向等要求,因此在完成目標捕獲后保持服務航天器姿態穩定是非常必要的。

Aghili和Hillenbrand[1-2]采用了基于視覺測量的慣性參數辨識算法。文獻[3-4]基于牛頓歐拉方程,采用外部推力作為激勵,對慣性參數辨識問題進行了研究。然而,精確測量發動機推力大小和方向是非常困難的。為了解決上述問題,研究人員提出了基于動量守恒方程對慣性參數進行辨識的方法,只需驅動機械臂向任意方向運動,獲取不同的測量信息,且只需要對速度進行測量即可。這樣一方面可以提高準確性,另一方面還能節省燃料。

Psiaki[5]提出了一種采用總體最小二乘法進行參數辨識的方法,該方法需要對姿態和姿態角速率進行測量,且已經用于實際航天器中。Thienel等[6]提出對航天器慣性參數全部分量進行辨識的自適應方法。Sanyal等[7]在剛體姿態角速度的自適應跟蹤問題中討論了參數辨識方法。Murotsu等[8]考慮了機械臂抓取未知目標后對未知目標進行慣性參數辨識的問題,分別給出了基于動量守恒和牛頓歐拉方程的辨識方法。Yoshida等[9]給出了基于動量守恒的自由漂浮空間機器人的慣性參數辨識方法,考慮了重力梯度力矩,辨識過程無需測量力矩和加速度信息,并在ETS-VII中得到了應用。郭琦等[10]基于動量守恒定律給出了雙臂4自由度空間機器人捕捉未知目標的參數辨識問題。田富洋等[11-12]基于動量守恒分別研究了空間單臂和雙臂機器人抓取目標時航天器本體和目標參數辨識問題,同時分析了避免方程組奇異的問題。劉宇等[13]考慮空間機器人參數誤差問題,基于動量守恒,基于偏差模型的最小二乘法和遺傳算法對動力學參數進行了辨識。何光彩等[14]利用前向神經網絡辨識了空間機器人本體質量。王洪柳[15]基于最小二乘在線估計目標參數,從線性方程組最小二乘角度分析辨識過程中的參數辨識誤差,同時分析了測量噪聲對參數辨識的擾動影響。上述基于動量守恒方法的辨識大多假設系統動量為零,與實際情況不符,動量不為零時算法失效,此外辨識過程中未考慮基座姿態穩定要求。

Ma等[16]考慮了動量不為零時的辨識問題,采用了分步法,但只是利用機械臂對航天器本體進行辨識。金磊等[17]在此基礎上做了改進,采用了增量形式的辨識方程,提出采用三組測量信息的辨識方法,但未考慮激勵問題,且利用測量信息樣本少,辨識結果受樣本選取方式和測量精度的影響較大。Nguyen等[18-19]提出一種空間冗余機械臂自適應零反作用控制與慣性參數辨識結合的方法,通過在設計控制器中考慮空間飛行器慣性矩陣的變化,達到同時辨識空間飛行器慣性矩陣的目的,但也是基于系統線動量為零的假設的。

本文針對已有研究存在的不足,提出一種同時實現動量未知情況下非合作目標參數辨識和基座姿態穩定控制的方法,給出了算法詳細流程。辨識方法基于機械臂動量守恒原理,并以自適應零反作用控制輸入為持續激勵,采用分步辨識的方法。為了消除未知項的影響,采用辨識方程的增量形式作為新的參數估計方程。算法同樣適用于動量已知的情況,具有廣泛適用性。

1 機械臂參數辨識方程

圖1中SI表示慣性坐標系,S0表示基座航天器本體坐標系,Si表示第i節臂的連體坐標系,原點位于關節處,x軸沿桿向,由鉸鏈i指向鉸鏈i+1,z軸為關節轉動軸,y軸右手定則確定。r0為基座質心位置矢量,ai為關節i到桿i質心位置矢量,pi為關節i的位置矢量,ki為關節旋轉軸指向矢量,φi為關節角速率,ωt為目標轉動角速度矢量。如無特殊說明,文中的所有矢量都需在慣性系中給出。

圖1 自由漂浮空間機械臂Fig.1 Free-floating space manipulator

本節給出基于動量守恒的參數辨識方程。首先給出如下假設:1)機械臂抓捕目標后末端執行器與目標剛性連接,末節臂與目標成為一個整體,即物體n,因此末節臂慣性參數發生變化;2)組合體除末節臂以外,其余各體的慣性參數都是已知的;3)所有速度信息可通過敏感器測量得到。待辨識參數為新的末節臂的質量mn,質心位置nan,以及關于末節臂質心的轉動慣量nIn。上標n表示在坐標系Sn中,下標n表示物體n。

機械臂抓捕目標前,機械臂系統和目標分別具有線動量Pm和Pt以及角動量Lm和Lt,系統總的線動量P=Pm+Pt,總的角動量L=Lm+Lt,由動量守恒得,抓捕目標后系統的線動量和角動量仍為P和L。抓捕后系統動量可分別為兩部分,分別為基座與前n-1節臂組成的子系統1的動量以及末節臂子系統2的動量。前者由于可通過測量信息直接獲得,也稱為可測量部分;后者由于慣性參數是通過估計方法得到的,因此稱為估計部分。本文擬采用遞推最小二乘方法進行參數辨識。因此,首先要將動量方程寫成關于未知參數的線性方程形式。根據上述動量分解方法,系統動量方程可重新寫為

(1)

(2)

位置矢量

(3)

(4)

(5)

(6)

其中

由于本文重點驗證特殊激勵下的辨識算法的可行性以及對任意運動狀態系統的適用性,因此式(6)中暫不引入測量誤差的影響。

類似地,將式(3)和式(4)代入式(2)中得:

(7)

(8)

其中

將式(8)寫成矩陣形式:

(9)

其中

式(6)、式(9)共同構成了慣性參數辨識的基本方程。接下來將以這兩個方程為基礎給出具體的辨識算法。

2 參數辨識算法

2.1增量形式估計方程

在前文的公式推導中已得到形如式(6)和式(9)的辨識方程。從式(6)中的測量矩陣H1中可以看到存在與系統總的線動量P相關的項。同樣,在式(9)中,測量矩陣H2中同樣包含線動量P,此外在測量值Z2中還存在系統角動量L。機械臂抓捕目標后系統總動量是機械臂動量和目標動量的總和。在實際任務中,Pm和Lm始終可以看做是已知量,對于非合作目標則無法得到Pt和Lt的精確信息,也就無法知道P和L的精確值。以往的研究中大多假設抓捕前后系統的總動量都為0或者為其他已知常數,然后再做進一步處理,而實際上,兩者并不能一概而論統一處理。

通過觀察式(6)等式兩邊項可以發現,方程中等式左邊所有量都是可測量的,右邊項測量矩陣H1中僅P為未知項,因此取式(6)的增量形式以消除該常值項

(10)

其中

符號Δ(·)表示(·)在時間t到t+Δt內的增量。

可以發現,通過選取增量形式的估計方程可以將常值項P消去,由于機械臂始終處于運動狀態,因此測量值也是時變的,從而得到了時變的測量矩陣ΔH1和測量向量ΔZ1,構成了新的參數估計方程。

接下來對式(9)做同樣處理可得

(11)

其中

從ΔH2和ΔZ2的表達式中可以看出,常值項L已經消去,但是P仍然保留在方程中無法消去。考慮到線動量僅與系統各體質量、質心位置和速度相關,與轉動慣量無關,因此這里采用分步辨識的方法。

分步辨識方法的基本思路是先利用式(10)對目標質量和質心位置進行辨識,待參數收斂或者達到給定收斂條件后,依據所得辨識結果和速度測量信息計算系統線動量估計值,然后依據式(11)開始執行轉動慣量值的辨識。參數估計方法可以采用文獻[20]中給出的遞推最小二乘法。

2.2辨識激勵問題

文獻[21]指出,只有在機械臂的輸入運動充分并且持續地給予激勵時,參數辨識算法才能收斂到它的精確值。Nguyen-Huynh T C和Sharf I[18-19]針對空間冗余機械臂抓捕空間非合作自旋目標任務,討論了減小抓捕目標后機械臂運動對基座姿態的干擾的問題,提出了一種自適應零反作用控制算法。該算法基于動量守恒和參數自適應問題中的遞推最小二乘法,實現了基座的零反作用控制,且證明了激勵的充分和持續性。本文采用這種激勵方式,控制律設計過程可參照原文獻。需要指出的是,激勵的本質是關節處的控制力矩輸入,采用的是PD控制方法,屬于反饋控制,反饋控制的狀態量是關節角和角速度,是動力學方程在系統激勵下的輸出,二者關系如圖2所示。

圖2 待辨識參數與辨識激勵的關系Fig.2 Parameters to be identified and control input

2.3參數辨識算法流程

至此,整個基于自適應零反控制的慣性參數辨識算法可描述為以下流程:

1)初始化t=0,給定初值X1(0),P(0)=p1I,令t=t+Δt;

2)以自適應控制為輸入,仿真計算t時刻關節角速率φ(t)和ω0(t)的值,令t=t+Δt;

3)重復步驟2)得到下一時刻測量值φ(t)和ω0(t);

4)以式(6)為估計方程,依據遞推最小二乘算法更新增益矩陣P1(t);

5)遞推計算下一時刻估計值X1(t);

6)令t=t+Δt,判斷是否滿足收斂條件,如滿足,轉到步驟7),否則返回步驟3);

7)將mn和nan的估值代入式(1)可得線動量P;

8)初始化初值X2(0),P(0)=p2I,t=t+Δt;

9)重復步驟3);

10)以式(9)為估計方程更新增益矩陣P2(t);

11)遞推計算X2(t);

12)令t=t+Δt,判斷是否滿足收斂條件,如滿足,算法結束,否則返回步驟9)。

3 仿真算例

圖3所示為空間9自由度冗余機械臂抓捕非合作自旋目標示意?;教炱髻|量mb=500 kg,轉動慣量Ib=diag([83.61 83.61 83.61])kg·m2,尺寸為1×1×1 m3,各桿完全相同,mi=5 kg,i=1,…,9,轉動慣量為Ii=diag([0.1 1.5 1.5])kg·m2,桿長0.2 m,目標質量為mt=100 kg,Itxx=Ityy=Itzz=5 kg·m2,Itxy=Ityz=-0.1 kg·m2,Itxz=0.1 kg·m2,尺寸為0.5×0.5×0.5m3。初始時刻基座角速度ω0=[0 0 0]T,線速度v0=[0 0 0]T,各關節轉角為0,速度都為0.01 rad/s。初始構型為,臂1、3、5、7、9的連體坐標軸與基座本體系平行,臂2、4、6、8連體系互相平行,且S2為S1繞x軸順時針轉過90°所得。

圖3 空間9自由度機械臂Fig.3 9 DOFs space manipulator

末節臂抓捕目標與目標形成組合體后慣性參數的理論值為質量m=105 kg,質心位置rc=[0.433 3 0 0] m,轉動慣量Ixx=5.1 kg·m2,Iyy=7.083 3 kg·m2,Izz=7.083 3 kg·m2,Ixy=Iyz=-0.1 kg·m2,Ixz=0.1 kg·m2,采用遞推最小二乘估計方法,分別選取初值X10=0,X20=0,選取P10=p1I,P20=p2I,p1=p2=1×105,得到的辨識結果如圖4~圖7所示。

圖4給出的是質量辨識誤差,從圖中可以看出辨識算法收斂速度非??欤蠹s在4 s左右參數估計結果收斂,且辨識精度很高,不考慮測量誤差情況下質量辨識誤差收斂到零。圖5為基于動量方程辨識得到的質心位置誤差,同樣可以看出辨識算法收斂快,辨識精度高。

圖6是基于角動量方程辨識得到的轉動慣量誤差曲線,從圖中可以看出,轉動慣量達到收斂所需時間更長,但也在25 s左右達到收斂,最后辨識誤差都收斂到了零。由此證明本文所給的基于遞推最小二乘和基于自適應零反作用控制的慣性參數辨識十分高效且具有高精度,適合實時在線操作。

圖4 質量辨識誤差Fig.4 Estimation error of mass

圖5 質心位置辨識誤差Fig.5 Estimation error of position center of mass

圖6 轉動慣量辨識誤差Fig.6 Estimation error of inertia tensor

為了驗證自適應零反作用控制算法,圖7給出了基座角速度變化曲線,從圖中可以看出,基座姿態角速度量級小于10-3,也即近似為零,這也證明了控制算法的有效性。

圖7 基座角速度Fig.7 Attitude angular velocity of base

4 結束語

本文研究了空間冗余機械臂抓捕非合作目標后對未知目標的慣性參數做在軌辨識的算法。將分步辨識法與增量形式的估計方程結合,解決了系統動量未知時的傳統基于動量的辨識算法失效的問題。此外,采用自適應零反作用控制輸入作為激勵,可以保證激勵充分而持久。數值仿真結果證明了辨識算法的有效性和快速性,因此適合實時在線操作。仿真結果同時也證明了激勵輸入下基座姿態響應幾乎為零,說明即使原機械臂在抓捕非合作目標后受到了強烈的干擾,但基座姿態依然能夠保持穩定,符合實際工程需求。

本文給出的算法是基于測量信息可測且精確的前提下給出的,后續研究將結合實際工程應用,考慮測量可行性和測量精度,開展更進一步研究。

References)

[1]AGHILI F,PARSA K. Motion and parameter estimation of space objects using laser-vision data[J]. Journal of Guidance, Control, and Dynamics,2009,32(2):537-549.

[2]HILLENBRAND U,LAMPARIELLO R. Motion and parameter estimation of a free-floating space object from range data for motion prediction[R]. ESA,2005:461-470.

[3]LAMPARIELLO R,HIRZINGER G. Modeling and experimental design for the on-orbit inertial identification of free-flying space robots[C]∥International Design Engineering Technical Conference & Computers and Information in Engineering Conference. Long Beach,2005.

[4]BERGMANN E,DZIELSKI J. Spacecraft mass property identification with torque-generating control[J]. Journal of Guidance,Control,and Dynamics,1990,13(1):99-103.

[5]PSIAKI M L. Estimation of a spacecraft′s attitude dynamics parameters by using flight data[J]. Journal of Guidance,Control,and Dynamics,2005,28(4):594-603.

[6]THIENEL J K,LUQUETTE,SANNER R M. Estimation of spacecraft inertia parameters[C]∥AIAA Guidance, Navigation and Control Conference and Exhibit. Honolulu,2008.

[7]SANYAL A K,CHELLAPPA M J. Globally convergent adaptive tracking of spacecraft angular velocity with inertia identification and adaptive linearization[C]∥Proceedings of 42nd IEEE Conference on Decision and Control. Maul:IEEE,2003:2704-2709.

[8]MUROTSU Y,SENDA K,OZAKI M. Parameter identification of unknown object handled by free-flying space robot[J]. Journal of Guidance,Control,and Dynamics,1994,17(3):488-494.

[9]YOSIDA K, ABIKO S. Inertia parameter identification for a free-flying space robot[C] ∥AIAA Guidance, Navigation and Control Conference and Exhibit. Monterey, California, USA, August 5-8, 2002.

[10]郭琦, 洪炳镕. 雙臂四自由度空間機器人捕捉未知目標的參數辨識[J].機器人, 2005,27(6):512-516.

GUO Q,HONG B R. Parameter identification of unknown object handled by a dual-arm four-degree-of-freedom space robot[J]. Robot, 2005,27(6):512-516(in Chinese).

[11]田富洋,吳洪濤,趙大旭,等. 在軌服務雙臂空間機器人的參數辨識[J]. 華南理工大學學報,2010,38(2):73-78.

TIAN F Y,WU H T, ZHAO D X,et al. Parameter identification of on-orbit-servicing dual-arm space robot[J]. Journal of South China University of Technology,2010,38(2):73-78(in Chinese).

[12]田富洋,吳洪濤,趙大旭,等. 在軌空間機器人的參數辨識[J].中國空間科學技術, 2010, 30(1):10-17.

TIAN F Y,WU H T, ZHAO D X,et al. Parameter identification of orbital free-floating space robot[J]. Chinese Space Science and Technology,2010,30(1):10-17(in Chinese).

[13]劉宇,夏丹,李瑰賢,等. 基于角動量守恒的空間機器人動力學參數辨識[J]. 宇航學報, 2010, 31(3):695-700.

LIU Y,XIA D,LI G X,et al. Dynamic parameter identification for a space robot based on angular momentum conservation[J]. Journal of Astronautics,2010,31(3):695-700(in Chinese).

[14]何光彩,洪炳镕,郭恒業. 基于參數辨識的冗余自由飛行空間機器人多臂協調運動規劃[J].宇航學報, 2000, 21(1):85-89.

HE G C,HONG B R,GUO H Y. Multi-arm coordinated motion planning of redundant free-flying space robot based on parameter identification[J]. Journal of Astronautics,2000,21(1):85-89(in Chinese).

[15]王洪柳. 空間機器人抓取過程中的目標參數辨識[D]. 哈爾濱:哈爾濱工業大學, 2011.

WANG H L. Parameter identification of target in space robot grasping process[D]. Harbin:Harbin Institute of Technology,2011(in Chinese).

[16]MA O,DANG H,PHAM K. On-orbit identification of inertia properties of spacecraft using a robotic arm[J]. Journal of Guidance,Control,and Dynamics,2008,31(6):1761-1771.

[17]金磊, 徐世杰. 空間機器人抓取未知目標的質量特性參數辨識[J].宇航學報, 2012, 33(11): 1570-1576.

JIN L,XU S J. Inertial parameter identification of unknown object captured by a space robot[J]. Journal of Astronautics,2012,33(11):1570-1576(in Chinese).

[18]NGUYEN-HUYNH T C. Adaptive reactionless motion for space manipulator when capturing an unknown tumbling target[C]∥IEEE International Conference on Robotics and Automation. Shanghai,2011:4202-4207.

[19]NGUYEN-HUYNH T C. Adaptive reactionless motion with joint limit avoidance for robotic capture of unknown target in Space[C].IEEE/RSJ International Conference on Intelligent Robots and Systems. Vilamoura,2012:1155-1160.

[20]劉勝, 張紅梅. 最優估計理論[M]. 北京:科學出版社,2011:68-70.

[21]ARMSTRONG B. On find exciting trajectories for identification experiments involving systems with nonlinear dynamics[J]. International Journal of Robotics Research,1989,8(6):28-48.

(編輯:高珍)

Inertia parameter identification of on-orbit uncooperative object

KONG Xianglong,LI Wenlong*,LI Wenfeng,ZHAO Yi

Shanghai Institute of Satellite Engineering,Shanghai 201109,China

Most algorithms of inertia parameters identification assumed the system momentum known and few attention was paid to the attitude stabilization of base. An incremental form of the identification equation was given based on the conservation principle of momentum and a fractional step method was proposed. Firstly, estimation equations about the mass and position of mass center were obtained and the unknown constant linear momentum was eliminated in incremental formulation. The estimate of linear momentum was calculated with the estimated parameters and then substituted into estimation equations about the inertia tensors. The incremental equations were obtained again to eliminate the unknown constant angular momentum. An adaptive reactionless controller was used to persistently excite the robot arm, which guaranteed the convergence of the algorithm as well as the stabilization of base attitude. Simulations results show that inertial parameters are accurately identified and the angular velocities of the base are less than 10-3, which can be viewed as nearly undisturbed.

uncooperative target;attitude stabilization;redundant degree of freedom;inertia parameter identification;reaction null-space;adaptive control

10.16708/j.cnki.1000-758X.2016.0042

2016-01-07;

2016-03-01;錄用日期:2016-05-11;

時間:2016-07-1213:26:45

http:∥www.cnki.net/kcms/detail/11.1859.V.20160712.1326.004.html

國家“863計劃”(2015AA6301)

孔祥龍(1986-),男,碩士研究生,kxl.longxiao@163.com

李文龍(1988-),男,博士研究生,liwenlongzacao@126.com,主要研究方向為衛星總體設計,在軌服務與維護技術

V43

A

http:∥zgkj.cast.cn

引用格式:孔祥龍,李文龍,李文峰,等.空間非合作目標慣性參數在軌辨識[J].中國空間科學技術, 2016,36(4):17-23.

KONGXL,LIWL,LIWF,etal.Inertiaparameteridentificationofon-orbituncooperativeobject[J].ChineseSpaceScienceandTechnology, 2016,36(4):17-23(inChinese).

猜你喜歡
機械測量
機械革命Code01
電腦報(2020年35期)2020-09-17 13:25:53
調試機械臂
當代工人(2020年8期)2020-05-25 09:07:38
把握四個“三” 測量變簡單
ikbc R300機械鍵盤
電腦報(2019年40期)2019-09-10 07:22:44
滑動摩擦力的測量和計算
滑動摩擦力的測量與計算
簡單機械
測量的樂趣
機械班長
測量
主站蜘蛛池模板: 波多野结衣国产精品| 在线五月婷婷| 精品伊人久久大香线蕉网站| 手机在线免费不卡一区二| 免费一级α片在线观看| 玖玖精品在线| 亚洲男人在线| 欧美日本在线观看| 亚洲中文字幕在线一区播放| 亚洲欧洲日产无码AV| 一级毛片网| 9啪在线视频| 老熟妇喷水一区二区三区| 中文字幕有乳无码| 亚洲欧洲日韩综合| 久久亚洲天堂| 亚卅精品无码久久毛片乌克兰| 国产不卡网| 激情爆乳一区二区| 亚洲欧美成人在线视频| 五月激情综合网| 国产原创自拍不卡第一页| 欧美成人国产| 亚洲一区国色天香| 毛片网站免费在线观看| 亚洲人成网站在线播放2019| 最新国产麻豆aⅴ精品无| 亚洲国产天堂久久九九九| 久久人与动人物A级毛片| 精品国产免费第一区二区三区日韩| 婷婷亚洲天堂| 成人免费视频一区| 又爽又黄又无遮挡网站| 亚洲一区网站| 91久久国产热精品免费| 五月激情婷婷综合| 再看日本中文字幕在线观看| 波多野结衣的av一区二区三区| 国产欧美日韩18| 欧美精品导航| 自拍欧美亚洲| 毛片基地视频| 91九色视频网| 国产视频大全| 婷婷午夜天| 精品国产电影久久九九| 国产精品蜜芽在线观看| 国产9191精品免费观看| 国产成人精品高清不卡在线| 九九视频免费在线观看| 最新午夜男女福利片视频| 波多野结衣视频一区二区| 91久草视频| 欧美精品另类| 久久亚洲国产最新网站| 国产a网站| 日韩精品亚洲一区中文字幕| 亚洲品质国产精品无码| 好久久免费视频高清| 99久久国产综合精品女同| 成人国产精品一级毛片天堂| 欧洲极品无码一区二区三区| 亚洲精品天堂在线观看| 欧美区一区二区三| 国产午夜福利在线小视频| 成人久久18免费网站| 国产精品大白天新婚身材| A级毛片无码久久精品免费| 黄色网站在线观看无码| 日本a∨在线观看| 国产啪在线| 国产精品一线天| 久久精品丝袜| 免费aa毛片| 国产成人综合日韩精品无码首页 | 91精品啪在线观看国产| 亚洲成人高清在线观看| 国产精品国产主播在线观看| 亚洲av色吊丝无码| 日本久久久久久免费网络| 毛片网站在线播放| 国产精品成人AⅤ在线一二三四|