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

空間失穩目標線陣成像畸變分析與三維重建

2020-03-13 08:32:20李榮華王振宇肖余之薛豪鵬
宇航學報 2020年2期
關鍵詞:方向測量

李榮華,王振宇,陳 鳳,肖余之,薛豪鵬

(1. 大連交通大學機械工程學院,大連 116028;2. 上海宇航系統工程研究所,上海 201109)

0 引 言

太空垃圾又稱空間碎片或軌道碎片,是指圍繞地球軌道的無用人造物體。太空垃圾小到由人造衛星碎片、漆片、粉塵,大到整個飛船殘骸構成。其中大型太空垃圾(衛星由于故障或任務結束而被放棄)不但占用其原在的有限軌道資源,還對其他航天器的安全造成不可估計的危害。為解決這一問題,各國陸續開展了對失效衛星的環繞檢測、運動跟蹤、交會對接、在軌維修等操作[1-4]。而實現這些目的,就必需獲取非合作目標的三維空間信息,進行三維模型重建。激光雷達由于其受光照和距離等因素影響小,抗干擾能力強等優點,是目前監測空間目標較為理想的工具,常見的成像方式為線陣式與面陣式[5]。

面陣式激光雷達[6]采用多束密集的激光束直接向各個方向漫射,最大優勢是能夠快速記錄整個場景,避免了掃描過程中目標或激光雷達移動帶來的各種問題。但對遠距離空間運動,面陣式成像分辨率低,對目標的測量精度仍然不能滿足測量要求。線陣式激光雷達在空間作業中,相較于面陣式負載較小,成像分辨率高,獲取的線陣數據同樣能夠用于運動參數估計及三維模型重建。對于靜態目標,已知目標與雷達相對運動信息,運用線陣式激光雷達獲取不同時刻線陣數據,可直接組合線陣信息得到測量目標的三維模型;但對于運動目標,目標位姿一直處于變化狀態,不同時刻獲取的是當前時刻目標掃描位置線陣數據,掃描數據存在畸變,所測量的運動目標輪廓與目標的實際外形相比一般會存在形變,直接組合線陣信息并不能準確地對被測目標進行三維重建。對于非平衡(失穩)狀態下,蔡銀橋等[7]對星載三維激光雷達圖像畸變進行研究,研究了在不同掃描周期內三維激光雷達圖像變形,得出了控制成像時間是降低圖像畸變的關鍵因素。針對線陣畸變對目標方位估計性能的影響,李建等[8]通過理論、仿真以及試驗數據處理等方面進行分析,并提出了一種初步的陣形校正方法。

解決運動目標線陣掃描中存在的畸變問題,主要分為合作目標與非合作目標兩種情況下進行分析[9-10]。對于合作目標,即目標運動參數已知或目標可添加標記,趙竹新等[11-12]給出了一種利用線陣相機立體交會測量彈丸設計,根據彈丸目標的外形特點設計了一種新的標記線方法,避免了線陣相機掃描速度不足而難以實現彈道同步測量的難題。高俊釵等[13]根據線陣相機圖像采集特點,以標準圓形作為標定目標確定了物像關系和校正插值系數,并設計線性插值算法進而能夠對任意形狀的目標進行成像校正。在合作模式下,被測目標運動估計精度通常能夠達到要求,但是對于無法進行標定的測量目標或者無法增加輔助設備的測量情況并不適用。對于非合作目標,王盈等[14]提出了一種適用于復雜空間目標的在軌激光雷達成像仿真流程,并強調了面積法在復雜目標成像可見判斷方面的適用性,適用于大表面目標測量。本課題組前期通過相鄰兩次掃描的畸變點云配準來獲取測量目標的瞬時運動參數,該方法僅使用單一載荷實現了非合作空間目標相對位姿的測量,但是僅能估算自旋狀態的空間失穩目標,對于復雜運動的空間失穩目標或者兩次獲取點云重合量達不到配準要求的情況有待進一步研究。

綜上所述,本研究以實現空間失穩目標三維重建為目的,針對空間失穩目標,因其非合作性、形狀多樣性、運動狀態不確定性,且線陣掃描式激光雷達測量動態目標發生畸變而無法精確恢復目標三維模型問題,首先對掃描畸變從測量角精度與測距精度兩大判定條件進行分析,然后根據測量工況對目標點云進行帶畸變配準。該方法不僅能夠有效解決空間失穩目標線陣成像畸變問題,而且可以實現較高精度的動態非合作目標三維重建,為線陣式測量技術進一步發展提供新的研究方向。

