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

地基GNSS接收偽距信號中的日重復項

2022-02-17 14:00:12張帝鑫陳界宏蘇小寧孫楊軼孟國杰葉大綱
導航定位學報 2022年1期

張帝鑫,陳界宏,蘇小寧,孫楊軼,林 愷,徐 銳,孟國杰,葉大綱

地基GNSS接收偽距信號中的日重復項

張帝鑫1,陳界宏1,蘇小寧2,孫楊軼1,林 愷1,徐 銳3,孟國杰4,葉大綱5

(1. 中國地質大學(武漢) 地球物理與空間信息學院,湖北 武漢 430074;2. 蘭州交通大學 測繪與地理信息學院,甘肅 蘭州 730070;3. 四川省地震局,四川 成都 610041;4. 中國地震局地震預測研究所 地震預測重點實驗室,北京 100036;5. 臺北大學 不動產與城鄉環境學系,臺灣 新北 23741)

多路徑誤差每天重復出現,是存在于偽距信號中的日重復項。偽距中的多路徑誤差可通過濾波的方法提取,但在多路徑誤差信號和噪聲信號頻率未知或者重疊的情況下,不能將二者完全分開。本研究嘗試提出一種新技術解決這個問題,將偽距中衛星與臺站的距離、鐘差、硬件延遲、電離層延遲、對流層干延遲項消除后,得到包含多路徑誤差、對流層濕延遲和觀測噪聲的殘差數據,再通過主成分分析的方法提取多路徑誤差。通過主成分分析法得到重構的第一主成分每天都重復出現,其亦為偽距信號中的日重復項,理論上是多路徑誤差。然而,日重復項的數值大小是前人研究的4到5倍。這表明偽距信號中還包含了除多路徑誤差以外的其他日重復項。通過分析17個臺站的結果,本研究發現上述的其他日重復項,在使用相同類型接收機主機的臺站之間的相關性更高。這反映了該日重復項與接收機主機型號有關,可能是與接收機主機相關的系統誤差。

多路徑誤差;日重復項;主成分分析法;接收機主機誤差

0 引言

全球衛星導航系統(global navigation satellite system, GNSS)衛星信號在向地表發射并傳播到地面接收機的過程中,會受到衛星鐘差、衛星硬件延遲、電離層延遲、對流層延遲、地面多路徑誤差、接收機鐘差、接收機硬件延遲等影響[1]。這些因素均會影響定位精度,需要采取一定的方法和措施加以消除或削弱。鐘差和硬件延遲能通過差分技術消除[2]。電離層延遲和對流層延遲通過建立模型消除[3-4]。然而,多路徑誤差不僅不能通過差分技術消除,而且無法建立適用于每個臺站的改正模型。

多路徑誤差是由被臺站附近的建筑物、植被或者地面等所反射的衛星信號進入接收機天線,對直達衛星的信號產生干擾而產生。由于衛星軌道和臺站周圍環境基本不變,多路徑誤差每天重復出現,其是存在于信號中的日重復項。多路徑誤差嚴重時會引起信號的失鎖,是GNSS觀測中的一個重要的誤差源[1]。事后數據處理是削弱多路徑誤差主要采取的方法。多路徑誤差根據研究對象可分為三類,分別是載波相位、偽距和定位結果。載波相位上的多路徑誤差可通過建立單個臺站的數學模型減弱[5],或者通過分析信號的信噪比進行削弱[6-7]。偽距中多路徑誤差,可通過偽距和載波相位的線性組合進行評估[8]。但是,由于評估值中還包含了觀測噪聲,還需進一步處理,如通過有限脈沖響應(finite impulse response, FIR)濾波[9]或者自適應濾波[10]。在利用載波相位和偽距信號定位之前,未削弱原始信號中多路徑誤差,導致定位結果中受多路徑誤差的影響。對此,文獻[11]提出小波濾波法,文獻[12]采用自適應小波變換方法。

