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

多傳感器空間配準的病態性及其穩健估計

2017-11-09 09:22:16潘江懷喬慧
兵工學報 2017年10期

潘江懷, 喬慧

(江蘇自動化研究所, 江蘇 連云港 222006)

多傳感器空間配準的病態性及其穩健估計

潘江懷, 喬慧

(江蘇自動化研究所, 江蘇 連云港 222006)

為提高多傳感器空間配準的精度和穩定性,以最小二乘空間配準算法為基礎,分析多傳感器空間配準病態性產生的原因,指出系統誤差估計Fisher信息矩陣的病態性是影響空間配準質量的主要因素。運用條件指標對病態性程度進行評價,并根據信息矩陣的奇異值呈階梯型分布特點,提出一種混合的奇異值修正穩健估計算法。該算法通過只對較小的條件指標進行修正,實現對系統誤差觀測方程中不確定成分的有效抑制和確定成分的保留,能較好解決多傳感器空間配準系統誤差實時估計問題;通過兩種不同場景下的系統誤差估計結果證明,該算法的使用可顯著地提高系統誤差參數估計的精度和穩定性。

信息處理技術; 空間配準; 病態性; 穩健估計; 奇異值分解; 多傳感器

0 引言

在多傳感器融合跟蹤系統中,傳感器探測的系統誤差是影響融合精度的最主要因素[1],系統誤差的存在將直接對數據融合中的跟蹤門計算、目標航跡狀態估計和協方差計算、航跡相關波門及準則產生不利影響[2]。若直接將各傳感器的探測數據進行融合,則將由于系統誤差的存在而使跟蹤結果惡化,導致目標一批變多批或航跡呈鋸齒、航跡斷航,甚至不如單傳感器的跟蹤效果[3]。空間配準是指消除探測源/探測系統所存在的系統誤差或將其控制在容許范圍內的處理過程。在多平臺多傳感器系統中,一般是指將分布在不同地理位置的傳感器目標探測信息轉換到一個統一坐標系或融合中心過程中,使由于傳感器探測、平臺定位不準以及坐標轉換等原因產生的系統誤差(配準誤差)最小化[4],以實現同一目標位置信息在不同傳感器中的迭合和正確交接,確保后續數據融合正確進行。因此,在多傳感器融合跟蹤系統中,必須對傳感器進行空間配準[5]。

在多傳感器系統誤差估計過程中,影響偏差可觀測性的因素有很多,偏差的種類、偏差之間相互耦合的程度、目標與傳感器之間的幾何關系等多種因素都會造成偏差估計可觀測性的變化,從而導致偏差估計可觀測度難以進行定量的描述[6]。偏差估計精度除了受到估計算法的影響,主要是由系統可觀測度決定,因此如何在未知可觀測度條件下定量地確定系統的配準效果就成為一個難題[7]。現有的空間配準算法如歐洲和北美洲廣泛使用SAGE系統監視系統中的實時質量控制(RTQC)法[4]、考慮傳感器測量噪聲的最小二乘(LS)配準法[5]及廣義LS(GLS)配準法[6]、精確極大似然(EML)配準法[7]、極大似然(MLR)配準法[8]及文獻[9-14]提出的空間配準算法,在可觀測度較高、目標狀態觀測完整、傳感器布站合理且累計了一定的觀測數據條件下,通常都可以獲得比較好的配準結果。當可觀測度較低時,偏差的估計值將發生巨大的波動[15-16]。這是因為觀測數據中緩變的觀測信息被劇烈變化的隨機噪聲湮沒,觀測數據的變化主要由于隨機誤差引起,導致了估計結果的波動,若直接根據這種估計值去修正目標航跡,則后續的融合結果將必然會受到質疑[17]。

本文基于LS空間配準算法,分析了LS系統誤差估計的病態性產生原因,并根據系統誤差估計的Fisher信息矩陣奇異值分布為階梯型分布特點,提出一種新的奇異值修正(CSVD)算法,對奇異特征值進行修正。該算法可兼顧系統誤差參數估計的分辨率與方差,實現了對系統誤差觀測方程中不確定成分的有效抑制和確定成分的有效保留,能有效解決多傳感器空間配準系統誤差實時估計問題,并從兩種不同的仿真場景進行了驗證比較。結果表明:當信息矩陣的條件數明顯分梯形分布時,該方法可顯著提高系統誤差估計結果的精度和穩定性。

1 多傳感器系統誤差的LS估計

