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

GPT/2模型用于GPS大氣可降水汽反演的精度分析

2016-04-11 01:20:46范士杰臧建飛劉焱雄秦學彬耿東哲
測繪工程 2016年3期

范士杰,臧建飛,劉焱雄,秦學彬,華 亮,耿東哲

(1.中國石油大學 地球科學與技術(shù)學院,山東 青島 266580;2.國家海洋局第一海洋研究所,山東 青島 266061;3.東方地球物理公司裝備服務(wù)處測量服務(wù)中心,天津 大港 300280)

?

GPT/2模型用于GPS大氣可降水汽反演的精度分析

范士杰1,臧建飛1,劉焱雄2,秦學彬3,華亮1,耿東哲1

(1.中國石油大學 地球科學與技術(shù)學院,山東 青島 266580;2.國家海洋局第一海洋研究所,山東 青島 266061;3.東方地球物理公司裝備服務(wù)處測量服務(wù)中心,天津 大港 300280)

摘要:氣象參數(shù)(溫度T、氣壓P)是GPS大氣可降水汽(PWV)反演中必不可少的數(shù)據(jù),也是PWV反演的重要誤差源之一。文中主要對GPT/2(GPT、GPT2)模型用于PWV反演的精度進行驗證和分析。基于非差精密單點定位(PPP)技術(shù),選取SuomiNet網(wǎng)9個測站的觀測數(shù)據(jù),借助研制的PPP軟件,分別采用GPT模型、改進的GPT2模型以及測站實測氣象數(shù)據(jù)進行大氣可降水汽(PWV)反演。以實測氣象數(shù)據(jù)處理結(jié)果為參考,對兩種模型解算的PWV進行了對比和精度分析。結(jié)果表明:改進的GPT2模型優(yōu)于GPT模型,尤其是當測站的高程較大時,GPT2模型的穩(wěn)定性更優(yōu)、適用性更廣;采用GPT2模型解算的PWV偏差均值小于±1.0 mm,精度(RMS)優(yōu)于±1.5 mm。在缺少實測氣象數(shù)據(jù)的情況下,利用GPT2模型數(shù)據(jù)仍然能夠取得較為理想的PWV反演結(jié)果。

關(guān)鍵詞:精密單點定位;氣象參數(shù);GPT模型;GPT2模型;大氣可降水汽

水汽是大氣中的活躍成分,大氣水汽的時空變化對于理解天氣過程、天氣預報和天氣系統(tǒng)演變等研究具有重要意義。利用精密單點定位(precise point positioning,PPP)技術(shù)進行大氣可降水汽(precipitable water vapor,PWV)反演,具有模型簡單、處理速度快、站間不相關(guān)、無需引入遠距離測站即可直接估計測站絕對時延等顯著優(yōu)勢[1-2]。然而要想獲得高精度的水汽產(chǎn)品,必須同步獲取測站處的氣象參數(shù)(溫度T、氣壓P)。對于沒有安裝氣象傳感器的測站,則只能通過其它手段獲取這些參數(shù)[3],但這勢必會影響PWV反演的精度。

目前,GNSS氣象學中最常用的氣象模型是GPT(global pressure and temperature)模型[4],它是基于9階9次球諧函數(shù)建立的經(jīng)驗模型,其時空分辨率較差(時間分辨率為1 d,空間分辨率為 20°×20°),只能近似反映溫度、氣壓的時空變化。Kouba分析了在高緯度地區(qū)和低高度角的情況下,GPT模型估計的氣壓誤差對PPP解算的高程方向的誤差影響[5];Steigenberger等分析了在不考慮大氣負荷改正的情況下,采用GPT模型估算的測站高程時間序列比采用歐洲中尺度數(shù)值預報模型數(shù)據(jù)計算的高程時間序列具有更高的重復率[6];于勝杰等分析了GPT模型計算的氣壓誤差對天頂總延遲ZTD(zenith tropospheric delay)估計的影響,發(fā)現(xiàn)在天氣不存在劇烈變化的情況下,由GPT模型得到的ZTD能夠滿足氣象業(yè)務(wù)的需求[7];陳澍分析了GPT模型用于PWV反演的可行性,發(fā)現(xiàn)該模型適合于無實測氣象數(shù)據(jù)且對水汽反演精度要求不高的測站[8]。GPT2模型是GPT模型的改進形式,空間分辨率提高至5°×5°,時間分辨率仍為1 d;該模型考慮了半年周期項影響,其氣壓的精度優(yōu)于GPT模型[3]。在PPP中引入GPT2模型,其單天解的坐標精度可達到mm級[9]。