偽距上的多路徑誤差需要通過偽距和載波相位的線性組合評估后,再通過濾波提取。在觀測噪聲的頻率和多路徑誤差的頻率未知或者重疊的情況下,濾波的方法不能將二者完全分離。為了解決這個問題,本研究根據多路徑誤差每天重復出現的特征,嘗試提出一種新方法。在將偽距中的衛星與臺站距離、鐘差、硬件延遲、電離層延遲、對流層干延遲項消除后,得到包含多路徑誤差、對流層濕延遲、觀測噪聲的殘差數據,再根據多路徑誤差重復出現的特征,通過主成分分析法提取多路徑誤差。

1 原始數據

本研究使用了17個位于中國香港衛星定位參考站網的地基GNSS臺站的觀測數據。17個臺站的位置如圖1所示。

圖1 臺站分布圖(倒三角形、正三角形、圓形,正方形均表示臺站位置)

根據臺站觀測文件提供的信息可知,4個臺站使用TRIMBLE NETR9型號的天寶接收機主機(簡稱天寶主機)和TRM59800.00 SCIS/SCIT型號的天寶接收天線(簡稱天寶天線),11個臺站使用LEICA GR50型號的徠卡接收機主機(簡稱徠卡主機)和LEIAR25.R4 LEIT型號的徠卡接收天線(簡稱徠卡天線),另外T430臺站使用LEICA GR50型號的徠卡主機和TRM59800.00 SCIT型號的天寶天線,HKKS臺站使用TRIMBLE NETR9型號的天寶主機和TRM59800.00 SCIT型號的天寶天線,在3月30號接收主機更換為LEICA GR50型號的徠卡主機。本研究使用2018年德國地球科學研究中心的快速精密星歷文件,用于計算衛星與臺站之間的幾何距離。

2 多路徑誤差的提取

GNSS衛星信號在傳播過程中會受到諸多因素影響。其偽距的觀測方程[1]表達式為

本研究從原始的偽距觀測數據中,扣除星站距離、電離層延遲和對流層靜力延遲,并通過去線性趨勢減去鐘差、硬件延遲,得到包含多路徑誤差、對流層濕延遲、觀測噪聲的殘差數據,即研究數據。此處以HKOH臺站接收北斗衛星導航系統(BeiDou navigation satellite system, BDS)C03同步衛星的B1頻段信號為例展示。圖2(a)為原始偽距觀測數據。在多種攝動因素的影響下,同步衛星會產生周期性振蕩,在赤道面上下來回運動[21]。因此在一天的時間尺度內,偽距數據呈余弦變化規律。圖2(b)為星站距離,其數值與偽距接近,變化趨勢與偽距一致。圖2(c)為B1頻段的電離層延遲量,其是通過雙頻載波相位計算,是相對值。圖2(d)為對流層靜力延遲,其數值在2.5~2.6 m之間,一天內變化幅度較小。圖2(e)為偽距扣除星站距離、電離層延遲、對流層靜力延遲、鐘差和硬件延遲后,得到的包含多路徑誤差、對流層濕延遲、觀測噪聲的殘差數據。

圖2 HKOH臺站與BDS C03同步衛星之間的物理量

本研究中,臺站的觀測數據和氣象數據,以及快速精密星歷文件均被使用。任何一種數據缺失或者不完整,均無法求得研究數據。以2018年1月為例,31 d中有4 d數據不完整,如圖3所示。

圖3 2018年1月的殘差數據

在圖3中,每天數據都進行了等間距(3 m)縱向平移。部分波動特征每天都存在,并且隨著日期增加會偏移。這些重復出現的波動特征可能是由多路徑誤差產生的。研究表明,BDS C03同步衛星軌道振蕩周期與地球自轉周期一致,約為23 h 56 min 4 s[22]。與衛星軌道位置緊密相關的多路徑誤差周期與地球自轉周期一致[23]。與前一天相比,每天的多路徑誤差提前4 min左右出現。為了提取多路徑誤差,需要將每天的數據中相同的波動特征移動至相同時刻。考慮到GNSS臺站的采樣間隔為30 s,以起始第一天的時間為基準,將之后第天的時間依次增加(-1)×4 min,而后截取所有天數共有的時間段內的數據用于提取多路徑。圖4是圖3的數據經過移動后并截取的結果,數據長度為2672個歷元。

