• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    GRACE重力反演中的軌道數(shù)值積分方法分析

    2017-03-31 02:21:59范東明谷延超
    宇航學(xué)報 2017年3期
    關(guān)鍵詞:歷元階數(shù)步長

    余 彪,范東明,游 為,谷延超,蘇 勇

    (1.西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院,成都611756;2.西南石油大學(xué)土木工程與建筑學(xué)院,成都610500)

    GRACE重力反演中的軌道數(shù)值積分方法分析

    余 彪1,范東明1,游 為1,谷延超1,蘇 勇2

    (1.西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院,成都611756;2.西南石油大學(xué)土木工程與建筑學(xué)院,成都610500)

    首先,對比分析Runge-Kutta積分法、Adams-Cowell聯(lián)合并行積分法、Gauss-Jackson積分法以及天體力學(xué)中常用的外插積分法用于GRACE衛(wèi)星軌道積分及變分方程解算時的優(yōu)缺點,得出Gauss-Jackson算法優(yōu)于其它算法并推薦了相對較優(yōu)的階數(shù)和步長,從算法的推導(dǎo)過程中分析出了Gauss-Jackson算法優(yōu)于其它算法的原因。其次,對比不同積分方法抵抗誤差的能力,結(jié)果表明這些方法都不能有效抵抗誤差的干擾;并分析初始狀態(tài)向量誤差和攝動力誤差對軌道積分的影響,結(jié)果表明,衛(wèi)星初始狀態(tài)向量中的速度誤差(0.1mm/s)對軌道積分的影響大于位置誤差(10mm)對軌道積分的影響,攝動力中的隨機(jī)誤差對軌道積分的影響較大且無規(guī)律可循。針對Gauss-Jackson算法中的不足,提出基于移動窗口的多項式內(nèi)插算法和Gauss-Jackson算法相結(jié)合的組合方法(改進(jìn)Gauss-Jackson算法),通過模擬和實測數(shù)據(jù)的計算表明該方法在保證積分精度的前提下不僅提高了積分效率而且可以得到任意時刻的積分值。

    軌道積分;積分精度;積分效率;誤差影響

    0 引 言

    在利用重力恢復(fù)與氣候?qū)嶒?Gravity recovery and climate experiment, GRACE)衛(wèi)星數(shù)據(jù)反演靜態(tài)或時變地球重力場模型以及下一代重力衛(wèi)星(如NGGM: Next generation gravity measurements和GRACE-FO:GRACE follow-on)反演重力場的模擬計算中,需要用數(shù)值積分方法計算衛(wèi)星參考軌道和動力學(xué)反演方法中的變分方程等。而高精度、高效率和高穩(wěn)定性的數(shù)值積分方法對于地球重力場模型的計算至關(guān)重要。

    目前,廣泛使用的數(shù)值積分方法主要有Runge-Kutta積分法(簡稱RK方法)、Adams-Cowell聯(lián)合并行積分法(簡稱AC方法)、Gauss-Jackson積分法(簡稱GJ方法)、天體力學(xué)中常用的外插積分法(簡稱EX方法)和Collocation積分法(配置法)[1-9]。不同的積分方法優(yōu)缺點不同[1-2],單步法(如RK方法)可以自起步積分,但是計算效率較低;多步法(如AC方法,GJ方法)雖然計算效率高于單步法,但是無法自起步;Collocation積分法一般用于受力相對簡單的高軌衛(wèi)星軌道積分,如國際上著名的GNSS定軌定位軟件Bernese就是采用這種方法,但低軌GRACE衛(wèi)星,受力比較復(fù)雜,實際計算時并不能很好地模型化衛(wèi)星受力情況,采用此種方法進(jìn)行軌道積分時積分精度不如其它方法,所以本文未對其進(jìn)行對比分析。即使同一種積分方法,采用不同的階數(shù)和步長進(jìn)行計算時也會得到不同精度的結(jié)果。在對衛(wèi)星進(jìn)行軌道積分時,衛(wèi)星初始狀態(tài)向量及攝動力中含有各種誤差,為了得到高精度和高穩(wěn)定性的積分結(jié)果需要分析誤差對軌道積分的影響。低軌衛(wèi)星軌道積分往往效率比較低,而且一般只能得到歷元為步長的整數(shù)倍處的軌道狀態(tài),若采用GPS解算的幾何軌道作為標(biāo)準(zhǔn)軌道進(jìn)行地球重力場模型反演,則可能需要將歷元為非步長整數(shù)倍處的軌道歸算到歷元為步長整數(shù)倍處的軌道。

    因此本文首先從軌道積分精度、積分效率及狀態(tài)轉(zhuǎn)移矩陣積分精度方面對比分析現(xiàn)有的常用數(shù)值積分方法,然后分析了初始狀態(tài)向量和攝動力誤差對軌道積分的影響。針對目前軌道積分算法的不足,基于本文的研究提出了改進(jìn)GJ算法,該算法采用適當(dāng)階數(shù)和較大步長的GJ算法,保證積分精度并結(jié)合基于移動窗口的多項式內(nèi)插算法計算出步長范圍內(nèi)任意歷元處的軌道。由于采用了較大的積分步長,使得積分效率相比于采用5s步長的積分效率成倍增加,而積分精度并不會降低;由于采用了基于移動窗口的多項式內(nèi)插算法可以內(nèi)插出任意時刻的積分軌道,降低了不規(guī)則歷元處軌道處理的難度。

    1 數(shù)值積分方法

    RK方法是最常用的單步積分法且存在多種不同算法[1],即便是同一階方法也存在不同的算法且算法之間計算精度相差較大,本文采用軌道積分精度較高的8階RK算法[5-6,8-9],8階RK積分方法每積分一步都需計算10次右函數(shù),對于GRACE衛(wèi)星實際觀測數(shù)據(jù),計算一次右函數(shù)(主要指衛(wèi)星所受的所有攝動力)大約需要耗時0.1s(在配有Intel Core四核處理器,主頻為2.50GHz的PC機(jī)上采用Fortran語言計算),這對于實際軌道積分計算效率低下。

    EX積分法是單步積分法[1,4],其計算過程分為兩步,首先以不同步長計算某一歷元處的近似值然后對近似值外插逐漸逼近軌道真值。通過計算發(fā)現(xiàn),對于GRACE衛(wèi)星軌道,三角EX法[1]比菱形EX法[4]具有更高的數(shù)值穩(wěn)定性,所以本文采用三角EX算法。

    AC方法為多步積分法,計算效率高于單步法,為了提高計算精度,本文在計算中結(jié)合預(yù)估校正迭代算法[1-3,5],并與單步法的計算效率進(jìn)行對比,本文在計算中采用8階RK方法提供初值。

    GJ算法也為多步積分法,為了起到更好的對比效果,GJ方法起步值的計算并未采用文獻(xiàn)[3]和[10]中介紹的Mid-Corrector公式,而采用與AC方法相同的起步算法(8階RK方法),計算中同樣結(jié)合預(yù)估校正迭代算法。

    2 不同方法軌道和狀態(tài)轉(zhuǎn)移矩陣積分對比

    通過比較不同方法的軌道和狀態(tài)轉(zhuǎn)移矩陣積分精度,得出最優(yōu)方法;比較同一方法不同階數(shù)和步長的計算結(jié)果得出計算精度較高的階數(shù)和步長范圍。目前GRACE重力反演是分弧段進(jìn)行處理的,弧段長度一般定為一天,但衛(wèi)星繞地球運行時受地球引力的作用是連續(xù)累積的,從理論上說采用連續(xù)的弧段而不進(jìn)行分弧段處理有益于地球重力場的反演;并且由于目前GRACE時變地球重力場模型的時間尺度一般為一個月,本文在積分時直接采用一個月作為積分的時間跨度。

    為了得到軌道的解析解,考慮衛(wèi)星只受地球中心引力的二體問題。根據(jù)開普勒軌道理論計算衛(wèi)星軌道解析解,將不同積分方法以不同的階數(shù)和步長計算積分軌道并與開普勒軌道比較得到積分誤差。采用表1中的GRACE衛(wèi)星軌道參數(shù)作為初始?xì)v元的開普勒軌道根數(shù),并以不同步長計算得到一個月的開普勒軌道。

    表1 GRACE衛(wèi)星初始軌道根數(shù)

    2.1 衛(wèi)星軌道積分誤差對比

    對于EX方法,外插計算一次算法階數(shù)增加2[1],所以該方法只有偶數(shù)階。當(dāng)選取的階數(shù)很高時,積分軌道仍收斂,但是該算法的積分效率會隨著階數(shù)的增加迅速降低(第2.2節(jié)詳細(xì)說明),考慮到實際計算效率本文只分析前20階的計算精度。圖1(b)表示不同階數(shù)和步長的EX方法積分軌道誤差。

    對于AC方法,滿足積分精度要求的階數(shù)為8階到13階。圖1(c)表示不同階數(shù)和步長的AC方法積分軌道誤差。從圖1(c)可以看出,該方法在步長小于15s時,采用任何一種階數(shù)進(jìn)行積分都無法滿足mm級的誤差要求(如常用的12階5s步長積分一個月的最大誤差為40.4mm)。

    對于GJ方法,滿足積分精度要求的階數(shù)為3階到12階。圖1(d)表示不同階數(shù)和步長的GJ方法積分軌道誤差。從圖1(d)可以看出,GJ方法積分誤差普遍比上述幾種方法小(可達(dá)兩個量級),而且該方法3到12階5s步長的積分誤差皆滿足要求。

    對比分析圖1的結(jié)果可以得出以下兩點結(jié)論:1)對于衛(wèi)星軌道積分并不是積分步長越短越好、積分階數(shù)越高越好。其可能的原因為,當(dāng)積分步長較小時,積分計算的歷元個數(shù)增多,則軌道積分累計誤差或計算機(jī)的累計舍入誤差會增大;當(dāng)積分步長較大時,用于內(nèi)插右函數(shù)值的弧段較長,則內(nèi)插出的右函數(shù)值與實際右函數(shù)值的差值較大[11],從而影響軌道積分的精度。當(dāng)積分的階數(shù)較大時,積分下一個歷元時所用到的前面歷元的個數(shù)增多,則用于內(nèi)插右函數(shù)的弧段增長,內(nèi)插出的右函數(shù)誤差可能增大[11];而當(dāng)階數(shù)較小時,則用于內(nèi)插右函數(shù)的多項式的階數(shù)小,右函數(shù)的插值余項變大[11],即右函數(shù)的近似誤差大。所以積分步長和積分階數(shù)應(yīng)根據(jù)精度、效率和觀測數(shù)據(jù)采樣要求合理選擇。2)從不同方法軌道積分誤差角度看,GJ方法優(yōu)于其它方法。通過對比分析各積分算法的推導(dǎo)過程,以GJ算法和AC算法為例來說明GJ算法優(yōu)于其它算法的原因。GJ算法與AC算法同為多步法,而且公式推導(dǎo)過程相似。兩套算法都是基于基本數(shù)學(xué)物理公式:

    (1)

    (2)

    之前可認(rèn)為都是一樣的。AC算法將式(2)化為:

    (3)

    GJ算法將式(2)轉(zhuǎn)化為:

    (rn+1-rn)-(rn-rn-1)=▽2rn+1=

    (4)

    2.2 衛(wèi)星軌道積分效率對比

    在衛(wèi)星軌道積分時,除了要考慮積分誤差,還應(yīng)該考慮積分效率,影響軌道積分效率的主要因素是右函數(shù)的計算次數(shù)。一般情況下,單步法的計算效率低于多步法,但為了提高多步法的計算精度,常常結(jié)合預(yù)估校正迭代算法,本節(jié)主要比較結(jié)合預(yù)估校正迭代算法的多步法和單步法的計算效率。

    對于8階RK方法,積分一步大約需要花費1s的時間。圖2(a)反應(yīng)了不同步長8階RK方法計算一個月的積分軌道時右函數(shù)計算次數(shù)。

    對于EX方法,該方法需要以不同的步長計算某一歷元處衛(wèi)星的狀態(tài)得到EX方法的初值,在計算初值的過程中需要多次計算右函數(shù)。圖2(b)反應(yīng)的是不同步長和階數(shù)的EX方法計算一個月的積分軌道時右函數(shù)計算次數(shù)。從圖2(b)可以看出,EX方法右函數(shù)的計算次數(shù)會隨著階數(shù)的增高成倍的增加,而且相比于同階的RK方法,EX方法右函數(shù)的計算次數(shù)明顯要更多。

    對于AC聯(lián)合并行預(yù)估校正迭代算法,雖然需要進(jìn)行迭代計算,但在對GRACE衛(wèi)星軌道的積分計算中發(fā)現(xiàn)不需要迭代幾次即可滿足精度要求,所以并不會增加多少額外的計算時間。圖2(c)反應(yīng)的是不同步長的偶數(shù)階AC方法計算一個月的積分軌道時右函數(shù)的計算次數(shù)。通過與圖2(a)和圖2(b)對比發(fā)現(xiàn),該方法的右函數(shù)計算次數(shù)明顯少于RK方法和EX方法。

    對于GJ方法,采用同樣的積分策略得到圖2(d)。從圖中可以看出,GJ方法也不需要迭代幾次就可滿足精度要求,且階數(shù)越高計算效率相對越高。對比圖2(c)與圖2(d)可知,AC方法與GJ方法的計算次數(shù)幾乎相等,且都小于RK方法和EX方法。

    通過分析比較不同方法用于GRACE衛(wèi)星軌道積分時右函數(shù)計算次數(shù),可以得出EX方法計算時右函數(shù)計算次數(shù)最多,所以積分效率最低,其次為RK方法,AC方法和GJ方法。

    2.3 狀態(tài)轉(zhuǎn)移矩陣積分誤差對比

    在用動力學(xué)法反演地球重力場模型時需要計算狀態(tài)轉(zhuǎn)移矩陣和參數(shù)敏感矩陣,一般通過對變分方程進(jìn)行數(shù)值積分實現(xiàn)[1,5-6]。積分過程中,首先要積分每個歷元處衛(wèi)星的狀態(tài)向量,再由衛(wèi)星狀態(tài)向量計算狀態(tài)轉(zhuǎn)移矩陣或參數(shù)敏感矩陣的導(dǎo)數(shù)(右函數(shù)),最后由右函數(shù)計算下一歷元的狀態(tài)轉(zhuǎn)移矩陣或參數(shù)敏感矩陣。對于二體問題,可通過解析法得到狀態(tài)轉(zhuǎn)移矩陣,文獻(xiàn)[12]提供了二體問題狀態(tài)轉(zhuǎn)移矩陣解析解算法。本文在二體問題下比較不同方法用于狀態(tài)轉(zhuǎn)移矩陣積分時的積分誤差,仍然采用表1中的GRACE衛(wèi)星軌道參數(shù)作為初始?xì)v元的開普勒軌道根數(shù)。由于在積分狀態(tài)轉(zhuǎn)移矩陣時需要積衛(wèi)星軌道,EX方法積分衛(wèi)星軌道效率低下,所以此處并未對EX方法進(jìn)行比較。狀態(tài)轉(zhuǎn)移矩陣積分效率主要受軌道積分的影響,所以軌道積分效率高的方法狀態(tài)轉(zhuǎn)移矩陣積分同樣效率高。表2為不同方法積分一個月的狀態(tài)轉(zhuǎn)移矩陣誤差(狀態(tài)轉(zhuǎn)移矩陣36個元素中積分誤差最大值),表中符號*表示積分發(fā)散。從表2可以看出,GJ方法積分狀態(tài)轉(zhuǎn)移矩陣誤差比AC方法小一個量級,與RK積分方法的積分誤差相當(dāng)。

    表2 不同參數(shù)下的各種方法狀態(tài)轉(zhuǎn)移矩陣積分誤差最大值

    基于第2.1~2.3節(jié)的分析考慮到GJ算法用于GRACE衛(wèi)星軌道和狀態(tài)轉(zhuǎn)移矩陣積分相對于其它三種積分方法有較大優(yōu)勢,表3給出了GJ算法積分誤差相對較小的階數(shù)以及對應(yīng)的步長,其它三種積分方法不推薦用于長弧段的GRACE衛(wèi)星軌道及狀態(tài)轉(zhuǎn)移矩陣積分。

    表3 GJ算法相對較優(yōu)階數(shù)及對應(yīng)步長

    3 誤差對軌道積分影響分析

    由于軌道積分需要初始時刻的衛(wèi)星狀態(tài)向量,而從GRACELevel1B的GNV1B文件或其它軌道文件中獲得的數(shù)據(jù)含有誤差,該誤差會影響積分軌道的計算。除了初始時刻的衛(wèi)星狀態(tài)含有誤差之外,積分計算時攝動力中含有的誤差也會對積分軌道產(chǎn)生影響。

    3.1 初始狀態(tài)誤差對軌道積分的影響

    為了分析衛(wèi)星初始時刻狀態(tài)向量誤差對軌道積分的影響,在軌道積分時分別僅對初始狀態(tài)向量中的位置和速度向量加入不同量級的誤差,分析該誤差對軌道積分的影響。通過采用不同的積分方法積分后發(fā)現(xiàn)初始狀態(tài)誤差對不同積分方法的影響幾乎相同,因此本文以GJ方法為例分析初始狀態(tài)誤差對軌道積分的影響。圖3(a)為僅對初始時刻衛(wèi)星位置向量加入不同量級誤差(100mm級誤差表示X、Y、Z方向各加100mm誤差,其它情形類似)后積分軌道誤差圖。圖3(b)為僅對衛(wèi)星初始時刻速度向量加入不同量級誤差(mm/s級誤差表示X、Y、Z方向各加1mm/s誤差,其它情形類似)后積分軌道誤差圖。為了分析初始狀態(tài)中位置和速度向量誤差對積分軌道的綜合影響,在初始狀態(tài)中共同加入不同量級的位置和速度誤差后得到積分軌道誤差圖3(c)。

    通過比較分析圖3的計算結(jié)果可以得出以下結(jié)論:初始狀態(tài)誤差對積分軌道的影響為隨著積分時間增加積分軌道誤差呈現(xiàn)線性累加特性。衛(wèi)星初始狀態(tài)誤差對軌道積分的影響遠(yuǎn)遠(yuǎn)大于積分算法本身的誤差。對比圖3(a)、3(b)和3(c)發(fā)現(xiàn),衛(wèi)星初始位置和速度向量中速度向量誤差(0.1mm/s)對軌道積分的影響大于位置向量誤差(10mm)對軌道積分的影響,若位置誤差為10mm量級,則當(dāng)速度誤差為0.01mm/s量級時,兩者對軌道積分的影響才能保持一致,而在當(dāng)前衛(wèi)星定軌精度下是很難達(dá)到的。其可能的原因為,軌道積分時速度為位置的右函數(shù),而速度的右函數(shù)為衛(wèi)星所受的力。位置誤差所引起的攝動力變化相對較小,從而引起的速度變化也相對較小,所以在隨后的積分中位置誤差相對較?。欢绻跏紶顟B(tài)中速度存在誤差,由于速度為位置的右函數(shù)則在對位置進(jìn)行數(shù)值積分時受該速度誤差影響較大,會使位置產(chǎn)生相對較大的偏差,具有較大偏差的位置對應(yīng)的攝動力偏差也相對較大,則在隨后的積分中位置誤差較大。

    3.2 攝動力中的隨機(jī)誤差對軌道積分的影響

    在對衛(wèi)星進(jìn)行軌道積分時,由于先驗力模型及非保守力測量數(shù)據(jù)的不準(zhǔn)確會對軌道積分產(chǎn)生較大的影響。非保守力數(shù)據(jù)的誤差大致為10-10ms-2量級,保守力中主要是先驗重力場模型的誤差較大,雖然在重力反演中會對先驗重力場模型和非保守力數(shù)據(jù)進(jìn)行校正,但如果誤差較大會影響參數(shù)的估計,本文取兩個不同的先驗重力場模型并取至不同的階次得到先驗重力場模型誤差引起的衛(wèi)星引力誤差量級。通過分別計算GGM05C模型(取至140階次)和EGM2008模型(取至70階次)對衛(wèi)星的引力并對引力的差做統(tǒng)計(時間為1d,間隔為5s)得到圖4(a)。根據(jù)圖4(a)的結(jié)果在二體問題下,對中心引力加上一個服從正態(tài)分布N(0, (5×10-7ms-2)2)的隨機(jī)誤差,圖4(b)為以GJ方法計算100次攝動力誤差對軌道積分結(jié)果影響圖。從圖4(b)可以看出,攝動力誤差對軌道積分的影響可超過20000m,大于初始狀態(tài)誤差對軌道積分的影響且無規(guī)律可循,影響也是隨機(jī)的?;诒疚牡姆治龅贸鋈糁苯硬捎?個月作為積分弧段進(jìn)行重力反演,則必須考慮如何減小該誤差的影響。

    4 改進(jìn)Gauss-Jackson方法

    雖然GJ算法具有較高的計算精度且有多種階數(shù)和步長可以選擇,但是該算法只能計算出歷元為步長整數(shù)倍處的積分軌道。由于GRACE LEVEL1B數(shù)據(jù)中給出的是5s間隔的標(biāo)準(zhǔn)軌道和星間距觀測數(shù)據(jù),所以一般采用5s步長進(jìn)行軌道積分,對于長弧段軌道積分積分效率較低?;谇拔牡姆治龅贸霾捎幂^大的步長進(jìn)行積分也可以達(dá)到較高的積分精度,而且提高了積分效率。為了得到步長內(nèi)任意點的積分軌道可采用適當(dāng)?shù)膬?nèi)插算法[13],只要步長的大小控制在一定的范圍內(nèi),并選擇合適的內(nèi)插算法,內(nèi)插的精度可以很高。內(nèi)插計算不僅可以得到5s步長整數(shù)倍處的軌道值,而且可以得到任意歷元處的軌道值,當(dāng)以GPS解算的幾何軌道作為標(biāo)準(zhǔn)軌道時可簡化幾何軌道的預(yù)處理。

    本文選擇8階30s作為積分階數(shù)和步長,在內(nèi)插算法上,采用10點9階基于移動窗口的多項式內(nèi)插算法[6,14],在起始段和末尾段采用同階拉格朗日算法[11],除了弧段兩端外,其它區(qū)域能保證為內(nèi)插窗口的中心區(qū)域。通過模擬數(shù)據(jù)計算得到該種內(nèi)插算法內(nèi)插出的軌道誤差在10-8~10-9m量級,起始段和末尾段內(nèi)插出的軌道誤差為10-6~10-7m量級,狀態(tài)轉(zhuǎn)移矩陣誤差在10-11量級,完全能滿足軌道積分和狀態(tài)轉(zhuǎn)移矩陣積分要求。內(nèi)插算法的計算時間相對于右函數(shù)的計算時間幾乎可以忽略不計,所以改進(jìn)GJ算法不僅保證了軌道積分精度,而且提高了軌道積分效率。圖5為二體問題下采用改進(jìn)GJ方法計算的積分軌道誤差圖,完全能滿足GRACE衛(wèi)星積分精度要求。

    為了分析改進(jìn)GJ方法在實際衛(wèi)星軌道積分中的性能,本文以GRACE Level1B數(shù)據(jù)以及相應(yīng)的力模型(見表4)進(jìn)行軌道積分試驗。計算時分別采用8階5s步長的GJ方法和改進(jìn)GJ方法進(jìn)行積分,并從積分精度和積分效率方面進(jìn)行比較。由于真實數(shù)據(jù)中存在數(shù)據(jù)間斷和粗差,需要對間斷數(shù)據(jù)進(jìn)行填充并剔除粗差。為了計算方便,此處利用一天的實際數(shù)據(jù)進(jìn)行軌道積分(未對初始狀態(tài)向量和力模型參數(shù)校正)并與JPL提供的軌道數(shù)據(jù)(GRACE GNV1B)對比。

    對比發(fā)現(xiàn),改進(jìn)GJ方法的最大積分誤差比8階5s步長GJ方法略小(200mm左右),并在計算效率上有了較大提高。采用8階5s步長GJ方法用時約57min,采用改進(jìn)GJ方法用時約18min,積分效率提高了將近3倍。

    表4 軌道積分計算所用的力模型

    5 結(jié) 論

    本文對比了多種軌道積分方法并分析了積分起始數(shù)據(jù)誤差對軌道積分結(jié)果的影響,提出了改進(jìn)GJ算法?;趯Ρ确治龅贸隽艘韵掠幸娴慕Y(jié)論:

    1)對于GRACE衛(wèi)星長弧段軌道積分,GJ算法由于在積分的過程中沒有直接采用前面歷元的積分軌道值使得其積分誤差比RK、EX和AC算法小(可達(dá)2個量級)。

    2)EX算法在積分階數(shù)很高的情況下仍可以保持較高的積分精度,但是對于低軌衛(wèi)星計算效率低下;AC算法和GJ算法同為多步法,但對于GRACE衛(wèi)星軌道積分AC算法不及GJ算法,不推薦在長弧段軌道積分中使用EX算法和AC算法,RK算法只建議用于多步法的起步值計算。

    3)對于GRACE衛(wèi)星軌道和狀態(tài)轉(zhuǎn)移矩陣積分,GJ算法較適宜的積分階數(shù)為8~10階,對應(yīng)的步長為5~40s和5~50s。

    4)衛(wèi)星初始狀態(tài)誤差和力模型誤差對軌道積分的影響遠(yuǎn)大于算法本身的誤差;初始狀態(tài)中的速度誤差(0.1mm/s)對軌道積分的影響比位置誤差(10mm)大將近1個量級,且該類誤差對軌道積分的影響呈現(xiàn)線性累加特性;力模型誤差對軌道積分的影響較大且無規(guī)律可循。

    5)GJ算法中取合適的較大積分步長(30s)并結(jié)合適宜階數(shù)(9階)的基于移動窗口的多項式內(nèi)插算法不僅可以大大提高積分效率而且可以得到任意時刻衛(wèi)星的積分軌道和狀態(tài)轉(zhuǎn)移矩陣,有利于GRACE衛(wèi)星不規(guī)則歷元的幾何軌道和星間觀測值聯(lián)合進(jìn)行衛(wèi)星重力反演。

    [1] Montenbruck O, Gill E. Satellite orbits: models, methods, and applications[M]. Berlin: Springer-Verlag, 2005.

    [2] Beutler G. Methods of celestial mechanics[M]. Berlin: Springer-Verlag, 2005.

    [3] Berry M M, Healy L M. Implementation of Gauss-Jackson integration for orbit propagation [J]. The Journal of the Astronautical Sciences, 2004, 52(3): 331-357.

    [4] Bulirsch R, Stoer J. Numerical treatment of ordinary differential equations by extrapolation methods [J].Numerische Mathematik, 1966, 8(1): 1-13.

    [5] 王正濤. 衛(wèi)星跟蹤衛(wèi)星測量確定地球重力場的理論與方法[D]. 武漢: 武漢大學(xué), 2005.[Wang Zheng-tao. Theory and methodology of earth gravity field recovery by satellite-to-satellite tracking data[D]. Wuhan: Wuhan University, 2005.]

    [6] 游為. 應(yīng)用低軌衛(wèi)星數(shù)據(jù)反演地球重力場模型的理論和方法[D]. 成都: 西南交通大學(xué), 2011. [You Wei. Theory and methodology of earth′s gravitational field model recovery by LEO data[D]. Chengdu: Southwest Jiaotong University, 2011.]

    [7] 龔曉穎, 蔡詩響, 李征航. 小步長低階Collocation方法軌道積分性能研究[J]. 宇航學(xué)報,2013, 34(5): 642-648. [Gong Xiao-ying, Cai Shi-xiang, Li Zheng-hang. Performance research of small-step length and low-order collocation method in orbit determination[J]. Journal of Astronautics, 2013, 34(5): 642-648.]

    [8] Xu G C. GPS theory, algorithms and applications[M]. Berlin: Springer-Verlag, 2007.

    [9] 劉大杰, 施一民, 過靜珺.全球定位系統(tǒng)(GPS)的理論與數(shù)據(jù)處理[M]. 上海: 同濟(jì)大學(xué)出版社, 1997: 257-263.

    [10] 羅志才,周浩, 鐘波, 等. Gauss-Jackson積分算法分析與驗證[J]. 武漢大學(xué)學(xué)報·信息科學(xué)版, 2013, 38(11):1364-1368.[Luo Zhi-cai, Zhou Hao, Zhong Bo, et al. Analysis and validation of Gauss-Jackson integral algorithm[J]. Geomatics and Information Science of Wuhan University, 2013, 38(11): 1364-1368.]

    [11] 李大美, 李素貞, 朱方生. 計算方法[M]. 武漢: 武漢大學(xué)出版社, 2012: 94-130.

    [12] 蔣虎, 黃珹, 胡小工. 二體問題狀態(tài)轉(zhuǎn)移矩陣的分析解和數(shù)值解比較[J]. 云南天文臺臺刊,1999, 4: 19-29. [Jiang Hu, Huang Cheng, Hu Xiao-gong. Comparison between analytic and numerical solutions of state transition matrices of two body problem[J]. Publications of the Yunnan Astronomical Observatory, 1999, 4: 19-29 .]

    [13] 鄭偉, 許厚澤, 鐘敏, 等. 不同插值法對下一代衛(wèi)星重力反演精度的影響[J]. 宇航學(xué)報, 2014, 35(3): 269-276. [Zheng Wei, Xu Hou-ze, Zhong Min, et al. Influences of different interpolation methods on accuracy of next-generation satellite gravity recovery[J]. Journal of Astronautics, 2014, 35(3): 269-276.]

    [14] 游為, 沈云中, 范東明, 等. 基于衛(wèi)星軌道擾動理論的重力反演方法[J]. 地球物理學(xué)報, 2010, 53(11): 2574-2581. [You Wei, Shen Yun-zhong, Fan Dong-ming, et al. The algorithm of earth′s gravitational field recovery based on satellite′s orbital perturbation theory[J]. Chinese Journal of Geophysics, 2010, 53(11): 2574-2581.]

    通信地址:四川省成都市郫都區(qū)西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院(611756)

    電話:18382157390

    E-mail:yubiao0221@qq.com

    游 為(1985-),男,講師,博士,主要從事衛(wèi)星重力測量及測量數(shù)據(jù)處理研究。本文通信作者。

    通信地址:四川省成都市郫都區(qū)西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院(611756)

    電話:18123216821

    E-mail: youwei1985@foxmail.com

    (編輯:牛苗苗)

    Analysis of Orbit Numerical Integration Methods in Earth′s Gravitational Field Recovery by GRACE

    YU Biao1, FAN Dong-ming1, YOU Wei1, GU Yan-chao1, SU Yong2

    (1.Faculty of Geosciences and Environmental Engineering, Southwest Jiao-tong University, Chengdu 611756, China; 2.School of Civil Engineering and Architecture, Southwest Petroleum University, Chengdu 610500, China)

    This paper assesses the existing numerical integration methods including the Runge-Kutta methods, the Adams-Cowell methods, the Gauss-Jackson methods and the extrapolation methods in the computation of the GRACE satellite orbits and state transition matrix. According to the results, we recommend the Gauss-Jackson methods with optimal parameters to integral orbit and also analyze the advantages of the Gauss-Jackson methods relative to the other methods. Then the resistance capacity of the random errors which are included in the satellite initial state for those methods is analyzed, and the numerical results indicate that all methods show the similar little resistance capacity. The integral orbits are more sensitive to the satellite initial state velocity errors (0.1mm/s) than the initial state position errors(10mm). Compared to the errors in the satellite initial state, the perturbation force model errors have a significant impact on the orbit integral precision. Finally, we proposes a modified method which combines the Gauss-Jackson algorithm and the moving-window polynomial interpolation algorithm to overcome the shortage of the large step size of the Gauss-Jackson method. The simulated and actual results show that the modified method can give the satellite position and state transition matrix at any time and has a significant improvement in computing efficiency while possessing high accuracy.

    Orbit integration; Integral precision; Integral efficiency; Error influence

    2016-10-24;

    2017-01-03

    國家自然科學(xué)基金(41574018,41404018)

    P228.1

    A

    1000-1328(2017)03-0253-09

    10.3873/j.issn.1000-1328.2017.03.005

    余 彪(1990-),男,碩士生,主要從事衛(wèi)星重力測量研究。

    猜你喜歡
    歷元階數(shù)步長
    基于Armijo搜索步長的BFGS與DFP擬牛頓法的比較研究
    關(guān)于無窮小階數(shù)的幾點注記
    確定有限級數(shù)解的階數(shù)上界的一種n階展開方法
    歷元間載波相位差分的GPS/BDS精密單點測速算法
    Recent advances of TCM treatment of childhood atopic dermatitis
    Clinical observation of Huatan Huoxue Formula in treating coronary heart disease with hyperlipidemia
    Mechanism of sex hormone level in biological clock disorder induced acne and analysis of TCM Pathogenesis
    基于逐維改進(jìn)的自適應(yīng)步長布谷鳥搜索算法
    一種新的多址信道有效階數(shù)估計算法*
    關(guān)于動態(tài)電路階數(shù)的討論
    波多野结衣高清作品| 国产av麻豆久久久久久久| 深爱激情五月婷婷| 宅男免费午夜| 在线天堂最新版资源| 天堂影院成人在线观看| 又粗又爽又猛毛片免费看| 日韩精品中文字幕看吧| 精品久久久久久久人妻蜜臀av| 色视频www国产| 国产在视频线在精品| 天堂网av新在线| 搞女人的毛片| 精品一区二区三区人妻视频| 亚洲精品成人久久久久久| 丰满人妻一区二区三区视频av | 国产亚洲精品久久久久久毛片| 精品电影一区二区在线| 国产精品久久久人人做人人爽| 长腿黑丝高跟| 色精品久久人妻99蜜桃| 日本撒尿小便嘘嘘汇集6| 亚洲久久久久久中文字幕| 窝窝影院91人妻| 婷婷亚洲欧美| 国产精品影院久久| 日韩欧美三级三区| 成人午夜高清在线视频| 无人区码免费观看不卡| 观看免费一级毛片| 色老头精品视频在线观看| 成人av一区二区三区在线看| 丰满的人妻完整版| 叶爱在线成人免费视频播放| 99热6这里只有精品| 免费看日本二区| 国产熟女xx| 国产精品久久久久久久电影 | 性色avwww在线观看| 欧美日韩亚洲国产一区二区在线观看| 成人精品一区二区免费| 亚洲成人免费电影在线观看| 搡老岳熟女国产| 亚洲18禁久久av| 久久久久久大精品| 在线观看午夜福利视频| 特大巨黑吊av在线直播| 亚洲国产欧美网| 亚洲人成网站在线播放欧美日韩| 18禁黄网站禁片免费观看直播| 五月玫瑰六月丁香| 中亚洲国语对白在线视频| 国产视频内射| 看黄色毛片网站| av福利片在线观看| 日韩欧美一区二区三区在线观看| 性色av乱码一区二区三区2| 夜夜看夜夜爽夜夜摸| av中文乱码字幕在线| 男女床上黄色一级片免费看| 亚洲精品影视一区二区三区av| 午夜激情欧美在线| 欧美另类亚洲清纯唯美| 特级一级黄色大片| 亚洲精品在线观看二区| 亚洲熟妇熟女久久| 成人无遮挡网站| 首页视频小说图片口味搜索| 色综合婷婷激情| 全区人妻精品视频| 久久精品国产亚洲av香蕉五月| 18禁在线播放成人免费| 国产成人aa在线观看| 一本精品99久久精品77| 国产乱人伦免费视频| 国产成人系列免费观看| 免费看美女性在线毛片视频| x7x7x7水蜜桃| 国产av麻豆久久久久久久| www日本在线高清视频| 亚洲精品日韩av片在线观看 | 亚洲成人久久爱视频| 免费在线观看日本一区| 国产精品久久久久久久电影 | 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 欧美另类亚洲清纯唯美| 天天一区二区日本电影三级| 国产男靠女视频免费网站| www.999成人在线观看| 男女视频在线观看网站免费| 精品午夜福利视频在线观看一区| 香蕉久久夜色| 久久久久亚洲av毛片大全| 国产高清视频在线观看网站| 午夜免费激情av| 欧美极品一区二区三区四区| 久久精品国产亚洲av涩爱 | 亚洲中文字幕一区二区三区有码在线看| 午夜福利欧美成人| 国产黄色小视频在线观看| 亚洲人成电影免费在线| 内地一区二区视频在线| 香蕉久久夜色| 亚洲真实伦在线观看| 亚洲国产欧洲综合997久久,| 成人午夜高清在线视频| 首页视频小说图片口味搜索| 免费观看精品视频网站| 亚洲精品在线美女| 久久久久久久久大av| 69av精品久久久久久| 国产97色在线日韩免费| 国产午夜精品论理片| 亚洲精品日韩av片在线观看 | 国产精品永久免费网站| 全区人妻精品视频| 久久精品国产清高在天天线| 色精品久久人妻99蜜桃| 欧美国产日韩亚洲一区| 在线a可以看的网站| 在线视频色国产色| 舔av片在线| 国产成人av教育| 精品一区二区三区av网在线观看| 精品不卡国产一区二区三区| 亚洲激情在线av| 99久久无色码亚洲精品果冻| 亚洲第一欧美日韩一区二区三区| 亚洲无线观看免费| 神马国产精品三级电影在线观看| 国产一区二区亚洲精品在线观看| 精品欧美国产一区二区三| 日本三级黄在线观看| 欧美激情久久久久久爽电影| 搡女人真爽免费视频火全软件 | 中文字幕精品亚洲无线码一区| xxx96com| 欧美日韩黄片免| 丰满人妻一区二区三区视频av | 看黄色毛片网站| 色哟哟哟哟哟哟| 免费电影在线观看免费观看| 欧美激情在线99| 在线观看日韩欧美| 亚洲欧美日韩高清在线视频| 1000部很黄的大片| 亚洲国产欧洲综合997久久,| 日韩高清综合在线| 成熟少妇高潮喷水视频| 免费av不卡在线播放| avwww免费| 欧美一区二区国产精品久久精品| 亚洲人成电影免费在线| 国产激情偷乱视频一区二区| 久久精品人妻少妇| 有码 亚洲区| 欧美色视频一区免费| 九色国产91popny在线| ponron亚洲| 国产真实伦视频高清在线观看 | 露出奶头的视频| 国产精品精品国产色婷婷| 51午夜福利影视在线观看| 啪啪无遮挡十八禁网站| 高清毛片免费观看视频网站| 日韩高清综合在线| 国产美女午夜福利| 女人被狂操c到高潮| 无人区码免费观看不卡| 又粗又爽又猛毛片免费看| 无遮挡黄片免费观看| 美女高潮喷水抽搐中文字幕| АⅤ资源中文在线天堂| 久久久久久国产a免费观看| 听说在线观看完整版免费高清| 在线观看美女被高潮喷水网站 | 亚洲熟妇熟女久久| 日韩高清综合在线| 午夜久久久久精精品| 国产乱人伦免费视频| 色噜噜av男人的天堂激情| 日韩av在线大香蕉| 国产在视频线在精品| 亚洲国产高清在线一区二区三| 搡老岳熟女国产| 18禁美女被吸乳视频| 久久精品国产99精品国产亚洲性色| 美女被艹到高潮喷水动态| 色哟哟哟哟哟哟| 亚洲国产精品sss在线观看| 午夜亚洲福利在线播放| 亚洲国产中文字幕在线视频| 亚洲成人免费电影在线观看| 国产欧美日韩一区二区三| 欧美黑人巨大hd| 亚洲欧美精品综合久久99| 18+在线观看网站| 午夜两性在线视频| 日韩欧美国产一区二区入口| 在线观看一区二区三区| 国产精品一及| 久久久精品欧美日韩精品| 日本 欧美在线| 91久久精品国产一区二区成人 | 久久午夜亚洲精品久久| 美女高潮的动态| 色尼玛亚洲综合影院| 日韩大尺度精品在线看网址| 国产免费男女视频| 夜夜看夜夜爽夜夜摸| 久久精品人妻少妇| 一a级毛片在线观看| 亚洲av电影在线进入| 在线观看免费午夜福利视频| 色综合婷婷激情| 人人妻人人澡欧美一区二区| 免费大片18禁| 午夜福利成人在线免费观看| 精品不卡国产一区二区三区| 床上黄色一级片| 国产97色在线日韩免费| 色老头精品视频在线观看| 国产精品香港三级国产av潘金莲| 日韩欧美精品免费久久 | 久久伊人香网站| 免费av不卡在线播放| 日韩欧美三级三区| 亚洲人与动物交配视频| 国产成+人综合+亚洲专区| 精华霜和精华液先用哪个| 欧美激情在线99| 精品久久久久久久久久久久久| 中文亚洲av片在线观看爽| 国产精品一区二区三区四区久久| 亚洲人与动物交配视频| 亚洲欧美日韩东京热| 午夜a级毛片| 最后的刺客免费高清国语| 高潮久久久久久久久久久不卡| 身体一侧抽搐| 欧美一区二区亚洲| 亚洲国产中文字幕在线视频| 精品一区二区三区视频在线观看免费| 久久久久久久午夜电影| 欧美一级毛片孕妇| 动漫黄色视频在线观看| 国内毛片毛片毛片毛片毛片| 国产探花在线观看一区二区| 久久国产乱子伦精品免费另类| 亚洲欧美日韩高清专用| 亚洲一区二区三区不卡视频| 久久中文看片网| www国产在线视频色| 欧美日韩综合久久久久久 | 成年版毛片免费区| 99热只有精品国产| 国产精华一区二区三区| 国产av麻豆久久久久久久| 亚洲 国产 在线| 19禁男女啪啪无遮挡网站| 成人av在线播放网站| 天堂av国产一区二区熟女人妻| 好男人电影高清在线观看| 久久精品综合一区二区三区| 国产麻豆成人av免费视频| 免费av观看视频| 亚洲av成人不卡在线观看播放网| 乱人视频在线观看| 日韩欧美免费精品| 麻豆国产av国片精品| 高清在线国产一区| 成人国产综合亚洲| 少妇人妻一区二区三区视频| 国产精品爽爽va在线观看网站| 成人av在线播放网站| 国产 一区 欧美 日韩| 欧美乱码精品一区二区三区| 亚洲精品在线观看二区| 国产一区二区亚洲精品在线观看| 女警被强在线播放| 一区二区三区国产精品乱码| 又粗又爽又猛毛片免费看| 午夜免费观看网址| 欧美在线一区亚洲| 成年免费大片在线观看| 国产欧美日韩精品一区二区| 久久久久久人人人人人| 日本精品一区二区三区蜜桃| 99riav亚洲国产免费| 亚洲人成网站在线播| 757午夜福利合集在线观看| 91久久精品电影网| 国产精品免费一区二区三区在线| 亚洲黑人精品在线| bbb黄色大片| 亚洲国产日韩欧美精品在线观看 | 欧美成狂野欧美在线观看| 欧美丝袜亚洲另类 | 国产伦精品一区二区三区四那| 亚洲一区二区三区色噜噜| 成年女人永久免费观看视频| 亚洲精品影视一区二区三区av| 制服人妻中文乱码| 亚洲成人久久性| 长腿黑丝高跟| 一个人观看的视频www高清免费观看| 好男人在线观看高清免费视频| 在线播放国产精品三级| 国产精品99久久久久久久久| 欧美又色又爽又黄视频| 哪里可以看免费的av片| 亚洲成av人片免费观看| 亚洲在线自拍视频| 一区二区三区激情视频| 久久精品国产亚洲av香蕉五月| 69av精品久久久久久| 成人亚洲精品av一区二区| 国产伦在线观看视频一区| 亚洲精品在线观看二区| 久久久久久国产a免费观看| 1024手机看黄色片| 欧美一级a爱片免费观看看| 精品一区二区三区视频在线观看免费| 91av网一区二区| 搡老妇女老女人老熟妇| 精品午夜福利视频在线观看一区| 男人舔奶头视频| 亚洲五月天丁香| 脱女人内裤的视频| 日本 av在线| 欧美日韩一级在线毛片| 国语自产精品视频在线第100页| 亚洲精品粉嫩美女一区| 欧洲精品卡2卡3卡4卡5卡区| 日韩免费av在线播放| av视频在线观看入口| 欧美一级a爱片免费观看看| 亚洲专区中文字幕在线| 免费看美女性在线毛片视频| 色综合婷婷激情| 免费电影在线观看免费观看| 夜夜夜夜夜久久久久| 久久香蕉精品热| 天天添夜夜摸| 啦啦啦韩国在线观看视频| 草草在线视频免费看| e午夜精品久久久久久久| 精品一区二区三区人妻视频| av视频在线观看入口| 亚洲va日本ⅴa欧美va伊人久久| 19禁男女啪啪无遮挡网站| 啦啦啦观看免费观看视频高清| 一a级毛片在线观看| 欧美日韩福利视频一区二区| 日韩免费av在线播放| 久久久久国内视频| 日本熟妇午夜| 久久久久久久精品吃奶| 亚洲国产日韩欧美精品在线观看 | 亚洲精品456在线播放app | 色综合婷婷激情| 国产不卡一卡二| 成年人黄色毛片网站| av专区在线播放| 男女午夜视频在线观看| 国产亚洲av嫩草精品影院| 老司机在亚洲福利影院| 男女视频在线观看网站免费| 欧美绝顶高潮抽搐喷水| 国产精品爽爽va在线观看网站| 久久久久久久久中文| 黄片小视频在线播放| 色综合欧美亚洲国产小说| 欧美av亚洲av综合av国产av| xxxwww97欧美| 一区二区三区高清视频在线| 观看美女的网站| 精品久久久久久久人妻蜜臀av| 99精品欧美一区二区三区四区| 91麻豆av在线| 国内揄拍国产精品人妻在线| 特大巨黑吊av在线直播| 亚洲欧美日韩无卡精品| 国产精品 国内视频| 香蕉久久夜色| 十八禁网站免费在线| 日本一本二区三区精品| 国产精品,欧美在线| 亚洲欧美精品综合久久99| 久久欧美精品欧美久久欧美| 99精品久久久久人妻精品| 精品不卡国产一区二区三区| 国产精品 国内视频| 99国产精品一区二区三区| 久久久久久大精品| 中文字幕久久专区| 在线免费观看的www视频| 国产精品久久电影中文字幕| 欧美av亚洲av综合av国产av| e午夜精品久久久久久久| 中文字幕av成人在线电影| 久久久色成人| 亚洲精品在线美女| 桃红色精品国产亚洲av| 亚洲专区国产一区二区| 嫩草影院入口| 欧美精品啪啪一区二区三区| av在线天堂中文字幕| 国产av麻豆久久久久久久| www.色视频.com| 两人在一起打扑克的视频| 三级国产精品欧美在线观看| 美女cb高潮喷水在线观看| 一个人看的www免费观看视频| 国产高清有码在线观看视频| 熟妇人妻久久中文字幕3abv| 99国产综合亚洲精品| 国产欧美日韩精品一区二区| 国产成年人精品一区二区| 亚洲成a人片在线一区二区| 久久久久久九九精品二区国产| 禁无遮挡网站| 亚洲国产高清在线一区二区三| 中国美女看黄片| 手机成人av网站| 51午夜福利影视在线观看| 韩国av一区二区三区四区| 琪琪午夜伦伦电影理论片6080| 精品国内亚洲2022精品成人| 波多野结衣高清作品| 九色国产91popny在线| 久久这里只有精品中国| 成熟少妇高潮喷水视频| 搡老岳熟女国产| 中文字幕精品亚洲无线码一区| 两个人的视频大全免费| 欧美成人性av电影在线观看| 两性午夜刺激爽爽歪歪视频在线观看| 噜噜噜噜噜久久久久久91| 国产精品久久电影中文字幕| 日韩欧美 国产精品| 性色avwww在线观看| av福利片在线观看| 99久久无色码亚洲精品果冻| 国产黄片美女视频| 国产精品 国内视频| 老司机在亚洲福利影院| 亚洲精华国产精华精| 黄色日韩在线| 欧美成人a在线观看| 婷婷精品国产亚洲av| 亚洲av第一区精品v没综合| 欧美日韩国产亚洲二区| eeuss影院久久| 国产精品久久视频播放| 18+在线观看网站| 欧美bdsm另类| 久久久色成人| 又黄又爽又免费观看的视频| 综合色av麻豆| 精品一区二区三区av网在线观看| 国产成人系列免费观看| 午夜激情欧美在线| 亚洲 欧美 日韩 在线 免费| 国产精品久久久久久人妻精品电影| 好看av亚洲va欧美ⅴa在| 久久久久久久午夜电影| 日韩大尺度精品在线看网址| 他把我摸到了高潮在线观看| 色av中文字幕| 美女大奶头视频| 亚洲欧美激情综合另类| 国产爱豆传媒在线观看| 男人和女人高潮做爰伦理| 国产视频一区二区在线看| 国产精品久久久久久人妻精品电影| 国产成+人综合+亚洲专区| 中文字幕熟女人妻在线| a级一级毛片免费在线观看| 国产欧美日韩精品亚洲av| 国产亚洲精品综合一区在线观看| 啪啪无遮挡十八禁网站| 免费在线观看影片大全网站| 毛片女人毛片| 动漫黄色视频在线观看| 日韩高清综合在线| 国语自产精品视频在线第100页| 欧美激情在线99| 亚洲国产精品合色在线| 69人妻影院| 欧美绝顶高潮抽搐喷水| 亚洲欧美一区二区三区黑人| 国产精品国产高清国产av| 免费人成视频x8x8入口观看| 欧美不卡视频在线免费观看| 老鸭窝网址在线观看| 免费无遮挡裸体视频| 亚洲av熟女| 欧美日韩瑟瑟在线播放| 一夜夜www| 特级一级黄色大片| 亚洲性夜色夜夜综合| 国产熟女xx| 亚洲性夜色夜夜综合| 国产亚洲精品综合一区在线观看| 亚洲国产欧美人成| 免费一级毛片在线播放高清视频| 美女黄网站色视频| 午夜精品一区二区三区免费看| 亚洲,欧美精品.| 一级毛片女人18水好多| 亚洲精品粉嫩美女一区| 亚洲av中文字字幕乱码综合| 嫩草影院精品99| 欧美乱码精品一区二区三区| 免费看美女性在线毛片视频| 日韩国内少妇激情av| 91在线精品国自产拍蜜月 | 成年女人永久免费观看视频| 日日摸夜夜添夜夜添小说| 久久久久久久午夜电影| 日韩欧美在线乱码| 极品教师在线免费播放| 免费看光身美女| 色播亚洲综合网| 搡老熟女国产l中国老女人| 在线观看美女被高潮喷水网站 | 日本三级黄在线观看| 国内揄拍国产精品人妻在线| 日韩欧美国产在线观看| 欧美黄色片欧美黄色片| 高潮久久久久久久久久久不卡| 无遮挡黄片免费观看| 精品午夜福利视频在线观看一区| 精品国产三级普通话版| 久久国产精品人妻蜜桃| 国产成人福利小说| 99久久精品一区二区三区| 亚洲狠狠婷婷综合久久图片| 波多野结衣巨乳人妻| 黄色视频,在线免费观看| 亚洲国产日韩欧美精品在线观看 | 欧美3d第一页| 亚洲美女黄片视频| 99精品久久久久人妻精品| 波多野结衣高清无吗| 最新在线观看一区二区三区| 美女免费视频网站| 欧美极品一区二区三区四区| 亚洲美女视频黄频| 国产色婷婷99| 国产一区二区三区在线臀色熟女| 一本综合久久免费| 一夜夜www| 悠悠久久av| 午夜激情欧美在线| 2021天堂中文幕一二区在线观| 午夜福利成人在线免费观看| 欧美午夜高清在线| 啦啦啦观看免费观看视频高清| 丰满乱子伦码专区| 国产伦在线观看视频一区| 国产高清有码在线观看视频| 伊人久久精品亚洲午夜| 欧美xxxx黑人xx丫x性爽| 欧美不卡视频在线免费观看| 给我免费播放毛片高清在线观看| 国产高清视频在线播放一区| 欧美丝袜亚洲另类 | 欧美色欧美亚洲另类二区| 亚洲国产精品999在线| 亚洲第一欧美日韩一区二区三区| 伊人久久大香线蕉亚洲五| 国产亚洲av嫩草精品影院| 99热精品在线国产| 夜夜躁狠狠躁天天躁| 欧美日本亚洲视频在线播放| 国产一区二区三区在线臀色熟女| 国产aⅴ精品一区二区三区波| 午夜免费观看网址| av黄色大香蕉| 久久99热这里只有精品18| 亚洲一区高清亚洲精品| 搡老岳熟女国产| 国产v大片淫在线免费观看| 欧美+亚洲+日韩+国产| 俺也久久电影网| 中文字幕高清在线视频| 国产一区二区在线av高清观看| 中文字幕人妻熟人妻熟丝袜美 | 叶爱在线成人免费视频播放| 人妻丰满熟妇av一区二区三区| 久久久久久久久大av| av黄色大香蕉| 亚洲,欧美精品.| 又爽又黄无遮挡网站| 亚洲真实伦在线观看| 国产中年淑女户外野战色| 欧美区成人在线视频| 欧美在线一区亚洲| 国产午夜精品论理片| 久久精品影院6| 男女那种视频在线观看| 久久精品国产亚洲av涩爱 | 特大巨黑吊av在线直播| 成人国产一区最新在线观看| 一本久久中文字幕| 国产精品99久久99久久久不卡| 国产成人福利小说| 国产精品精品国产色婷婷| 日韩大尺度精品在线看网址| 国产真人三级小视频在线观看| 美女 人体艺术 gogo| 每晚都被弄得嗷嗷叫到高潮| 男人和女人高潮做爰伦理|