2012年全國(guó)數(shù)學(xué)建模競(jìng)賽優(yōu)秀選靳高曹

1、 填寫）參賽（由第九屆“杯”數(shù)學(xué)建模競(jìng)賽學(xué)校中國(guó)海洋大學(xué)參賽隊(duì)號(hào)104230061.隊(duì)員2.曹安州3.參賽（由 填寫）第九屆“杯”數(shù)學(xué)建模競(jìng)賽題 目 基于伴隨同化方法的空間飛行器軌道估計(jì)與誤差分析摘要：利用觀測(cè)獲得的空間飛行器的軌道參數(shù)，結(jié)合動(dòng)力學(xué)方程對(duì)空間飛行器的軌道進(jìn)行估計(jì)，能夠?yàn)轱w行器類別、飛行意圖的判斷提供信息基礎(chǔ)。問題一對(duì)觀測(cè)的三維位置估計(jì)是基于簡(jiǎn)化運(yùn)動(dòng)方程，利用觀測(cè)零時(shí)刻的位置和速度作為初始條件對(duì)運(yùn)動(dòng)方程數(shù)值求解，得到觀測(cè)的運(yùn)動(dòng)軌跡，并得出觀測(cè)運(yùn)行周期。通過與中低軌近圓軌道運(yùn)行軌跡和周期的對(duì)比，驗(yàn)證了計(jì)算結(jié)果的準(zhǔn)確性。問題二中 06 號(hào)和 09 號(hào)觀測(cè)對(duì) 0 號(hào)空間飛行器的觀測(cè)數(shù)據(jù)

2、時(shí)間不同步而且數(shù)據(jù)中包含了隨機(jī)誤差（白噪聲），故首先利用低通濾波器對(duì)觀測(cè)數(shù)據(jù)進(jìn)行濾波處理，再將觀測(cè)數(shù)據(jù)插值為時(shí)間的序1列。然后將觀測(cè)坐標(biāo)系下的觀測(cè)數(shù)據(jù)轉(zhuǎn)化到基礎(chǔ)坐標(biāo)系中，根據(jù)逐點(diǎn)交匯定位的思路，計(jì)算出空間飛行器的所在位置。在此基礎(chǔ)上，分別利用伴隨同化方法和數(shù)據(jù)擬合方法給出了空間飛行器的位置和速度，并給出位置的估計(jì)殘差。兩種方法得到的結(jié)果均與觀測(cè)數(shù)據(jù)吻合較好，使用伴隨同化方法的誤差比使用數(shù)據(jù)擬合方法的小一個(gè)量級(jí)。問題三文基礎(chǔ)上考慮系統(tǒng)誤差對(duì)觀測(cè)數(shù)據(jù)產(chǎn)生的影響，根據(jù)對(duì)系統(tǒng)誤差的認(rèn)識(shí)和模型假設(shè)，認(rèn)為逐點(diǎn)交匯定位方法可以對(duì)系統(tǒng)誤差進(jìn)行估計(jì)，推導(dǎo)出了觀測(cè)數(shù)據(jù)和系統(tǒng)誤差之間的轉(zhuǎn)化關(guān)系。利用最的d 、d 和

3、d ）。小二乘原理，結(jié)合觀測(cè)數(shù)據(jù)得到系統(tǒng)誤差（兩顆利用系統(tǒng)誤差對(duì)原觀測(cè)數(shù)據(jù)進(jìn)行校正后，采用伴隨同化方法給出了空間飛行器的位置和速度估計(jì)，并計(jì)算了位置的估計(jì)殘差。2基于伴隨同化方法的空間飛行器軌道估計(jì)與誤差分析1. 問題重述空間信息對(duì)抗是現(xiàn)代信息中重要的和保障，是防御體系的重要組成部分, 目前已成為國(guó)際前沿研究熱點(diǎn)問題之一。利用中低軌近圓軌道衛(wèi)星探測(cè)空間飛行器發(fā)射與軌道參數(shù)是實(shí)現(xiàn)對(duì)飛行器和作出反應(yīng)的第一步。第一問要求基于基礎(chǔ)坐標(biāo)系下觀測(cè)主動(dòng)段的簡(jiǎn)化運(yùn)動(dòng)方程和 09 號(hào)在零時(shí)刻的位置和速度信息計(jì)算其在 50.0s、100.0s、150.0s、200.0s、250.0s五個(gè)時(shí)刻的三維位置。第二問要求

4、利用 06 號(hào)觀測(cè)和 09 號(hào)觀測(cè)對(duì) 0 號(hào)空間飛行器的仿真觀測(cè)數(shù)據(jù)，在空間飛行器的簡(jiǎn)化運(yùn)動(dòng)方程框架下建立模型給出 0號(hào)空間飛行器的軌道估計(jì)，并給出從 50s 到 170s 間隔 10 s 的采樣時(shí)間點(diǎn)的位置和速度，同時(shí)給出估計(jì)殘差，還要描繪三個(gè)位置分量和三個(gè)速度分量對(duì)時(shí)間的變化曲線。第三問在第二問的基礎(chǔ)上同時(shí)考慮系統(tǒng)誤差，要求在完成第二問問題的同時(shí)找到對(duì)系統(tǒng)誤差進(jìn)行正確估計(jì)的方案及其估計(jì)結(jié)果。第四問探討只有 09 號(hào)觀測(cè)對(duì) 01 號(hào)空間飛行器進(jìn)行軌道估計(jì)，完成第三問問題，并進(jìn)一步考慮在同時(shí)有多顆觀測(cè)觀測(cè)多個(gè)空間飛行器的情況下的聯(lián)合系統(tǒng)誤差估計(jì)問題。2. 問題分析問題一：根據(jù)理論力學(xué)中的兩體問