對于不同平臺的n部三坐標雷達而言,需要估計的系統誤差參數有3n個之多,(Δri,Δβi,Δεi)分別表示為第i個傳感器的距離系統誤差、方位角系統誤差、俯仰角系統誤差,i=1,…,n. 為簡化推導,下面以兩部傳感器為例對空間配準的算法進行分析。如圖1所示,假設兩傳感器的大地坐標為(L1,B1,H1)和(L2,B2,H2),兩傳感器在地心地固(ECEF)坐標系中的位置為(xp1,yp1,zp1)、(xp2,yp2,zp2),k時刻兩傳感器對目標距離、方位角和俯仰角的探測為(r1(k),β1(k),ε1(k))、(r2(k),β2(k),ε2(k)),轉換到兩傳感器地理坐標系下的位置為(x1(k),y1(k),z1(k))、(x2(k),y2(k),z2(k)). 以下省略時刻k.

圖1 有系統誤差下的多平臺傳感器探測Fig.1 Multi-sensor detection with bias

若沒有任何誤差,兩傳感器同一時刻對同一目標的探測在ECEF坐標系中應該迭合,即

(1)

式中:TR1和TR2分別為兩傳感器地理坐標系到ECEF坐標系的旋轉矩陣;(xi,t,yi,t,zi,t)為目標在傳感器i(i=1,2)地理坐標系下的真值:

(2)

(3)

兩傳感器對同一目標的探測在ECEF坐標系中的坐標之差如下:

(4)

簡記為

Z=Ab,

(5)

式中:A=[TR1Jf1,-TR2Jf2];b=[Δr1,Δβ1,Δε1,Δr2,Δβ2,Δε2]T.

當有k(k>2)個量測值時,(5)式的LS解為

=(ATA)-1ATZ.

(6)

定義系統誤差參數b的Fisher信息矩陣為

J=ATA.

(7)

2 LS系統誤差估計的病態性

LS估計系統誤差的條件是估計參數b的信息矩陣J可逆,即A列滿秩。當A的列線性相關或弱相關時,信息矩陣是病態的。

2.1 常見的病態場景

1) 傳感器平臺相距較近。當兩傳感器相距較近時,兩傳感器對同一目標的探測距離、方位角、俯仰角近似相等,即有(r1(k),β1(k),ε1(k))≈(r2(k),β2(k),ε2(k)),而 (7) 式中矩陣A的系數完全由傳感器對目標探測的距離、方位角、俯仰角所決定,將導致矩陣A的列向量具有相關性,從而使信息矩陣J是病態的。

2) 觀測數據密集情況。當目標相對于傳感器平臺靜止不動或者活動范圍很小時,傳感器在觀測時間內對目標的距離、方位角、俯仰角探測近似相等,即(r1(1),β1(1),ε1(1))≈…≈(r1(k),β1(k),ε1(k)),此時也會導致矩陣A的列向量具有相關性,從而使信息矩陣J是病態的。

2.2 病態性的評估

一般地,通過條件數對方程病態性進行評估。對于(5)式,小擾動引起的估計相對系統誤差有如下關系:

(8)

當A非奇異時,A條件數定義如下:

(9)

式中:λl為矩陣A的第l個特征值。

可見,系數矩陣A的最大最小特征值越分散,cond{A}越大,解向量b的相對誤差越大。條件數是一個相對數,其值的大小表示了最小特征根相對于最大特征根“小”的程度,即條件數從某種程度上度量了病態性的嚴重程度。但條件數不能判定觀測矩陣A的數據列之間存在幾個相關關系,也不能確切斷定出每個相關性到底存在于哪些數據列之間,這些信息對于全面把握病態性產生的機理和實質,以及判定哪個系統誤差參數的可觀測性高低具有重要的參考價值。

定義條件指標:

(10)

稱ηj法矩陣的第j個條件指標,λ1為最大的特征值,顯然ηj≥1.

條件指標直接反映了具體的系統誤差參數項估計受測量噪聲干擾的影響程度。若出現高條件指標,則矩陣A的該列元有微小擾動時,該系統誤差參數項的估計結果將引起相當大的變化,從而觀測矩陣A的數據列之間存在相關關系。大量模擬研究表明[18-19],若病態性很弱,則條件指標小于100;若病態性較強,則條件指標在100~1 000之間;若病態性很嚴重,則條件指標在1 000以上。

3 多傳感器系統誤差CSVD穩健估計算法