目前關(guān)于GPT2模型對PWV反演的精度分析以及GPT/2兩種模型用于PWV反演的精度比較的文獻較少。基于非差PPP技術(shù),利用SuomiNet網(wǎng)9個測站的GPS觀測數(shù)據(jù)和同步氣象數(shù)據(jù),借助研制的PPP軟件,本文對GPT2模型用于GPS PWV反演的精度進行驗證,對GPT/2兩種模型用于PWV反演的精度進行了比較和分析。

1數(shù)學模型

1.1基于PPP的PWV估計

本文精密單點定位(PPP)采用傳統(tǒng)的雙頻消電離層組合觀測模型,其相位觀測值的觀測方程為

(1)

式中:Lc為雙頻消電離層組合相位觀測值;p為站星幾何距離;c為光速;dtr為接收機鐘差;dts為衛(wèi)星鐘差;λ1為L1載波波長;Nc為對應(yīng)的相位模糊度參數(shù);天頂干延遲ZHD(zenith hydrostatic delay)由經(jīng)驗模型計算(如Saastamoinen模型);天頂濕延遲ZWD(zenith wet delay)作為未知參數(shù),與接收機坐標等參數(shù)一起進行平差計算;Mdry和Mwet分別表示干、濕延遲的映射函數(shù)(如GMF模型[10]);M為大氣水平梯度的映射函數(shù);GN和GE分別為大氣水平梯度改正的北向和東向分量;θ為衛(wèi)星方位角;E為衛(wèi)星高度角;ε為包含觀測噪聲、多路徑效應(yīng)以及未被模型化的各項誤差影響。

大氣可降水汽(PWV)和PPP解算得到的天頂濕延遲(ZWD)密切相關(guān),二者之間的物理轉(zhuǎn)換關(guān)系可表示如下[11]:

PWV=Π·ZWD,

(2)

(3)

式中:Π為轉(zhuǎn)換系數(shù);ρ為水的密度;RV為水汽氣體常數(shù),一般取為461.495 J/(Kg×K);ω為水汽分子與干空氣分子的摩爾質(zhì)量比,其值為0.622;k1,k2,k3為氣體常數(shù);Tm為大氣加權(quán)平均溫度,可通過本地的探空數(shù)據(jù)建立的本地化模型得到[12-13],也可采用已有的經(jīng)驗模型計算,如Mendes給出的全球Tm模型為[14]

Tm=50.4+0.789TS.

(4)

式中,TS為測站的溫度(K)。

從上述PWV的估計過程可以看出,測站氣壓PS與ZHD的計算有關(guān),因此氣壓誤差直接影響著ZHD的計算精度,從而影響到ZWD和PWV的估計精度;而測站溫度TS與大氣加權(quán)平均溫度的計算有關(guān),溫度誤差影響著大氣加權(quán)平均溫度的計算精度,進而影響到轉(zhuǎn)換因子和PWV估值的精度。研究發(fā)現(xiàn),9 hPa的氣壓誤差引起的ZHD誤差約為20.5 mm,從而引起的PWV誤差為3.3 mm[7];而當Tm產(chǎn)生2.5K的誤差時,所引起的PWV誤差約為0.2mm[2]。因此,高精度的氣象數(shù)據(jù)對于PWV精確反演至關(guān)重要。

1.2GPT/2模型

1.2.1GPT模型

GPT模型是目前空間大地測量中比較常用的氣象模型,它是基于歐洲中尺度數(shù)值預報模型ERA40的3年(1999—2002年)全球氣壓、氣溫的月平均資料建立的經(jīng)驗模型[4]。

(5)

(6)

式中:P為測站大氣壓;h為測站大地高;T為測站溫度;dT/dh為溫度在高程方向上的衰減率;Pr為平均海平面大氣壓;hr為平均海平面高。

平均海平面大氣壓和溫度可通過下式計算。

(7)

(8)

式中:a為平均海平面大氣壓或溫度;a0為平均海平面大氣壓或溫度的均值;A為平均海平面大氣壓或溫度的年變化振幅,其可通過式(8)計算;doy為年積日;Pnm為勒讓德多項式;φ為測站大地緯度;λ為測站大地經(jīng)度;Anm和Bnm為球諧系數(shù)。

利用GPT模型計算測站氣壓和溫度的方法如下:首先輸入測站的位置和年積日,利用式(7)和式(8)計算平均海平面的大氣壓和溫度;然后利用式(5)和式(6)將平均海平面的大氣壓和溫度直接改正到測站高程處,從而得到測站的氣壓和溫度值。