5、題的解，觀測(cè)的運(yùn)動(dòng)軌跡是橢圓?；诹銜r(shí)刻的位置和于基礎(chǔ)坐標(biāo)系下觀測(cè)主動(dòng)段的簡(jiǎn)化運(yùn)動(dòng)方程和 09 號(hào)速度信息，采用常微分方程初值問題的數(shù)值解法求解 09 號(hào)觀測(cè)的位置。問題二：由測(cè)數(shù)據(jù)中包含隨機(jī)誤差，首先進(jìn)行濾波消除該誤差。其次要解決雙星逐點(diǎn)交匯定位，其中包含觀測(cè)數(shù)據(jù)從觀測(cè)坐標(biāo)系向基礎(chǔ)坐標(biāo)系的轉(zhuǎn)化和如何確定空間飛行器的空間位置。最后根據(jù)已有的觀測(cè)數(shù)據(jù)選擇適當(dāng)?shù)哪Ｐ湍M空間飛行器的運(yùn)動(dòng)軌跡。應(yīng)當(dāng)注意兩顆二時(shí)，應(yīng)當(dāng)在濾波后將觀測(cè)數(shù)據(jù)插值為時(shí)間觀測(cè)時(shí)間不同步，在解決問題的序列。問題三：首先考慮系統(tǒng)誤差對(duì)觀測(cè)的影響，推導(dǎo)出觀測(cè)值、精確值與d 、d 、d 的轉(zhuǎn)化關(guān)系。進(jìn)而利用多組觀測(cè)數(shù)據(jù)進(jìn)行最小二乘擬合，

6、求出d 、d 和d 。對(duì)觀測(cè)數(shù)據(jù)進(jìn)行校正后，重復(fù)問題二的過程。問題四：針對(duì)定位問題，假設(shè)軌道的數(shù)學(xué)表達(dá)形式，采用待定系數(shù)法和最小二乘法確定待定系數(shù)，從而確定空間飛行器的運(yùn)動(dòng)參數(shù)。對(duì)于多顆觀測(cè)多個(gè)空間飛行器而言，可以借鑒問題三中的方法計(jì)算最小二乘意義下的系統(tǒng)誤差。33. 符號(hào)說明與3.1 符號(hào)說明r(t) ：空間飛行器在基礎(chǔ)坐標(biāo)系下的位置矢量；r (t) ： r(t) 對(duì)時(shí)間t 的二階導(dǎo)數(shù)，即加速度；G ：地球引力常數(shù)（ G 3.9860051014m3 / s2 ）；mmv (t) ：相對(duì)于火箭尾部噴口的噴射速度；rm(t) ：瞬時(shí)質(zhì)量， m (t) 表示質(zhì)量變化率；oldold、： 濾 波前

7、 觀測(cè) 對(duì) 于空 間飛 行器的 無量 綱化 簡(jiǎn)數(shù)據(jù) （old zsold ysoldxs ;oldxsold）； newnew、：濾波后觀測(cè) 對(duì)于空間飛行 器 的無量綱化簡(jiǎn)數(shù)據(jù) new ysnew zsnewxs;newxs）；new ：兩條異面直線 MA 與 NB 的夾角；1 ：06 號(hào)觀測(cè)線的夾角；與 09 號(hào)觀測(cè)的連線同 06 號(hào)觀測(cè)與空間飛行器的連2 ：06 號(hào)觀測(cè)線的夾角。與 09 號(hào)觀測(cè)的連線同 09 號(hào)觀測(cè)與空間飛行器的連3.2白噪聲：平均值函數(shù)為 0、自相關(guān)函數(shù)為 函數(shù)的隨機(jī)過程；估計(jì)殘差：觀測(cè)數(shù)據(jù)與模擬值之間的差；伴隨同化方法：伴隨同化方法是變分原理和最優(yōu)控制論相結(jié)合的所要解

8、決的實(shí)際問題作為條件最小值問題來解決；法,它將4. 模型假設(shè)(1)在存在各種誤差的條件下，每?jī)深w觀測(cè)與空間飛行器之間的連線，不一定在同一平面內(nèi)，假設(shè)兩條異面直線的公垂線段的中點(diǎn)為空間飛行器所在的位置；(2)假設(shè)隨機(jī)誤差為直接疊加在觀測(cè)數(shù)據(jù)上的白噪聲；(3)假設(shè)系統(tǒng)誤差只考慮與相關(guān)的誤差，即不同觀測(cè)4的系統(tǒng)誤差相互沒有關(guān)聯(lián)，同一觀測(cè)對(duì)不同空間飛行器的系統(tǒng)誤差是一樣的；(4)假設(shè)去掉隨機(jī)誤差的觀測(cè)在逐點(diǎn)交匯定位空間飛行器時(shí)，異面直線公垂線的長(zhǎng)度可以表征兩顆系統(tǒng)誤差。5. 模型建立與求解5.1 問題一求解Gm r t 3的運(yùn)動(dòng)方程rt r t 對(duì)應(yīng)的運(yùn)由理論力學(xué)的知識(shí)可知，觀測(cè)動(dòng)應(yīng)當(dāng)為橢圓運(yùn)動(dòng)，而地

9、心則是橢圓的其中一個(gè)焦點(diǎn)1。采用常微分方程初值問題的數(shù)值解法2可以得到 09 號(hào)觀測(cè)的位置。將觀測(cè)的簡(jiǎn)化運(yùn)動(dòng)方程分解到 x 、 y 、 z 三個(gè)方向：d 2 x dt 2 d 2 y dt 2 d 2 zdt 2G mxr3G myr3G（5-1） mzr3r x2 y2 z2采用如下的差分格式進(jìn)行離散：M n1 2M n M n1G m M ，其中， r x y z， M222nnnnn代表 x ，rn 3ty 或 z 。由于該差分方程是一個(gè)三層時(shí)間模式，所以需要用初始時(shí)刻的位置和速度先計(jì)算下一時(shí)刻的位置，才能完成迭代計(jì)算。即 x1 x0 u0ty y v 100t（5-2）z1 z0 w0

