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

地球自轉(zhuǎn)參數(shù)的RFFLS短期預(yù)報(bào)算法研究

2017-10-12 03:21:11韓恒星黨亞民許長(zhǎng)輝
測(cè)繪通報(bào) 2017年9期
關(guān)鍵詞:模型

韓恒星,黨亞民,許長(zhǎng)輝

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

地球自轉(zhuǎn)參數(shù)的RFFLS短期預(yù)報(bào)算法研究

韓恒星1,2,黨亞民1,2,許長(zhǎng)輝2

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

地球自轉(zhuǎn)參數(shù)(ERP)是實(shí)現(xiàn)地心天球坐標(biāo)系(geocentric celestial reference system,GCRS)與國(guó)際地球坐標(biāo)系(international terrestrial reference system,ITRS)相互轉(zhuǎn)換的必要參數(shù),是國(guó)際GNSS服務(wù)組織(IGS)和國(guó)際GNSS監(jiān)測(cè)評(píng)估系統(tǒng)(iGMAS)分析中心的重要產(chǎn)品。本文針對(duì)最小二乘地球自轉(zhuǎn)參數(shù)預(yù)測(cè)算法會(huì)造成數(shù)據(jù)飽和以及新舊數(shù)據(jù)在數(shù)據(jù)處理及預(yù)報(bào)中被同等對(duì)待等問(wèn)題,將遺忘因子引入最小二乘預(yù)測(cè)算法,進(jìn)而提高ERP預(yù)報(bào)精度。遺忘因子遞推最小二乘算法能防止數(shù)據(jù)飽和,降低舊數(shù)據(jù)的影響,加強(qiáng)新數(shù)據(jù)的作用,降低在求解擬合參數(shù)時(shí)出現(xiàn)秩虧矩陣求逆的幾率,提高預(yù)報(bào)精度。本文詳細(xì)推導(dǎo)了遺忘因子遞推最小二乘表達(dá)式,探究了最佳遺忘因子,并通過(guò)ERP試驗(yàn)將該方法和原最小二乘的試驗(yàn)結(jié)果及LS-AR模型的預(yù)報(bào)結(jié)果作對(duì)比,發(fā)現(xiàn)僅用遺忘因子最小二乘模型預(yù)測(cè)就可以達(dá)到與LS-AR組合模型預(yù)測(cè)相當(dāng)?shù)木取?/p>

地球自轉(zhuǎn)參數(shù);預(yù)報(bào);遺忘因子;遞推最小二乘

Abstract: The earth’s rotation parameters(ERP)is an essential parameter for the conversion between the geocentric celestial reference system(GCRS)and the international terrestrial reference system(ITRS).It is the important products of the International GNSS Service Organization(IGS)and the International GNSS Monitoring and Evaluation System Analysis Center(iGMAS).In this paper,the least squares earth rotation parameter prediction algorithm will cause the data saturation and the old and new data in the data processing and forecasting are equally treated and so on,the forgetting factor into the least squares prediction algorithm,and thus improve the accuracy of ERP forecast.The forgetting factor recursive least squares algorithm can prevent data saturation,reduce the influence of old data,strengthen the function of new data,reduce the probability of inversion matrix of rank loss matrix when solving the fitting parameters,and improve the prediction accuracy.In this paper,the least squares expression of forgetting factor is deduced in detail,and the best forgetting factor is explored.The experimental results of this method and the results of LS-AR model are compared with those of LS-AR model.The least squares model prediction of the forgetting factor can achieve the same accuracy as the LS-AR model.

Keywords: earth rotation parameters;forecasting;forgetting factor;the recursive least squares