1.2.2GPT2模型

GPT2模型是GPT模型的發(fā)展,它是基于10 a(2001—2010年)歐洲中尺度數(shù)值預報模型的再分析資料建立的經(jīng)驗模型[3]。它將全球地表劃分為5°×5°的空間格網(wǎng),除了給出這些格網(wǎng)點的氣壓、溫度值外,還提供了它們的溫度垂直衰減率、水汽壓等,并考慮了這些參數(shù)的年變化振幅和半年變化振幅,具體公式如下[9]:

(9)

式中:X表示格網(wǎng)點的氣壓、溫度等參數(shù)值;A0為對應(yīng)參數(shù)的平均值;A1,B1為對應(yīng)參數(shù)的年變化振幅;A2,B2為對應(yīng)參數(shù)的半年變化振幅,這些參數(shù)可從模型外的一個格網(wǎng)文件中獲得;doy為年積日。

采用GPT2模型計算測站的氣壓和溫度的方法如下:首先輸入測站的位置和年積日,利用式(9)計算測站周圍距離最近的4個格網(wǎng)點的氣壓和溫度;然后根據(jù)各個格網(wǎng)點與測站的高程之差,將各格網(wǎng)點的氣壓和溫度改正到測站高程處,其中格網(wǎng)點的溫度垂直衰減率用式(9)計算,而氣壓改正則是根據(jù)格網(wǎng)點的虛溫以指數(shù)形式[3]完成;最后采用雙線性內(nèi)插得到測站處的氣壓、溫度值。

2試驗驗證和結(jié)果分析

2.1數(shù)據(jù)來源

本文選取位于北美的SuomiNet網(wǎng)sa01、sa06、sa10、sa11、sc02、sc04等9個GPS站點,下載2011-07-02至2011-07-07(年積日183~188)共計6 d的GPS觀測數(shù)據(jù)和同步氣象觀測數(shù)據(jù)。GPS數(shù)據(jù)采樣率為30 s;而不同站點的氣象數(shù)據(jù)的采樣率有所不同,均在1~60 s之間。在本文的后續(xù)數(shù)據(jù)處理中采用就近的原則選取與GPS數(shù)據(jù)采樣率相對應(yīng)的氣象參數(shù)。GPS站點的位置分布如圖1所示。

圖1 GPS站點分布示意圖

2.2試驗方案設(shè)計

借助研制的PPP軟件,分別采用實測氣象數(shù)據(jù)(方案1)、GPT模型(方案2)和改進的GPT2模型(方案3)數(shù)據(jù)進行GPS PWV反演。PPP數(shù)據(jù)處理的基本策略[15]如下:取衛(wèi)星截止高度角為10°,采用IGS事后精密星歷(時間間隔為15 min)和鐘差(時間間隔為30 s)產(chǎn)品;天頂干延遲(ZHD)采用Saastamoinen模型計算;天頂濕延遲(ZWD)作為未知參數(shù),采用分段參數(shù)估計方法,其間隔取為30 min,并采用隨機游走過程模擬其動態(tài)變化;映射函數(shù)選取目前最常用的GMF模型[10];同時顧及水汽分布的各向異性,引入大氣水平梯度改正參數(shù);考慮地球自轉(zhuǎn)改正、相對論效應(yīng)、天線相位中心改正等誤差影響;最后利用序貫最小二乘方法,對接收機坐標、接收機鐘差、相位模糊度、ZWD等參數(shù)進行平差計算。大氣加權(quán)平均溫度則采用Mendes模型。

在上述PWV信息提取中,3種試驗方案所采用的模型和數(shù)據(jù)處理策略相同,不同之處僅在于氣象數(shù)據(jù)的來源不同。

2.3結(jié)果分析

2.3.1ZWD結(jié)果

以sa06站為例,以方案1為標準,將方案2、方案3利用模型計算的氣壓值和ZWD估值,分別與方案1的實測氣壓數(shù)據(jù)和ZWD估值進行比較和求差,可得到GPT和GPT2模型的氣壓誤差為ΔP1、ΔP2,ZWD偏差為ΔZWD1、ΔZWD2。兩種模型的氣壓誤差和ZWD偏差時間序列以及ZWD偏差隨氣壓誤差的變化如圖2所示。從圖2可以看出,兩種模型的氣壓誤差ΔP1和ΔP2以及ZWD偏差ΔZWD1和ΔZWD2在整個試驗期間的變化基本一致,互差很小,說明利用GPT和GPT2模型計算的氣壓值以及ZWD估值的精度相當。同時可以看出,ZWD偏差與氣壓誤差的變化密切相關(guān),隨著氣壓誤差的增大而減小,而隨著氣壓誤差的減小而增大。這主要是因為氣壓誤差直接影響到ZHD的計算,且ZHD與氣壓值成正比。因此,在天頂總延遲量保持不變的情況下,當ZHD誤差增大時,則ZWD誤差一定減小;反之亦然。