由系統誤差參數b的Fisher信息矩陣定義可知,J是一個實對稱矩陣,根據矩陣奇異值分解定理,對于任意實對稱矩陣Jm×m,存在正交矩陣Vm×m使得:

J=VSVT,

(11)

J的Mooer-Penrose廣義逆為

J+=VS-1VT,

(12)

系統誤差參數估計的方差為

(13)

向量形式:

(14)

TSVD法就是將較小的奇異值舍棄,即Σ=diag(λ1,λ2,…,λr)舍棄較小的,使其為0,以減少參數解的方差、增加解算的穩定性。若令D=diag(λ1,λ2,…,λr1) (r1

(15)

(15)式向量形式表示為

(16)

在具體采用TSVD法時,主要需要考慮按什么原則去截斷。

(17)

參數估計的向量形式為

(18)

參數估計的方差為

(19)

從(18)式和(19)式可以看出,適當選擇α可以減少估計的方差、增加估計的穩定性。但上述修正方法存在如下缺陷:1)當奇異值分布呈均勻下降型時,MSVD法對不可靠的成分進行了抑制,它不但修正了較小的奇異值,也修正了相對較大的奇異值,從而使模型中的確定成分發生了畸變,降低了分辨率;2)TSVD法將小于方差門限值的奇異值全部截掉,將引起分辨率的明顯變壞,因為保留的和截掉的成分差別不明顯。

在進行奇異值分解之后,相對大的奇異值及相應的左右特征向量表示模型參數中比較肯定的和可靠的部分,很小的奇異值和相應的特征向量則表示不可靠的部分。TSVD法在處理病態問題的奇異值時,刪除了小的奇異值及對應的特征向量,實際上刪掉了模型參數中不可靠的部分,以減少解的方差。但是這種做法同時也可能嚴重損害了解估計的分辨率。

對于多傳感器空間配準問題,同時對傳感器的測距和測向系統誤差進行估計時,其觀測矩陣的特征值的大小順序呈階梯狀態的分布,如(20)式所示:

η1≤η2≤η3<η4≤η5≤η6.

(20)

理想的做法是修正奇異值較小的,而對較大的奇異值不進行修正或給予不同于較小的奇異值的修正。本文綜合了TSVD法和MSVD法的優點,提出了一種新的CSVD方案,即對大的奇異值盡量趨近于用TSVD法進行修正,對小的奇異值趨近于用MSVD法進行修正。具體做法為:對超過門限值t的奇異值保持不變,對小于門限值t的奇異值進行修正,修正公式為

(21)

相應的修改后的奇異值為

Σδ=diag(δ1,δ2,…,δr).

(22)

相應的配準偏差估計結果為

(23)

1)配準估計結果的分辨率和方差。偏差估計結果(23)式的分辨率為

(24)

當λp≥t時,(23)式解估計的方差為

(25)

由方差公式(24)式可以看出,與LS法估計的方差相比,新的解估計的方差減小了。本文的方法兼具了MSVD法和TSVD法的優點,通過降低分辨率減小了解的方差,解估計的質量得到了提高。

2)最小奇異值門限t的確定。判斷病態觀測方程式的Fisher信息矩陣ATA的條件數大于103,可以得到

(26)

設奇異值λq(q=1,2,…,r)滿足

(27)

則選擇最小奇異值門限t為

t=λq.

(28)

4 仿真驗證

4.1 可觀測性較好的情況

目標1的初始位置[121.05°,26.6°,2 km],運動速度[121.2°,26.5°,0 km];平臺1初始位置[121.2°,26.5°,0 km],運動速度[0 m/s,0 m/s,0 m/s];平臺2初始位置[121.4°,26.5°,0 km],運動速度[0 m/s,0 m/s,0 m/s];傳感器1探測系統誤差:Δr1=-500 m,Δβ1=0.9°,Δε1=0.5°,噪聲σr1=100 m,σβ1=0.3°,σε1=0.3°;傳感器2探測系統誤差:Δr2=400 m,Δβ2=0.9° ,Δε2=-0.5°,噪聲σr2=100 m,σβ2=0.3° ,σε2=0.3°.

仿真場景如圖2所示,仿真結果如圖3~圖9所示。

圖2 仿真場景示意圖Fig.2 Schematic diagram of simulation scene

圖3 條件指標變化情況Fig.3 Changes of condition indexes

圖4 CSVD修正之后條件指標變化情況Fig.4 Changes of modified condition indexes