地球不僅反映了地球的整體運(yùn)動(dòng)狀態(tài),同時(shí)表征固體地球與地核、地幔、海洋及大氣在各種時(shí)空尺度上的相互作用過(guò)程,它們?cè)跊](méi)有外力的情況下構(gòu)成了復(fù)雜的地球動(dòng)力學(xué)系統(tǒng)。地球自轉(zhuǎn)變化可以用地球定向參數(shù)(EOP)來(lái)描述,地球定向參數(shù)包括地球自轉(zhuǎn)參數(shù)、歲差、章動(dòng),由于當(dāng)前的IAU2000/IAU2006歲差、章動(dòng)模型與實(shí)際的觀測(cè)符合很好,對(duì)于大多數(shù)應(yīng)用歲差、章動(dòng)資料采用理論模型計(jì)算值就可以滿足其要求,因此地球自轉(zhuǎn)參數(shù)就成為進(jìn)一步研究的重點(diǎn)主題。地球自轉(zhuǎn)參數(shù)(ERP)包括極移參數(shù)、日長(zhǎng)變化、UT1-UTC,極移參數(shù)包括X、Y兩個(gè)方向上的極移分量。ERP是衛(wèi)星精密定軌、深空探測(cè)、導(dǎo)航等領(lǐng)域的重要基礎(chǔ)數(shù)據(jù),由于空間測(cè)地技術(shù)獲得的ERP存在幾天到兩個(gè)星期的滯后性,空間導(dǎo)航等領(lǐng)域?qū)RP高精度預(yù)報(bào)值的需求與日俱增,ERP高精度預(yù)報(bào)已成為亟待解決的問(wèn)題。國(guó)際上常用的地球自轉(zhuǎn)參數(shù)預(yù)報(bào)方法有最小二乘、卡爾曼濾波、神經(jīng)網(wǎng)絡(luò)方法、協(xié)方差法等,使得地球自轉(zhuǎn)參數(shù)預(yù)報(bào)的精度逐漸提高。隨著航空航天技術(shù)和大地測(cè)量技術(shù)的發(fā)展,地球自轉(zhuǎn)參數(shù)的預(yù)報(bào)方法得到重大改進(jìn),預(yù)報(bào)精度也達(dá)到更高的水平。2012年徐君毅運(yùn)用基于截距修正的GM(1.1)模型進(jìn)行極移預(yù)測(cè)[1];2011年張昊等利用最小二乘與差分自回歸滑動(dòng)模型進(jìn)行極移的短期預(yù)報(bào)[2];2012年,許雪晴對(duì)地球定向參數(shù)預(yù)報(bào)高精度方法作了研究[3],并運(yùn)用AR模型間隔方式和迭代方式預(yù)報(bào)地球自轉(zhuǎn)參數(shù)[4];2012年王小輝等將經(jīng)驗(yàn)?zāi)J椒纸膺\(yùn)用在極移超短期預(yù)報(bào)中[5];2012年嚴(yán)鳳等對(duì)殘差序列進(jìn)行差分運(yùn)算之后,再采用自回歸模型進(jìn)行預(yù)報(bào)[6];2016年雷雨對(duì)地球自轉(zhuǎn)參數(shù)的高精度預(yù)報(bào)方法進(jìn)行了進(jìn)一步研究[7]。本文借鑒前輩的成果與經(jīng)驗(yàn),將遺忘因子遞推最小二乘應(yīng)用于地球自轉(zhuǎn)參數(shù)的預(yù)報(bào)中,探究最佳預(yù)報(bào)遺忘因子,進(jìn)一步提高預(yù)報(bào)穩(wěn)定性和預(yù)報(bào)精度,并將該改進(jìn)的模型投入實(shí)際應(yīng)用。

1 LS-AR模型

現(xiàn)代研究表明,極移的時(shí)間尺度變化主要有:長(zhǎng)趨勢(shì)項(xiàng)、Chandler項(xiàng)、周年項(xiàng)、半年項(xiàng)及高頻極移。其中高頻極移變化較大且無(wú)一定規(guī)律性,因此在構(gòu)造最小二乘擬合模型時(shí),包含的固定線性項(xiàng)和周期項(xiàng)是:線性趨勢(shì)項(xiàng)、Chandler項(xiàng)、周年項(xiàng)及半周年項(xiàng),通過(guò)對(duì)一段時(shí)間的極移序列進(jìn)行最小二乘擬合并求解模型參數(shù),以達(dá)到預(yù)測(cè)的目的。其模型公式為

X(t)=a0+a1t+a2cos(2πt/T1)+a3sin(2πt/T1)+a4cos(2πt/T2)+a5sin(2πt/T2)+a6cos(2πt/T3)+a7sin(2πt/T3)

Y(t)=b0+b1t+b2cos(2πt/T1)+b3sin(2πt/T1)+b4cos(2πt/T2)+b5sin(2πt/T2)+b6cos(2πt/T3)+b7sin(2πt/T3)