1 線陣激光雷達畸變分析

1.1 畸變產生原因

線陣掃描型激光雷達在探測失穩目標過程中,會發生時空非定常點云數據畸變,所掃描的運動目標數據與目標的實際模型相比會存在差異。這種差異是由于線陣雷達的掃描等效速度與目標運動速度不一致所導致的。所謂的線陣雷達的掃描等效速度是指線陣雷達單個像素的寬度與掃描頻率的乘積,比如傳感器的單個像素寬度為0.001 mm,掃描頻率為100k 線/s,則掃描等效速度為1 m/s。當雷達的掃描等效速度與目標的運動速度相等時,目標的線陣成像與目標的實際模型一致;當雷達的掃描等效速度大于目標的實際運動速度時,目標的線陣成像與其實際目標模型相比顯得被拉伸了,相反,當雷達的掃描等效速度小于目標的運動速度時,目標的線陣成像與其實際目標模型相比顯得被壓縮了。運動目標線陣成像的這種特點,從成像效果看,是一種成像失真。如圖1所示,假設目標總共需要線陣雷達單程掃描8次,奇數次雷達掃描方向與目標運動方向一致,雷達的掃描等效速度小于目標的實際運動速度,目標表面同向壓縮;偶數次雷達掃描方向與目標運動方向相反,雷達的掃描等效速度大于目標的實際運動速度,目標表面反向拉伸。最終導致線陣雷達掃描運動目標后產生不容忽視的畸變誤差。

圖1 目標表面掃描畸變產生示意圖

1.2 畸變容忍度臨界值

因為線陣掃描雷達受掃描棱鏡加工精度影響,角分辨率大小會出現周期性變化,所以即使目標是靜止狀態,線陣雷達測量的點云也會出現畸變。由此可見,當目標受運動影響出現的畸變大小如果不大于雷達自身測量精度,則恢復畸變的過程將失去意義;反之,如果運動畸變大于雷達測量畸變,則有理由認定為目標失穩,需要畸變補償。此外,由于受雷達角分辨率精度影響的畸變程度無法通過疊加來判斷,而且考慮到為了提高畸變恢復效果及重建精度,以雷達測角精度為參考的失穩臨界值向下降低一個數量級。

以雷達測量精度為衡量基準計算目標畸變容忍度臨界值,即當目標由于運動導致的數據畸變程度超過了雷達自身測量精度的10%(降低一個數量級)時的目標運動速率認為是畸變容忍度臨界速率。現從雷達測角精度與測距精度兩方面對畸變容忍度臨界值進行判定。

1.2.1角精度判定條件

當畸變容忍角度為r時,根據雷達指標:雷達環繞測量分辨率為m×n=400×400,線掃描寬度為Ls=50,單方向掃描次數為m/Ls=8,單次掃描時間為t=100 ms=0.1 s。可以推算出相鄰點掃描所花時間t′=3.125×10-5s:

(1)

重建所能容忍的相鄰點間畸變角度臨界值r為激光雷達測角精度(30″),即r=30″,則目標失穩臨界自旋速率ωr=26.6667(°)/s:

(2)

目標畸變容忍度臨界自旋速率為ω=26.6667(°)/s;即當目標自旋速率大于該值時,即使忽略雷達測角精度,章動的復合運動的角重建誤差也會超過指標。

1.2.2測距精度判定條件

假定已知信息:目標自旋角速率為ωm=12(°)/s,章動角速率為ωz=2.5(°)/s,雷達視場角為20°×20°,目標與雷達工作距離為L=80 m,成像分辨率為m×n=400×400,目標主體參數為2.000 m×3.500 m×2.500 m,雷達測距精度為dm=0.05 m。目標章動方向和自旋運動方向在某一時刻同向時,目標表面一點達到最大的運動角速率ωmax=14.5(°)/s:

ωmax=ωm+ωz

(3)

考慮最大畸變情況,對下面兩種情況分別計算。

1)目標運動方向與掃描方向平行時

目標運動方向與掃描方向平行時,雷達掃描運動示意圖如圖2所示,根據雷達掃描距離L=80 m,雷達視場角為20°×20°,可得X軸單方向雷達最大掃描距離為LX=28.212 m(Y軸方向同理):

