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

        ?

        利用矩諧分析的局部大地水準(zhǔn)面確定

        2018-08-31 06:12:50路媛琦章傳銀
        測(cè)繪通報(bào) 2018年8期
        關(guān)鍵詞:水準(zhǔn)面階數(shù)重力

        路媛琦,章傳銀,蔣 濤

        (1. 青島市勘察測(cè)繪研究院,山東 青島 266032; 2. 中國(guó)測(cè)繪科學(xué)研究院,北京 100039; 3. 山東科技大學(xué),山東 青島 266590)

        大地水準(zhǔn)面是一個(gè)與靜止的平均海水面重合并延伸到大陸內(nèi)部的、包圍整個(gè)地球的、封閉的重力位水準(zhǔn)面,可以作為定義全球或區(qū)域高程系統(tǒng)的基準(zhǔn)面[1]。確定大地水準(zhǔn)面是建立高程基準(zhǔn)統(tǒng)一的有效途徑,并可結(jié)合GNSS實(shí)現(xiàn)大地高到正高的轉(zhuǎn)換,實(shí)現(xiàn)衛(wèi)星三維空間定位。

        全球大地水準(zhǔn)面的確定通常采用球諧展開(kāi)模型,但確定局部大地水準(zhǔn)面所利用的局部重力數(shù)據(jù)無(wú)法滿足球諧函數(shù)的正交性,球諧展開(kāi)模型并不適用于表達(dá)區(qū)域大地水準(zhǔn)面。區(qū)域大地水準(zhǔn)面確定的常用方法一種是基于Stokes或Molodensky理論求解大地測(cè)量邊值問(wèn)題,我國(guó)第二代似大地水準(zhǔn)面模型(CQG2000)即是采用此類(lèi)傳統(tǒng)方法構(gòu)建。另一種被許多學(xué)者應(yīng)用的解析法為球冠諧分析;李建成等首次提出將球冠諧分析應(yīng)用于區(qū)域重力場(chǎng)建模,證明了球冠諧展開(kāi)表達(dá)區(qū)域重力場(chǎng)的可行性。球冠諧分析在大區(qū)域范圍內(nèi)建立大地水準(zhǔn)面是一種理想方法,但應(yīng)用于較小范圍的大地水準(zhǔn)面確定時(shí),存在數(shù)值計(jì)算上的困難,精度較差。

        確定局部大地水準(zhǔn)面的譜方法除球冠諧分析外,還有一種方法是矩諧分析。矩諧分析在地磁學(xué)研究中應(yīng)用廣泛;Alldredge最早使用矩諧分析用于區(qū)域地磁場(chǎng)建模。1992年邊少鋒[1]將矩諧分析引入地球重力場(chǎng)逼近,由于矩諧分析的基函數(shù)是三角函數(shù),球冠諧分析的基函數(shù)是非整階締合Legendre函數(shù),展開(kāi)到較高階次時(shí)前者的計(jì)算遠(yuǎn)比后者穩(wěn)定、快捷。

        本文利用EGM-2008重力位模型,加入標(biāo)準(zhǔn)差為2 mGal的高斯白噪聲,模擬陸地重力及飛行高度為4 km的航空重力值;并在研究區(qū)域內(nèi)建立直角坐標(biāo)系,求解地球引力位的Laplace方程,得到由矩諧系數(shù)表達(dá)的重力場(chǎng)參數(shù)表達(dá)式;將模擬重力數(shù)值代入矩諧展開(kāi)式中求得矩諧系數(shù)、殘余重力擾動(dòng)及殘余大地水準(zhǔn)面;最后恢復(fù)大地水準(zhǔn)面并與真實(shí)值進(jìn)行對(duì)比,以此說(shuō)明矩諧分析確定大地水準(zhǔn)面的可靠性、準(zhǔn)確性。

        1 矩諧分析

        已知地球外部空間某一區(qū)域或地球表面的重力觀測(cè)值(重力異?;蛑亓_動(dòng)),求定該區(qū)域內(nèi)地球重力場(chǎng)的矩諧系數(shù),這一過(guò)程叫作矩諧分析[2]。

        1.1 矩諧展開(kāi)模型

        矩諧分析(rectangular harmonic analysis,RHA)的基本原理是對(duì)地球某一地區(qū)構(gòu)建區(qū)域大小的矩形面,以此近似球面;以矩形面的中心點(diǎn)為原點(diǎn)建立局部直角坐標(biāo)系(如圖1所示);并在其中求解Laplace方程,在選定一定的截?cái)嗨胶螅?chǎng)和引力位是包含有限個(gè)待定系數(shù)的已知函數(shù),可以用所研究區(qū)域內(nèi)一組觀測(cè)點(diǎn)上的引力場(chǎng)值確定這些系數(shù),并用這些系數(shù)表示不同波長(zhǎng)的重力場(chǎng)信息,便于求解大地水準(zhǔn)面[3-6]。

        地球引力位滿足Laplace方程[6-7]

        (1)

        Laplace方程在直角坐標(biāo)系中表示為

        (2)

        采用分離變量法,引力位可寫(xiě)為

        V(x,y,z)=f(x)g(y)h(z)

        (3)

        將式(3)代入式(2)得局部直角坐標(biāo)系中Laplace方程的解[3]為

        (4)

        式中,Cnm為矩諧系數(shù),即引力位系數(shù);ψ(x),ψ(y)為基函數(shù),可分別表示為

        (5)

        其中

        (6)

        τnm定義為

        (7)

        式中,Dx、Dy分別表示矩形計(jì)算區(qū)域在東西(x)和南北(y)方向上的距離。

        圖1 矩諧分析使用的局部直角坐標(biāo)系

        實(shí)用上矩諧級(jí)數(shù)Cnm不可能展開(kāi)到無(wú)窮階數(shù),應(yīng)截?cái)嘀聊硞€(gè)最高階數(shù)N和次數(shù)M,則扣除參考橢球的正常引位之后,擾動(dòng)位可展開(kāi)為如下有限階次的矩諧級(jí)數(shù)

        (8)

        由Molodensky邊值條件得重力異常與擾動(dòng)位滿足如下近似關(guān)系

        (9)

        于是聯(lián)合式(8)和式(9)可推出重力異常的矩諧展開(kāi)式為

        (10)

        式中,r為地心距離。

        同理,重力擾動(dòng)、大地水準(zhǔn)面差距的矩諧展開(kāi)式分別為

        (11)

        (12)

        式中,γ0為正常重力。

        1.2 截?cái)嚯A數(shù)

        在實(shí)際應(yīng)用中,無(wú)論球諧級(jí)數(shù)還是矩諧級(jí)數(shù)都不可能展開(kāi)到無(wú)窮多次,因此對(duì)級(jí)數(shù)表達(dá)式選擇適當(dāng)?shù)慕財(cái)嚯A數(shù)非常重要,過(guò)低的截?cái)嚯A數(shù)會(huì)丟失有價(jià)值的信息,而過(guò)高的截?cái)嚯A數(shù)會(huì)導(dǎo)致計(jì)算量增加,從而使結(jié)果不穩(wěn)定,出現(xiàn)所謂的“龍格”現(xiàn)象[8-9]。因此需要根據(jù)數(shù)據(jù)的數(shù)量和精度要求[10],確定模型的截?cái)嚯A數(shù)Nmax。最優(yōu)截?cái)嚯A數(shù)的確定仍舊困難,通常得到最高截?cái)嚯A數(shù)后不斷進(jìn)行嘗試性的計(jì)算比較,以獲取最優(yōu)截?cái)嚯A數(shù)。

        最高截?cái)嚯A數(shù)公式[8]為

        最佳截?cái)嚯A數(shù)公式[11]為

        1.3 邊界效應(yīng)

        對(duì)于區(qū)域大地水準(zhǔn)面而言,其研究范圍的局部性限制了重力場(chǎng)信號(hào)的完全周期性;但矩諧展開(kāi)模型采用的周期性函數(shù)Fourier級(jí)數(shù)成立的前提是假定待求信號(hào)在計(jì)算區(qū)域內(nèi)為周期函數(shù)[12]。式(8)在計(jì)算區(qū)域內(nèi)收斂于真實(shí)地球擾動(dòng)位,在計(jì)算區(qū)域邊界處近似等于兩邊界擾動(dòng)位的平均值,因此在邊界處會(huì)產(chǎn)生Gibbs振蕩現(xiàn)象。因此可擴(kuò)展所求區(qū)域,以此降低Gibbs震蕩現(xiàn)象的影響,使得計(jì)算區(qū)域的重力場(chǎng)信號(hào)具有周期性[13-15]。

        舉例說(shuō)明,若研究范圍大小的矩形長(zhǎng)、寬分別為為Dx、Dy,數(shù)據(jù)點(diǎn)范圍長(zhǎng)、寬同為dx、dy;可將計(jì)算區(qū)域的平面范圍擴(kuò)展為Dx+dx、Dy+dy,此時(shí)待估重力場(chǎng)信號(hào)在區(qū)間D+d上滿足周期性條件,此時(shí)式(7)為

        當(dāng)擴(kuò)展參數(shù)增大到一定大小時(shí)可以有效降低邊緣效應(yīng),但過(guò)大的擴(kuò)展參數(shù)會(huì)使求解過(guò)程存在病態(tài),同樣影響帶球擾動(dòng)位系數(shù)的估算精度[3]。

        2 數(shù)值計(jì)算與分析

        2.1 基于地面重力數(shù)據(jù)的算例與分析

        為驗(yàn)證矩諧分析基于陸地重力試驗(yàn)數(shù)據(jù)的可靠性,設(shè)計(jì)如下試驗(yàn):獲取緯度為[29°,34°]、經(jīng)度為[111°,116°]圍成的5°×5°區(qū)域,Dx=504.37 km,Dy=601.34 km,格網(wǎng)間隔為2.5′×2.5′,共計(jì)14 400個(gè)計(jì)算點(diǎn)。利用EGM2008重力位模型的2~2190階系數(shù)計(jì)算得到重力異常的模擬觀測(cè)值,并加入數(shù)學(xué)期望為0、標(biāo)準(zhǔn)差為2 mGal的高斯白噪聲;移去基于GOCO05S模型的2~200階系數(shù)計(jì)算的參考重力異常,得到殘余重力異常;根據(jù)殘余重力異常的矩諧展開(kāi)式得到矩諧系數(shù);根據(jù)大地水準(zhǔn)面差距的矩諧表達(dá)式求得殘余大地水準(zhǔn)面,恢復(fù)參考大地水準(zhǔn)面得到最終的大地水準(zhǔn)面;最后利用EGM2008模型2~2190階位系數(shù)計(jì)算的大地水準(zhǔn)面高作為真實(shí)值,比較基于矩諧分析的結(jié)果與真實(shí)值間差距,說(shuō)明矩諧分析的可靠性與精度。

        根據(jù)上文所述可知,矩諧分析受周期延拓邊界效應(yīng)的制約,并成為影響其精度的關(guān)鍵因素。本試驗(yàn)利用擴(kuò)展參數(shù)法以減小邊界效應(yīng)的影響,圖2中虛線范圍為3°×3°中心區(qū)域。由于在不同截?cái)嚯A數(shù)下大地水準(zhǔn)面誤差變化規(guī)律相同,因此以截?cái)嚯A數(shù)30為例進(jìn)行邊界效應(yīng)的研究,并繪制誤差分布圖。

        選取擴(kuò)展參數(shù)為截?cái)嚯A數(shù)30、擴(kuò)展參數(shù)為50 km 繪制大地水準(zhǔn)面對(duì)比圖,如圖3、圖4所示。

        2.2 基于航空重力數(shù)據(jù)的算例與分析

        由于航空重力測(cè)量獲取的是飛行高度處的重力擾動(dòng)值,但所需值為大地水準(zhǔn)面上的擾動(dòng)值,而將飛行高度處的數(shù)值向下延拓至大地水準(zhǔn)面的過(guò)程中,噪聲會(huì)被放大,重力場(chǎng)信息因此失真,這是一個(gè)病態(tài)求解的過(guò)程。但航空重力測(cè)量對(duì)于獲取高頻重力場(chǎng)信息有不可替代的作用,因此如何能在延拓過(guò)程中將重力場(chǎng)信息高度還原,獲得穩(wěn)定的重力場(chǎng)解是航空重力數(shù)據(jù)處理的難點(diǎn)問(wèn)題。

        為驗(yàn)證矩諧分析的全面性、可靠性,本文基于航空重力數(shù)據(jù)設(shè)計(jì)如下試驗(yàn)。獲取緯度為[29°,34°]、經(jīng)度為[111°,116°]圍成的5°×5°區(qū)域,Dx=504.37 km,Dy=601.34 km,格網(wǎng)間隔為2.5′×2.5′,共計(jì)14 400個(gè)計(jì)算點(diǎn)。利用EGM2008重力位模型的2~2190階系數(shù)分別計(jì)算飛行高度為4 km的重力擾動(dòng)的模擬觀測(cè)值,并加入數(shù)學(xué)期望為0、標(biāo)準(zhǔn)差為2 mGal的高斯白噪聲;移去基于GOCO05S模型的2~200階系數(shù)計(jì)算的參考重力擾動(dòng),得到殘余重力擾動(dòng);根據(jù)殘余重力擾動(dòng)的矩諧展開(kāi)式得到矩諧系數(shù);根據(jù)大地水準(zhǔn)面差距的矩諧表達(dá)式和重力擾動(dòng)的矩諧展開(kāi)式求得殘余大地水準(zhǔn)面;最后恢復(fù)參考大地水準(zhǔn)面及參考重力擾動(dòng)得到最終的大地水準(zhǔn)面;最后利用EGM2008模型2~2190階位系數(shù)計(jì)算的大地水準(zhǔn)面作為真實(shí)值,比較基于矩諧分析的計(jì)算結(jié)果與真實(shí)值間的差距。

        圖2 截?cái)嚯A數(shù)30時(shí)取不同擴(kuò)展參數(shù)矩諧分析所得大地水準(zhǔn)面的誤差

        圖3 真實(shí)大地水準(zhǔn)面

        圖4 由重力異常值所得大地水準(zhǔn)面

        與重力異常處理方法相同,航空重力值也選用擴(kuò)展參數(shù)法以減小邊界效應(yīng)的影響。由于不同截?cái)嚯A數(shù)時(shí)的大地水準(zhǔn)面誤差變化規(guī)律相同,因此選取截?cái)嚯A數(shù)30為例繪制誤差分布圖,如圖5所示。

        圖5 截?cái)嚯A數(shù)30時(shí)取不同擴(kuò)展參數(shù)矩諧分析所得大地水準(zhǔn)面的誤差

        選取擴(kuò)展參數(shù)為截?cái)嚯A數(shù)30、擴(kuò)展參數(shù)為50 km 繪制大地水準(zhǔn)面對(duì)比圖,如圖6、圖7所示。

        圖6 真實(shí)大地水準(zhǔn)面值

        圖7 由重力擾動(dòng)值所得大地水準(zhǔn)面

        3 分析與結(jié)論

        如圖2、圖5所示,當(dāng)截?cái)嚯A數(shù)固定時(shí),針對(duì)每一截?cái)嚯A數(shù),分別選取0、30、50、70 km的擴(kuò)展參數(shù)計(jì)算其大地水準(zhǔn)面的誤差分布。由于各截?cái)嚯A數(shù)的變化規(guī)律一致,因此選取截?cái)嚯A數(shù)30為案例進(jìn)行分析。表1、表2分別給出了不同研究范圍內(nèi)取不同擴(kuò)展參數(shù)時(shí)矩諧分析所得的大地水準(zhǔn)面的誤差統(tǒng)計(jì)信息。

        表15°×5°范圍時(shí)取不同擴(kuò)展參數(shù)時(shí)矩諧分析所得大地水準(zhǔn)面的誤差統(tǒng)計(jì)m

        范圍擴(kuò)展參數(shù)Δ最小值最大值平均值均方差5°×5°0-0.3380.311-0.0150.06730-0.3590.299-0.0150.06750-0.5460.273-0.0150.06370-0.5540.269-0.0140.063

        對(duì)于全圖5°×5°區(qū)域范圍,在4個(gè)角點(diǎn)處誤差出現(xiàn)極大值,邊界效應(yīng)明顯,對(duì)于本研究區(qū)域尤以左上角部分最為嚴(yán)重,誤差最大值達(dá)到3 dm;在擴(kuò)展參數(shù)為0、30 km時(shí)均方誤差基本不變;當(dāng)擴(kuò)展參數(shù)不斷增大時(shí),整體誤差逐步降低,當(dāng)擴(kuò)展參數(shù)為50 km 時(shí),均方誤差產(chǎn)生小型跳變,最大值降低到2.7 cm,均方誤差為6.3 cm;繼續(xù)增大擴(kuò)展參數(shù),達(dá)到70 km時(shí),并不能減小均方誤差,且會(huì)增長(zhǎng)計(jì)算時(shí)間,影響計(jì)算效率。

        表23°×3°范圍時(shí)取不同擴(kuò)展參數(shù)時(shí)矩諧分析所得大地水準(zhǔn)面的誤差統(tǒng)計(jì)m

        范圍擴(kuò)展參數(shù)Δ最小值最大值平均值均方差3°×3°0-0.1570.1190.0150.04830-0.1510.1160.0150.04750-0.1510.1150.0150.04270-0.1510.1150.0140.042

        可以明顯看出圖示中心3°×3°區(qū)域在沒(méi)有加入擴(kuò)展參數(shù)時(shí)相較于周邊區(qū)域誤差較小,但依然受到長(zhǎng)波系統(tǒng)誤差的影響;當(dāng)加入擴(kuò)展參數(shù)時(shí),隨著擴(kuò)展參數(shù)的增大,中心3°區(qū)域大地水準(zhǔn)面的誤差逐步明顯降低,與5°區(qū)域變化趨勢(shì)一致;在擴(kuò)展參數(shù)增大到50 km時(shí)誤差值發(fā)生一個(gè)小型跳變,降低到4.2 cm;繼續(xù)增大擴(kuò)展參數(shù)(Δ=70 km)均方誤差沒(méi)有明顯變化,并不能提高大地水準(zhǔn)面的精度。

        按同樣的方式對(duì)航空重力數(shù)據(jù)進(jìn)行統(tǒng)計(jì),誤差統(tǒng)計(jì)見(jiàn)表3、表4。在不同范圍的不同擴(kuò)展參數(shù)下其變化趨勢(shì)與重力異常數(shù)據(jù)一致。

        綜上所述,邊界效用確實(shí)影響矩諧分析的精度,但采用擴(kuò)展參數(shù)方法后,精度上升明顯。擴(kuò)展參數(shù)為50 km時(shí)矩諧分析得到的大地水準(zhǔn)面的精度最高,擴(kuò)展參數(shù)過(guò)小(如0、30 km)對(duì)精度的提高無(wú)法達(dá)到最佳,當(dāng)擴(kuò)展參數(shù)過(guò)大時(shí)(如70 km),并不能繼續(xù)提高精度,相反會(huì)大大增加解算時(shí)間以及復(fù)雜程度;因此實(shí)際計(jì)算中宜選取擴(kuò)展參數(shù)為50 km進(jìn)行解算。

        表35°×5°范圍時(shí)取不同擴(kuò)展參數(shù)時(shí)矩諧分析所得大地水準(zhǔn)面的誤差統(tǒng)計(jì)m

        范圍擴(kuò)展參數(shù)Δ最小值最大值平均值均方差5°×5°0-0.3340.273-0.0140.06830-0.3540.300-0.0150.06650-0.5460.275-0.0150.06370-0.5490.271-0.0140.062

        表43°×3°范圍時(shí)取不同擴(kuò)展參數(shù)時(shí)矩諧分析所得大地水準(zhǔn)面的誤差統(tǒng)計(jì)m

        范圍擴(kuò)展參數(shù)Δ最小值最大值平均值均方差3°×3°0-0.1580.104-0.0220.04630-0.1560.107-0.0200.04450-0.1530.109-0.0180.04370-0.1540.110-0.0170.043

        限制矩諧分析精度的另一大難題是截?cái)嚯A數(shù)的選取。為獲得針對(duì)矩諧分析的適當(dāng)?shù)慕財(cái)嚯A數(shù),避免選取截?cái)嚯A數(shù)較大時(shí)產(chǎn)生“龍格”現(xiàn)象,截?cái)嚯A數(shù)較小時(shí)無(wú)法獲取有效的數(shù)據(jù)值,本文針對(duì)擴(kuò)展參數(shù)為50的不同截?cái)嚯A數(shù)的矩諧展開(kāi)結(jié)果進(jìn)行分析,誤差分布如圖8所示。數(shù)值統(tǒng)計(jì)分析見(jiàn)表5。

        表5取不同展開(kāi)階數(shù)時(shí)矩諧分析所得大地水準(zhǔn)面的誤差統(tǒng)計(jì)m

        范圍截?cái)嚯A數(shù)(N=M)最小值最大值平均值均方差3°×3°20-0.1560.119-0.0190.05925-0.1520.117-0.0160.04830-0.1530.119-0.0180.04335-0.1510.115-0.0160.04340-0.1520.116-0.0160.043

        可以看出,當(dāng)截?cái)嚯A數(shù)由20增加至40時(shí)誤差逐步減小,均方誤差值由5.9 cm減小至4.3 cm,誤差減小了27%;其中尤以截?cái)嚯A數(shù)從20增加至25時(shí)變化最為明顯,均方誤差減少1.1 cm;當(dāng)截?cái)嚯A數(shù)達(dá)到30時(shí),誤差變化趨于穩(wěn)定,繼續(xù)增大截?cái)嚯A數(shù)并不能減小誤差。相較于截?cái)嗟?5階時(shí),截?cái)嚯A數(shù)為30階時(shí)誤差分布更為均勻,最大值與最小值差距較小。而截?cái)嚯A數(shù)達(dá)到40時(shí)雖然均方誤差值不變,但計(jì)算時(shí)間大大增加,不利于進(jìn)行大量數(shù)據(jù)點(diǎn)的解算。故在實(shí)際計(jì)算中可選取截?cái)嚯A數(shù)30進(jìn)行矩諧分析。

        圖8 選取擴(kuò)展參數(shù)為50 km時(shí)不同截?cái)嚯A數(shù)下矩諧分析所得大地水準(zhǔn)面的誤差

        4 結(jié) 語(yǔ)

        針對(duì)球諧分析及球冠諧分析在求解小區(qū)域大地水準(zhǔn)面時(shí)的解算難點(diǎn),本文利用更穩(wěn)當(dāng)?shù)木刂C分析求解大地水準(zhǔn)面。為驗(yàn)證理論的可靠有效性,設(shè)計(jì)完成了基于EGM2008、基于陸地重力數(shù)據(jù)和航空重力數(shù)據(jù)的試驗(yàn)。試驗(yàn)結(jié)果表明,利用矩諧分析求解大地水準(zhǔn)面,基于陸地觀測(cè)數(shù)據(jù)精度可達(dá)到4 cm,基于航空重力觀測(cè)數(shù)據(jù)精度可達(dá)到4 cm。相比于應(yīng)用廣泛的stokes理論求解大地測(cè)量邊值問(wèn)題所得結(jié)果,精度提高明顯。因此,矩諧分析是確定區(qū)域大地水準(zhǔn)面的理想方法,值得應(yīng)用推廣。

        猜你喜歡
        水準(zhǔn)面階數(shù)重力
        瘋狂過(guò)山車(chē)——重力是什么
        關(guān)于無(wú)窮小階數(shù)的幾點(diǎn)注記
        確定有限級(jí)數(shù)解的階數(shù)上界的一種n階展開(kāi)方法
        仰斜式重力擋土墻穩(wěn)定計(jì)算復(fù)核
        一張紙的承重力有多大?
        GPS似大地水準(zhǔn)面精化及精度分析
        一種新的多址信道有效階數(shù)估計(jì)算法*
        關(guān)于動(dòng)態(tài)電路階數(shù)的討論
        重力異常向上延拓中Poisson積分離散化方法比較
        顧及完全球面布格異常梯度項(xiàng)改正的我國(guó)似大地水準(zhǔn)面精化
        在线观看网址你懂的| 成人网站在线进入爽爽爽| 欧美成人精品午夜免费影视| 精品久久久久久久久午夜福利| 国产乱人伦AⅤ在线麻豆A| 日本一区二区三区综合视频| 女人18毛片a级毛片| 亚洲欧洲无码av不卡在线| 国产69口爆吞精在线视频喝尿| 亚洲大胆美女人体一二三区| 亚洲中文字幕日产无码| 日韩人妻无码一区二区三区| 欧美日韩国产在线观看免费| 日本成年少妇人妻中文字幕| 欧美高清视频手机在在线| 性激烈的欧美三级视频| 亚洲中文欧美日韩在线人| 亚洲小少妇一区二区三区| 99噜噜噜在线播放| 韩国无码av片在线观看网站 | 成人免费无码大片a毛片软件| āV第三区亚洲狠狠婷婷综合久久| 免费一区二区三区av| 亚洲精品国产一二三区| 97伦伦午夜电影理伦片| 亚洲精品成人国产av| 成人男性视频在线观看| 国产精品国三级国产av| 无码欧亚熟妇人妻AV在线外遇 | 丰满人妻无奈张开双腿av| 免费a级毛片无码a∨免费| 日本成人精品在线播放| 亚洲精品乱码8久久久久久日本 | 国产女人18毛片水真多18精品| 久久久久久久人妻无码中文字幕爆| 久久熟女五十路| 日韩美女av一区二区| 久久精品无码一区二区三区免费 | 亚洲免费无毛av一区二区三区| 一区二区三区日韩亚洲中文视频| 精品日产卡一卡二卡国色天香|