式中,a0為常數(shù)項(xiàng);t為UTC時(shí);a1為線性趨勢(shì)項(xiàng)的系數(shù);a2、a3為半年項(xiàng)的系數(shù);a4、a5為周年項(xiàng)的系數(shù);a6、a7為錢德勒周期項(xiàng)的系數(shù)。T1為半年項(xiàng)的周期;T2為周年項(xiàng)的周期;T3為錢德勒項(xiàng)的周期。

AR模型表示平穩(wěn)隨機(jī)序列xt(t=1,2,…,n)在t時(shí)刻以前的規(guī)律性變化和t時(shí)刻白噪聲的關(guān)系,其數(shù)學(xué)模型為

式中,φ1、φ2、…、φP為模型參數(shù);at為白噪聲;p為模型階數(shù);at~N(0,σ),σ為白噪聲的方差。式(3)稱為p階自回歸模型,簡(jiǎn)記為AR(P)。

從最小二乘求擬合參數(shù)的矩陣求解過(guò)程中可以發(fā)現(xiàn)最小二乘有以下不足:

(1) LS模型進(jìn)行ERP預(yù)報(bào)時(shí)發(fā)現(xiàn)在反復(fù)求解擬合參數(shù)過(guò)程中對(duì)矩陣求逆時(shí)會(huì)產(chǎn)生秩虧。

(2) 每有一次新數(shù)據(jù)的加入都要重新求逆矩陣,工作量大,且矩陣秩虧的頻率會(huì)增加,影響預(yù)報(bào)進(jìn)程。

(3) 在預(yù)報(bào)過(guò)程中新舊數(shù)據(jù)同等對(duì)待,然而經(jīng)試驗(yàn)發(fā)現(xiàn)通過(guò)對(duì)新舊數(shù)據(jù)設(shè)置不同的權(quán)重得出的預(yù)報(bào)結(jié)果是有很大的差異的。

根據(jù)最小二乘這些不足,引入遺忘因子遞推最小二乘來(lái)解決以上問(wèn)題。

2 遺忘因子遞推最小二乘

通常情況下,最小二乘準(zhǔn)則意義下的動(dòng)態(tài)模型為

U(k)=V(k)B(k)+e(k)

式中,V(k)為因子數(shù)m的前k個(gè)因子向量X(i),i=1,2,…,k組成的矩陣;B(k)為前k個(gè)模型系數(shù)組成的矩陣;U(k)為前k個(gè)預(yù)報(bào)向量Z(k),i=1,2,…,k,組成的矩陣或向量,e(k)為前k個(gè)隨機(jī)擾動(dòng)組成的矩陣或者向量。在因子矩陣和預(yù)報(bào)向量矩陣中,加上衰減因子β(0<β≤1)后,因子矩陣和預(yù)報(bào)向量矩陣記作

V*(k)=[βk-1X(1),βk-2X(2),…,βX(k-1),X(k)]T

U*(k)=[βk-1Z(1),βk-2Z(2),…,βZ(k-1),Z(k)]T

通常定義

P(k)=[V*(k)TV*(k)]-1

稱之為逆相關(guān)矩陣。如果因子數(shù)為m,那么它是一個(gè)m×m的方陣。由于最小二乘方法很容易得到B(k)的估計(jì)

接著可以得到

式中,μ=β2,0<μ≤1,μ稱為遺忘因子,其大小可表征遺忘速度。當(dāng)遺忘因子μ=1時(shí),考慮遺忘因子線性自適應(yīng)建模算法就是普通的線性自適應(yīng)建模算法。顯然可以得到

P-1(k)=μP-1(k-1)+X(k)TX(k)

然后有

通常令

Γ(k)=P(k)X(k)T

因此,上式可以寫(xiě)成

另外利用矩陣反演公式改寫(xiě)成

最終可得到

Γ(k)=P(k-1)X(k)T[X(k)P(k-1)×X(k)T+μI]-1

式中,P(k)為系數(shù)估計(jì)誤差方陣,是一對(duì)稱矩陣;Γ(k)為卡爾曼增益矩陣。

3 試驗(yàn)及試驗(yàn)結(jié)果分析

本文試驗(yàn)數(shù)據(jù)皆來(lái)自北斗分析中心(CGS),選取2017年儒略日為57 835—57 847之間的地球自轉(zhuǎn)參數(shù)超快速部分作為試驗(yàn)基礎(chǔ)數(shù)據(jù),試驗(yàn)以極移和日長(zhǎng)變化3項(xiàng)為主,利用遺忘因子最小二乘遞推方法進(jìn)行ERP預(yù)報(bào)。