多路徑誤差每天都會重復出現,而對流層濕延遲和觀測噪聲具有隨機性。以連續若干天包含多路徑誤差、對流層濕延遲、觀測噪聲的殘差數據作為若干個原始指標,通過主成分分析法,可提取每天重復出現的部分,即多路徑誤差。主成分分析法是利用降維的思想,在損失最少信息的前提下把多個指標轉為幾個綜合指標的多元統計方法。其本質是利用原始變量的線性組合形成幾個綜合指標(主成分),在保存原始變量主要信息的前提下起到降維與簡化的作用。每一個主成分都是各原始變量的線性組合,并且每個主成分之間互不相關。在信號處理上,主成分分析法用于原始數據降噪。

圖4 2018年1月經時間移動并截取的殘差的數據

則矩陣的協方差矩陣為

可求得協方差矩陣的特征方程為

第主成分的貢獻率表示為

根據累計貢獻率選擇前個主成分,將原來的個變量轉變為個變量。在數據處理過程中,需

對前個主成分進行重構,得到一個近似的矩陣為

3 試驗結果

對HKOH臺站27 d的數據做主成分分析,其各個主成分的貢獻率如表1所示。

表1 主成分的貢獻率

從表1可以看出:第一至第五主成分的貢獻率分別為68.71 %、7.63 %、4.81 %、3.96 %和2.23 %;第6至第20主成分的貢獻率從2.17 %減小至0.15 %。判斷哪些主成分包含多路徑誤差的成分,需要根據重構后的各個主成分所體現特征辨別。圖5(a)和圖5(b)分別是通過縱向等間距(2 m)移動后的重構的第一主成分和第二主成分。

圖5中重構的第一主成分基本一致且重復出現,而重構的第二主成分在27 d里存在明顯差異。根據貢獻度大小可知,重構的其他主成分每天的差異更大,此處不做展示。因此,通過主成分分析法,從包含多路徑誤差、對流層濕延遲、觀測噪聲的殘差數據中提取的第一主成分,是偽距信號中的日重復項,每天重復出現,理論上為多路徑誤差。

圖5 一月份重構的主成分

本研究對2018年2月的數據進行與2018年1月數據相同的處理,得到2018年2月的第一主成分貢獻率為70.11 %,并且將第一主成分進行重構。圖6為1月和2月重構的第一主成分,圖6中每天的數據都經過縱向等間距(2 m)移動,其中2月份只有為20 d。隨著天數的增加,日重復項(重構的第一主成分)的波動特征每天較前一天提前出現,并且從1月持續到2月。若重構的第一主成分是真實多路徑誤差,這表明利用本研究的方法提取多路徑誤差是可行的。

圖6 1月份和2月份的重構的第一主成分

4 結果分析

為了驗證這個方法的可靠性,將結果與用前人的研究方法得到的結果比較。前人通過偽距和載波相位的線性組合,來評估偽距的多路徑誤差[8,23]。該方法忽略載波相位的多路徑誤差和觀測噪聲,利用線性因子的組合將星站距離等項消除,得到偽距中的多路徑誤差和觀測噪聲。此處將前人方法得到的偽距的多路徑誤差和觀測噪聲稱評估值。與評估值相比,通過主成分分析法得到的重構的第一主成分不含觀測噪聲,僅包含多路徑誤差。然而,重構的第一主成分分布在-2 m到3m之間,而評估值分布在-0.5 m到0.5 m之間(圖7)。重構的第一主成分是評估值的4到5倍。這表明在重構的第一主成分中,除了多路徑誤差這個日重復項,還包含了其他日重復項。

圖7 比較兩種方法的結果