圖2 方案2和方案3的ZWD偏差及氣壓誤差的時間序列(sa06站點)

表1給出了所有測站方案2和方案3估算的ZWD偏差和氣壓誤差的統(tǒng)計結(jié)果。由表1可以看出,對于大多數(shù)測站,GPT和GPT2模型計算的氣壓值的精度相當,與實測氣象數(shù)據(jù)的偏差均值在2.5 hPa以內(nèi),RMS優(yōu)于3.5 hPa;由此造成利用兩種模型解算的ZWD估值的偏差均值小于5.5 mm,RMS優(yōu)于7.5 mm,兩種模型解算的ZWD估值的精度也基本一致。

表1 方案2和方案3的ZWD偏差及氣壓誤差的統(tǒng)計結(jié)果

但是,對于sa11和sa46站點,方案3的統(tǒng)計結(jié)果要明顯優(yōu)于方案2的結(jié)果。分析兩種方案的氣壓誤差會發(fā)現(xiàn),sa11和sa46站點GPT模型的氣壓誤差的偏差均值分別達到8.0 hPa和-4.0 hPa,RMS達到8.0 hPa和4.8 hPa,明顯高于其他站點;由此導致利用GPT模型解算的ZWD估值的偏差均值分別達到-17.7 mm和8.8 mm,RMS達到17.9 mm和10.5 mm,同樣明顯超過其他站點的誤差水平。分析其原因可能是因為sa11和sa46站點的高程較大(sa11、sa46的大地高分別為2 211 m、732 m),采用GPT模型時,由式(5)直接將平均海平面大氣壓改正到測站高程處,有可能產(chǎn)生較大的誤差;而改進的GPT2模型先將測站周圍4個格網(wǎng)點的氣壓值改正到測站同一高程面上,然后采用雙線性內(nèi)插的方法得到測站點的氣壓值,因此效果會更好一些。

2.3.2PWV結(jié)果

仍以sa06站點為例,將方案2、方案3利用模型計算的溫度值和PWV估值,與方案1的實測溫度數(shù)據(jù)和PWV估值進行比較并求差,可得到GPT和GPT2模型的溫度誤差為ΔT1、ΔT2,PWV偏差為ΔPWV1、ΔPWV2。兩種模型的氣壓誤差、溫度誤差和PWV偏差的時間序列以及PWV偏差隨氣壓誤差、溫度誤差的變化如圖3所示。從圖3中可以看出,兩種方案的PWV偏差的變化趨勢基本一致,互差較小;PWV偏差與氣壓誤差和溫度誤差的變化密切相關(guān),隨著氣壓誤差的增加而減小,而隨著溫度誤差的增加而增大。PWV由ZWD轉(zhuǎn)化而來,因此氣壓誤差對PWV的影響與氣壓誤差對ZWD的影響相同;溫度誤差對PWV的影響是通過對轉(zhuǎn)換系數(shù)Π中加權(quán)平均溫度的影響完成的,因此PWV偏差與溫度誤差的變化趨勢基本一致。

圖3 方案2和方案3的PWV偏差、氣壓誤差及溫度誤差的時間序列(sa06站點)

將所有測站方案2和方案3的PWV偏差、溫度誤差結(jié)果進行統(tǒng)計,如表2所示。從表2可以看出,GPT模型與實測氣象數(shù)據(jù)的整體偏差均值為2.8 ℃,RMS為7.0 ℃;而GPT2模型與實測氣象數(shù)據(jù)的整體偏差均值為1.0 ℃,RMS為5.5 ℃。因此,GPT2模型計算的溫度值的精度略優(yōu)于GPT模型。

但是,除sa11和sa46站點之外,兩種模型解算的PWV估值的精度相當,其偏差均值都小于1.2 mm,RMS均優(yōu)于1.5 mm。而對于sa11和sa46站點,GPT2模型解算的結(jié)果明顯優(yōu)于GPT模型;GPT2模型的統(tǒng)計結(jié)果與其他站點相似,而GPT模型的統(tǒng)計結(jié)果則出現(xiàn)了異常,其原因可能與氣壓誤差和ZWD偏差的異常情況類似。