經(jīng)過(guò)大量試驗(yàn)找出了適合預(yù)報(bào)的最佳遺忘因子,利用該最佳遺忘因子進(jìn)行極移兩個(gè)方向(PMX和PMY)及日長(zhǎng)變化的預(yù)報(bào)。試驗(yàn)結(jié)果及分析如圖1所示。

圖1 最佳遺忘因子試驗(yàn)探究

一般遺忘因子的取值必須選擇接近于1的整數(shù),一般不小于0.9。如果是線性系統(tǒng)且時(shí)變特征明顯,該因子的選擇應(yīng)該是更大點(diǎn)的值。由于遺忘因子選取的不確定性,進(jìn)行了因子探究試驗(yàn),從圖1中可以看出,遺忘因子為0.95時(shí)得到了精度最高的預(yù)報(bào)結(jié)果,也就是說(shuō)最佳的遺忘因子為0.95。

利用探究出的最佳遺忘因子,確定了預(yù)報(bào)最終的遺忘因子遞推最小二乘模型,將北斗分析中心的數(shù)據(jù)利用最小二乘模型(LS)、最小二乘和自回歸組合模型(LS-AR)、遺忘因子遞推最小二乘模型(RFFLS)分別進(jìn)行試驗(yàn)預(yù)報(bào),將預(yù)報(bào)結(jié)果進(jìn)行分析。

從圖2中就可以看出,在日長(zhǎng)變化的預(yù)報(bào)中LS模型的預(yù)報(bào)結(jié)果稍有起伏,當(dāng)加入AR模型以后精度沒(méi)有大的改善但是穩(wěn)定性有所增強(qiáng),應(yīng)用遺忘因子最小二乘模型的預(yù)報(bào)精度中期精度稍差,不如LS-AR模型的預(yù)報(bào)精度。

圖2 日長(zhǎng)變化不同模型預(yù)報(bào)結(jié)果對(duì)比

圖3 極移X方向不同模型預(yù)報(bào)結(jié)果對(duì)比

從圖3中可以看出,在極移X方向上,LS的預(yù)報(bào)結(jié)果上下起伏較大,經(jīng)過(guò)AR模型進(jìn)行殘差擬合預(yù)報(bào)后,發(fā)現(xiàn)LS-AR模型得出的結(jié)果在穩(wěn)定性和精度上明顯優(yōu)于LS模型,遺忘因子最小二乘模型的預(yù)報(bào)精度前期的精度稍差,但優(yōu)于LS預(yù)報(bào)模型的精度,在后期的精度與LS-AR模型精度相當(dāng),甚至優(yōu)于LS-AR模型。

圖4 極移Y方向不同模型預(yù)報(bào)結(jié)果對(duì)比

從圖4中可以看出,在極移Y方向上,LS單獨(dú)的預(yù)報(bào)精度變化大,經(jīng)過(guò)AR模型進(jìn)行殘差擬合預(yù)報(bào)后,LS-AR模型得出的結(jié)果在穩(wěn)定性和精度上優(yōu)于LS模型,遺忘因子最小二乘模型的預(yù)報(bào)精度個(gè)別點(diǎn)精度稍差,但優(yōu)于LS預(yù)報(bào)模型的精度,總體的預(yù)報(bào)精度與LS-AR模型精度相當(dāng)。

4 結(jié) 論

經(jīng)過(guò)以上預(yù)報(bào)試驗(yàn)與分析得出以下結(jié)論:

(1) 在ERP預(yù)報(bào)中LS模型單獨(dú)預(yù)報(bào)的結(jié)果精度稍差且預(yù)報(bào)結(jié)果不穩(wěn)定,當(dāng)加入了AR模型進(jìn)行殘差擬合預(yù)報(bào)后,精度和穩(wěn)定性都有不同程度的提高,在極移的預(yù)報(bào)上精度提高最為明顯。

(2) 將遺忘因子遞推最小二乘應(yīng)用于ERP預(yù)報(bào)時(shí),在極移預(yù)報(bào)中的作用較為明顯,與LS模型單獨(dú)進(jìn)行預(yù)報(bào)精度有一定程度的提高,與LS-AR組合模型的預(yù)報(bào)精度相當(dāng),在日長(zhǎng)變化的預(yù)報(bào)中作用不明顯。