由于對流層延遲和觀測噪聲的隨機性大,可排除這二者是日重復項的可能。在包含多路徑誤差、對流層濕延遲、觀測噪聲的殘差數據中,還可能包含之前扣除掉的物理量的殘留項。根據重構的第一主成分的數值變化幅度以及在2個月內會重復出現的特征,可排除鐘差、硬件延遲、對流層延遲和電離層延遲的殘差項的影響。而星站距離的殘留項有可能與日重復項有關。其中,星站距離的殘留項包含衛星和臺站天線的相位中心偏差和相位中心變化,以及臺站和衛星的坐標與真實坐標之間的偏差。其中,天線相位中心偏差屬于系統偏差,可當作固定值[24]。衛星和臺站天線的相位中心變化不會超過10 cm[25]。臺站坐標與真實值的偏差為常值,而通過插值后的衛星坐標與真實坐標的偏差,有可能是日重復項。

如果重構的第一主成分中的日重復項主要是由衛星坐標與真實坐標的偏差引起,則在中國香港地區的17個臺站都應有相似的情況。由于大部分臺站缺少氣象觀測數據,以下研究將偽距減去星站距離、電離層延遲、鐘差和硬件延遲,再減去多路徑誤差和觀測噪聲的評估值,從而得到包含對流層延遲的殘差數據。圖8展示了HKOH臺站2018年1月的包含對流層延遲的殘差數據(圖8中每天數據經過縱向3 m等間距平移)。

圖8 2018年1月的包含對流層延遲的殘差數據

重構的第一主成分中除多路徑誤差之外的日重復項亦存在于該殘差數據中。每天的殘差數據的變化幅度在4 m左右,而對流層延遲一天變化在分米級。因此,殘差值的大部分為該日重復項。依舊以BDS C03同步衛星的B1頻段偽距信號作為研究對象,對2018年每天任意兩個臺站的包含對流層濕延遲的殘差數據做皮爾遜(Pearson)相關性分析[26]。在數據齊全的情況下,任意兩個臺站每天都可計算出一個相關系數。表2為17個臺站中任意兩個臺站在2018年內相關系數的中位數,用于反映不同臺站之間殘差數據中日重復項的相關程度。

使用同為天寶的主機和天線的HKCL、HKTK、HKLM和HKQT臺站之間的相關系數基本大于0.40(表2中淺灰色部分)。使用同為徠卡的主機和天線的11個臺站,它們相互之間的相關系數基本分布在0.80到0.92之間(表2中深灰色部分)。然而,使用同為徠卡的主機和天線的臺站與使用同為天寶的主機和天線的臺站之間的相關系數基本分布在0.20左右。這表明包含對流層延遲的殘差數據中的日重復項在使用相同類型的主機和天線的臺站之間的相關系數,大于使用不同類型的主機和天線的臺站之間的相關系數。另外,在圖9(a)中,使用徠卡主機和天寶天線的T430臺站與使用同為徠卡的主機和天線的臺站的相關系數的中位數為0.87,而T430臺站與使用同為天寶的主機和天線的臺站的中位數為0.23。在圖9(b)中,HKKS臺站于年積日(day of year, DOY)第149天將天寶主機更換為徠卡主機。在更換主機前,HKKS臺站與其他使用天寶主機的臺站的相關系數的中位數為0.84,更換后相關系數降低為0.17。而與使用徠卡主機的臺站的相關系數的中位數在更換接收機后明顯變大,從原來的0.24變為0.93。上述結果表明,使用相同類型的主機的臺站的日重復項相關性更高,且與接收天線無關。

表2 任意兩個臺站在2018年內相關系數的中位數

注:表身中省略另一半重復的數值。

圖9 與其他16個臺站的相關系數

如果日重復項是衛星坐標與真實坐標的偏差,則17個臺站之間的相關系數應處于同一水平,而不是使用相同類型的主機的臺站的相關性更高。故日重復項并非由衛星坐標與真實坐標的偏差引起的。基于使用相同類型的主機的臺站的日重復項相關性更高的情況,本研究推斷包含對流層延遲的殘差數據中的日重復項與接收機主機類型有關。由于每個接收機主機的鐘差是獨立的,故不能引起上述結果。其可能是由于不同類型的接收機主機在處理信號的過程中的差異導致,為與接收機主機相關的系統誤差。