10、t4 2G運(yùn)動(dòng)（近似圓周運(yùn)動(dòng)）的方程式： m r，結(jié)合根據(jù)物理學(xué)中關(guān)于r2T 2初始時(shí)刻 09 號(hào)觀測(cè)檢驗(yàn)數(shù)值解的準(zhǔn)確性，步。的位置，估計(jì) 09 號(hào)觀測(cè)的周期大致為 9000s。為了模擬了 50000s，時(shí)間步長(zhǎng)取 0.1s，共模擬了 500000圖 1 和圖 2 分別給出了數(shù)值解的結(jié)果。通過圖 1 和圖 2 可以看出，數(shù)值求解計(jì)算得到的軌跡是一個(gè)橢圓，與理論結(jié)果吻合，且數(shù)值計(jì)算得到的 09 號(hào)觀5測(cè)在基礎(chǔ)坐標(biāo)系中的位置具有明顯的周期特性，周期9028s 左右，與估的位置。09 號(hào)觀測(cè)計(jì)結(jié)果基本吻合，因此數(shù)值求解結(jié)果可以代表 09 號(hào)觀測(cè)在 50.0s、100.0s、150.0s、200.0s

11、、250.0s 的三維位置見表 1。圖 1數(shù)值求解得到的 09 號(hào)觀測(cè)的運(yùn)動(dòng)軌跡圖 2 數(shù)值求解得到的 09 號(hào)觀測(cè)在基礎(chǔ)坐標(biāo)系中的坐標(biāo)隨時(shí)間的變化6表 109 號(hào)觀測(cè)在基礎(chǔ)坐標(biāo)系下的三維坐標(biāo)時(shí)間（s）X（m）Y（m）Z（m）50.0.0.0.0.01.77327E+061.50108E+061.22716E+069.51807E+056.75353E+058.16133E+068.12670E+068.08263E+068.02918E+067.96642E+064.51704E4.68502E4.84755E5.00446E5.15556E+065.2 問題二的建模與求解5.2.1濾波消除

12、白噪聲本題假設(shè)隨機(jī)誤差為直接疊加在觀測(cè)數(shù)據(jù)上的白噪聲，白噪聲指的是一個(gè)白色隨機(jī)過程 t ，它的平均值函數(shù)與自相關(guān)函數(shù)滿足以下條件： t E t 0R t1,t2 E t1 t2 在使用觀測(cè)數(shù)據(jù)之前，需要消除白噪聲對(duì)實(shí)際數(shù)據(jù)的影響。（5-3）（5-4）使用不同的低通濾波器3對(duì)觀測(cè)數(shù)據(jù)進(jìn)行濾波分析，利用得到的新數(shù)據(jù)與觀測(cè)數(shù)據(jù)作差，差的平均值越接近 0，認(rèn)為越符合白噪聲的特點(diǎn)，把此時(shí)的差作為白噪聲。不同低通濾波器差別體現(xiàn)在權(quán)重系數(shù) 的選取，分別選取 0.8, 0.85, 0.9, 0.95 對(duì)給出數(shù)據(jù)進(jìn)行濾波處理，濾波后的數(shù)據(jù)new 和 new 與原始數(shù)據(jù)中old 和 old 之間符合如下關(guān)系式：i

13、 1 i 1 i 1 （5-5）newold2i 1 （5-6）newold2其中，i 2,., n 1。經(jīng)過比較，當(dāng) 0.95 時(shí)濾波效果最好，數(shù)據(jù) 和 的誤差平均值分別為-7.79E-09 和-3.31E-09。白噪聲變化曲線如下圖（清晰起見，僅列出數(shù)據(jù) 上施加的白噪聲，數(shù)據(jù) 有類似結(jié)果）。7圖 309 號(hào)觀測(cè)觀測(cè)值 的白噪聲5.2.2插值觀測(cè)數(shù)據(jù)的時(shí)間序列并非從 50s 開始，而且兩個(gè)觀測(cè)的觀測(cè)時(shí)間序列不同步，為方便對(duì)數(shù)據(jù)處理，需要插值為時(shí)間的序列。09 號(hào)觀測(cè)對(duì) 0號(hào)空間飛行器的觀測(cè)數(shù)據(jù) 和 隨時(shí)間的變化見圖 4。圖 409 號(hào)觀測(cè)觀測(cè)值（ 和 ）隨時(shí)間變化曲線從圖 4 中可以看出，觀

14、測(cè)數(shù)據(jù) 和 均隨時(shí)間緩慢變化，故此處可以使用分段線性插值方法2將觀測(cè)數(shù)據(jù)插值為從 50s 開始的、與給定觀測(cè)時(shí)間間隔相同8（0.2s）的時(shí)間序列。5.2.3坐標(biāo)系轉(zhuǎn)換將觀測(cè)坐標(biāo)系中的觀測(cè)轉(zhuǎn)化到基礎(chǔ)坐標(biāo)系中是實(shí)現(xiàn)雙星逐點(diǎn)交匯定位的基礎(chǔ)4。這一部分給出了坐標(biāo)轉(zhuǎn)換的基本思路及其推導(dǎo)公式。如圖 5 所示， Oc - XcYc Zc 表示基礎(chǔ)坐標(biāo)系，Os - XsYs Zs 表示觀測(cè)坐標(biāo)系，Q表示觀測(cè)坐標(biāo)系中的任意一點(diǎn)。矢量Os Q 可以分解到Os X s 、OsYs 和Os Zs 方向，其分量分別為Qx 、Qy 和Qz 。由于Os xs 也即OcOs的方向是確定的，并且Os zs 向北， Os ys

15、向東，因此可以將Qx 、Qy 和Qz 在基礎(chǔ)坐標(biāo)系中分解，不妨設(shè)： Qx a1 i b1 j c1 kQy a2 i b2 j c2 k（5-7）Qz a i b j c k其中， i 、 j 和k 分別是基礎(chǔ)坐標(biāo)系的333向量。所以有：OsQ a1 a2 a3 i b1 b2 b3 j c1 c2 c3 k（5-8）再結(jié)合Os 在基礎(chǔ)坐標(biāo)系中的坐標(biāo)，就可以給出Q 點(diǎn)在基礎(chǔ)坐標(biāo)系中的坐標(biāo)。圖 5觀測(cè)坐標(biāo)系向基礎(chǔ)坐標(biāo)系轉(zhuǎn)化示意圖設(shè)Os 在基礎(chǔ)坐標(biāo)系中的坐標(biāo)為x0 , y0 , z0 ，觀測(cè)到的空間飛行器的參數(shù)為和 。 和 提供了空間飛行器在觀測(cè)坐標(biāo)系中的方向信息，只要再知道該方向上的任意一點(diǎn)，結(jié)