圖2 平行時雷達掃描運動示意圖

(4)

X軸單方向掃描次數為m/Ls=8,則單次掃描距離L′X=3.527 m:

L′X=LX/(m/Ls)

(5)

目標主體參數中寬度3.500 m與高度2.500 m均小于3.527 m,因此雷達掃面目標主體掃描一次或兩次。當雷達掃面目標主體掃描兩次時,用時最長,且畸變最大。最大用時tmax=0.014 s:

tmax=(1+3.500/28.212)t/8

(6)

圖3 目標運動半徑

如圖3所示,目標運動半徑r為:

(7)

圖4 目標表面點畸變量示意圖

圖4中,目標表面P點為真實位置,P′為畸變位置。畸變量Ldis=0.0071 m:

Ldis=2rsin(tmaxωmax/2)

(8)

計算顯示:在畸變最大的情況下,表面一點P的畸變量為0.0071 m,而雷達測距精度小于等于0.05 m,因此畸變可忽略。

2)當目標運動方向與掃描方向垂直時

當目標運動方向與掃描方向垂直時雷達掃描運動示意圖如圖5所示。

圖5 垂直時雷達掃描運動示意圖

同理當目標運動方向與掃描方向平行時,可得目標運動方向與掃描方向垂直狀態畸變量Ldis=0.0071 m:

Ldis=2rsin(tmaxωmax/2)

(9)

計算顯示:在畸變最大的情況下,表面一點P的畸變量為0.0071 m,而雷達測距精度小于等于0.05 m,因此畸變可忽略。

3)根據測距精度計算最大運動速度容忍度ωr

根據測距精度0.05 m,畸變影響在10%獲得失穩速率臨界:

ωr=2arcsin(0.05×10%/2/r)/tmax=10.15(°)/s

(10)

若章動為2.5(°)/s,自旋速率最大容忍度為7.65(°)/s。

1.3 臨界值判定

綜合兩種臨界值判定條件:根據角精度判定條件,目標畸變容忍度臨界自旋速率ω=26.6667(°)/s;根據測距精度判定條件,當目標運動方向與掃描方向平行與垂直兩種狀態下,畸變量均低于測距精度一個數量級,且根據測距精度,目標畸變容忍度臨界自旋速率ωr=10.15(°)/s。在本測量系統已知數據下,結合費馬引理,目標畸變容忍度臨界自旋速率取極大值,即ωr=26.6667(°)/s。當目標自旋速率大于等于ωr時,初次配準時需要對點云數據進行畸變恢復;當目標自旋速率小于ωr時,初次配準時對點云的畸變恢復將失去意義。

2 帶畸變點云配準

2.1 坐標系建立

目標坐標系坐標原點設為待測目標幾何中心位置,ZO軸方向設為待測目標自旋軸方向,YO軸方向設為太陽翼中軸線方向,基于右手法則確定XO軸方向。雷達坐標系以雷達中心為坐標原點,環繞檢測半徑設為80 m,三軸方向與初始獲取目標反饋信息的方向一致。世界坐標系與雷達坐標系重合,掃描工況笛卡爾坐標系如圖6所示。

圖6 掃描工況坐標系

2.2 點云配準實質

點云配準是通過匹配具有重疊部分的數據集,將不同視場下得到的三維數據變換到同一坐標系下,得到信息更加豐富完整的目標三維點云數據,計算出點云間的歐拉變換關系,即旋轉變換矩陣R和平移向量T。

現采用線陣掃描型激光雷達進行N次掃描,得到N個視場點云數據,假設兩組相鄰點云數據Ni和Ni+1是對同一目標在不同成像視角得到的,配準目的便是通過平移和旋轉變換將Ni+1變換到Ni坐標系下,定義為N′i+1:

N′i+1=R(α,β,γ)×Ni+1+T(tX,tY,tZ)

(11)

其中,R(α,β,γ)為旋轉矩陣;T(tX,tY,tZ)為平移矢量;α,β,γ分別為繞X軸、Y軸、Z軸逆時針旋轉的歐拉角;tX,tY,tZ為三個坐標軸正方向的平移量。旋轉矩陣R(α,β,γ)可以表示為:

R(α,β,γ)=

(12)

通過選取合適的變換角度和平移量,將Ni+1變換到Ni坐標系下,進行坐標系的統一處理,最后得到的N′i+1與Ni即為所需要的配準結果。

