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

        ?

        顧及雙程聲徑的常梯度聲線跟蹤水下定位算法

        2022-03-07 12:03:36閆鳳池王振杰聶志喜李偉嘉
        測繪學報 2022年1期
        關鍵詞:信號

        閆鳳池,王振杰,趙 爽,聶志喜,孫 振,李偉嘉

        中國石油大學(華東),山東 青島 266580

        海洋是人類資源環(huán)境可持續(xù)發(fā)展的希望[1]。海洋導航定位技術是海洋開發(fā)和海洋高技術應用的基本前提,是進行海洋活動的重要基礎[2]。由于電磁波在海水中衰減速度過快,無法滿足水下定位的需要[3],而聲波因其良好的傳播特性,可以進行遠距離、大深度的傳播,得到廣泛的研究和應用,成為水下定位的主要觀測手段[4]。

        海洋中,聲速受海水溫度、鹽度、密度、海流等各種復雜海洋環(huán)境的影響,給高精度水下聲學定位帶來了極大挑戰(zhàn)[5]。水下聲學定位中最為顯著的誤差是聲線彎曲誤差和聲速誤差,為了消除或減弱聲線彎曲誤差的影響,需要進行聲速改正[6]。目前,常用的水下定位聲速改正方法有加權平均聲速法、等效聲速法、聲線跟蹤法等[7]。其中,加權平均聲速法最為簡單,但定位精度受限[8];等效聲速法可以有效地提高定位的精度,但無法完全消除聲線折射效應的影響[9];聲線跟蹤法對聲信號進行了聲線彎曲改正,計算精度最高[10]。

        為削弱聲速誤差等系統(tǒng)誤差對定位的影響,可以采用差分定位的方法,包括空間差分和歷元間差分兩種模式??臻g差分類似于GNSS差分定位技術[11],削弱觀測值中的共性誤差,但也抵消了非差觀測值中的垂向信息,降低了豎直方向的定位精度。文獻[12]提出歷元間差分方法,通過對相鄰觀測歷元之間的差分,削弱與時間相關聯(lián)的系統(tǒng)誤差影響。文獻[13]將抗差估計與差分定位結合,提出了基于抗差估計的歷元間單差和雙差定位模式,不僅削弱時間關聯(lián)誤差的影響,同時有效控制了異常誤差的影響。

        文獻[14]顧及聲速相關系統(tǒng)誤差的相似性,提出采用差分組合觀測值方法。假設局部范圍內(nèi)聲速穩(wěn)定,聲速對測距影響近似相同,借助距離差解算海底應答器坐標,雖然削弱聲信號傳播時間上的系統(tǒng)誤差影響,但該方法垂直定位結果不穩(wěn)定。文獻[15]提出利用最小二乘反演法確定海底應答器陣列位置,通過不同應答器之間的基線長和相對深度差對應答器坐標進行約束。但距離計算采用平均聲速法,且模擬試驗中測量船始終位于應答器陣列中心位置上方海面,因此有一定局限性。以上方法雖然在一定程度上提高了水下定位的精度,但由于未能考慮到聲信號往返的實際原理,因此定位精度有限。

        考慮到聲速誤差不僅與復雜的海洋環(huán)境有關,也受到觀測儀器自身的限制[16]。當利用觀測區(qū)域的聲速剖面進行聲線跟蹤時,理論上可以完全消除聲線彎曲影響[17],但考慮到水下聲信號傳播的實際情況及船載換能器在聲信號傳播過程中的位移,需要構建更加合理的定位模型[18]。

        目前,水下聲學定位方法通常采用基于聲信號的單向傳播路徑(單程聲徑),忽略了測量船持續(xù)走航引起的換能器在聲信號傳播過程中的位置移動,進而會降低應答器定位精度[19]。針對上述問題,本文結合水下聲學定位中聲信號實際傳播過程,構建了顧及聲信號雙向傳播路徑(雙程聲徑)的定位模型,發(fā)展了基于雙程聲徑的分層常梯度聲線跟蹤定位算法。

        1 單程聲徑水下定位方法

        單程聲徑水下定位的原理是將船載換能器測量的聲信號往返時間取二分之一作為其單程傳播時間,結合已知聲速剖面,獲得船載換能器與海底應答器之間的斜距,然后利用多個斜距建立誤差方程組,通過最小二乘求解非線性方程組獲得海底應答器位置[20]。單程聲徑水下定位原理如圖1所示。

        圖1 單程聲徑水下定位原理Fig.1 Principle of high precision underwater positioning with single-trip acoustic path

        1.1 單程聲徑水下定位模型

        在單程聲徑水下定位過程中,假設海底應答器坐標為Xo(x0,y0,z0),通過船載GNSS接收機獲取測量船的位置信息[21],如圖1所示。結合測量船的姿態(tài)傳感器測得的姿態(tài)信息,計算出船載換能器各歷元的坐標Xi(xi,yi,zi),觀測歷元i=1,2,…,n。船載換能器測量的聲信號往返時間間隔2ti,則單向傳播時間為ti。在已知聲速剖面的情況下,通過單程聲徑聲線跟蹤算法計算船載換能器到海底應答器之間的幾何距離ρi,具體計算方法見1.2。

        單程聲徑水下定位模型為[22]

        ρi=f(Xi,Xo)+δρdi+εi

        (1)

        式中,f(Xi,Xo)是船載換能器到海底應答器的真實幾何距離;δρdi是與聲速和測時有關的系統(tǒng)誤差;εi是隨機誤差。

        船載換能器到海底應答器的真實幾何距離f(Xi,Xo)表示為

        (2)

        將式(1)線性化,可得

        (3)

        在忽略各類系統(tǒng)誤差的情況下,單程聲徑水下定位觀測方程為

        L=AdX+Δ

        (4)

        式中,L是n維的觀測向量;A是n×3的系數(shù)矩陣;dX是海底應答器坐標改正數(shù);Δ是由多個歷元的εi+biεXi組成的矩陣。

        觀測向量

        通過最小二乘求解海底應答器坐標改正數(shù)為

        dX=(ATPA)-1ATPL

        (5)

        式中,P為觀測值的權矩陣,可基于聲線入射角確定或根據(jù)驗后殘差進行確定[12]。

        觀測值改正數(shù)

        V=AdX-L

        (6)

        分別計算單位權中誤差和協(xié)方差陣

        (7)

        式中,r為多余觀測數(shù);n為總觀測數(shù);t為必要觀測數(shù),一般為3;QXX為平差參數(shù)的協(xié)因數(shù)陣。

        1.2 單程聲徑聲線跟蹤算法

        聲線跟蹤的理論基礎是聲速分層假設,即任何復雜的聲速剖面結構,都可以近似為由多層具有簡單結構的聲速層組成。此假設用每個獨立層內(nèi)的聲速折線分布來代替整個聲速剖面的連續(xù)變化。目前,常用的聲速分層形式有常聲速分層和常聲速梯度分層兩種,其中,常聲速分層認為每層內(nèi)聲速不變,聲線沿直線傳播;常聲速梯度分層認為層內(nèi)聲速線性變化,聲線沿曲線傳播。本文采用的是常梯度聲線跟蹤法。假設將聲信號經(jīng)過的水柱分為N個等梯度層,層內(nèi)聲速的變化是常梯度的,即認為層內(nèi)聲速線性變化[24]。通過逐層跟蹤的方法計算聲信號的水平位移、豎直位移和傳播時間,如圖2所示。

        圖2 單程聲徑聲線跟蹤算法Fig.2 Acoustic ray-tracing algorithm of single-trip acoustic path

        在水深zi對應的i層,用ci表示聲信號的傳播速度。第i層內(nèi)的聲速梯度gi可以用式(8)表示為[25]

        (8)

        式中,Δzi為第i層的水層厚度,ci+1為第i+1層聲速。

        在第i層內(nèi),聲線曲率處處相等,聲線在層內(nèi)的軌跡為一段圓弧。該圓弧的曲率半徑為[26]

        (9)

        式中,p是Snell常數(shù)。

        由Snell定律和常梯度聲速假設[27],聲信號在第i層內(nèi)的水平距離Δyi,豎直距離Δzi和時間Δti為

        (10)

        式中,θi和θi+1分別為聲信號的入射角和出射角。

        對完整層進行聲線跟蹤后,需要依據(jù)聲信號剩余傳播時間進行剩余層的計算。假設剩余時間為tr,則

        (11)

        聲信號在剩余層的水平位移和豎直位移為

        (12)

        聲信號的總水平位移和總豎直位移為

        (13)

        改正后幾何距離為

        (14)

        將ρi代入式(1),構建單程聲徑聲線跟蹤的觀測方程,通過式(5)求解海底應答器坐標改正數(shù),通過迭代計算完成單程聲徑水下定位[28]。

        該方法得到的單程聲徑幾何距離與真實聲線基本一致,有效地提高了傳統(tǒng)的水下聲學定位方法的定位精度,這也是目前水下定位中最常用的方法之一。但單程聲徑模型存在局限性,即忽略了換能器在聲信號往返過程的位置變化,直接利用換能器聲信號發(fā)射或接收時刻的坐標進行解算,獲得海底應答器的坐標,從而降低了海底應答器位置解算精度。

        2 雙程聲徑水下高精度定位

        測量船在海面持續(xù)走航過程中,船載換能器在聲信號發(fā)射和接收時的位置會發(fā)生變化,故聲信號往、返的路徑存在差異。若忽略換能器的位置變化,只考慮聲線單向傳播過程,將會在海底應答器定位解算時引入系統(tǒng)性偏差,降低水下定位的精度。針對單程聲徑水下定位方法中存在的問題,為削弱換能器位置變化對水下定位的影響,有必要研究更符合水下聲學定位實際的模型和解算方法。如圖3所示,本文利用船載換能器聲信號發(fā)射和接收時刻的位置計算聲信號傳播的雙程聲徑幾何距離,構建雙程聲徑水下高精度定位模型,通過解算,實現(xiàn)海底應答器的高精度定位。

        圖3 雙程聲徑水下高精度定位原理Fig.3 Principle of high precision underwater acoustic positioning with round-trip acoustic path

        2.1 雙程聲徑水下高精度定位模型

        (15)

        將式(15)線性化后,得

        bidX+δρdi+εi+niεXi

        (16)

        在忽略各類系統(tǒng)誤差的情況下,雙程聲徑水下高精度定位觀測方程為

        L=BdX+Δ

        (17)

        式中,L是n維的觀測向量;B表示n×3的系數(shù)矩陣;dX是海底應答器三維未知坐標改正數(shù);Δ是隨機誤差向量。

        將觀測方程展開可得

        通過最小二乘求解海底應答器坐標改正數(shù)

        dX=(BTPB)-1BTPL

        (18)

        觀測值改正數(shù)

        V=BdX-L

        (19)

        式中,P為觀測值的權矩陣,可基于聲線入射角確定或根據(jù)驗后殘差進行確定。

        分別計算單位權中誤差和協(xié)方差陣

        (20)

        式中,r為多余觀測數(shù);n為總觀測數(shù);t為必要觀測數(shù),一般為3;QXX為平差參數(shù)的協(xié)因數(shù)陣。

        2.2 雙程聲徑聲線跟蹤算法

        顧及聲信號的實際往返傳播路徑,基于分層常梯度聲線跟蹤算法進行雙程聲徑聲線跟蹤,如圖4所示,結合聲信號往返時間,計算聲信號往返的幾何距離,實現(xiàn)海底應答器位置的準確解算。

        圖4 雙程聲徑聲線跟蹤算法Fig.4 Acoustic ray-tracing algorithm of round-trip acoustic path

        如圖4所示,在進行雙程聲徑聲線跟蹤前,需要已知聲速剖面,根據(jù)聲速剖面計算各層的聲速梯度gi,由于聲信號往返時間較短,故認為聲速剖面結構未發(fā)生變化,已知聲信號自上向下的傳播時間和自下向上傳播的傳播時間之和等于總傳播時間,采用相同的分層結果對往返的聲線進行跟蹤,具體過程可以分為以下3步。

        第1步:聲信號自上向下傳播。設在聲信號發(fā)射時刻的聲速剖面的第i層內(nèi),聲線曲率處處相等,聲線在層內(nèi)的軌跡為一段圓弧,該圓弧的曲率半徑為Rsi。

        由Snell定律和常梯度聲速假設,聲信號在第i層內(nèi)傳播的水平距離,豎直距離和時間為

        (21)

        第2步:聲信號自下向上傳播。設在聲信號接收時刻的聲速剖面的第i層內(nèi),聲線曲率處處相等,聲線在層內(nèi)的軌跡為一段圓弧,該圓弧的曲率半徑為Rri,由于聲線自下向上進行追蹤,故每層的聲速梯度為-gi。

        由Snell定律和常聲速梯度假設,聲信號在第i層內(nèi)傳播的水平距離,豎直距離和時間為

        (22)

        L=B′dX+Δ′

        (23)

        通過最小二乘求解海底應答器坐標改正數(shù)

        dX=(B′TPB′)-1B′TPL′

        (24)

        得改正后的海底應答器坐標

        Xnew=Xo+dX

        (25)

        將改正后的海底應答器坐標Xnew代入第1步,再次進行雙程聲徑聲線跟蹤,并進行收斂判斷,直至滿足收斂閾值要求,則雙程聲線跟蹤完成。利用雙程聲徑聲線跟蹤建立傳播時間與往返距離的關系,利用雙程聲線跟蹤獲得的聲信號往返幾何距離構建雙程聲徑觀測方程,通過式(24)解算海底應答器的坐標改正數(shù),完成雙程聲徑水下高精度定位。

        與單程聲徑聲線跟蹤相比,雙程聲徑聲線跟蹤顧及了測量船在走航中換能器的位置變化,通過將聲線跟蹤過程分解為自上向下的去程和自下向上的返程,并結合聲學觀測時間,實現(xiàn)聲信號往返雙程聲徑的準確跟蹤,該算法相比單程聲徑聲線跟蹤定位更符合水下聲信號傳播實際。

        3 試驗與分析

        2019年7月,筆者所在團隊在我國南海海域開展了海洋大地測量基準與海洋導航新技術深海綜合試驗,試驗區(qū)域平均水深約3000 m。測量船搭載GNSS接收機、高精度姿態(tài)傳感器、聲速剖面儀和海面長基線定位系統(tǒng)等設備。本文中的往返時間是已知觀測值,通過測量船的相關儀器設備觀測得到。為了減小聲速剖面時空差異對定位的影響,試驗采用走航觀測前后的聲速剖面平均值進行定位解算。分別采用圓走航觀測和交叉十字觀測。

        首先,測量船以圓形航跡對海底基準點進行了連續(xù)觀測,圓走航半徑約為1.5倍水深,測量船航速設定約為1.00 m/s,即1.95 kn,聲信號往返過程中測量船平均水平位移為7.85 m,觀測歷元共1080個,走航航跡及應答器位置示意圖如圖5所示。在測區(qū)內(nèi)進行聲速測量,繪制聲速剖面如圖6所示。經(jīng)過粗差剔除、坐標系統(tǒng)轉換及姿態(tài)改正等數(shù)據(jù)預處理后,分別采用單程聲徑聲線跟蹤和雙程聲徑聲線跟蹤方法進行解算,觀測值權矩陣基于聲線入射角確定,圓走航觀測下個各歷元的聲線入射角如圖7所示,海底應答器解算坐標及其精度統(tǒng)計見表1,精度指標對應于協(xié)方差陣Dxx的對角線上的元素。

        圖5 圓走航觀測航跡與海底應答器位置 Fig.5 Circularsurveying track and the seafloor transponder

        圖6 聲速剖面 圖7 圓走航聲線入射角Fig.6 Sound velocityFig.7 Incidence beam angles profile of circular surveying

        表1 圓走航觀測下海底應答器定位結果統(tǒng)計Tab.1 Statistical table of positioning results for seafloor transponder by the circular surveying track m

        由表1可以看出,單程聲徑聲線跟蹤與雙程聲徑聲線跟蹤水平方向互差為0.161 m,豎直方向互差0.034 m,即在圓走航觀測模式下,換能器位置變化對海底應答器定位的影響主要體現(xiàn)在水平方向。比較單程聲徑聲線跟蹤和雙程聲徑聲線跟蹤的定位內(nèi)符合精度,可以發(fā)現(xiàn)雙程聲徑無論在水平方向還是豎直方向的內(nèi)符合精度均明顯優(yōu)于單程聲徑。

        圓走航觀測模式下單程聲徑聲線跟蹤和雙程聲徑聲線跟蹤的殘差值,如圖8所示。可以看出兩種方法的殘差分布較為接近,波動范圍為-1.500 m~1.500 m,均呈現(xiàn)出正態(tài)隨機分布的特點。分析圓走航觀測殘差的統(tǒng)計特征,包括平均值(mean)、最大值(max)、標準差(std)、均方根(RMS)見表2。

        圖8 圓走航觀測殘差分布直方圖Fig.8 Histogram of residual distribution in circular surveying track

        表2 圓走航觀測殘差統(tǒng)計結果Tab.2 The residual statistical results of circular surveying track m

        由表2可以看出,雙程聲徑聲線跟蹤在平均值、最大值、標準差、均方根上均略優(yōu)于單程聲徑聲線跟蹤。由圖7可以看出,圓走航觀測模式下,聲線入射角波動較小,單程聲徑和雙程聲徑的變化基本可以忽略??紤]到圓走航觀測結構對稱,在平差過程中部分系統(tǒng)誤差被消除,故雙程聲徑聲線跟蹤的優(yōu)勢體現(xiàn)不明顯,需要進一步考慮交叉十字觀測結構對定位的影響。

        測量船沿交叉十字航跡對多個海底應答器進行連續(xù)、同步觀測,測量船航跡和海底應答器位置如圖9所示,測量船的航速設定約為1.83 m/s,即3.56 kn,聲信號往返過程中測量船平均水平位移為11.63 m。5個應答器對應觀測歷元分別為1795、1827、1821、1838和1735個。以5號應答器為例,根據(jù)其概略坐標計算不同觀測歷元下該應答器對應的聲線入射角變化,如圖10所示。不同觀測歷元下中心點對應的聲線入射角最小0.65°,最大45.29°,平均聲線入射角29.26°。可以發(fā)現(xiàn),當測量船沿交叉十字航跡進行測量時,聲線入射角波動較大,觀測值權矩陣基于聲線入射角確定。圖11為十字交叉觀測5號應答殘差分布。

        圖9 交叉十字觀測航跡與海底應答器Fig.9 Coss tracking lines and the seafloor transponders

        圖10 各歷元下5號應答器對應的聲線入射角Fig.10 Incidence beam angles of No.5 transponder under each epoch

        圖11 十字交叉觀測5號應答器殘差分布直方圖Fig.11 Histogram of residual distribution with No.5 transponder in cross tracking lines

        經(jīng)過粗差剔除、坐標系統(tǒng)轉換及姿態(tài)改正等數(shù)據(jù)預處理后,分別采用單程聲徑聲線跟蹤和雙程聲徑聲線跟蹤方法進行解算,各應答器解算坐標及其精度統(tǒng)計見表3、表4。

        表3 交叉十字觀測定位解算結果統(tǒng)計Tab.3 Statistical table of cross tracks location m

        對表3和表4試驗結果分析:

        表4 十字交叉觀測5號應答器殘差統(tǒng)計結果Tab.4 The residual statistical results of No.5 transponder in cross tracking lines m

        (1) 從海底應答器坐標解算結果分析,單程聲徑聲線跟蹤與雙程聲徑聲線跟蹤水平方向互差最大為0.285 m,最小為0.166 m;豎直方向互差最大為0.145 m,最小為0.057 m。不難發(fā)現(xiàn)在交叉十字觀測模式下,換能器位置變化對海底應答器定位的影響依舊主要體現(xiàn)在水平方向,且較圓走航兩種方法互差進一步加大。對比5號應答器在圓走航和十字交叉觀測下的互差,單程聲徑聲線跟蹤在兩種不同觀測航跡下水平方向互差為0.372 m,豎直方向互差為0.533 m;雙程聲徑聲線跟蹤在兩種不同觀測航跡下水平方向互差為0.306 m,豎直方向互差為0.335 m??梢钥闯鲭p程聲徑聲線跟蹤在不同觀測航跡下的互差明顯小于單程聲徑聲線跟蹤,即該方法具有較高的定位精度。

        (2) 從定位精度分析,各海底應答器雙程聲徑聲線跟蹤的內(nèi)符合精度均優(yōu)于單程聲徑聲線跟蹤,且非中心位置的4個海底應答器坐標的內(nèi)符合精度基本一致。對于中心位置的應答器,雙程聲徑內(nèi)符合精度較單程聲徑提高了近一半。結果表明,雙程聲線跟蹤可以有效提高水下定位精度,尤其在非對稱觀測條件下,效果更加明顯。對比不同觀測航跡的內(nèi)符合精度,可以發(fā)現(xiàn),圓走航模式下的內(nèi)符合精度優(yōu)于十字交叉走航的內(nèi)符合精度,證明了觀測航跡的結構對于定位的精度有較為明顯的影響。

        (3) 為進一步評價定位結果的可靠性,以5號海底應答器為例,采用不同觀測方式下5號海底應答器的計算結果均值為參考。由于海底應答器坐標的真值無法獲取,通常是采用多次觀測的平均值作為最或然值,即不同航跡下單程、雙程聲線跟蹤坐標解的平均值。因為不同觀測方案是在不同日期進行的,互為獨立觀測,可以進行相互檢核。分析得到單程聲徑聲線跟蹤與參考值水平方向差為0.281 m,豎直方向差為0.299 m,雙程聲徑聲線跟蹤與參考值水平方向差為0.141 m,豎直方向差為0.135 m。從殘差平均值來看,單程聲徑為0.219 m,雙程聲徑為0.154 m;從殘差均方根來看,單程聲徑為2.561 m,雙程聲徑為1.038 m。可以看出,當采用雙程聲徑聲線跟蹤時,殘差均方根等明顯變小,較單程聲徑聲線跟蹤有更好的定位效果。如果進一步優(yōu)化函數(shù)模型,改化算法,有望獲得更優(yōu)的坐標精度,優(yōu)化處理模型和改化算法涉及面較多,將另文討論。

        4 結 論

        在水下高精度定位中,聲速誤差是影響聲學定位精度的重要原因,采用聲線跟蹤的方法可以有效地削弱聲線折射效應的影響,但基于單程聲徑水下定位模型,沒有考慮聲信號往返過程中船載換能器的位置變化,導致單程聲徑聲線跟蹤算法存在建模誤差,降低了水下聲學定位精度。本文算法,考慮了水下定位過程中聲信號傳播的實際情況,利用測量船持續(xù)走航引起的換能器位置變化信息,通過建立雙程聲徑水下定位模型,實現(xiàn)水下聲學高精度定位,得到了以下結論:

        (1) 在圓走航觀測情況下,測量船沿圓航跡行進,換能器與應答器之間的往返聲徑變化較小,雙程聲徑聲線跟蹤對海底應答器定位精度有一定的提升,主要體現(xiàn)在水平方向。

        (2) 在交叉十字觀測情況下,測量船沿直線航跡行進,換能器與應答器之間的往返聲徑不斷變化,需要考慮換能器位置變化對海底應答器定位的影響。通過對雙程聲徑和單程聲徑兩種解算方法的海底應答器坐標解殘差的對比,可以發(fā)現(xiàn)雙程聲徑聲線跟蹤能有效提高海底應答器的定位精度,尤其在非對稱觀測條件下,改善效果更加明顯。

        本文在解算過程中忽略了系統(tǒng)誤差的影響,下一步將針對水下聲學定位的系統(tǒng)誤差,考慮采用半?yún)?shù)法等最優(yōu)參數(shù)估計方法,研究系統(tǒng)誤差處理策略,完善雙程聲徑聲線跟蹤算法,以期實現(xiàn)更高精度的水下聲學定位。

        猜你喜歡
        信號
        信號
        鴨綠江(2021年35期)2021-04-19 12:24:18
        完形填空二則
        7個信號,警惕寶寶要感冒
        媽媽寶寶(2019年10期)2019-10-26 02:45:34
        孩子停止長個的信號
        《鐵道通信信號》訂閱單
        基于FPGA的多功能信號發(fā)生器的設計
        電子制作(2018年11期)2018-08-04 03:25:42
        基于Arduino的聯(lián)鎖信號控制接口研究
        《鐵道通信信號》訂閱單
        基于LabVIEW的力加載信號采集與PID控制
        Kisspeptin/GPR54信號通路促使性早熟形成的作用觀察
        久久久国产精品va麻豆| 日日骚一区二区三区中文字幕| 国产熟妇一区二区三区网站| 91精品国产综合久久久密臀九色| 久久青青草原精品国产app| 免费无码成人av在线播放不卡| www.五月激情| 96中文字幕一区二区| 脱了老师内裤猛烈进入| 好大好深好猛好爽视频免费| 在线免费欧美| 青青草是针对华人绿色超碰| 三级全黄裸体| 免费人成无码大片在线观看| 久久亚洲成a人片| 少妇精品偷拍高潮少妇在线观看| 欧美a级在线现免费观看| 久久久久久国产精品美女| 精品国产高清a毛片| 青青草手机视频免费在线播放| 国内免费高清在线观看| 亚洲精品自产拍在线观看| 亚洲成a人片在线观看中| 国产精品综合一区久久| 小蜜被两老头吸奶头在线观看| 精品十八禁免费观看| 狼人狠狠干首页综合网| 2021国产精品视频网站| 国产专区国产av| 国模少妇无码一区二区三区 | av中国av一区二区三区av | 国产又黄又猛又粗又爽的a片动漫| 久99久精品免费视频热77| 在线观看在线观看一区二区三区| 亚洲色大成网站www永久网站| 日韩AV不卡一区二区三区无码| 久久伊人精品只有这里有| 久久伊人这里都是精品| 人妻少妇av中文字幕乱码| 久久99精品这里精品动漫6| 国产片在线一区二区三区|