16、合的位置（即Os ），就可以唯一確定出9和空間飛行器連線在基礎(chǔ)坐標(biāo)系中的精確表達(dá)式。當(dāng)存在兩顆及以上的點(diǎn)定位方法確定空間飛行器在基礎(chǔ)坐標(biāo)系中的位置。時(shí)，就可以利用逐設(shè)某點(diǎn)在觀測(cè)坐標(biāo)系中的坐標(biāo)為 x1, y1, z1 ，在基礎(chǔ)坐標(biāo)系中的坐標(biāo)為x2 , y2 , z2 。下面給出坐標(biāo)轉(zhuǎn)化的推導(dǎo)過程。向量Os xs 在基礎(chǔ)坐標(biāo)系中的法向量為：xyzv 0 0 0,（5-9）1x2 y2 z2x2 y2 z2x2 y2 z2000000000由于Os zs 指向北，所以 Os zs 在平面Os Oc zc 內(nèi)。從而平面 OsOc zc 的法向量即為Os ys 的方向向量，因此Os ys 在基礎(chǔ)坐標(biāo)系中

17、的法向量為：yxv 0 0,0（5-10）2x2 y 2x2 y 20000由于Os xs 和Os ys 已經(jīng)確定，所以O(shè)s zs 在基礎(chǔ)坐標(biāo)系中的法向量為：x z y z x2 y2v 0 0 0 0 00,3222x z y z x yx z y z x yx z y z x y2 22 2222 22 2222 22 2220 00 0000 00 0000 00 000（5-11）于是觀測(cè)坐標(biāo)系中的點(diǎn)x1, y1, z1 在基礎(chǔ)坐標(biāo)系中的坐標(biāo)為：x2 , y2 , z2 x1v1 y1v2 y2 v3 x0 , y0 , z0 （5-12）5.2.4雙星逐點(diǎn)交匯定位由于實(shí)際觀測(cè)中存在各

18、種誤差，每?jī)深w觀測(cè)與空間飛行器的連線不一定在同一個(gè)平面內(nèi)。在這種情況下，認(rèn)為兩條異面直線的公垂線段的中點(diǎn)為空間飛行器所在的位置，公垂線段的長(zhǎng)度可以表征兩顆的觀測(cè)誤差。如圖 6 所示，M 和 N 分別表示兩顆觀測(cè)的位置，直線MA 和 NB 分別表示從兩顆上觀測(cè)到的空間飛行器與所在的直線， AB 為兩直線的公垂線段， AB 的中點(diǎn)即空間飛行器的實(shí)際位置。10圖 6空間內(nèi)異面直線的及其公垂線示意圖下面給出空間飛行器實(shí)際位置的推導(dǎo)公式。M 和 N 的位置分別為x1, y1, z1 和x2 , y2 , z2 ， MA 和 NB 的設(shè)兩顆觀測(cè)向量分別為1, 1,1 和2 , 2 , 2 （可以通過 5.

19、2.3 得到）。則有 x x y y z z arccos121121121 （5-13）2 1x x y y 2 z z2212121 x x y y z z arccos212212212 （5-14）x x 2 y y 2 z z22121212 arccos12 12 1 2 由空間異面直線的定理可知：（5-15） cos1 cos2 cos x x y y z z2222MA（5-16）sin 2 212121NB cos2 cos1 cos x x y y 2 z z2（5-17）sin 2 212121進(jìn)而可以得到 Ax3 , y3 , z3 和 Bx4 , y4 , z4 點(diǎn)的

20、坐標(biāo)：x3 , y3 , z3 x1, y1, z1 MA 1, 1,1（5-18）x4 , y4 , z4 x2 , y2 , z2 NB 2, 2, 2 （5-19）則空間飛行器在基礎(chǔ)坐標(biāo)系中的位置x, y, z為：x, y, z , , 1x , y , z x , y, z MA , , NB （5-20）11 1222111222211采用上述方法處理 06 號(hào)和 09 號(hào)觀測(cè)對(duì) 0 號(hào)空間飛行器的觀測(cè)資料，可以得到 0 號(hào)空間飛行器在基礎(chǔ)坐標(biāo)系中的軌跡，分別如圖 7 和圖 8 所示。圖 70 號(hào)空間飛行器在基礎(chǔ)坐標(biāo)系中的軌跡圖 80 號(hào)空間飛行器三維坐標(biāo)隨時(shí)間的變化5.2.5伴隨同

21、化方法確定空間飛行器的軌道伴隨同化方法是建立在嚴(yán)格的數(shù)學(xué)基礎(chǔ)之上的一種有效的變分同化技術(shù)，它將變分原理和最優(yōu)控制理論相結(jié)合5。伴隨同化方法可以用來解決許多不同類型的實(shí)際問題，其中包括：對(duì)模型中的參數(shù)進(jìn)行優(yōu)化，對(duì)初始條件或邊界條件進(jìn)行優(yōu)化，對(duì)定常流或環(huán)流問題以及表面熱通量問題進(jìn)行研究等。伴隨同化方12法在本文中的應(yīng)用是以空間飛行器在基礎(chǔ)坐標(biāo)系下的簡(jiǎn)化運(yùn)動(dòng)方程作為約束條件，由觀測(cè)到的空間飛行器的位置，得到空間飛行器的軌道參數(shù)，進(jìn)而得到軌道信息。伴隨同化方法中的 Lagrange 乘子法是建立在嚴(yán)格的數(shù)學(xué)基礎(chǔ)之上的法，它將所要解決的實(shí)際問題作為條件極值問題來求解。在此過程中將動(dòng)力學(xué)方以及定解條件作為