圖5 兩傳感器距離系統誤差估計結果Fig.5 Estimated range deviations of two sensors

圖6 兩傳感器方位角系統誤差估計結果Fig.6 Estimated azimuth deviations of two sensors

圖7 兩傳感器俯仰角系統誤差估計結果Fig.7 Estimated elevation deviations of two sensors

圖8 CSVD實時配準前后兩傳感器探測結果Fig.8 Detection results of two sensors before and after real time CSVD registration

圖9 實時配準前后兩傳感器誤差比較Fig.9 Comparison of detection errors of two sensors before and after real time CSVD registration

4.2 可觀測性較差的情況

目標1初始位置[120.95°,26.6°,2 km],運動速度[0 m/s,100 m/s,0 m/s];平臺1初始位置[121.2°,26.5°,0 km],運動速度[0 m/s,0 m/s,0 m/s];平臺2初始位置[121.4°,26.5°,0 km],運動速度[0 m/s,0 m/s,0 m/s];傳感器1探測系統誤差:Δr1=-500 m,Δβ1=0.9°,Δε1=0.5°,噪聲σr1=100 m,σβ1=0.3°,σε1=0.3°;傳感器2探測系統誤差:Δr2=400 m, Δβ2=0.9°,Δε2=-0.5°,噪聲σr2=100 m,σβ2=0.3°,σε2=0.3°.

仿真場景如圖10所示,仿真結果如圖11~圖16所示。

圖10 仿真場景示意圖Fig.10 Schematic diagram of simulation scene

圖11 條件指標變化情況Fig.11 Change of condition indexes

圖12 CSVD修正之后條件指標變化情況Fig.12 Changes of modified condition indexes

圖13 兩傳感器距離系統誤差估計結果Fig.13 Estimated range deviations of two sensors

圖14 兩傳感器方位角系統誤差估計結果Fig.14 Estimated azimuth deviations of two sensors

圖15 兩傳感器俯仰角系統誤差估計結果Fig.15 Estimated elevation deviations of two sensors

圖16 CSVD實時配準前后兩傳感器探測結果Fig.16 Detection results of two sensors before and after real time CSVD registration

圖17 實時配準前后兩傳感器誤差比較Fig.17 Comparison of detection errors of two sensors before and after real time CSVD registration

對比圖3和圖17中兩種場景的條件指標可以看出:1)第1種場景的條件指標明顯比第2種場景的要低;2)兩個場景之間一開始均存在3個變量之間的相關性,開始階段均只能估計出相對的系統誤差;3)隨著數據的積累,第1種場景的可觀測性明顯變好,當數據量超過100時所有的待估計參數之間是完全獨立的,此時使用經典的算法均可得到有效的結果。

5 結論

在多傳感器空間配準中,信息矩陣的病態性是影響系統誤差參數估計質量的主要因素之一。本文針對LS法系統誤差估計過程中由于信息矩陣的病態性導致估計結果無效的情況,分析LS法系統誤差估計的病態性產生原因,并根據Fisher信息矩陣奇異值的分布為階梯型分布特點,提出一種新的CSVD方案對奇異特征值進行修正,該算法兼顧參數估計的分辨率與方差,實現了對觀測方程中不確定成分的有效抑制和確定成分的有效保留,有效解決了多傳感器空間配準系統誤差實時估計問題,并從兩種不同的仿真場景進行了驗證比較,結果表明:當信息矩陣的條件數明顯分梯形分布時,該方法可顯著提高系統誤差估計結果的精度和穩定性。

References)

[1] Bo Y, Chen Z, Yin M, et al. Improved different dimen-sional sensors combined space registration algorithm[J]. Mathematical Problems in Engineering, 2015(1):1-9.

[2] 程志剛. 多雷達組網系統空間配準方法研究[J]. 艦船科學技術, 2015, 37(12):131-134.

CHENG Zhi-gang. Research on radar space registration method for multi-radar network system[J]. Ship Science and Technology, 2015, 37(12):131-134.(in Chinese)

[3] 胡勤振, 蘇洪濤, 劉子威. 配準誤差下的多基地雷達目標檢測算法[J]. 電子與信息學報, 2017,39(1):87-94.

HU Qin-zhen, SU Hong-tao, LIU Zi-wei. Target detection algorithm for multistatic radar with registration errors[J]. Journal of Electronics & Information Technology, 2017,39(1):87-94.(in Chinese)