2.3 運動特性初值

1)統計配準參數與對應時間數組

通過上一步配準,獲取目標相鄰掃描時間內運動的旋轉矩陣與平移向量,以初始時刻目標所在位置為位姿起始點,對應位姿增量疊加可以獲得目標位姿-時間函數。通過配準獲得參數如式(13)所示:

(13)

2)計算第i次旋轉平移運動增量

旋轉平移矩陣形式為(旋轉次序依次為:Z-Y-X):

(14)

(15)

旋轉運動增量:

(16)

平移運動增量:

(17)

式中:sX與cX代表繞X軸旋轉的正弦與余弦值;Y軸、Z軸同理。

2.4 迭代配準

實現點云精配準現今主流的核心算法是ICP算法,該算法精度很大程度上取決于預先給定的位置姿態初值,即需要先進行點云粗配準,航天器一般具有較為明顯的平面幾何特征,現采用RANSAC算法提取目標點云特征點進行粗配準[15]。

1)點云粗配準

設由線陣雷達獲取的點云數據集為P,從中隨機選取三個點:P1(x1,y1,z1);P2(x2,y2,z2);P3(x3,y3,z3)。然后根據三點位置確定平面S:z=Ax+By+C,參數A,B,C由式(18)確定:

(18)

統計平面S上的點的個數,設定平面厚度為2σ,計算P中任何一點Pi到平面S的距離di,di由下式計算得到:

(19)

統計di<σ的點的個數,記為S的得分數。重復上述步驟K次,選擇出得分最高的平面S*,其中K的計算公式為:

(20)

式(20)中m為點云中總的點數,n為點云中特征點點數,由于m和n值都很大,采用近似計算,由下式計算K值:

(21)

式中:ε為位于平面S*之外點所占比例的值(估計值),Ф為經過K次采樣之后最后被選中的概率。

如圖7所示,運用RANSAC算法提取相鄰點云平面幾何特征,圖7(a)為相鄰點云數據顯示,圖7(b)為兩組點云相似特征點連線,據此進行點云粗配準。

圖7 RANSAC算法

2)點云精配準

精細配準是一種更深層次的配準方法。經過前一步粗配準,得到了變換估計值。將此值作為初始值,在經過不斷收斂與迭代的精細配準后,達到更加精準的效果。ICP算法(最鄰近點迭代算法)最早由文獻[16-17]提出,文獻[16]采用的是待配準表面上點的配準方法,文獻[17]采用曲面點集與切平面的配準方式。本文處理的點云數據形式為點集,而后者所采用的ICP算法對每個點的法向量進行計算,計算量較大,因此本文主要采用前者的待配準表面上點的配準方法。

算法基本原理:兩點云P與Q,對于點云P中的每個點pi,搜索其在點云Q上的最近點作為與之相對應的點qi,然后依據對應關系求解使得式(22)中目標函數達到最小的剛體變換,即旋轉矩陣R和平移向量T,并將該變換作用于源點云,迭代進行這一過程直到滿足某一設定的收斂準則。

(22)

迭代終止條件:當目標相對于真值位姿精度不大于1°(3σ,三軸)時迭代終止。

3 三維重建

3.1 重建視場站位分布

假定失穩衛星最大自旋速度約為12(°)/s,定義旋轉軸過模型質心,方向沿Z軸方向,并且目標在短時間內為勻速運動。現對Y軸方向采用最大自旋速度進行高壓仿真,即三個方向速率分別為:vX=0(°)/s,vY=12(°)/s,vZ=0(°)/s。環繞目標沿Y軸按照每12度1幀提取點云可視化數據信息,如圖8所示。

圖8 多視場點云數據獲取

3.2 畸變恢復

以雷達第一次掃描的初始時刻,目標所在空間位置建立目標坐標系,如果想把每一條畸變線信息恢復到目標形貌真實位置(x0,y0,z0),需要計算該條畸變點云線相對初始時刻所在坐標系的平移量和旋轉分量,經過前面的步驟已經獲取了目標天體的運動參數模型,可以通過對不同時刻采集的線信息進行逆運算以把多組線信息整合到同一時刻的空間位置,ti代表統計雷達在獲取每一條線信息相對初始時刻所經過的掃描時間。則每一條雷達掃描線點云的逆向數據恢復模型為:

(23)

式中:

(24)

(25)

通過對多組相鄰畸變點云數據的配準,可以獲取沿時間延長的多組有限個相鄰掃描時間內的相對位姿增量,如圖9所示。

圖9 目標瞬時運動狀態

3.3 自修正式重建

本研究所采用的目標模型由衛星主體及太陽能帆板構成。在激光雷達環繞檢測過程中,兩側的太陽能帆板會因主體及自身遮擋而導致點云關鍵部位特征突變,從而導致相鄰點云配準時產生極大躍性偏差,對重建結果造成極大干擾。運用自修正原理,通過去除運動參數中存在的跳躍點,對環繞檢測一周數據進行擬合,恢復躍性點處數據,并據此進行三維重建。

失穩目標在雷達相鄰的兩次掃描時間內能夠保證有效重合掃描區域,基于此,對目標相鄰掃描點云進行拼接;從而在目標旋轉周期內,重復掃描、配準過程完成目標表面點云的三維重建。過程如下:1)獲取第j次掃描點云,并作為重建起始,經上一步實現點云畸變信息恢復;2)獲取第j+1次掃描點云,同樣恢復畸變信息,與第i次點云進行配準,舍去重合區域;3)繼續獲取下一時刻掃描點云,重復上述過程;4)以當前時刻掃描點云與第j次點云配準重合度作為終止重建條件:當兩處點云存在重合區域或者重合度大于設定閾值,認為目標完成一次周期自轉,即完成重建工作。依據配準擬合結果,將各旋轉角度對應結果還原到目標初始坐標系下,點云數據依次疊加,得到目標初始三維重建模型。

4 仿真校驗

為了證明算法的有效性,進而驗證本文提出的帶畸變自修正式三維重建方法的可行性,現應用VS2013結合點云庫(Point cloud library,PCL)對上述思想逐一進行仿真試驗與校驗。(文中算法結果均在CPU為Intel(R)Core(TM)i5-3210M CPU,2.50 GHz,內存為4 GB的PC機上運行所得,操作系統為Windows 7)

4.1 獲取目標衛星三維點云模型

現對目標進行仿真重建,如圖10所示,外形尺寸為16000 mm×6000 mm×3000 mm,本體尺寸為3500 mm×2000 mm×2500 mm。

圖10 目標點云模型

4.2 配準

按照每12度1幀,逆時針環繞目標沿Y軸提取點云數據,可獲得N個視場點云數據:

(26)

其中第1組0°為基準初始位置,其后每隔12°取樣一次,最后一組為環繞檢測360°后得到的可視化點云數據。32組目標環繞點云數據組如圖11所示。

圖11 32組目標環繞點云

按照上面提出的帶畸變自修正式三維重建方法,針對獲取的32組點云數據進行校驗。

1)第i組與第i+1組相鄰視角點云數據ICP配準。

2)獲取該組旋轉矩陣Ri與平移向量Ti。

3)將第i+1組點云根據旋轉矩陣Ri與平移向量Ti,配準恢復到第i組點云狀態下,相鄰組點云配準如圖12所示。

圖12 相鄰點云配準

4)將32組點云數據相鄰兩兩進行配準,獲取目標運動增量、解算目標瞬時運動參數。

5)然后對目標運動參數進行自修正式處理,去除特征突變干擾點。

6)依次恢復點云到初始幾何位置,完成初始目標三維重建。

7)通過系統進一步迭代,位姿增量自收斂,直至最終滿足重建精度要求。

4.3 配準結果分析

將32組點云數據,按照相鄰組旋轉配準,按照每12度1幀取樣,每組旋轉角度真值為12°,結果如表1所示。

表1 相鄰組旋轉配準結果

如圖13所示,實線為32組點云數據第一次配準結果,虛線為真值,點劃線為擬合結果。其中當旋轉角度為72°及252°時,模型單側帆板在視角可視范圍內消失,使得ICP配準時產生明顯偏差。根據第2節闡述的方法,剔除兩組較大偏差點,將剩余30組配準數據進行擬合,得到擬合結果。

圖13 相鄰點云第1次配準

依據配準擬合結果,將各旋轉角度對應結果還原到目標初始坐標系下,點云數據依次疊加,得到目標初始三維重建模型,第1次配準耗時為18337 ms。

4.4 迭代結束