5 結束語

本研究對中國香港地區GNSS臺站接收的BDS C03同步衛星B1頻段偽距信號進行處理,得到包含多路徑誤差、對流層濕延遲、觀測噪聲的殘差數據,再利用主成分分析法從數據中提取每天重復出現的多路徑誤差。從HKOH臺站2018年1月和2月的數據中,提取的第一主成分在重構后每天重復出現,為偽距中的日重復項,符合多路徑誤差的特征。

然而,當將重構的第一主成分與前人研究比較時發現,重構的第一主成分的大小是前人研究的4~5倍。這可能是由除了多路徑誤差之外的日重復項引起。本研究通過數值變化及重復出現的特征,排除了對流層延遲、電離層延遲、鐘差和硬件延遲為該日重復項的可能性,并對星站距離殘留項進行檢驗。本研究通過對中國香港地區的17個臺站的殘差數據做相關性分析,驗證星站距離殘留項對17個臺站的影響是否相似。結果發現,使用同一類型的接收機主機的臺站之間的日重復項呈現高度相關,而使用不同類型的接收機主機的臺站之間的日重復項的相關系數比前者小,而且該結果與接收天線的類型無關。相關性分析的結果表明,日重復項與星站距離的殘留項無關,而是與接收機主機類型有關。該日重復項可能是與接收機主機相關的系統誤差。

致謝:感謝中國香港衛星定位參考站網和德國地球科學研究中心提供的數據。

[1] H?KANSSON M, JENSEN A B O, HOREMUZ M, et al. Review of code and phase biases in multi-GNSS positioning[J]. GPS Solutions, 2017, 21(3):849-860.

[2] 李政航, 黃勁松. GPS測量與數據處理[M]. 武漢: 武漢大學出版社, 2016:79-82.

[3] KLOBUCHAR J A. Ionospheric time-delay algorithm for single frequency GPS users[J]. IEEE Transactions on Aerospace and Electronic, 1987, 23(3): 325-331.

[4] HOPFIELD H S. Two-quartic tropospheric refractivity profile for correcting satellite data[J]. Journal of Geophysical Research, 1969, 74(18):4487-4499.

[5] GEORGIADOU Y, KLEUSBER A. On carrier signal multipath effects in relative GPS positioning[J]. Manuscripta Geodaetica, 1998, 13(3):172-179.

[6] 張波, 黃勁松, 蘇林. 利用信噪比削弱GPS多路徑效應的研究[J]. 測繪科學, 2003, 28(3):32-35.

[8] KLEUSBERG A, TEUNISSEN P. GPS for geodesy[M]. New York: Springer-Verlag, 1996:156-160.

[9] HAN S W, RIZOS C. GPS multipath mitgation using FIR filters[J]. Survey Review, 2000, 35(277):487-498.

[10] GE L, HAN S, RIZOS C. Multipath mitigation of continuous GPS measurements using an adaptive filter[J]. GPS Solutions, 2000, 4(2):19-30.

[11] 黃丁發,丁曉利, 陳永奇, 等. GPS多路徑效應影響與結構振動的小波濾波篩分研究[J].測繪學報, 2001, 30(2):36-41.

[12] ZHONG P, DING X L, ZHENG D W, et al. Adaptive wavelet transform based on cross-validation method and its application to GPS multipath mitigation[J]. GPS Solutions, 2008, 12(2):109-117.

[13] 王青平, 關玉梅, 王紫燕, 等. GPS衛星軌道三維坐標插值算法比較[J]. 地球物理學進展, 2014, 29(2):573-579.

[14] 周佩元, 杜蘭, 路余, 等. 多星定軌條件下北斗衛星鐘差的周期性變化[J].測繪學報, 2015,44(12):1299-1306.