(3) 日長(zhǎng)變化和極移對(duì)不同預(yù)報(bào)模型的敏感性不同,在進(jìn)行極移和日長(zhǎng)變化預(yù)報(bào)過(guò)程中,選取預(yù)報(bào)模型時(shí)應(yīng)根據(jù)不同的預(yù)報(bào)類型、不同特性加以區(qū)分。

[1] 徐君毅.基于截距修正的GM(1.1)模型在極移預(yù)測(cè)中的應(yīng)用[J].大地測(cè)量與地球動(dòng)力學(xué),2010,30(1):88-91.

[2] 張昊,王琪潔,朱建軍.加權(quán)最小二乘與AR組合模型在極移預(yù)測(cè)中的應(yīng)用研究[J].天文學(xué)進(jìn)展,2011,29(3):343-352.

[3] 許雪晴.地球定向參數(shù)高精度預(yù)報(bào)方法研究[D].上海:中國(guó)科學(xué)院上海天文臺(tái),2012.

[4] 許雪晴,周永宏.AR模型間隔方式和迭代方式預(yù)報(bào)地球自轉(zhuǎn)參數(shù)對(duì)比[J].中國(guó)科學(xué)院上海天文臺(tái)刊,2012(33):20-27.

[5] 王小輝,王琪潔,劉建.經(jīng)驗(yàn)?zāi)J椒纸庠跇O移超短期預(yù)報(bào)中的應(yīng)用[J].天文學(xué)報(bào),2012,53(3):519-526.

[6] 嚴(yán)鳳,姚宜斌,地球自轉(zhuǎn)參數(shù)短期預(yù)報(bào)方法及實(shí)現(xiàn)[J].大地測(cè)量與地球動(dòng)力學(xué),2012(4):71-75.

[7] 雷雨.地球自轉(zhuǎn)參數(shù)的高精度預(yù)報(bào)方法研究[D].北京:中國(guó)科學(xué)院大學(xué),2016.

[8] SU X,LIU L,HOUTSE H,et al.Long-term Polar Motion Prediction Using Normal Time-frequency Transform[J].Journal of Geodesy,2013,88(2):145-155.

[9] 張昊,王琪潔,朱建軍,等.對(duì)錢德勒參數(shù)進(jìn)行時(shí)變修正的CLS-AR模型在極移預(yù)測(cè)中的應(yīng)用[J].武漢大學(xué)學(xué)報(bào)(信息科學(xué)版),2012,37(3):286-289.

[10] 張龍平,黨亞民,成英燕,等.北斗GEO/IGSO/MEO衛(wèi)星定軌地面站構(gòu)型影響分析及其優(yōu)化[J].測(cè)繪學(xué)報(bào),2016,45(S2):82-92.

[11] 王宇譜,呂志平.小波神經(jīng)網(wǎng)絡(luò)日常預(yù)報(bào)算法研究[J].大地測(cè)量與地球動(dòng)力學(xué),2012,32(1):127-131.

[12] 王琪潔.神經(jīng)網(wǎng)絡(luò)的地球自轉(zhuǎn)變化預(yù)報(bào)[D].上海:中國(guó)科學(xué)院上海天文臺(tái),2007.

[13] 王琪潔,廖德春,周永宏,等.日長(zhǎng)變化預(yù)報(bào)中BP神經(jīng)網(wǎng)絡(luò)拓?fù)浣Y(jié)構(gòu)的選擇[J].中國(guó)科學(xué)院上海天文臺(tái)刊,2008,49(1):93-100.

[14] 黨亞民,秘金鐘,成英燕.全球?qū)Ш叫l(wèi)星系統(tǒng)原理與應(yīng)用[M].北京:測(cè)繪出版社,2007.

[15] SCHUH H,ULRICH M,EGGER D,et al.Prediction of Earth Orientation Parameters by Artificial Neural Networks[J].Journal of Geodesy,2002,76(5):247-258.

[16] KOSEK W,MCCARTHY D D,LUZUM B J.El Nio Impact on Polar Motion Prediction Errors[J].Studia Geophysica et Geodaetica,2001,45(4):347-361.

Short-termForecastingofEarthRotationParameterBasedonForgettingFactorRecursiveLeastSquares

HAN Hengxing1.2,DANG Yamin1,2,XU Changhui2

(1. Shandong University of Science and Technology,Qingdao 266590,China; 2. Chinese Academy of Surveying & Mapping, Beijing 100830,China)

