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

        ?

        GPT2w模型在南極地區(qū)精度分析

        2018-10-26 08:09:32姚宜斌單路路王澤民
        測繪學(xué)報 2018年10期
        關(guān)鍵詞:模型

        孔 建,姚宜斌,單路路,王澤民

        1. 武漢大學(xué)中國南極測繪研究中心,湖北 武漢 430079; 2. 武漢大學(xué)測繪學(xué)院,湖北 武漢 430079; 3. 武漢大學(xué)地球空間環(huán)境與大地測量教育部重點實驗室,湖北 武漢 430079

        GNSS定位過程中由于非色散效應(yīng)影響,對流層延遲不能采用頻率組合的方式消除,因此其延遲量是影響定位精度的主要誤差源之一[1]。天頂方向?qū)α鲗友舆t影響一般能達(dá)到2~3 m左右,當(dāng)高度角在10°以下時,傾斜對流層誤差影響可達(dá)20 m以上[2-3]。目前常采用模型改正的方法來削弱對流層延遲影響,即將信號傳播路徑上的對流層延遲通過映射函數(shù)投影到天頂方向,然后對天頂對流層延遲(zenith tropospheric delay,ZTD)進(jìn)行建模[4-5]。

        傳統(tǒng)的對流層氣象參數(shù)模型如文獻(xiàn)[6—8]提出的模型改正精度可達(dá)厘米級,但模型計算要求輸入實測氣象參數(shù)。由于實測氣象參數(shù)難以獲得,又發(fā)展出一系列不依賴實測氣象參數(shù)的對流層經(jīng)驗?zāi)P汀N墨I(xiàn)[9—10]建立了UNB系列模型,UNB3m模型精度已達(dá)到4.9 cm;文獻(xiàn)[11]指出,歐盟星基廣域增強(qiáng)系統(tǒng)(The European Geostationary Navigation Overlay Service,EGNOS)采用EGNOS模型,平均精度約為5 cm;文獻(xiàn)[12]構(gòu)建了IGGtrop模型,雖然精度有所改進(jìn),但是參數(shù)數(shù)量太多;文獻(xiàn)[13]建立了TropGrid模型,以格網(wǎng)形式提供氣象元素的時序變化特征參數(shù);而后文獻(xiàn)[14]對該模型進(jìn)行了簡化,提出了TropGrid2格網(wǎng)模型;文獻(xiàn)[15]建立了基于球諧函數(shù)表達(dá)的全球非氣象參數(shù)模型GZTD,計算簡便且精度較高。文獻(xiàn)[16]在文獻(xiàn)[17]提出的GPT2模型的基礎(chǔ)上,增加了大氣加權(quán)平均溫度Tm和水汽壓垂直梯度λ兩個參數(shù),建立了GPT2w模型,可同時獲得計算ZHD和ZWD所需的氣象元素,是目前精度最高、使用較為廣泛的對流層延遲模型[18]。GPT2w模型采用ECMWF(European Center for Medium-range Weather Forecast)提供的2001—2010年共10年的高質(zhì)量的ERA-37數(shù)據(jù)進(jìn)行建模,氣象參數(shù)可由式(1)[18]計算得到

        (1)

        式中,r(t)為氣象參數(shù);DOY為年積日;A0為平均值;(A1,B1)為年周期振幅;(A2,B2)為半年周期振幅。

        式(1)得到的氣象參數(shù)是相對地面測站高度的,其他高度上是通過地面數(shù)值進(jìn)行高度歸算得到的,如式(2)所示[19](http:∥ggosatm.hg.tuwien.ac.at/DELAY/SOURCE/GPT2w/)

        (2)

        式中,T0、P0、e0分別為格網(wǎng)點處溫度、氣壓、水汽壓;dT為溫度遞減率;dh為高度;T、P、e分別為格網(wǎng)點對應(yīng)高度為dh處的溫度、氣壓、水汽壓;gm為重力,gm=9.806 65 m/s2;dMtr為干空氣的摩爾質(zhì)量;dMtr=28.965×10-3kg/mol;Rg為氣體常數(shù),Rg=8.314 3 J/K/mol;O為比濕;λ為水汽壓垂直梯度。

        對流層延遲由天頂靜力學(xué)延遲(zenith hydrostatic delay,ZHD)和天頂濕延遲(zenith wet delay,ZWD)組成,一般而言,ZHD可利用Saastamoinen模型[20]精確計算得出

        (3)

        式中,P、φ、H分別代表格網(wǎng)點處總氣壓、緯度和正高。

        而GPT2w模型計算ZWD可由加權(quán)平均溫度、水汽壓及其垂直梯度[16]通過式(4)得到

        (4)

        進(jìn)而可得到對流層延遲

        ZTD=ZHD+ZWD

        (5)

        從式(1)可以看到,GPT2w模型氣象參數(shù)變化周期項僅考慮到年周期項和半年周期項,沒有考慮更短周期變化,如季節(jié)周期,日周期項。文獻(xiàn)[21]指出,GPT2w模型提供的氣壓、溫度、水汽壓在全球范圍內(nèi)均具有很高的精度和穩(wěn)定性,但在南北極地區(qū)精度仍然較低。文獻(xiàn)[16]利用全球341個IGS站對GPT2w模型進(jìn)行了精度評估,發(fā)現(xiàn)其在全球范圍具有很好的精度,但是統(tǒng)計過程中南極區(qū)域IGS測站分布較少且不均勻。文獻(xiàn)[22—24]對比了目前的對流層延遲模型,指出GPT2w模型較其他模型精度高,與IGS站得到的對流層延遲最為接近。文獻(xiàn)[25]指出,GPT2w模型得到的氣壓等參數(shù)由于沒有考慮其日周期項,在寒冷季節(jié)存在較大誤差,可能更適用于熱帶地區(qū)。

        南極地區(qū)地理位置特殊,氣象條件不同于地球其他區(qū)域,目前專門研究GPT2w在南極區(qū)域精度評價的文章較少,而隨著后勤保障的巨大改善,在南極地區(qū)開展的大地測量科學(xué)研究日趨增多,因此有必要對GPT2w模型在南極地區(qū)的精度進(jìn)行評估。另外,由于GPT2w模型沒有考慮日周期項,本文利用中國第33次南極科考實測探空數(shù)據(jù)進(jìn)一步檢驗了模型在短時間尺度下的精度狀況。

        1 數(shù)據(jù)來源及數(shù)據(jù)處理方法

        無線電探空數(shù)據(jù)可以提供垂直氣象元素(如氣壓、溫度、水汽壓等)的剖面信息,且垂直分辨率和精度較高,經(jīng)常被作為參考(真值),用來評估其他方式獲取的氣象元素的精度[26-27]。本文選取兩種類型的探空數(shù)據(jù)來對GPT2w模型提供的氣壓、溫度、水汽壓進(jìn)行精度評估:一是全球站點無線電探空資料數(shù)據(jù)集(the integrated global radiosonde archive,IGRA)中南極探空站數(shù)據(jù),二是中國第33次南極科考期間施放的探空氣球?qū)崪y數(shù)據(jù)。

        美國國家氣候數(shù)據(jù)中心(NCDC)研制的IGRA包含全球1538個探空站點自1938年至今的多層次、多要素逐日探空資料。該數(shù)據(jù)是目前全球范圍內(nèi)時空密度最高、資料最為完整的一套高質(zhì)量探空站點資料數(shù)據(jù)集[28]。探空站施放探空氣球的時間為世界時每天的00時和12時,但由于天氣、設(shè)備故障等原因,會導(dǎo)致數(shù)據(jù)缺失。

        IGRA探空站能直接提供氣壓、溫度、水汽壓等氣象參數(shù)。南極科考期間施放的探空氣球數(shù)據(jù)只能提供氣壓和溫度,不能直接獲得水汽壓,可通過露點溫度對水汽壓進(jìn)行求解[29],具體公式如下

        (6)

        式中,a、b為常數(shù):溫度高于0℃時,a=17.26,b=237.29;溫度低于-40℃時,a=21.87,b=265.49;若溫度介于兩者之間,a、b用線性內(nèi)插獲得。

        為了對GPT2w模型在高度上的精度進(jìn)行評估,本文按照氣球上升高度進(jìn)行分層,以探究模型隨高度的精度變化。由于對流層及平流層沒有明顯的分層特性,結(jié)合數(shù)據(jù)的垂直分布特點采用不均勻分層[30],共分為9層,分層高度見表1,其中高度表示離地面高度。

        表1 高度分層信息

        采用偏差(Bias)和均方根誤差(RMS)作為精度評定標(biāo)準(zhǔn),分別對分層數(shù)據(jù)進(jìn)行評估。Bias和RMS分別由式(7)和式(8)計算得到

        (7)

        (8)

        式中,xmodel,i表示模型計算值;xreal,i表示探空數(shù)據(jù)氣象元素(本文視為真值);n為數(shù)據(jù)總個數(shù)。

        2 GPT2w模型南極地區(qū)精度檢驗及分析

        2.1 利用IGRA探空站數(shù)據(jù)進(jìn)行模型氣象參數(shù)精度分析

        取南極地區(qū)探空站2016年全年數(shù)據(jù)作為參考值,共選取14個具有可用數(shù)據(jù)的探空站,與GPT2w模型計算的氣壓、溫度、水汽壓進(jìn)行分層對比。探空站坐標(biāo)信息見表2。由于南極特殊的地理位置,全年僅分為寒、暖兩季,4月—10月是寒季,11月—次年3月是暖季。取1月和7月分別作為暖季、寒季代表月份,對模型精度的季節(jié)性特征進(jìn)行初步分析,并分層分季節(jié)統(tǒng)計Bias和RMS。

        表2 探空站坐標(biāo)信息

        表3和表4為GPT2w模型氣壓、溫度、水汽壓的分層分季節(jié)Bias和RMS統(tǒng)計結(jié)果。從表中可以看出,隨著層數(shù)(高度)的增加,三者的Bias和RMS均呈現(xiàn)增大的趨勢;當(dāng)?shù)竭_(dá)第7、8層左右,隨著高度增加,氣壓的Bias有減小趨勢,其Bias最大出現(xiàn)在第8層(10~15 km)左右,數(shù)值為7月的41 mb左右;溫度的Bias逐漸減小為負(fù)數(shù),但其絕對值仍遠(yuǎn)大于地面附近數(shù)值,其Bias絕對值最大出現(xiàn)在第9層(15 km以上),數(shù)值為7月的-45 K左右;水汽壓的Bias和RMS也有增加趨勢,但是沒有氣壓和溫度變化顯著,其Bias最大出現(xiàn)在第9層(15 km以上),數(shù)值為7月的38 mb。對比發(fā)現(xiàn),在全年時間尺度上,GPT2w在地面附近精度較高,但是隨高度增加,精度大致呈下降趨勢。GPT2w模型一般提供地面高度的氣象參數(shù)初值,其他高度氣象參數(shù)的計算都是從地面進(jìn)行歸算的,隨著高度增加,誤差會累積,從而導(dǎo)致精度降低。

        表3 1月GPT2w模型氣壓、溫度、水汽壓的精度統(tǒng)計

        表4 7月GPT2w模型氣壓、溫度、水汽壓的精度統(tǒng)計

        對比1月與7月統(tǒng)計結(jié)果,氣壓、溫度和水汽壓精度整體1月精度數(shù)值小于7月,說明1月精度較7月高,即GPT2w模型精度在南極區(qū)域呈現(xiàn)出季節(jié)性特征。

        2.2 利用南極探空氣球?qū)崪y數(shù)據(jù)進(jìn)行模型氣象參數(shù)精度分析

        本文選取2016年11月到2017年4月第33次南極考察期間施放的環(huán)南極18顆探空氣球獲得的實測探空氣象數(shù)據(jù),通過高度分層對比其與GPT2w模型計算的氣溫、氣壓、水汽壓,對模型在南極地區(qū)的短時間尺度精度進(jìn)行進(jìn)一步分析。由于探空氣球是在航行過程中施放的,每個特定位置施放一個,達(dá)到既定高度氣球會自行爆炸,因此每顆氣球只能記錄當(dāng)天的一組數(shù)據(jù),共有18組數(shù)據(jù),分布在不同時間和地點。表5為施放探空氣球的坐標(biāo)信息。

        表5 探空氣球坐標(biāo)信息

        2.2.1 利用ECMWF氣壓分層數(shù)據(jù)驗證氣球?qū)崪y數(shù)據(jù)

        歐洲中期天氣預(yù)報中心(European Centre for Medium-Range Weather Forecasts,ECMWF)是全球數(shù)值天氣預(yù)報組織,可提供高精度的ERA-Interim產(chǎn)品,在全球范圍內(nèi)包括南北極地區(qū)都具有很高的精度[20]。本文選取ECMWF氣壓分層數(shù)據(jù),與探空氣球?qū)崪y氣壓數(shù)據(jù)進(jìn)行對比,來驗證氣球氣象數(shù)據(jù)的可靠性,以15185729號探空氣球數(shù)據(jù)為例,選取與其時間及地點對應(yīng)的ECMWF氣壓分層數(shù)據(jù),兩者隨高度的變化如圖1所示??梢钥吹?,兩者符合較好,說明南極探空氣球?qū)崪y數(shù)據(jù)具有較高的可靠性,可以用來對GPT2w模型精度進(jìn)行評估。

        2.2.2 利用探空氣球?qū)崪y數(shù)據(jù)對模型氣象參數(shù)進(jìn)行精度分析

        表6為Bias和RMS統(tǒng)計結(jié)果。可以看出,對于氣壓和溫度而言,其趨勢與探空站統(tǒng)計結(jié)果類似,表現(xiàn)為地面精度最高,隨高度增加,精度下降。當(dāng)?shù)竭_(dá)第7、8層左右時,氣壓的Bias有減小趨勢,Bias最大出現(xiàn)在第8層(10~15 km),數(shù)值為31 mb左右;溫度的Bias逐漸減小為負(fù)數(shù),但其絕對值仍遠(yuǎn)大于地面附近數(shù)值,Bias絕對值最大出現(xiàn)在第9層(15 km以上),數(shù)值為-60 K左右;而水汽壓的Bias和RMS都比較小,在0附近波動,保持在基本穩(wěn)定狀態(tài),這是由于實測數(shù)據(jù)與GPT2w模型的水汽壓均是利用指數(shù)模型計算得來的,兩者之間差異較小。

        圖1 ECMWF氣壓分層數(shù)據(jù)與氣球?qū)崪y氣壓數(shù)據(jù)對比Fig.1 Comparison between ECMWF pressure data and measured pressure data from sounding balloons

        表6 利用實測數(shù)據(jù)對GPT2w模型氣壓、溫度、水汽壓的精度統(tǒng)計

        與2.1節(jié)相比,雖然GPT2w沒有考慮日周期項,但是溫度和氣壓的天內(nèi)精度與月統(tǒng)計精度在地面附近精度相當(dāng)。但是垂直高度上,天內(nèi)精度與月統(tǒng)計精度相比有所下降,特別是與1月精度對比,隨著高度的增加精度在下降,可見未考慮日周期變化項對模型精度有所影響。

        圖2、圖3、圖4是GPT2w模型分別與探空站1月、7月及氣球?qū)崪y數(shù)據(jù)對比得到的氣壓、溫度、水汽壓的Bias,圖中給出了各層的Bias分布。對于氣壓和溫度,三者的趨勢有一致性,表現(xiàn)為地面附近精度較高,隨著高度增加精度下降;本文使用的實測數(shù)據(jù)采集于2016年11月至2017年4月,對應(yīng)于2.1節(jié)統(tǒng)計的暖季區(qū)間,整體而言,GPT2w模型與實測數(shù)據(jù)的精度對比結(jié)果差于探空站1月對比精度,而優(yōu)于7月對比精度,與2.1節(jié)得到的暖季精度略優(yōu)于寒季精度的結(jié)果相符。而與探空站水汽壓增幅較大不同,氣球?qū)崪y數(shù)據(jù)的水汽壓較小且波動不大。

        圖2 GPT2w模型與探空站1月、7月及氣球?qū)崪y氣壓數(shù)據(jù)對比的BiasFig.2 Comparison of pressure between GPT2w and the radiosonde stations in January and July and sounding balloons,respectively

        圖3 GPT2w模型與探空站1月、7月及氣球?qū)崪y溫度數(shù)據(jù)對比的BiasFig.3 Comparison of temperature between GPT2w and the radiosonde stations in January and July and sounding balloons,respectively

        2.3 利用探空數(shù)據(jù)進(jìn)行模型對流層延遲精度分析

        2.3.1 利用探空數(shù)據(jù)計算對流層延遲的方法及原理

        探空數(shù)據(jù)計算對流層延遲是根據(jù)其定義進(jìn)行求解的。GPS衛(wèi)星信號在對流層中的延遲為[31]

        (9)

        圖4 GPT2w模型與探空站1月、7月及氣球?qū)崪y水汽壓數(shù)據(jù)對比的BiasFig.4 Comparison of water vapor pressure between GPT2w and the radiosonde stations in January and July and sounding balloons,respectively

        因此,ΔL可以進(jìn)一步表示為

        (10)

        (11)

        式中,N=(n-1)·106表示大氣折射率。

        大氣折射率計算公式為[32]

        Nd+Nw

        (12)

        (13)

        1.75×10-4t2+1.44×10-6t3)

        (14)

        將對流層延遲繼續(xù)沿高度進(jìn)行積分,則得到GPS信號穿過對流層天頂方向的總延遲

        (15)

        2.3.2 利用探空數(shù)據(jù)對模型對流層延遲進(jìn)行精度分析

        利用上述計算方法,分別計算GPT2w模型與探空數(shù)據(jù)計算得到的對流層延遲,結(jié)果見表7、表8。

        表7 GPT2w模型與探空站計算結(jié)果對比

        數(shù)據(jù)對比時段:2017-01-16

        表8 GPT2w模型與實測探空數(shù)據(jù)對比

        由兩表可以看出,由GPT2w模型及探空數(shù)據(jù)計算得到的ZTD差距約為幾個厘米,這與其他學(xué)者驗證得到的GPT2w模型對流層延遲在全球范圍內(nèi)精度為厘米級的結(jié)論是一致的[16]。同時可以看出,ZWD數(shù)值較小,僅為幾個厘米,一般來說,在大氣中,ZHD占90%左右,ZWD占10%左右,而南極處于高緯度地區(qū),水汽含量很少,故ZWD相對來說更小。且兩者ZWD相差最多在氣球編號為15188344處,約為5.7 cm,數(shù)值較小。

        在探空數(shù)據(jù)計算結(jié)果作為真值的情況下,統(tǒng)計模型ZWD和ZTD的Bias及RMS,結(jié)果見表9??梢钥闯?,GPT2w得到的對流層延遲與探空數(shù)據(jù)計算得到的對流層延遲較為接近,Bias及RMS均處于厘米級,說明該模型可用于一般定位等對流層延遲的求解。但在精密定位中,仍需對該模型進(jìn)行精化。

        表9 ZWD及ZTD的Bias、RMS統(tǒng)計

        3 結(jié)果討論及分析

        3.1 結(jié)果分析

        一般認(rèn)為,探空氣球在上升過程中,只有高度發(fā)生變化,經(jīng)緯度不發(fā)生變化。

        (1) 由上述分析可知,在南極地區(qū),Bias和RMS隨著高度的增加而發(fā)生變化,GPT2w模型在地面附近精度最高;用ECMWF氣壓分層數(shù)據(jù)對實測氣壓數(shù)據(jù)進(jìn)行檢驗,兩者符合較好,說明了實測數(shù)據(jù)具有可靠的測量精度,Bias和RMS隨高度變差主要是由于GPT2w模型本身誤差造成的。

        (2) GPT2w模型中,氣壓和溫度的計算均與高度有關(guān),氣壓采用高度的指數(shù)函數(shù)擬合,溫度采用高度的線性函數(shù)擬合得到,而水汽壓的計算與氣壓相關(guān)。因此,高度的變化,會導(dǎo)致模型精度的變化。GPT2w模型一般用來提供地面高度的氣象參數(shù)初值,其他高度的參數(shù)值模型內(nèi)部進(jìn)行了函數(shù)改正,即高度歸算,而這種改正本身存在一定誤差。且會隨著高度增加,空氣成分及其他因素的變化,進(jìn)一步導(dǎo)致模型精度下降,從而導(dǎo)致模型值與實測值出現(xiàn)較大不符。

        3.2 殘差擬合

        GPT2w模型計算的氣象參數(shù)值與實測數(shù)據(jù)在高度上存在較大的差異,本節(jié)對兩者之間的差異作進(jìn)一步的分析。以15187632號探空氣球數(shù)據(jù)為例,圖5給出的是實測氣壓與GPT2w模型計算氣壓隨時間變化對比圖。從圖中可以看出,在地面附近,兩者相差不大,隨著高度的增加,兩者之間的差異在增大,在最高處略有下降。從兩者的殘差圖(圖6)也可以看出這一現(xiàn)象。圖6中開始時段異常是由于探空氣球是在走航期間施放的,開始時段氣球運動處于不穩(wěn)定的上下浮動狀態(tài)。對兩者殘差進(jìn)行初步擬合,如圖7所示,可以看到,用三次多項式就可以較好的擬合兩者之間的差異。通過對溫度和水汽壓進(jìn)行類似計算,可以得到相同的結(jié)論,如圖8、圖9所示。

        因此,下一步可以對GPT2w模型在高度上與實測數(shù)據(jù)的差異進(jìn)行模型擬合改正,進(jìn)而對該模型在南極地區(qū)乃至全球垂直方向上的精度進(jìn)行改進(jìn)。

        圖5 氣球?qū)崪y氣壓數(shù)據(jù)與模型計算氣壓對比Fig.5 Comparison of pressure between results calculated by GPT2w and results measured by sounding balloons

        圖6 氣球?qū)崪y氣壓數(shù)據(jù)與模型計算氣壓殘差Fig.6 Residual of pressure between results calculated by GPT2w and results measured by sounding balloons

        圖7 氣壓殘差三階多項式擬合 Fig.7 Three-order polynomial fitting of pressure’s residual

        圖8 溫度殘差三階多項式擬合Fig.8 Three-order polynomial fitting of temperature’s residual

        圖9 水汽壓殘差三階多項式擬合Fig.9 Three-order polynomial fitting of water vapor pressure’s residual

        4 結(jié) 論

        本文通過探空站和探空氣球提供的氣象數(shù)據(jù)對GPT2w模型在南極區(qū)域計算的氣壓、溫度、水汽壓進(jìn)行精度檢驗,結(jié)果表明,該模型在地面附近精度可達(dá)到與全球其他區(qū)域精度一致的水平,但隨著高度的增加,Bias和RMS有不同程度的增加,精度有所下降。由模型計算得到的對流層延遲可達(dá)到厘米級精度,可用于一般定位等方面的對流層延遲改正。得到如下結(jié)論:

        (1) GPT2w一般用于地面氣象參數(shù)的計算,本文分析顯示,在地面附近,模型在南極地區(qū)精度與全球其他位置精度一致;且在南極區(qū)域,模型精度表現(xiàn)出季節(jié)性。

        (2) 通過與實測數(shù)據(jù)對比發(fā)現(xiàn),GPT2w模型雖然沒有考慮日周期項,但是對于南極區(qū)域而言,在地面高度上,計算的氣象參數(shù)天內(nèi)精度與月平均精度相當(dāng)。但是隨著高度的增加,精度較1月精度有所下滑,而略優(yōu)于7月精度。

        (3) 隨著高度的增加,模型值與探空數(shù)據(jù)差異增加,除了氣球升空過程中產(chǎn)生漂移等原因外,主要原因是模型本身計算氣象元素時所使用的高度歸算函數(shù)存在一定誤差。同時可以發(fā)現(xiàn),殘差規(guī)律性較強(qiáng),因此,可以在考慮季節(jié)性因素的基礎(chǔ)上,對其GPT2w模型加以改進(jìn),使其在垂直高度上精度得以改善。

        (4) 總體而言,由GPT2w模型計算的ZTD在南極地區(qū)有較高的精度,Bias和RMS均在厘米級,與全球其他地區(qū)得到的ZTD精度較為一致。

        猜你喜歡
        模型
        一半模型
        一種去中心化的域名服務(wù)本地化模型
        適用于BDS-3 PPP的隨機(jī)模型
        提煉模型 突破難點
        函數(shù)模型及應(yīng)用
        p150Glued在帕金森病模型中的表達(dá)及分布
        函數(shù)模型及應(yīng)用
        重要模型『一線三等角』
        重尾非線性自回歸模型自加權(quán)M-估計的漸近分布
        3D打印中的模型分割與打包
        国产午夜福利精品| 亚洲乱码日产精品bd| 亚洲成人电影在线观看精品国产 | 精品中文字幕久久久久久| 中文字幕人妻少妇精品| 久久精品色福利熟妇丰满人妻91| 午夜福利理论片在线观看| 国产av日韩a∨亚洲av电影| 亚洲av福利天堂在线观看| 日本中文字幕精品久久| 久久久亚洲精品一区二区三区| 狠狠色婷婷久久一区二区| 亚洲精品国产精品av| 国产在线一区二区av| 国产一区二区三区免费观看在线| 亚洲国际无码中文字幕| 亚洲综合精品在线观看中文字幕| 亚洲国产av一区二区四季| 97精品超碰一区二区三区| 日本国产视频| 国产精品国产三级在线专区| 男女男精品视频网站免费看| 亚洲av无码av制服另类专区 | 久久亚洲国产欧洲精品一| 亚洲白嫩少妇在线喷水| 成人免费看aa片| 色94色欧美sute亚洲线路二| 日本肥老熟妇在线观看| 电驱蚊液可以插一晚上吗| 99精品欧美一区二区三区| 97SE亚洲国产综合自在线不卡| 凹凸世界视频a一二三| 亚洲av永久无码天堂网小说区 | 亚洲av乱码国产精品观看麻豆| 国产亚洲欧洲aⅴ综合一区| 18禁美女裸身无遮挡免费网站| 尤物蜜芽福利国产污在线观看 | 甲状腺囊实性结节三级| 尤物蜜桃视频一区二区三区| 女人脱了内裤趴开腿让男躁| 欧美日韩综合网在线观看|