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

        ?

        航空重力測量厄特弗斯改正公式注記

        2015-01-11 02:10:24黃謨濤寧津生歐陽永忠陸秀平翟國君鄧凱亮吳太旗
        測繪學報 2015年1期
        關鍵詞:弗斯橢球重力

        黃謨濤,寧津生,歐陽永忠,,陸秀平,翟國君,鄧凱亮,吳太旗

        1.海軍海洋測繪研究所,天津300061;2.武漢大學測繪學院,湖北 武漢430079

        1 引 言

        厄特弗斯(E?tv?s)改正是航空重力測量最重要的改正項之一。它是由于測量載體相對地球運動,使重力測量儀器受到附加的離心力作用而增加的一項改正,是車載、船載、機載等運動模式下的重力測量都必須顧及的共同改正項[1-10]。由于車載和船載測量速度較低,一般使用球近似公式計算厄特弗斯改正就能滿足實用要求[4,10]。因此,對于車載和船載重力測量,厄特弗斯改正計算公式在使用上并不存在異議。但在航空重力測量中,由于飛機速度可達300~500km/h甚至更高,厄特弗斯改正計算公式的近似程度將直接影響到重力測量成果的質量,因此必須對此問題給予高度重視。由于歷史上多種主客觀的原因,在使用航空重力測量厄特弗斯改正公式問題上國內外至今仍未取得完全一致,近期出版的著作和發(fā)表的論文甚至作業(yè)標準仍在引用早期的近似公式[11-21]。我國航空重力測量作業(yè)部門在使用厄特弗斯改正公式時也出現比較明顯的引用錯誤,在一定程度上影響了航空重力測量成果的質量。針對此情況,本文全面分析比較了不同時期、不同形式厄特弗斯改正公式在數值上的差異,指出了我國目前使用近似公式存在的問題和統(tǒng)一使用嚴密公式的必要性,供實際作業(yè)部門參考。

        2 不同形式厄特弗斯改正公式來源與分析

        2.1 Thompson(1960)公式

        在運動載體上實施重力測量會產生厄特弗斯效應,最早是由匈牙利科學家E?tv?s于1919年發(fā)現并加以推證的。E?tv?s當年給出的球近似改正公式在船載重力測量中一直沿用至今,其公式形式為[4]

        式中,δgE為厄特弗斯改正數;ω為地球自轉角速度;V為載體運動速度;α為載體運動方位角;φ為測點大地緯度;R為地球平均半徑。1958年以后,國際上陸續(xù)開展了航空重力測量試驗,文獻[5]最先推出針對航空重力測量的厄特弗斯改正公式,其形式為

        式中,Rφ為地球橢球在緯度φ處的向徑;Vφ為由地球自轉引起的地球表面上緯度φ處的切線速度;V為飛機速度在地球表面的投影;Ve為V的東向分量;Vn為V的北向分量;h代表飛機相對平均海面高度。將Vφ=Rφωcosφ和Ve=Vsinα代入式(2)得不難看出,Thompson公式仍為球近似公式,

        與傳統(tǒng)使用的海面船測重力改正公式(1)相比,其前面增加了乘系數(1+h/Rφ),它的意義相當于把飛機速度在地球表面上的投影值V轉換為飛行高度上的速度。在Thompson公式發(fā)表后的一個時期內,國際上諸多文獻都一致將該公式作為標準公式使用[6-9]。我國早期的文獻[1—2,19],包括近期出版的文獻[20],也都一直引用該公式討論航空重力測量數據歸算問題。這里需要指出的是,Thompson公式中的V值代表的是飛機速度在地球表面的投影,即地面(或地表)速度(ground speed)[5],而不是通常理解的飛行高度上的速度Vh。我國相關引用文獻幾乎都忽略了以上兩種速度之間的差異[1-2,19-20],有的文獻將地面速度理解為飛機相對地球的速度[2]。實際上,在球近似下,V和Vh之間存在如下簡單的聯(lián)系

        按照原Thompson公式的推導思路,不難推出使用飛行高度上的速度Vh表示的另一形式的Thompson公式

        由于飛行高度h相對Rφ是個很小的量,因此V和Vh之間的差異也不大,但在高精度測量要求中必須對它們加以區(qū)別。至于Thompson當年為什么要用地面速度而不直接用飛行高空速度,是由當時的測速技術條件所決定的,因為當時比較成熟的測速手段是地面攝影(ground photograph)和多普勒雷達(Doppler radar)技術[22],兩種技術都能直接給出飛機的地面速度。

        2.2 Harlan(1968)公式

        針對Thompson公式球近似可能帶來不可忽視的誤差影響,1968年,Harlan從幾何學出發(fā),推出了兩組分別對應于地表速度和飛行高空速度、顧及地球橢球扁率的厄特弗斯改正公式[11],其形式為:

        當V代表地表速度時,有

        當V代表飛行高空速度時,有

        式(6)—(9)中,a代表地球橢球長半軸;f為橢球扁率;其他符號意義同前。

        Harlan公式發(fā)表后,在國內外學術界得到了普遍認可并在實踐中得到了廣泛應用[3,11-18,21]。但在這些大量的引用文獻中,不乏有由于作者可能的疏忽而造成的引用錯誤。其中,文獻[15]在引用時將式(7)右端第2項h/a前的“+”誤寫為“-”號;文獻[16]引用的公式為

        文中同時注明公式中的V值代表地表速度,且

        顯然,式(10)存在多處錯誤:①右端第1項分母中的括號部分應該位于分子;②公式中的V值代表的應該是飛行高空速度;③f代表的應該是橢球扁率而不是式(11)所示的表達式。另外,文獻[17]在引用式(9)時,也誤將式中的f值理解為

        令人更為遺憾的是,我國學者和作業(yè)規(guī)程在引用Harlan公式時也出現了理解上的偏差,文獻[3]和文獻[21]的引用公式為

        顯然,式(13)只是式(7)的簡單變形體,兩者應當是等價的,即式(13)中出現的地心向徑r應該是a。上述兩個文獻之所以都把a寫成r,可能是受到文獻[11]中譯本出現排版錯誤的誤導[23]。另外,式(13)所對應的V值代表的應該是地表速度,而文獻[3]和[21]都誤將其標明為飛行高空速度。美國Micro-g LaCoste生產的TAGS型航空重力儀在數據處理軟件中,也同樣使用式(7)作為航空重力測量厄特弗斯改正計算模型[24],筆者經比對分析發(fā)現,軟件中同樣存在直接使用飛行高空速度替代地表速度的問題。

        上述引用錯誤不僅在理論上是不嚴密的,在實際應用中也會引起較大的數值偏差。因此,應當及時加以糾正。

        2.3 嚴密計算公式

        由前面的分析得知,最早發(fā)表的Thompson公式只是一個球近似公式,后來發(fā)表的Harlan公式也只是一個顧及橢球扁率一階項影響的近似公式。實際上,前面提及的一些重要文獻在推演航空重力測量數學模型過程中,已經給出了厄特弗斯改正的嚴密計算公式[3,14,17,25],其形式為

        式中,V代表飛行高空速度;N和M分別代表地球橢球卯酉圈和子午圈曲率半徑,其計算式分別為

        式中,e為橢球第一偏心率,e2=2f-f2;其他符號含義同前。

        式(14)的具體推導過程見文獻[3]。由式(14)知,當近似取N≈M≈R(即球近似)時,式(14)就轉變?yōu)榍懊娴?Thompson公式,即式(5),如果近似取

        將其代入式(14),作冪級數展開并取至橢球扁率一階項,則不難得到以飛行高空速度表示的Harlan公式,即式(8)和式(9),其推導過程從略。

        假如用Vg表示飛行高空速度V在地球表面的投影,Vge和Vgn分別表示Ve和Vn的投影,根據關系式

        將其代入式(14),則可得到以地表速度表示的厄特弗斯改正嚴密計算式

        同樣,如果將近似式(17)和式(18)代入式(21),則可得到以地表速度表示、顧及橢球扁率一階項的Harlan公式,即式(6)和式(7)。

        至此,筆者已經理清了Thompson公式、Harlan公式與嚴密計算公式之間的關系,同時給出了從嚴密公式推導Harlan公式的技術途徑,其結果與Harlan從幾何學角度出發(fā)推導結果是等價的[22]。但這里必須指出的是,前面所指的嚴密公式其實也只是相對地球視為旋轉橢球而言的,即統(tǒng)一把橢球法線作為重力加速度的投影線處理。而實際上重力傳感器感應的應該是重力加速度在重力垂線方向上的投影,法線和垂線的差異即垂線偏差對厄特弗斯改正的影響大小,取決于測量區(qū)域重力場的復雜程度。文獻[22]認為,由于垂線偏差量值一般不超過30″,因此該項影響可以忽略不計。文獻[26]則通過估算認為,當垂線偏差變化率超過1″/km時,該項影響可達幾個mGal(1mGal=1×10-5m/s2)甚至幾十個mGal。由于篇幅所限,本文暫不涉及此問題的討論。

        在實際應用中,除了部分學者仍習慣于引用Harlan公式外,個別學者和機構在引用嚴密公式(14)時,往往又作了近似處理。文獻[27]和[28]把N和M都近似取為地球平均半徑R,即回歸到簡單的Thompson公式;文獻[29]則是把(N+h)和(M+h)都近似取為地面計算點的地心向徑r,即使用式(22)和式(23)

        為了估算曲率半徑近似誤差對厄特弗斯改正計算結果的影響,取厄特弗斯改正球近似公式的第2項進行討論

        對上式求微分并寫成中誤差形式得

        式中,mR代表球半徑誤差;mΔE代表由mR引起的厄特弗斯改正誤差。取R=6371km;V=500km/h,可得到表1的估算結果。

        表1 曲率半徑誤差對厄特弗斯改正計算結果的影響Tab.1 Error effect of curvature radius to the E?tv?s corrections

        由表1計算結果可以看出,在高精度測量要求條件下,大于5km的曲率半徑誤差對厄特弗斯改正計算結果的影響都是不可忽略的。而由橢球大地測量學得知,從赤道到兩極,卯酉圈曲率半徑N的變化幅度接近22km,子午圈曲率半徑M的變化幅度達64km,地球橢球平均半徑R與N的互差范圍為7~28km,R與M的互差范圍為-36~28km。由此可見,任何對曲率半徑作近似處理的做法都是不可取的。

        3 不同形式厄特弗斯改正公式數值差異分析

        3.1 模擬數值比對分析

        為了準確掌握不同形式厄特弗斯改正公式在數值上的差異大小,這里選擇嚴密公式(14)作為基準,分別計算 Thompson公式(3)和式(5)、Harlan公式(7)和式(9)、我國作業(yè)規(guī)范采用的公式(13)及俄羅斯生產設備GT-1A采用的公式(22)與基準值的差異。其中,為了說明引用錯誤可能帶來的誤差大小,各計算式中的V統(tǒng)一直接采用飛行高度上的速度,并取飛行高度h=3000m,V=500km/h;同時采用2000國家大地坐標系(CGCS2000)定義的地球橢球參數[30],即

        不同緯度測點相對于不同航向角上的計算比對結果分列于表2—表4。

        表2 φ=0°時的計算比對結果Tab.2 Comparison results from different formulas onφ=0° mGal

        表3 φ=45°時的計算比對結果Tab.3 Comparison results from different formulas onφ=45° mGal

        表4 φ=89°時的計算比對結果Tab.4 Comparison results from different formulas onφ=89° mGal

        由前面的分析得知,式(5)、式(9)和式(22)與嚴密公式(14)的差異主要源于曲率半徑的近似處理。式(3)和式(5)、式(7)和式(9)原本是等價的,但在實際引用中,往往出現飛行高空速度和地表速度混淆使用的情況,由此引起了額外的誤差。由表2—表4的計算結果可以看出,基于球近似的 Thompson 公 式(3)和 式(5),誤 差 量 值 接近2mGal;基于曲率半徑顧及橢球扁率一階項影響的 Harlan公式(7)和式(9),誤差量值不超過0.1mGal,但當出現飛行高空速度和地表速度混淆使用時,其帶來的計算誤差可達1 mGal;式(13)沒有嚴格區(qū)分飛行高空速度和地表速度,同時誤用地心向徑r代替地球橢球長半軸a,由此引起的計算誤差超過1mGal;式(22)忽略飛行高度同時對曲率半徑作近似處理,其誤差量值超過2mGal。

        對于區(qū)域性的航空重力測量,由于在單一測線上的飛行速度、航向和航高都基本保持穩(wěn)定,測點緯度變化范圍也相當有限。因此,以上由于各種原因引起的計算偏差一般都具有系統(tǒng)性誤差特征,它們對測量成果質量的影響不可忽視。此外,目前國際上航空重力測量的精度已達1~3mGal[14,18,27],如果不顧及以上計算偏差的影響,勢必限制測量成果精度的進一步改善。因為在測量學科領域,通常要求理論模型的計算精度應該比儀器觀測量精度高出1個數量級甚至更高。表2—表4的計算結果對應于飛行速度V=500km/h,如果飛行速度降低到300km/h甚至更小,則表中所列各類誤差都將相應減小。

        3.2 實測數據比對分析

        2012年,筆者所在單位在南海某海域開展了多種形式的航空重力測量飛行試驗,其中包括使用多臺套重力儀在東—西、南—北、西北—東南等特定方向上開展重復測線飛行試驗。測線長度東—西向約240km,南—北向約200km,西北—東南向約350km,測量期間飛行高度約1500m,飛行速度約400km/h。對重復測線試驗數據進行分析比對后發(fā)現,所有同向飛行的重復測線觀測數據內符合精度都很高,與儀器廠家標稱的測量精度基本相符。但對向(或稱正反向)飛行數據的比對情況有所不同,除南—北向外,東—西向和西北—東南向重復測線的對向觀測數據比對結果都顯示,兩者之間存在比較明顯的系統(tǒng)性偏差。經進一步分析研究后確認,除了某些環(huán)節(jié)上的測量環(huán)境效應改正還不夠完善外,厄特弗斯改正公式使用不當也是造成出現上述現象的主要原因。為了突出說明后者的影響,這里單獨列出兩型儀器(美國公司生產的TAGS型和LCRⅡ型)分別使用Harlan公式(7)和我國作業(yè)規(guī)范公式(13),計算厄特弗斯改正的對向飛行數據比對結果(統(tǒng)一采用飛行高度上的速度)。其中,每一方向都給出3組計算結果,分別對應于正向和反向飛行使用近似厄特弗斯改正公式帶來的誤差統(tǒng)計,以及正向與反向互差帶來的累積誤差統(tǒng)計,具體見表5和表6。

        表5 使用公式(7)時的實測數據比對結果Tab.5 Comparisons from practical measurements when using formulas(7)mGal

        從表5和表6的計算結果可以看出,盡管這里的試驗飛行速度和高度相對于前面的模擬數據都有所降低,但無論是使用公式(7)還是公式(13),由于公式使用不當引起的系統(tǒng)性偏差仍達0.5mGal,重復測線正反向互差最大可達0.8 mGal。顯然,與儀器廠家標稱的測量精度(優(yōu)于1 mGal)相比較[24],以上誤差影響量值大小都是不可忽略的。

        表6 使用公式(13)時的實測數據比對結果Tab.6 Comparisons from practical measurements when using formulas(13)mGal

        由于南—北向測線正向與反向計算誤差符號相同,同時量值較小,求互差時有相互抵消作用,因此在此方向上由公式使用不當引起的計算偏差幾乎可以忽略不計。

        文獻[3]介紹了使用LCRⅠ型重力儀在我國山西大同地區(qū)的試驗情況,其中測量飛行平均高度為3400m,飛行速度為360km/h。就該試驗而言,如果使用公式(7)計算厄特弗斯改正數,那么東—西方向測線正向飛行的計算偏差為-0.76mGal,反向飛行的計算偏差為0.43mGal,重復測線互差值的計算偏差為-1.19mGal;如果使用作業(yè)規(guī)范公式(13)進行計算,那么東—西方向測線正向飛行的計算偏差為-0.89mGal,反向飛行的計算偏差為0.30mGal,重復測線互差值的計算偏差為-1.19mGal。以上誤差量值比表5和表6的比對結果還要大一些。

        4 結論與建議

        綜合前面的分析、推演和比對結果可以看出:

        (1)盡管航空重力測量厄特弗斯改正的嚴密計算公式并不復雜,但國內外機構和學者在使用厄特弗斯改正公式問題上至今還缺乏統(tǒng)一。

        (2)與嚴密計算公式相比較,不同時期發(fā)表的不同形式的厄特弗斯改正公式都存在一定大小的系統(tǒng)性偏差。雖然Harlan(1968)公式逼近度較高,但經常出現使用不當的情況,由此引起的計算偏差不可忽略。

        (3)我國航空重力測量作業(yè)規(guī)范采用的厄特弗斯改正公式存在比較明顯的引用錯誤,建議有關部門盡快對其進行更正,統(tǒng)一采用厄特弗斯改正嚴密計算公式。

        [1] FANG Jun.Gravimetry and Figure of the Earth[M].Beijing:Science Press,1965.(方俊.重力測量與地球形狀學[M].北京:科學出版社,1965.)

        [2] GUAN Zelin,NING Jinsheng.Figure of the Earth and the External Gravitational Field[M].Beijing:Surveying and Mapping Press,1981.(管澤霖,寧津生.地球形狀及外部重力場[M].北京:測繪出版社,1981.)

        [3] SUN Zhongmiao.Theory,Methods and Applications of Airborne Gravimetry[D].Zhengzhou:Information Engineering University,2004.(孫中苗.航空重力測量理論、方法及應用研究[D].鄭州:信息工程大學,2004.)

        [4] HUANG Motao,ZHAI Guojun,GUAN Zheng,et al.The Determination and Application of Marine Gravity Field[M].Beijing:Surveying and Mapping Press,2005.(黃謨濤,翟國君,管錚,等.海洋重力場測定及其應用[M].北京:測繪出版社,2005.)

        [5] THOMPSON L G D,LACOSTE L J B.Aerial Gravity Measurements[J].Journal of Geophysical Research,1960,65(1):305-322.

        [6] NETTLETON L L,LACOSTE L J B,HARRISON J C.Tests of an Airborne Gravity Meter[J].Geophysics,1960,25(1):181-202.

        [7] NETTLETON L L,LACOSTE L J B,GLICKEN M.Quantitative Evaluation of Precision of Airborne Gravity Meter[J].Journal of Geophysical Research,1962,67(11):4395-4410.

        [8] HARRISON J C.The Measurement of Gravity[J].Proceedings of the IREE,1962,50(11):2302-2312.

        [9] GLICKEN M.Eotvos Corrections for a Moving Gravity Meter[J].Geophysics,1962,27(4):531-533.

        [10] DEHLINGER P.Marine Gravity[M].New York:Elservier Scientific Publishing Company,1978.

        [11] TORGE W.Gravimetry[M].Berlin:Walter de Gruyter,1989.

        [12] BELL R E,BERNARD J C,ROBERT W S.Airborne Gravimetry from a Small Twin Engine Aircraft over the Long Island Sound[J].Geophysics,1991,56(9):1486-1493.

        [13] FORSBERG R,OLESEN A V,KELLER K.Airborne Gravity Survey of the North Greenland Shelf 1998[M].Copenhagen:[s.n.],1999.

        [14] OLESEN A V.Improved Airborne Scalar Gravimetry for Regional Gravity Field Mapping and Geoid Determination[D].Copenhagen:University of Copenhagen,2003.

        [15] RIEDEL S.Airborne-based Geophysical Investigation in Dronning Maud Land,Antarctica[D].Bremen:University Bremen,2008.

        [16] NEUMEYER J,SCHAFER U,KREMER J,et al.Derivation of Gravity Anomalies from Airborne Gravimeter and IMU Recording—Validation with Regional Analytic Models U-sing Ground and Satellite Gravity Data[J].Journal of Geodynamics,2009,47:191-200.

        [17] ALBERTS B.Regional Gravity Field Modeling Using Airborne Gravimetry Data[C]∥Publications on Geodesy 70, Netherlands Geodetic Commission. Delft:[s.n.],2009.

        [18] JEKELI C.Moving-Base Gravimetry[C]∥Supplemental Notes to Gravimetic Geodesy, GS776.Columbus:[s.n.],2012.

        [19] ZHANG Shanyan.On the E?tv?s Problem for Aerial Gravity Measurement[J].Acta Geodaetica et Geophysica,1982,4:97-104.(張善言.航空重力測量的厄特弗斯問題[J].測量與地球物理集刊,1982,4:97-104.)

        [20] LüZhiping,QIAO Shubo.Foundation of Geodesy[M].Beijing:Surveying and Mapping Press,2010.(呂志平,喬書波.大地測量學基礎[M].北京:測繪出版社,2010.)

        [21] Headquarters of General Equipment.GJB 6561-2008.Rules for Operations of Airborne Gravimetry[S].Beijing:Military Standard Press of the Headquarters of General Equipment,2008.(總裝備部.GJB 6561-2008.航空重力測量作業(yè)規(guī)范[S].北京:總裝備部軍標出版發(fā)行部,2008.)

        [22] HARLAN R B.E?tv?s Corrections for Airborne Gravimetry[J].Journal of Geophysical Research,1968,73(14):4675-4679.

        [23] TORGE W.Gravimetry[M].Beijing:Earthquake Press,1993.(TORGE W.重力測量學[M].北京:地震出版社,1993.)

        [24] LIANG Xinghui,LIU Lintao,WU Pengfei,et al.Study on CHZ Gravimeter Applied in Airborne Gravimetry Involving Error Spectrum Characteristic [J]. Acta Geodaetica et Cartographica Sinica,2013,42(5):633-639.(梁星輝,柳林濤,吳鵬飛,等.顧及誤差頻譜特性的CHZ重力儀航空應用研究[J].測繪學報,2013,42(5):633-639.)

        [25] WALL R E.Airborne Gravimetry Erros Associated with Geoidal Undulations[J].Journal Geophysical Research,1971,76(29):7293-7295.

        [26] HWANG C,HSIAO Y S,SHIH H C,et al.Geodetic and Geophysical Results from a Taiwan Airborne Gravity Survey:Data Reduction and Accuracy Assessment[J].Journal Geophysical Research,2007,112(B10),doi:10.1029/2005JB004220.

        [27] SHIH H C.Multiple-Altitude Airborne Gravity Surveys and Applications to Geoid Determination and Kuroshio Current[D].Taiwan:National Chiao Tung University,2010.

        [28] Joint-Stock Company.Gravimeter Training Abstract for Airborne Gravimeter Model GT-1A[C]∥Proceedings of Scientific and Technological Enterprise.Moscow:[s.n.],2008.

        [29] Headquarters of General Equipment.GJB 6304-2008.China Geodetic System 2000[S].Beijing:Military Standard Press of the Headquarters of General Equipment,2008.(總裝備部.GJB 6304-2008.2000中國大地測量系統(tǒng)[S].北京:總裝備部軍標出版發(fā)行部,2008.)

        猜你喜歡
        弗斯橢球重力
        瘋狂過山車——重力是什么
        科學大眾(2022年23期)2023-01-30 07:04:16
        獨立坐標系橢球變換與坐標換算
        橢球槽宏程序編制及其Vericut仿真
        智能制造(2021年4期)2021-11-04 08:54:44
        總得有人去擦星星
        愛心樹(上)
        西西弗斯的神話
        航空世界(2018年12期)2018-07-16 08:34:32
        閣樓上的光
        橢球精加工軌跡及程序設計
        基于外定界橢球集員估計的純方位目標跟蹤
        仰斜式重力擋土墻穩(wěn)定計算復核
        中文字幕国产精品专区| 日本乱子人伦在线视频| 亚洲成a人片在线网站| 亚洲av天堂久久精品| 精品国产精品三级在线专区| 猫咪av成人永久网站在线观看| 丰满爆乳无码一区二区三区| 精品一区二区三区在线观看l| 成人自拍偷拍视频在线观看| 日韩夜夜高潮夜夜爽无码| 国产精品成人一区二区三区| 欧美日韩激情在线一区二区| 国产日本精品一区二区免费 | 色偷偷久久久精品亚洲| 麻豆久久久9性大片| 91精品全国免费观看青青| 日韩一区中文字幕在线| 国产精品久久久久久妇女| 男男车车的车车网站w98免费| 免费黄色福利| 少妇被猛烈进入中文字幕| 日本真人做爰免费视频120秒| 中文字幕亚洲无线码| 久久亚洲一级av一片| 就爱射视频在线视频在线| 亚洲第一页综合图片自拍| 中文字幕少妇AV| 亚洲精品av一区二区日韩| 国产成人无码a在线观看不卡| 国产乱人伦av在线无码| 国产精品女丝袜白丝袜 | 成人无码av免费网站| 欧美巨大xxxx做受l| 国产极品美女到高潮视频| 中文字幕av永久免费在线| 狼人香蕉香蕉在线28 - 百度| 精品免费人伦一区二区三区蜜桃| 国产少妇露脸精品自拍网站| 成人麻豆日韩在无码视频| 一本色道久久综合狠狠躁| 永久免费在线观看蜜桃视频|