22、約束條件，目標(biāo)是根據(jù)具體問題而設(shè)計(jì)的代價(jià)函數(shù)達(dá)到極小6。設(shè)依賴于時(shí)間變量的動(dòng)力學(xué)方x）為 F ( x, c)(5-21)t這里向量 x 的分量是隨時(shí)間變化的量（例如位移、速度等）， t 表示時(shí)間，向量c 表示模型中的參數(shù)， F 是線性或非線性算子。模型方程(5-21)不封閉，還需要相應(yīng)的定解條件才能確定其唯一解。 y 表示控制變量，它包括初始條件、邊界條件、模型中的參數(shù)等。y 一旦確定，模型(5-21)的唯一解 x( y) 也就隨之確定。當(dāng)然控制變量必須屬于一個(gè)可容許的控制集合 yad 。 yad 可根據(jù)有關(guān)控制變量的物理學(xué)性質(zhì)以及歷史經(jīng)驗(yàn)來確定。在解決問題的過程中，研究者的目的是得到模型方程

23、(5-21)的接近觀測(cè)結(jié)果的解 x 。這里的所謂接近用代價(jià)函數(shù)J (度量模型方程(5-21)的解與觀測(cè)結(jié)果之間的距離，通常定義為差的平方和)來定義。觀測(cè)結(jié)果包括位移、速度、高度等。*這樣，變分問題可簡(jiǎn)述為問題(P)：尋求屬于 yad 的 y 使得代價(jià)函數(shù) J 達(dá)到最小,y* 表示最優(yōu)的 y 。以上所述問題是以動(dòng)力學(xué)方程(5-21)作為約束的約束極值問題。此問題可轉(zhuǎn)化為無約束極值問題來求解。當(dāng)把約束條件作為強(qiáng)約束條件來對(duì)待時(shí)，經(jīng)典的Lagrange（條件極值）乘子法為函數(shù)L( x, , y) 。提供了堅(jiān)實(shí)的理論基礎(chǔ)。在此構(gòu)造 LagrangeL( x, , y) J ( x, y) , G(x;

24、 c) ， G( x; c) x F ( x; c) ，(5-22)t這里是 Lagrange 乘子，是定義在G( x;c) 所屬 Hilbert 空間的內(nèi)積。這樣約束極值問題(P)就被轉(zhuǎn)化成關(guān)于和 x , y 的無約束問題。在此需說明的是這里的約束條件并不是固定不變的，而是在同化過程中不斷得到調(diào)整，使得模擬結(jié)果盡可能地接近觀測(cè)結(jié)果。在強(qiáng)約束條件G( x;c) 0 的限制下，確定泛函 J ( x, y) 的駐點(diǎn)與確定 Lagrange 函數(shù)關(guān)于變量 x ,和 y 的駐點(diǎn)是等價(jià)的。表示 Lagrange 函數(shù)的駐點(diǎn)的方稱為約束最小值問題(P)的 Euler-Lagrange 方。Euler-La

25、grange最優(yōu)條件(最優(yōu)的 x*, *, y* )將由下述方確定：13L (x*, *, y*) 0 .(5-23)L (x*, , y.(5-24)xLy, ,*) 0 .(5-25)(方程(5-23)就是原來的模型方程。方程(5-24)中的算子是方程(5-23)中的算子的伴隨算子，通常稱方程(5-24)為方程(5-23)的伴隨方程。模型方程隨時(shí)間正向傳播信息，而伴隨方程隨時(shí)間反向數(shù)關(guān)于控制變量的梯度。信息。至于方程(5-25)，給出了 Lagrange 函Lagrange 乘子法，代價(jià)函數(shù)關(guān)于任何控制變量（初始條件、邊界條件、方程中的參數(shù)）的梯度在理論上都能夠表達(dá)出來，這決定了利用 La

26、grange 乘子法可。而且，在 Lagrange 乘子法中，模擬值與觀測(cè)值的差作以解決許多方面為驅(qū)動(dòng)伴隨方程的外力，很自然地加到了伴隨方程上7。Lagrange 乘子法與經(jīng)典的 Lagrange 條件極值問題的不同有以下幾點(diǎn)：經(jīng)典的 Lagrange 條件極值問題的約束是一個(gè)或者數(shù)個(gè)，而 Lagrange乘子法的約束每個(gè)時(shí)間步都存在一個(gè)；Lagrange 乘子法中的伴隨方程需要求解而不能簡(jiǎn)單消掉；Lagrange 乘子法需要求解控制方程和伴隨方程，以得到需要優(yōu)化參數(shù)的梯度；Lagrange 乘子法反復(fù)利用（3）中的梯度，選取適當(dāng)?shù)南陆邓惴?，不斷?yōu)化參數(shù)。假設(shè)火箭產(chǎn)生的推力加速度項(xiàng)在基礎(chǔ)坐標(biāo)系

27、中分解為 Fx 、 Fy 和 Fz ，且 Fx 、Fy 和 Fz 只與時(shí)間t 有關(guān)（在方程中作為隨時(shí)間變化的控制參數(shù)），所以空間飛行器在基礎(chǔ)坐標(biāo)系下的簡(jiǎn)化運(yùn)動(dòng)方程可以寫成如下的形式： d 2xG mr3xFx dt 2d 2 yG dt 2 d 2z y Fy mr3（5-26）G2 mr3z2 Fz dt 22通常代價(jià)函數(shù)可以構(gòu)造為： t2J (x x)2 ( y y)2 (z z)2 dt（5-27）2 t1其中 為常數(shù)（在本題中取為1）， x 、 x 、 z 代表觀測(cè)值， x 、 y 、 z 代表數(shù)14值模擬結(jié)果。利用 Lagrange 乘子法構(gòu)造 Lagrange 函數(shù)，如下所示：t2