[4] 史江林, 周學平, 劉長海,等. 基于EM算法的聯合數據關聯與空間配準[J]. 現代雷達, 2016, 38(1):51-55.

SHI Jiang-lin, ZHOU Xue-ping, LIU Chang-hai, et al. Joint data association and spatial registration based on EM[J]. Modern Radar, 2016, 38(1):51-55.(in Chinese)

[5] 方峰, 蔡遠利. 機動目標多傳感器組網空間配準方法[J]. 固體火箭技術, 2016, 39(4):574-579.

FANG Feng, CAI Yuan-li. Multi-sensor space registration for maneuvering target tracking [J]. Journal of Solid Rocket Technology, 2016, 39(4):574-579. (in Chinese)

[6] Okello N N, Chall S L.A joint radar registration and track-to-track fusion for distributed trackers[J]. IEEE Transactions on Aerospace and Electronic Systems, 2004, 40(3):808-823.

[7] Herman S M, Poore A B.Nonlinear least-squares estimation for radar and navigation biases[J]. Proceedings of SPIE, Signal and Data Processing of Small Targets, 2006, 6236:1701-1717.

[8] Okello N N, Ristic B. Maximum likelihood registration for multiple dissimilar radars[J]. IEEE Transactions on Aerospace and Electronic Systems, 2003, 39(3):1235-1240.

[9] 潘江懷, 李洪梅, 何佳洲.三維傳感器組網偏差估計方法[J].火力與指揮控制, 2008, 4(4):104-106.

PAN Jiang-huai, LI Hong-mei, HE Jia-zhou.Bias estimation in 3D sensor networks[J]. Fire Control and Command Control, 2008, 4(4):104-106.(in Chinese)

[10] 宋強, 崔亞奇, 何友.反饋式多目標多傳感器系統誤差融合估計技術[J]. 宇航學報, 2011, 32(1):115-122.

SONG Qiang, CUI Ya-qi, HE You.Multi-target multi-sensor systematic error fusion estimation technique with feedback[J]. Journal of Astronautics, 2011, 32(1):115-122 .(in Chinese)

[11] 王成飛, 王航宇, 石章松, 等. 基于高精度導航設備的海上多傳感器配準算法[J].控制理論與應用, 2011, 28(4):497-503.

WANC Cheng-fei, WANG Hang-yu, SHI Zhang-song, et al. Maritime sensor registration algorithm based on high-precision navigation equipment[J]. Control Theory and Applications, 2011, 28(4):497-503.(in Chinese)

[12] 崔亞奇, 熊偉, 何友.基于MLR的機動平臺傳感器誤差配準算法[J].航空學報, 2012, 33(1):118-128.

CUI Ya-qi, XIONG Wei, HE You.Mobile platform sensor registration algorithm based on MLR[J]. Acta Aeronautica et Astronautica Sinica, 2012, 33(1):118-128.(in Chinese)

[13] 連峰, 韓崇昭, 彭一峰, 等.基于廣義似然比的自適應在線配準算法[J]. 控制與決策, 2009, 24(1):23-28.

LIAN Feng, HAN Chong-zhao, PENG Yi-feng, et al.Adaptive on-line registration algorithm based on GLR[J]. Control and Decision, 2009, 24(1):23-28.(in Chinese)

[14] 吳衛華, 江晶. WGS-84坐標系下多空基無源傳感器最大似然配準[J]. 系統工程與電子技術, 2015, 37(2): 304-309.

WU Wei-hua, JIANG Jing. Maximun likelihood registration for passive sensors of multiple airborne platforms in WGS-84[J]. Systems Engineering and Electronics, 2015, 37(2):304-309.(in Chinese)

[15] Huang D, Leung H.A pseudo-measurement approach to simultaneous registration and track fusion[J]. IEEE Transactions on Aerospace and Electronic System, 2012, 48(3):2315-2331.

[16] Steafano F, Alfonso F, Fulvio G, et al.Least squares estimation and Cramer-Rao type lower bounds for relative sensor registration process[J]. IEEE Transactions on Signal Processing, 2011, 59(3):1075-1087.

[17] Boer Y,Driessen J N. Multi-target particle filter track before detect application[J]. IEEE Proceedings Radar, Sonar and Navigation, 2004, 151(6):351-357.

[18] Ong H T.Sensor registration using airlanes maximum likelihood solution[J]. Proceedingd of SPIE, Signal and Data Processing of Small Targets, 2003,5204:390-401.