表2 方案2和方案3的PWV偏差及溫度誤差的統(tǒng)計結(jié)果

3結(jié)束語

氣象參數(shù)(溫度T、氣壓P)是GPS大氣可降水汽(PWV)反演中必不可少的數(shù)據(jù),也是PWV反演的重要誤差源之一。本文采用非差PPP技術(shù),利用SuomiNet網(wǎng)的GPS觀測數(shù)據(jù)和同步氣象數(shù)據(jù),對GPT2模型用于PWV反演的精度進行了驗證,對GPT/2兩種模型計算的氣壓、溫度參數(shù)以及PWV反演的精度進行了比較和分析。結(jié)果表明:①以實測氣象數(shù)據(jù)為參考值,改進的GPT2模型計算的溫度參數(shù)的精度略優(yōu)于GPT模型;②除sa11、sa46等高程較大的站點之外,兩種模型計算的氣壓參數(shù)以及ZWD、PWV估值的精度均基本一致;但是當測站的高程較大時,GPT模型的氣壓值、溫度值、ZWD及PWV估值均出現(xiàn)異常情況,而GPT2模型仍然保持正常,其穩(wěn)定性更優(yōu)、適用性更廣;③采用GPT2模型解算的PWV偏差均值小于±1.0 mm,精度(RMS)優(yōu)于1.5 mm。因此,在缺少實測氣象數(shù)據(jù)的情況下,采用GPT2模型進行PWV反演仍能得到比較理想的結(jié)果。

參考文獻:

[1]何錫揚,張小紅,李星星,等.PPP 估計天頂對流層延遲方法與結(jié)果分析[J].測繪信息與工程,2010(1):3-5.

[2]張小紅,何錫揚,郭博峰,等.基于GPS非差觀測值估計大氣可降水量[J].武漢大學學報(信息科學版),2010,35(7):806-810.

[3]LAGLER K,SCHINDELEGGER M,B?HM J,et al.GPT2:Empirical slant delay model for radio space geodetic techniques[J].Geophysical Research Letters,2013,40(6):1069-1073.

[4]B?HM J,HEINKELMANN R,SCHUH H.Short note:a global model of pressure and temperature for geodetic applications[J].Journal of Geodesy,2007,81(10):679-683.

[5]KOUBA J.Testing of global pressure/temperature (GPT) model and global mapping function (GMF) in GPS analyses[J].Journal of Geodesy,2009,83(3-4):199-208.

[6]STEIGENBERGER P,BOEHM J,TESMER V.Comparison of GMF/GPT with VMF1/ECMWF and implications for atmospheric loading[J].Journal of Geodesy,2009,83(10):943-951.

[7]于勝杰,萬蓉,付志康.氣壓對GPS大氣可降水量解算的影響分析[J].大地測量與地球動力學,2013,33(2):87-90.

[8]陳澍.利用參考站網(wǎng)絡(luò)數(shù)據(jù)建立精確大氣層析模型研究[D].成都:西南交通大學,2011.

[9]CHEN W,GAO C,PAN S.Assessment of GPT2 Empirical Troposphere Model and Application Analysis in Precise Point Positioning[C]//China Satellite Navigation Conference (CSNC) 2014 Proceedings:Volume II.Springer Berlin Heidelberg,2014:451-463.

[10] BOEHM J,NIELL A,TREGONING P,et al.Global Mapping Function(GMF):A new empirical mapping function based on numerical weather model data[J].Geophysical Research Letters,2006,33:L07304.

[11] ASKNE J,NORDIUS H.Estimation of tropospheric delay for microwaves from surface weather data[J].Radio Science,1987,22(3):379-386.

[12] 朱爽.北京地區(qū)地基GPS加權(quán)平均溫度計算本地化模型研究[J].測繪工程,2014,23(4):28-32.

[13] 陳宏斌,熊永良,陳志勝,等.垂直不均勻分層的地基GPS層析水汽研究[J].測繪工程,2015,24(5):11-14.

[14] MENDES V B,PRATES G,SANTOS L,et al.An evaluation of the accuracy of models for the determination of the weighted mean temperature of the atmosphere[C]//Proceedings of ION.2000:433-438.

[15] 范士杰,劉焱雄,高興國,等.海上動態(tài)GPS大氣可降水量信息反演[J].中國石油大學學報(自然科學版),2012,36(3):84-87,92.