28、 d 2x d 2 y d 2zGGGL J x Fx y Fy z Fz dtmmm dt 2r3 dt 2r3 dt 2r3t1 （5-28）（5-26）及其對(duì)應(yīng)的初邊值條件下，代價(jià)函要求：的目的是尋求在滿足方數(shù)達(dá)到極小。為求得極小值，L 0 ， L 0 ， L 0（5-29）L 0 ， L 0 ， L 0（5-30）xzyLFxLLFz 0 ， 0 ， 0（5-31）Fy根據(jù)（5-30），能夠得到伴隨方程：d x )2y(z 2y2z223x()Gm(x x)0 dt 2r5 y2) d x(z2x2z 223y)G ( y y)0（5-32） dt 2mr5z )2d (x y22 y

29、 223z(x) (z )zG0mdt 2r5利用梯度下降算法，對(duì)控制參數(shù)進(jìn)行修正，根據(jù)上述推導(dǎo)得到的校正關(guān)系表達(dá)式為： FFx xF Fy （5-33）y FFz z其中 、 Fy 和 是優(yōu)化前的值， FFxFz、 F和 F是優(yōu)化后的值。xyz15方程（5-26）的差分格式為： x ynFxnxyn（5-34） znznn3t 2nrn n 2 n 2n 2伴隨方程（5-32）的差分格式為：n 1 2n n 1t 21 2 n nt 2 nn 1n 1n2t 2（5-35）利用伴隨同化方法的工作流程如下：給 Fx 、 Fy 和 Fz 賦值（初始猜測(cè)）。積分正向方程（5-34）。計(jì)算表征模擬結(jié)

30、果與觀測(cè)結(jié)果差的平方和的代價(jià)函數(shù)： t2x( J 2 y )2yx2 t1其中 為常數(shù)， x 、 x 、 z 代表觀測(cè)值， x 、 y 、 z 代表相應(yīng)的模擬結(jié)果，當(dāng)代價(jià)函數(shù)下降到一定程度或者其變化率很小時(shí)，進(jìn)入步驟(8)；如果不滿足要求，執(zhí)行步驟(4)。反向積分伴隨方程。計(jì)算代價(jià)函數(shù)關(guān)于 Fx 、 Fy 和 Fz 的梯度。校正 Fx 、 Fy 和 Fz 。返回步驟（1），利用改進(jìn)后的 Fx 、 Fy 和 Fz 開始下一次校正。結(jié)束對(duì) Fx 、 Fy 和 Fz 的校正，得到與觀測(cè)值充分靠近的數(shù)值模擬結(jié)果。利用伴隨同化方法得到從 50 s 到 170s 間隔 10s 的采樣點(diǎn)的位置和速度及位置的

31、估計(jì)殘差，如表 2-4。圖 7 和圖 8 分別表示 0 號(hào)空間飛行器在采樣時(shí)間范圍內(nèi)的位置和速度曲線。16表 20 號(hào)空間飛行器的位置模擬值Time（s）X（m）Y（m）Z（m）5060708090-1.11202E+06-1.12054E+06-1.13074E+06-1.14259E+06-1.15630E+06-1.17195E+06-1.18968E+06-1.20962E+06-1.23194E+06-1.25681E+06-1.28429E+06-1.31487E+06-1.34803E+066.20027E+066.20787E+066.21595E+066.22459E+066

32、.23364E+066.24312E+066.25299E+066.26322E+066.27378E+066.28469E+066.29591E+066.30747E+066.31933E+061.13305E+061.14326E+061.15539E+061.16960E+061.18602E1.20482E1.22614E1.25015E1.27705E1.30704E1.34027E1.37718E1.41750E+06表 30 號(hào)空間飛行器的速度模擬值u （m/s）v （m/s）w （m/s）Time（s）V （m/s）5060708090-6.96939E+02-9.53603E

33、+02-1.09724E+03-1.27859E+03-1.46847E+03-1.66910E+03-1.88317E+03-2.11213E+03-2.36035E+037 E+03- 2.89380 E+03- 3.22708 E+03- 3.36689 E+038.17465E+027.66832E+028.42609E+028.84760E+029.27313E+029.68954E+021.00522E+031.03991E+031.07414E+031.10637 E+031.13713 E+031.17957 E+031.15585 E+039.40316E+021.11394

34、E+031.31747E+031.53111E+031.76059E+032.00587E+032.26632E+032.54464E+032.84481E+033.16143 E+033.49719 E+033.90042 E+034.08279 E+031.42764E+031.65477E+031.91041E+032.18217E+032.47305E2.78357E3.11336E3.46665E3.84941E4.25126E4.68062E5.19592E5.42017E+0317表 40 號(hào)空間飛行器的位置估計(jì)殘差Time（s）X（m）Y（m）Z（m）5060708090-4.

35、14400E-01 1.72815E+01-2.56023E+01 1.33891E+011.38297E+01-1.07119E+01 1.46210E-011.22921E+011.84681E+00-2.10827E+01 5.31117E+01-7.14229E+014.15916E+02-2.81300E-02 1.35528E+01-2.04622E+01 5.99139E+006.07051E+00-3.74736E+00 6.14745E+00-5.42108E+00 2.17123E+013.11782E+005.64645E+004.98970E-01-2.18496E+0

36、1-4.11980E-01 5.68476E+003.25836E+001.08001E+013.87634E-2.52564E-1.05900E-01 3.62984E7.56858E1.00273E-4.12539E 3.05004E-2.90209E+02120圖 90 號(hào)空間飛行器在采樣時(shí)間內(nèi)的位置曲線18圖 100 號(hào)空間飛行器在采樣時(shí)間內(nèi)的速度模擬值5.2.6擬合方法確定空間飛行器的飛行軌跡本題中給出了觀測(cè)在 600 個(gè)時(shí)間采樣點(diǎn)處對(duì)空間飛行器的觀測(cè)數(shù)據(jù)，考慮使用擬合方法2得出空間飛行器的飛行軌跡。根據(jù)上文的逐點(diǎn)交匯定位思路，結(jié)合觀測(cè)對(duì)空間飛行器的觀測(cè)數(shù)據(jù) 和 ，可以推導(dǎo)出空間飛行

37、器在上述時(shí)間采樣點(diǎn)處的空間位置觀測(cè)數(shù)據(jù)X0 ,Y0 , Z 0。分別對(duì)三個(gè)空間位置坐標(biāo)分量進(jìn)行多項(xiàng)式擬合可以得到空間位置分量隨時(shí)間變化的函數(shù) X t、Y t和Zt，由此可以對(duì)空間飛行器在 50.0s 到 170.0s采樣時(shí)間點(diǎn)的位置進(jìn)行估計(jì)。根據(jù)軌道參數(shù)的匹配原理8，空間飛行器的速度分量變化分別是空間位置變化函數(shù)對(duì)時(shí)間的導(dǎo)數(shù)，即Vx t X t，Vy t Y t ，Vz t Zt，通過以上方程可以對(duì)空間飛行器在各采樣點(diǎn)處的速度給出估計(jì)。將簡(jiǎn)化運(yùn)動(dòng)方程表示成坐標(biāo)分量形式，以 x 方向?yàn)槔?，原運(yùn)動(dòng)方程可表示為如下等式（方便起見，令運(yùn)動(dòng)方程右端第二項(xiàng)vt m t F ）：tmt Gma t X t