[15] WILSON B D, MANNUCCI A J. Instrumental biases in ionospheric measurements derived from GPS data[C]//The Institute of Navigation (ION).Proceedings of the 6th International Technical Meeting of the Satellite Division of the Institute of Navigation. Salt Lake City: ION, 1993:1343-1351[2021-01-18].

[16] XUE J, SONG S, ZHU W. Estimation of differential code biases for Beidou navigation system using multi-GNSS observations: how stable are the differential satellite and receiver code biases?[J]. Journal of Geodynamics, 2016, 90(4):309-321.

[17] ELGERED G, DAVIS J L, HERRING T A, et al. Geodesy by radio interferometry: water vapor radiometry for estimation of the wet delay[J]. Journal of Geophysical Research, 1991, 96(B4): 6541-6555.

[18] 曲建光, 魏旭東, 王澤民, 等. 在高海拔地區Saastamoinen與Hopfield模型推算水汽含量差異的研究[J]. 武漢大學學報(信息科學版), 2003, 28(4):397-399.

[19] SAASTAMOINEN J. Atmospheric correction for troposphere and stratosphere in radio ranging of satellites[C]// American Geophysical Union (AGU). Proceedings of the Use of Artificial Satellites for Geodesy. Washington: AGU, 1972, 15:247-251[2021-01-18].

[20] BOEHM J, NIELL A, TREGONING P, et al.Global mapping function (GMF): a new empirical mapping function based on numerical weather model data[EB/OL].[2021-01-18].http://rses.anu.edu.au/geodynamics/gps/papers/boehms_ GMF.pdf.

[21] 華愛輝. 機動情況下的GEO衛星定軌方法研究[D]. 西安: 中國科學院研究生院(國家授時中心), 2008.

[22] 馮帥, 李亮, 程春, 等. 北斗空間信號誤差統計特征及描述方法研究[J].導航定位與授時, 2019, 6(3):95-102.

[23] WANG G X, ZHAO Q L, HU Z G, et al. Multipath analysis of code measurements for BeiDou geostationary satellites[J]. GPS Solutions, 2015, 19(1): 129-139.

[24] GU D F, LAI Y W, LIU J H, et al. Spaceborne GPS receiver antenna phase center offset and variation estimation for the Shiyan 3 satellite[J]. Chinese Journal of Aeronautics2016, 29(5): 1335-1344.

[25] SCHMID R, ROTHACHER M, THALLER D, et al. Absolute phase center corrections of satellite and receiver antennas[J]. GPS Solutions, 2005, 9(4): 283-293.

[26] BENESTY J, CHEN J, HUANG Y. Pearson correlation coefficient[M]//BENESTY J, KELLERMANN W. Noise reduction in speech processing. Berlin: Springer Nature Switzerland AG, 2009: 37-38[2021-01-18].

Daily repetition variations in pesudorange signals from ground-based GNSS receivers

ZHANG Dixin1, CHEN Chiehhung1, SU Xiaoning2, SUN Yangyi1, LIN Kai1, XU Rui3, MENG Guojie4, YEH Takang5

(1.Institute of Geophysics and Geomatics, China University of Geosciences (Wuhan), Wuhan 430074, China;2.Faculty of Geomatics, Lanzhou Jiaotong University, Lanzhou 730070, China;3.Earthquake Administration of Sichuan Province, Chengdu 610041, China;4.Institute of Earthquake Forecasting, Key Laboratory of Earthquake Forecast, China Earthquake Administration, Beijing 100036, China;5.Department of Real Estate and Built Environment, Taipei University, New Taipei 23741, China)

Multipath error are daily repetition variations, which appear in pesudorange every day. A filtering algorithm is a typical and efficient strategy to extract the multipath errors from the pesudorange signal. However, the straightforward method gets struggling if the frequency of multipath errors and noises are unknown or mixed up. This study attempts to propose a new technology to figure out this problem. We extract the multipath error from residual data by applying the principal component analysis. It should be noted that residual data containing multipath errors, tropospheric wet delay and noises, are obtained by pesudorange eliminating the distance between a satellite and a station, clock errors, hardware delay, ionospheric delay and tropospheric hydrostatic delay. Reconstructed first principal components obtained by principal component analysis have periodic characteristics called daily repetition variations, and they are multipath error theoretically. However, its value is 4 to 5 times bigger than previous studies. It indicates that there are other daily repetition variations besides multipath errors in pesudorange signals. We analyzed the data from 17 stations, and results show that the correlations of this daily repetition variations are higher among stations using the same type of receiver host. It suggests that the daily repetition variations in the pseudorange signals are associated with the type of receiver host, and it probably is a system error related to the receiver host.