[19] Watson G A, Rice T R.Sensor alignment and compensation for composite tracking[J]. Proceedings of SPIE, Signal and Data Processing of Small Targets, 2002, 4728: 354-367.

Ill-conditionAnalysisandRobustEstimationforMulti-sensorSpatialRegistration

PAN Jiang-huai, QIAO Hui
(Jiangsu Automation Research Institute, Lianyungang 222006,Jiangsu, China)

In order to improve the accuracy and stability of multi-sensor spatial registration, the causes of the ill condition of the multi sensor spatial registration are analyzed based on LS estimation for multi-sensor spatial registration. The analysis result shows that the ill-condition of the Fisher information matrix is the key factor which influences the quality of system error parameter estimation. The conditional index of information matrix is used to measure the ill-condition of spatial registration. A new algorithm of compound modified singular value decomposition (CSVD) for robust estimation is proposed according to the characteristics of trapezia distribution of the information matrix singular values. The algorithm can effectively suppress uncertainty components and retain deterministic components in the observation equation by modifing only minor conditional index, and can effectively solve the multi sensor real time estimation of spatial registration. Two different scenarios of the estimatied results show that CSVD method can be used to significantly improve the estimation accuracy and stability of system error parameters.

information processing technology; spatial registration; ill-condition; robust estimation; singular value decomposition; multi-sensor

2016-03-02

國家“973”計劃項目(613101)

潘江懷(1982—),男,高級工程師。E-mail:panjianghuai@163.com

TN953+.6

A

1000-1093(2017)10-1965-09

10.3969/j.issn.1000-1093.2017.10.012

主站蜘蛛池模板: 国产麻豆福利av在线播放| 国产麻豆另类AV| 日本伊人色综合网| 国产天天色| 伊人久热这里只有精品视频99| 国产乱人伦AV在线A| 亚洲人成网站色7777| 婷婷久久综合九色综合88| 亚洲男人天堂网址| 午夜性爽视频男人的天堂| 亚洲中文字幕无码mv| 女人天堂av免费| 最新国语自产精品视频在| 国产农村1级毛片| AⅤ色综合久久天堂AV色综合| 亚洲精品国产精品乱码不卞| 欧美成人午夜影院| 波多野结衣无码中文字幕在线观看一区二区 | 狼友av永久网站免费观看| 亚洲AⅤ综合在线欧美一区| 九九热免费在线视频| 欧美国产综合视频| 日韩中文精品亚洲第三区| 久久久精品久久久久三级| 青青久久91| 久久亚洲黄色视频| 国产一级α片| 网友自拍视频精品区| 日韩欧美国产精品| 四虎永久在线视频| 欧美日韩国产在线人| 久久这里只有精品国产99| 久久精品无码国产一区二区三区| 人妻丰满熟妇AV无码区| 人妻丰满熟妇αv无码| 中文字幕av无码不卡免费| av性天堂网| 国产精品女熟高潮视频| 狠狠做深爱婷婷久久一区| 中文字幕色在线| 国内精品久久人妻无码大片高| 男人天堂亚洲天堂| аv天堂最新中文在线| 国产香蕉国产精品偷在线观看| 久久精品国产免费观看频道 | 亚洲三级视频在线观看| 噜噜噜综合亚洲| 国产视频一区二区在线观看| 日韩在线观看网站| 香蕉eeww99国产精选播放| 国产美女无遮挡免费视频网站| 欧美笫一页| 欧美日韩中文国产va另类| 视频二区欧美| 免费黄色国产视频| 婷婷激情亚洲| 亚洲人成日本在线观看| 91精品国产综合久久不国产大片| 亚洲日韩精品欧美中文字幕| 久久无码免费束人妻| 亚洲欧美综合另类图片小说区| 一本久道久久综合多人| 毛片大全免费观看| 一级毛片在线免费看| 亚洲美女一区二区三区| 欧美日韩第三页| 亚洲日韩Av中文字幕无码| 一级毛片免费播放视频| 日韩国产黄色网站| 69国产精品视频免费| 五月六月伊人狠狠丁香网| 爆乳熟妇一区二区三区| 日韩 欧美 国产 精品 综合| 精品无码国产自产野外拍在线| 亚洲91精品视频| 久热这里只有精品6| 色噜噜狠狠狠综合曰曰曰| a毛片免费在线观看| 中文国产成人精品久久| 久草视频福利在线观看| 国产精品一区二区不卡的视频| 国产精品自拍露脸视频|