本文獲取的目標模型,在進行第6次迭代后,單次耗時為6779 ms,目標單位時間內運動增量均位于12°±1°內;在現有硬件配置環境下,配準總耗時為72467 ms,符合三維重建技術指標。第6次配準結果如表2所示,擬合結果及目標最終三維重建模型如圖14~圖15所示。

5 結 論

本研究以實現空間失穩目標三維重建為目的,解決目標失穩及線陣雷達掃描動態目標等多種因素所導致的目標三維模型無法精確恢復問題。首先根據角精度與測距精度兩大判定條件計算畸變容忍度最大臨界值,其次對目標點云進行自修正式迭代配準,直至滿足三維重建要求。針對本文獲取的目標模型,三維重建精度可達到±1°以內,重建總耗時在現有配置環境下為72467 ms。該方法不僅能夠在理論上指導三維目標的點云數據配準、非合作目標的三維姿態測量和目標三維重建等工作,還能在實際應用中為線陣式測量技術的進一步發展提供新的研究方向。

表2 相鄰組旋轉配準結果

圖14 相鄰點云第6次配準

圖15 目標最終三維重建模型

猜你喜歡
方向測量
2022年組稿方向
計算機應用(2022年2期)2022-03-01 12:33:42
2022年組稿方向
計算機應用(2022年1期)2022-02-26 06:57:42
2021年組稿方向
計算機應用(2021年4期)2021-04-20 14:06:36
2021年組稿方向
計算機應用(2021年3期)2021-03-18 13:44:48
2021年組稿方向
計算機應用(2021年1期)2021-01-21 03:22:38
把握四個“三” 測量變簡單
滑動摩擦力的測量和計算
滑動摩擦力的測量與計算
測量的樂趣
測量
主站蜘蛛池模板: 九九九九热精品视频| 在线色综合| 久久久久无码国产精品不卡| 国产精品视频系列专区| 黄色在线网| 在线观看91香蕉国产免费| 久久久久亚洲精品无码网站| 亚洲中文在线看视频一区| 视频一区视频二区日韩专区| 久久久噜噜噜| 内射人妻无套中出无码| 亚洲欧美激情小说另类| 久久精品国产999大香线焦| 69av在线| 国产精品密蕾丝视频| jizz在线观看| 激情无码字幕综合| 99久久性生片| 97在线视频免费观看| 欧美激情视频一区| 亚洲国产清纯| 国产亚洲精| 欧美成人亚洲综合精品欧美激情 | 制服无码网站| 国产精欧美一区二区三区| 青青青国产在线播放| 中文无码毛片又爽又刺激| 久久综合一个色综合网| 国产成人欧美| 免费国产福利| 免费在线看黄网址| 亚洲美女一级毛片| 国产超薄肉色丝袜网站| 亚洲人成网站18禁动漫无码| a亚洲视频| 福利在线一区| 久久精品国产免费观看频道| 激情六月丁香婷婷| 性做久久久久久久免费看| 97国产精品视频人人做人人爱| 伊人成人在线视频| 午夜不卡视频| 国产综合色在线视频播放线视 | 中文字幕 欧美日韩| 欧美一级高清视频在线播放| 国产精品人成在线播放| 亚洲综合亚洲国产尤物| 亚洲狠狠婷婷综合久久久久| 亚洲欧美日韩精品专区| 国产麻豆精品久久一二三| 国产欧美中文字幕| 国产精品网址你懂的| 第九色区aⅴ天堂久久香| 国产乱视频网站| AV天堂资源福利在线观看| 久久黄色毛片| 毛片一区二区在线看| 免费激情网站| 午夜老司机永久免费看片| 国产精品蜜芽在线观看| 亚洲欧美另类久久久精品播放的| 99久久国产综合精品2023| 一级毛片高清| 国产91精品调教在线播放| 999国内精品久久免费视频| 亚洲AV无码乱码在线观看裸奔| 中文纯内无码H| 亚洲第一在线播放| 色男人的天堂久久综合| 一区二区午夜| 日韩精品亚洲一区中文字幕| av在线无码浏览| 成人在线亚洲| 免费欧美一级| 色噜噜综合网| 在线精品亚洲一区二区古装| 波多野结衣无码视频在线观看| 国产一二三区视频| 欧洲熟妇精品视频| 毛片免费高清免费| 亚洲中文字幕日产无码2021 | 亚洲swag精品自拍一区|