multipath error; daily repetition variations; principal component analysis; receiver host

P228

A

2095-4999(2022)01-0076-09

張帝鑫,陳界宏,蘇小寧,等. 地基GNSS接收偽距信號中的日重復項[J]. 導航定位學報, 2022, 10(1): 76-84.(ZHANG Dixin, CHEN Chiehhung, SU Xiaoning, et al. Daily repetition variations in pesudorange signals from ground-based GNSS receivers[J]. Journal of Navigation and Positioning, 2022, 10(1): 76-84.)

10.16547/j.cnki.10-1096.20220111.

2021-03-01

國家自然科學基金委聯合基金項目(U2039205)。

張帝鑫(1995—),男,廣東惠州人,碩士研究生,研究方向為地基GNSS接收機的多路徑誤差。

陳界宏(1977—),男,臺灣臺北人,博士,教授,研究方向為震前多物理參量耦合。

主站蜘蛛池模板: 成人韩免费网站| 欧美精品xx| 欧美日韩免费观看| 国产精品吹潮在线观看中文| 久久九九热视频| 亚洲综合精品第一页| 欧美福利在线| 国产欧美精品专区一区二区| 色婷婷丁香| 九色视频线上播放| 中文纯内无码H| 无码中文字幕加勒比高清| 99热国产这里只有精品无卡顿"| 欧美在线视频不卡第一页| 欧美午夜精品| 亚洲精品va| 一级爆乳无码av| 一级爱做片免费观看久久| 538国产视频| 一本久道久综合久久鬼色| 欧美日韩久久综合| 欧美日本激情| 国产女人18毛片水真多1| www.亚洲一区二区三区| 日韩精品高清自在线| 欧美午夜小视频| 成人午夜免费观看| 91毛片网| 又爽又黄又无遮挡网站| 欧美国产在线看| 尤物成AV人片在线观看| 国产91精选在线观看| 99伊人精品| 欧美精品成人一区二区在线观看| 国产精品极品美女自在线网站| 国产精品福利社| 久久婷婷六月| AV熟女乱| 亚洲天堂在线免费| 激情午夜婷婷| 日本一区中文字幕最新在线| 免费观看亚洲人成网站| AV不卡无码免费一区二区三区| 国产十八禁在线观看免费| 久久久久中文字幕精品视频| 亚洲精品午夜天堂网页| 中国一级特黄视频| 亚洲另类色| 超碰aⅴ人人做人人爽欧美 | 久久伊伊香蕉综合精品| 欧美成人综合视频| 亚洲人成电影在线播放| 国产精品亚洲а∨天堂免下载| 91精品啪在线观看国产60岁| 粉嫩国产白浆在线观看| 午夜精品福利影院| 精品一区二区三区视频免费观看| 无码专区第一页| 丝袜久久剧情精品国产| 日韩午夜伦| 国产精品丝袜视频| 成人午夜久久| 国产肉感大码AV无码| 亚洲国产成人自拍| 国产91色在线| 国产精品嫩草影院视频| 国产新AV天堂| 国产精品久久久久鬼色| 久久婷婷六月| 亚洲黄色激情网站| 国产欧美中文字幕| 毛片网站在线播放| 91精品啪在线观看国产| 在线a网站| 中文字幕永久视频| 亚洲欧洲天堂色AV| 免费又爽又刺激高潮网址| 91www在线观看| 国产91高跟丝袜| 色窝窝免费一区二区三区| 国产成人高清亚洲一区久久| 欧美成人一级|