38、F t（5-36）xx3X t2 Y t2 Z t2其中， ax t X t ，那么 Fx t 的函數(shù)表達(dá)式可以求得，同樣可以求出 Fy t 和 可以確定，這對(duì)研究不同空F t，則對(duì)應(yīng)數(shù)據(jù)擬合結(jié)果的vt和mt關(guān)系式Ftz19間飛行器之間的差異有重要意義。根據(jù) 06 號(hào)和 09 號(hào)觀測(cè)對(duì) 0 號(hào)空間飛行器的觀測(cè)數(shù)據(jù)，由上述方法擬合得到空間飛行器的位置變化函數(shù)分別為：X t 11.2330t 2 536.781t 112750 Y t 1.88221t 2 583.230t 6166030 Z t 13.6602t 2 676.400t 1135180對(duì)比觀測(cè)數(shù)據(jù)變化趨勢(shì)與擬合方程可以得到如下圖像

39、。（5-37）圖 110 號(hào)空間飛行器在采樣時(shí)間內(nèi)的位置曲線擬合結(jié)果從圖像中可以看到，擬合結(jié)果與觀測(cè)結(jié)果幾乎完全吻合，二者之間的平均絕對(duì)誤差與平均相對(duì)誤差見表 5。表 5擬合結(jié)果的誤差統(tǒng)計(jì)表分量XYZ平均絕對(duì)誤差（102）平均相對(duì)誤差（10-4）6.976725.792951.008210.1613128.566426.89172通過上表統(tǒng)計(jì)的誤差可以看出，擬合結(jié)果的絕對(duì)誤差在百米量級(jí)，相對(duì)測(cè)數(shù)據(jù)的誤差在萬分之一量級(jí)。計(jì)算 0 號(hào)空間飛行器在各個(gè)采樣時(shí)間點(diǎn)的位置，具體結(jié)果見表 6。20表 60 號(hào)空間飛行器在采樣時(shí)間點(diǎn)的位置擬合值Time（s）X（m）Y（m）Z（m）5060708090-1.

40、11399E+06-1.12098E+06-1.13022E+06-1.14170E+06-1.15543E+06-1.17140E+06-1.18962E+06-1.21009E+06-1.23281E+06-1.25777E+06-1.28498E+06-1.31443E+06-1.34613E+066.19990E+066.20780E+066.21608E+066.22473E+066.23377E+066.24317E+066.25296E+066.26312E+066.27366E+066.28457E+066.29586E+066.30753E+066.31957E+061.13

41、551E+061.14378E+061.15477E+061.16850E+061.18496E1.20415E1.22607E1.25072E1.27811E1.30823E1.34108E1.37666E1.41498E+06分析在各個(gè)采樣時(shí)間點(diǎn)的擬合結(jié)果和對(duì)應(yīng)觀測(cè)數(shù)據(jù)可以給出各個(gè)時(shí)間點(diǎn)的位置估計(jì)殘差 X（m）、Y（m）和 Z（m），具體數(shù)值見下表：表 70 號(hào)空間飛行器的位置擬合值殘差Time（s）X（m）Y（m）Z（m）5060708090-1.97917 E+03-4.30955 E+024.96079 E+028.99390 E+028.82551 E+025.39132 E+02

42、5.85129 E+01-4.59417 E+02-8.61536 E+02-9.79739 E+02-6.29416 E+023.70579 E+022.30985 E+03-3.69853 E+02-6.05458 E+011.12606 E+021.52668 E+021.31800 E+024.84698 E+01-2.67036 E+01-9.98992 E+01-1.00952 E+02-1.09291 E+02-3.70661 E+016.59240 E+012.20154 E+022.45956 E+035.22928 E+02-6.14181 E+02-1.09027 E+0

43、3-1.06178 E-6.97452 E-6.97895 E5.77507 E1.07006 E1.19917 E7.68599 E-4.83890 E-2.81579 E+03同樣的可以給出采樣時(shí)間點(diǎn)的速度變化，圖 12 反映了三個(gè)速度分量在給定時(shí)間范圍內(nèi)的變化趨勢(shì)。21圖 120 號(hào)空間飛行器在采樣時(shí)間內(nèi)的速度曲線很顯然，三個(gè)方向的速度分量隨時(shí)間變化的曲線都近似直線，其方程分別為：ut 22.4660t 536.781vt 3.764424 583.230wt 27.3204t 676.400（5-38）根據(jù)速度分量的方程可以得到空間飛行器在采樣時(shí)間點(diǎn)處的速度，具體數(shù)值見表 8。表 80

44、 號(hào)空間飛行器的速度擬合值u （m/s）v （m/s）w （m/s）Time（s）V （m/s）5060708090-5.86518E+2-8.11178E+2-1.03584E+3-1.26050E+3-1.48516E+3-1.70982E+3-1.93448E+3-2.15914E+3-2.38380E+3-2.60846E+3-2.83312E+3-3.05778E+3-3.28246E+37.71451E+28.09095E+28.46739E+28.84383E+29.22027E+29.59671E+29.97315E+21.03496E+31.07260E+31.11025E+

45、31.14789E+31.18554E+31.22318E+36.89624E+29.62829E+21.23603E+31.50924E+31.78244E+32.05565E+32.32885E+32.60206E+32.87526E+33.14847E+33.42167E+33.69488E+33.96808E+31.18942E+031.49656E+031.82146E+032.15611E+032.49658E2.84080E3.18754E3.53606E3.88589E4.23669E4.58825E4.94041E5.29304E+03225.2.7兩種方法的比較通過伴隨同化