P228

A

0494-0911(2017)09-0015-04

2017-05-09;

2017-07-19

國(guó)家自然科學(xué)基金(41474011);國(guó)家重點(diǎn)研發(fā)計(jì)劃(2016YFB0501405);公益性行業(yè)專項(xiàng)(B1503);國(guó)家基礎(chǔ)測(cè)繪科技項(xiàng)目(2017KJ0205);中國(guó)第二代衛(wèi)星導(dǎo)航系統(tǒng)重大專項(xiàng)(GFZX0301040308-06)

韓恒星(1990—),男,碩士生,主要研究方向?yàn)榈厍蜃赞D(zhuǎn)參數(shù)(ERP)計(jì)算與預(yù)報(bào)。E-mail:496426305@qq.com

韓恒星,黨亞民,許長(zhǎng)輝.地球自轉(zhuǎn)參數(shù)的RFFLS短期預(yù)報(bào)算法研究[J].測(cè)繪通報(bào),2017(9):15-18.

10.13474/j.cnki.11-2246.2017.0278.

猜你喜歡
模型
一半模型
一種去中心化的域名服務(wù)本地化模型
適用于BDS-3 PPP的隨機(jī)模型
提煉模型 突破難點(diǎn)
函數(shù)模型及應(yīng)用
p150Glued在帕金森病模型中的表達(dá)及分布
函數(shù)模型及應(yīng)用
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
3D打印中的模型分割與打包
主站蜘蛛池模板: AV在线天堂进入| 一级毛片在线播放| 亚洲六月丁香六月婷婷蜜芽| 91精品啪在线观看国产60岁| 国产粉嫩粉嫩的18在线播放91| 色天天综合久久久久综合片| 很黄的网站在线观看| 2021精品国产自在现线看| 国产成年无码AⅤ片在线| 亚洲欧美综合另类图片小说区| 玖玖精品视频在线观看| 久久网综合| 精品第一国产综合精品Aⅴ| 免费视频在线2021入口| 久久婷婷人人澡人人爱91| 国产欧美成人不卡视频| 亚洲欧美一区二区三区麻豆| 亚洲大尺度在线| 日韩不卡高清视频| 色综合成人| 666精品国产精品亚洲| 伊人久久大香线蕉综合影视| 亚洲精品无码av中文字幕| 蜜芽国产尤物av尤物在线看| 国产h视频在线观看视频| 国产波多野结衣中文在线播放| 国产性爱网站| 91亚瑟视频| 欧美在线天堂| 强乱中文字幕在线播放不卡| 国产真实乱子伦精品视手机观看| 欧美日韩亚洲国产| 婷婷丁香在线观看| 国产成在线观看免费视频| 一本久道久久综合多人| 久久精品人人做人人爽电影蜜月| 亚洲中文字幕久久精品无码一区| 99精品国产高清一区二区| 成人av专区精品无码国产| 91区国产福利在线观看午夜| 9丨情侣偷在线精品国产| 一本大道视频精品人妻| 少妇精品在线| 国产精品污污在线观看网站| 国产肉感大码AV无码| 72种姿势欧美久久久久大黄蕉| 在线欧美日韩国产| 国产麻豆精品久久一二三| 国产幂在线无码精品| 婷婷午夜天| 久久久久人妻精品一区三寸蜜桃| 高潮毛片免费观看| 99免费视频观看| 在线看片中文字幕| 麻豆AV网站免费进入| 最新国产你懂的在线网址| 欧美成人一级| 在线免费不卡视频| 成人一级免费视频| 成人国产精品一级毛片天堂| 日本欧美中文字幕精品亚洲| 人妻91无码色偷偷色噜噜噜| 国产微拍精品| 天堂成人在线视频| 中文字幕 日韩 欧美| 黄色成年视频| 国产国模一区二区三区四区| 国产后式a一视频| 日韩 欧美 小说 综合网 另类| 久久精品一卡日本电影 | www.av男人.com| 性喷潮久久久久久久久| 性欧美在线| 亚洲区欧美区| 久久国语对白| 日韩欧美中文字幕在线韩免费| 日韩精品专区免费无码aⅴ| 欧美国产视频| 色婷婷综合激情视频免费看| 亚洲国产日韩一区| 人妻一区二区三区无码精品一区| 一本无码在线观看|