[責任編輯:劉文霞]

Accuracy analysis on GPS precipitable water vapor inversion usingGPT/2 models

FAN Shijie1,ZANG Jianfei1,LIU Yanxiong2,QIN Xuebin3,HUA Liang1,GENG Dongzhe1

(1.School of Geosciences,China University of Petroleum,Qingdao 266580,China;2.The First Institute of Oceanography,State Oceanic Administration,Qingdao 266061,China;3.Survey Center of BGP Equipment,CNPC,Tianjin 300280,China)

Abstract:Meteorological parameters (pressure,temperature) are essential to GPS precipitable water vapor inversion.The feasibilityof GPT/2(GPT、GPT2) model used in GPS precipitable water vapor (PWV) inversion is discussed,and the accuracy of PWV is analyzed in this paper.Based on precise point positioning (PPP) technique,using the GPS observations of SuomiNet sites,the precipitable water vapor (PWV) inversion experiment is carried out with measured meteorological data,GPTmodel and GPT2 model, respectively.The accuracy of PWV derived from GPT and GPT2 model is analyzed, compared with measured meteorological data.The results show that improved GPT2 model is more stable and applicablethan GPT model when the station height is bigger.The mean deviation of PWV with GPT2 modelis is less than ±1.0 mm and the RMS is better than ±1.5 mm. In the absence of measured meteorological data,the PWV inversion using GPT2 model can still get an ideal result.

Key words:precise point positioning;meteorological parameter;GPT model;GPT2 model;precipitable water vapor

中圖分類號:P228.4

文獻標識碼:A

文章編號:1006-7949(2016)03-0001-05

作者簡介:范士杰(1969-),男,副教授,博士.

基金項目:國家自然科學基金面上項目(41274011,41374044)

收稿日期:2015-01-02;修回日期:2015-06-08

主站蜘蛛池模板: 国产午夜一级毛片| 国产亚洲精久久久久久久91| 日本在线亚洲| 国产精品55夜色66夜色| 777午夜精品电影免费看| 国产精品专区第1页| 美女毛片在线| 国产精品妖精视频| 成人综合久久综合| 国产永久免费视频m3u8| 亚洲精品第五页| 欧美日韩另类国产| 色九九视频| 国产农村精品一级毛片视频| 在线观看无码av免费不卡网站| 久久99热这里只有精品免费看 | 伊人AV天堂| 久久精品娱乐亚洲领先| 动漫精品中文字幕无码| 国产日韩丝袜一二三区| 在线观看免费人成视频色快速| 不卡无码h在线观看| 亚洲无码日韩一区| 国产日本欧美在线观看| 国产精品无码作爱| 热思思久久免费视频| 欧美日韩免费| 国产夜色视频| 91久久偷偷做嫩草影院| 久久情精品国产品免费| 久久久久亚洲av成人网人人软件| 欧美日韩一区二区三区四区在线观看 | 亚洲欧洲免费视频| 国模极品一区二区三区| 国产簧片免费在线播放| 久久激情影院| 国产尤物在线播放| 精品人妻一区二区三区蜜桃AⅤ| 国产天天色| 香蕉eeww99国产精选播放| 一区二区在线视频免费观看| 久久性妇女精品免费| 日韩无码黄色| 国产免费网址| 亚洲精品国产综合99久久夜夜嗨| 国内精品自在自线视频香蕉| 免费在线a视频| 国产特级毛片| 亚洲无码高清免费视频亚洲| 精品福利视频导航| 久久久久国产一级毛片高清板| 国产日韩欧美精品区性色| 亚洲精品欧美日韩在线| 欧美无专区| 成人蜜桃网| 国产第一页第二页| 成人在线观看不卡| 久久久精品无码一区二区三区| 67194成是人免费无码| 国产精品九九视频| 青青青国产精品国产精品美女| 国产激情第一页| 日韩亚洲综合在线| 波多野结衣中文字幕一区二区| 欧美精品亚洲二区| 黄片在线永久| 婷婷色一区二区三区| 91精品国产一区自在线拍| 亚洲系列无码专区偷窥无码| 丝袜无码一区二区三区| 欧美成人免费| 国产成人综合久久精品尤物| 国产欧美日韩va| 精品人妻系列无码专区久久| 国产人人乐人人爱| 亚洲欧美国产视频| 亚洲国产精品无码AV| 91无码国产视频| 国产三级视频网站| 久久午夜影院| 亚洲人妖在线| 久久黄色免费电影|