46、方法和數(shù)據(jù)擬合的方法都實(shí)現(xiàn)了對(duì)空間飛行器位置和速度的估計(jì)，兩個(gè)方法得到的最終結(jié)果十分接近，而且結(jié)果與觀測(cè)之間的相對(duì)誤差都小于萬分之一。但是這兩個(gè)方法一種是數(shù)值模式計(jì)算，另一種則是數(shù)據(jù)擬合，二者各有優(yōu)勢(shì)，下面對(duì)兩種方法得到的結(jié)果進(jìn)行比較。下圖反映了伴隨同化方法和數(shù)據(jù)擬合方法分別得到的 0 號(hào)空間飛行器的運(yùn)行速度比較。為方便比較，根據(jù)軌道參數(shù)的匹配原理，結(jié)合觀測(cè)數(shù)據(jù)可以得到空間飛行器的速度，將該速度作為兩種方法所得結(jié)果的速度參考值，三者在同一坐標(biāo)系下進(jìn)行比較。圖 13兩種方法得到的速度比較從圖 13 可以看到，兩種方法得到的速度均與利用觀測(cè)數(shù)據(jù)生成的速度接近，下面分別給出兩種方法結(jié)果的誤差，分析伴

47、隨同化方法與數(shù)據(jù)擬合方法結(jié)果的精度。表 9兩種方法速度結(jié)果的誤差平均絕對(duì)誤差(m/s)平均相對(duì)誤差(%)速度3.39348E+014.31213E+014.24387E+012.0764.5022.081u vw伴隨同化方法7.58509E+014.38668E+019.28405E+014.7574.6254.853u vw擬合方法由上表可知，兩種方法得到速度的結(jié)果誤差都在 100m/s 以下，相對(duì)測(cè)生成的速度數(shù)據(jù)相對(duì)誤不超過 5%，伴隨同化方法得到的各個(gè)速度分量的誤小于數(shù)據(jù)擬合方法得到的誤差。因此，在第 3 和第 4 題的求解過程中，都采用伴隨同化方法求解。235.3 問題三的建模與求解采

48、用 5.2.1 節(jié)中的濾波方案可以消除隨機(jī)誤差。因此，在本節(jié)的中，所用到的數(shù)據(jù)是消除了隨機(jī)誤差的數(shù)據(jù)。設(shè)包含系統(tǒng)誤差的觀測(cè)為1, 1 ，即觀測(cè)實(shí)際觀測(cè)到的數(shù)據(jù)；不包含系統(tǒng)誤差的觀測(cè)為0 , 0 ，最終用于求解空間飛行器軌道的數(shù)據(jù)?？紤]到系統(tǒng)誤差由兩個(gè)平移誤差d 和d ，一個(gè)旋轉(zhuǎn)誤差d組成， 0 , 0 和1, 1 應(yīng)當(dāng)滿足如下的關(guān)系：0 1 cos d 1 sin d d（5-39） sin d cos d d 011又題目中，d 、d 和d 是三個(gè)常值小量，因此可以認(rèn)為sin d 和cosd 滿足： cosd 1， sin d d 。于是公式（5-34）可以簡(jiǎn)化為：0 1 1d d（5-40

49、） d d 011根據(jù)假設(shè)，可以計(jì)算出在1, 1 的觀測(cè)下計(jì)算得到的空間飛行器在基礎(chǔ)坐標(biāo)系中的位置，不妨設(shè)為 Px , y , z 。同一時(shí)刻觀測(cè)的位置為O x , y, z ，ppps000那么可以得到O P x x , y y , z z ，利用 5.2.3 節(jié)中的觀測(cè)坐標(biāo)系的單sp0p0p0位方向向量v1 、v2 和v3 ，可以得到Os P 在觀測(cè)坐標(biāo)系中三個(gè)方向上的分量長(zhǎng)度，即 x0 s 、 y0 s 和 z0s ： x0s Os P v1 Os P v2 y0s（5-41） z O P v0ss3進(jìn)而可以得到：0 y0sx0s0s（5-42） zx 00s聯(lián)立（5-40）、（5-41

50、）和（5-42），就可以得到系統(tǒng)誤差d 、d 、d 與已知量x , y, z 和 , 的關(guān)系。根據(jù)已有的觀測(cè)資料，利用最小二乘原理 ，200011可以求得d 、d 和d 的大小。表 10 給出了在現(xiàn)有的觀測(cè)下確定的 06 和 09號(hào)的系統(tǒng)誤差。24表 1006 和 09 號(hào)的系統(tǒng)誤差ddd06092.89070E-04-3.50795E-043.23265E-044.69959E-04-3.45788E-044.68311 E-04利用公式（5-35）去除觀測(cè)中的系統(tǒng)誤差后，采用 5.2.5 節(jié)中的方法重新計(jì)算空間飛行器的位置和速度，結(jié)果如表 11表 13、圖 14、圖 15 所示。表 110

51、 號(hào)空間飛行器的位置模擬值Time（s）X（m）Y（m）Z（m）5060708090-1.11261E+06-1.12110E+06-1.13128E+06-1.14311E+06-1.15680E+06-1.17244E+06-1.19015E+06-1.21008E+06-1.23238E+06-1.25723E+06-1.28471E+06-1.31526E+06-1.34850E+066.19825E+066.20588E+066.21398E+066.22265E+066.23173E+066.24123E+066.25113E+066.26137E+066.27196E+066.28289E+066.29411E+066.30571E+066.31750E+061.13401E+061.14421E+061.15633E+061.17054E+061.18695E1.20574E1.22706E1.25106E1.27795E1.30795E1.34116E1.37808E1.41831E+06表 120 號(hào)空間飛行器的速度模擬值Time（s）u （m/s）v （m/s）w （m/s）V （m/s）5060708090-6.94405E+02-9.51274E+02-1.09502E+03-1.27650E+03-1.46652