空中飛行器無源定位

1、評 分 專 用 頁編號：學(xué)員評閱記錄：評閱隊號評分備注教員評閱記錄：空中飛行器無源定位摘 要目標(biāo)定位技術(shù)是導(dǎo)航與制導(dǎo)技術(shù)的重要基礎(chǔ)。在現(xiàn)有的導(dǎo)航與制導(dǎo)技術(shù)中，衛(wèi)星定位技術(shù)是精度最高的，也是較為理想的導(dǎo)航與制導(dǎo)技術(shù)。本文研究了利用同步衛(wèi)星確定空中飛行器的位置參數(shù)及如何提高定位精度及定位效率的優(yōu)化選擇問題。針對問題一，根據(jù)測向陣列方向和地球同步衛(wèi)星夾角與飛行器位置的關(guān)系，建立了基于最小二乘法的測量、計算誤差平方和目標(biāo)函數(shù)模型，運(yùn)用最小二乘法去逼近真實(shí)參數(shù)，當(dāng)目標(biāo)函數(shù)最小時，所得數(shù)據(jù)為最接近真實(shí)參數(shù)的數(shù)據(jù)。在此基礎(chǔ)上，構(gòu)造了以余弦夾角法分析的聚類模型，剔除誤差較大的數(shù)據(jù)，以剩余數(shù)據(jù)為有效解，得到了更

2、加接近真實(shí)參數(shù)的數(shù)據(jù)，是模型進(jìn)一步得到優(yōu)化。針對問題二，運(yùn)用問題一構(gòu)建的數(shù)學(xué)模型，求得五個不同時刻飛行器的位置坐標(biāo)，接著采用最小二乘法曲線擬合的方法擬合出值、值、值與時間的函數(shù)關(guān)系并計算其擬合曲線的皮爾遜相關(guān)系數(shù)，驗(yàn)證擬合曲線可靠性，進(jìn)而解得時飛行器的空間位置坐標(biāo)。然后將不同時刻的擬合值與計算所得值之間的空間距離作為原始非負(fù)時間序列，建立模型，對距離進(jìn)行灰色預(yù)測，得到其距離誤差的預(yù)測值，最后得出時飛行器的距離誤差值作為其預(yù)測的可靠度。針對問題三，本文以幾何精度因子作為衛(wèi)星分布對定位精度影響的指標(biāo)。通過分析衛(wèi)星數(shù)目的影響,得到了隨衛(wèi)星數(shù)目增加單調(diào)遞減,但遞減幅度逐漸變小的變化規(guī)律。綜合考慮衛(wèi)星定

3、位精度和定位效率認(rèn)為6星組合方案最為適宜。計算所有的6星組合方案，找出其中最小的衛(wèi)星組合作為最終的優(yōu)選方案。在測量角度存在誤差限的情況下，我們對附表中的數(shù)據(jù)采用加入方差為正態(tài)分布的方法，然后進(jìn)行500次仿真實(shí)驗(yàn)，觀察加入方向誤差后產(chǎn)生的定位誤差。統(tǒng)計定位誤差，發(fā)現(xiàn)加入噪聲后有的定位結(jié)果誤差在以內(nèi)，即定位精度可以認(rèn)為。在本文的最后，針對每個問題對其結(jié)果進(jìn)行了分析、對每個問題解決方法的優(yōu)缺點(diǎn)進(jìn)行了分析，并提出了相應(yīng)的改進(jìn)方案。本文的特色在于運(yùn)用灰色預(yù)測模型對飛行器位置的定位距離誤差進(jìn)行科學(xué)的分析，將誤差的變化考慮進(jìn)模型，進(jìn)而提高了預(yù)測的精度。關(guān)鍵詞：無源定位 最小二乘法 灰色模型 幾何精度因子1

4、問題的重述目標(biāo)定位技術(shù)是導(dǎo)航與制導(dǎo)技術(shù)的重要基礎(chǔ)。在現(xiàn)有的導(dǎo)航與制導(dǎo)技術(shù)中，衛(wèi)星定位技術(shù)是精度最高的，也是較為理想的導(dǎo)航與制導(dǎo)技術(shù)。目前，較為成熟的衛(wèi)星導(dǎo)航系統(tǒng)有GPS系統(tǒng)、Galileo系統(tǒng)等。衛(wèi)星定位的基本原理是目標(biāo)接收機(jī)通過接收多顆衛(wèi)星的信號測量出目標(biāo)距各衛(wèi)星的距離（偽距），再通過一定的計算確定出目標(biāo)的位置。 對于空中飛行器，在其飛行過程中很容易接收到太空衛(wèi)星的信號?，F(xiàn)在考慮通過測量飛行器與地球同步衛(wèi)星的方向角來實(shí)現(xiàn)空中飛行器的自定位。在球心坐標(biāo)系下，空中飛行器P 的空間坐標(biāo)記為(x,y,z)，不妨設(shè)它同時能接收到N 顆同步衛(wèi)星的信號，其N 顆同步衛(wèi)星Xi的空間坐標(biāo)分別記為(xi,yi,

5、zi)(i=1,2,N)。為了方便檢測與同步衛(wèi)星的方向角，在空中飛行器上固定安裝了兩個相互垂直的測向陣列，它們的指向分別為d1(d1x,d1y,d1z)和d2(d2x,d2y,d2z)。地球同步衛(wèi)星Xi 與空中飛行器P的位置關(guān)系示意圖如圖所示，i，i分別表示空中飛行器P的測向陣列方向d1，d2與地球同步衛(wèi)星Xi(i=1,2,N)的夾角?，F(xiàn)在請你們建立數(shù)學(xué)模型研究解決下面的問題：（1）通過測量空中飛行器測向陣列方向d1和d2與多顆地球同步衛(wèi)星的夾角i和i，建立空中飛行器定位的數(shù)學(xué)模型；對于附表1所給出的9顆同步衛(wèi)星的數(shù)據(jù)，試確定空中飛行器P的位置參數(shù)。（2）在某些特殊情況下，空中飛行器能直接檢測

6、到的同步衛(wèi)星數(shù)量較少，可以利用空中飛行器在勻速飛行過程中多次檢測的結(jié)果來實(shí)現(xiàn)定位。針對這種情況，試建立空中飛行器定位的數(shù)學(xué)模型；對附表2中給出的3顆同步衛(wèi)星的檢測數(shù)據(jù)，確定空中飛行器P在第70秒時的位置參數(shù)，并分析其可靠性。（3）當(dāng)可用同步衛(wèi)星數(shù)量較多時，為了提高定位精度和定位效率，需要對可用的同步衛(wèi)星進(jìn)行一定的優(yōu)選。試研究具體的優(yōu)選策略，并通過仿真，分析在檢測方向角誤差限為0.1°時空中飛行器的定位方法和精度。2 模型的假設(shè)（1）假設(shè)附表中所有的數(shù)據(jù)都是真實(shí)可靠的；（2）假設(shè)飛行器是一個質(zhì)點(diǎn)，不考慮飛機(jī)的飛行姿態(tài)；（3）衛(wèi)星信號強(qiáng)度都足夠的強(qiáng)；（4）不考慮地球的實(shí)際曲率變化，認(rèn)為地

7、球是一個均勻球體。（5）衛(wèi)星環(huán)繞運(yùn)動以地心為中心的圓；（6）固定在飛行器上的兩個垂直測向陣列為單位向量；3 符號說明符號含義Xi(xi,yi,zi)同步衛(wèi)星的位置P(x,y,z)飛行器的位置PXi飛行器指向同步衛(wèi)星方向的矢量PXi與測向陣列d1(d1x,d1y,d1z)的夾角PXi與測向陣列d2(d2x,d2y,d2z)的夾角v(vx,vy,vz)飛行器的速度矢量j每個點(diǎn)到平均位置的距離與實(shí)際測量值之間的誤差與實(shí)際測量值之間的誤差d1d2測向陣列相互垂直關(guān)系確定的方程的解的誤差d1測向陣列模為1確定的方程的解的誤差d2測向陣列模為1確定的方程的解的誤差4 問題一模型的建立與求解4.1問題一的分

8、析同步衛(wèi)星處于地球赤道平面上，圍繞地球自轉(zhuǎn)軸旋轉(zhuǎn)，并且相對地球上一點(diǎn)靜止，因此可建立以地心為原點(diǎn)的三維坐標(biāo)系，衛(wèi)星位置參數(shù)可確定。分析測向陣列方向和地球同步衛(wèi)星夾角與飛行器位置的關(guān)系，2個測向陣列以及飛行器位置參數(shù)屬于未知參數(shù)，即至少需要9個包含上述參數(shù)的方程才能確定一個飛行器位置參數(shù)以及2個測向陣列。由于方向角以及測向陣列存在誤差，故可求得多組未知參數(shù)解，因此，精度難以滿足實(shí)際需求。因此，利用方向角測量、計算誤差，以及飛行器高度為約束條件，得到使得計算誤差平方和最小的目標(biāo)函數(shù)，將9顆衛(wèi)星同時代入關(guān)系式，利用最小二乘法，當(dāng)目標(biāo)函數(shù)最小時，可得到最逼近真實(shí)參數(shù)的參數(shù)解。在上述方法考慮了所有飛行器

9、參數(shù)解，即所有解都為有效解，并對最終結(jié)果產(chǎn)生影響。因此可進(jìn)一步優(yōu)化，可考慮按夾角余弦法聚類，將相近的參數(shù)解聚成一類，將個體或者數(shù)量明顯較少的類別剔除，剩下類別的參數(shù)解作為有效解，并以坐標(biāo)均值求解最終參數(shù)解。4.2基于最小二乘法的模型建立建立以地球中心為原點(diǎn)，以赤道平面為基本平面，X軸在基本平面內(nèi)由地心向外指向經(jīng)線與赤道的交點(diǎn)，Z軸過地心且垂直基本平面，Y軸與X，Z軸組成右手系的三維直角坐標(biāo)系(見圖 1 )：圖 1 球心坐標(biāo)系在此坐標(biāo)系下，同步衛(wèi)星的位置是，飛行器的位置是，如圖 2 所示：圖 2 飛行器與衛(wèi)星位置關(guān)系圖建立同步衛(wèi)星位置參數(shù)模型： （地球同步衛(wèi)星所處經(jīng)度角）將附表一的數(shù)據(jù)代入求得同

10、步衛(wèi)星在坐標(biāo)系中的位置參數(shù)如下：表 1 同步衛(wèi)星的位置參數(shù)衛(wèi)星編號xyz11020140914027354216003-144213962404-241853454105-271043230106-303322929207-332282596008-403241232809-4175658680測向陣列方向矢量為和，矢量PXi和d1、d2的夾角分別為i、i，則利用矢量點(diǎn)乘關(guān)系利用一顆衛(wèi)星可以得到如下非線性方程組：由于兩個測向陣列本身相互垂直且模為1，但考慮到測量、計算過程中存在誤差，因此可得式中，、分別為計算誤差。由上面的分析可知，方向角測量、計算等過程中誤差，故可得、為計算誤差,聯(lián)立9顆衛(wèi)星

11、和測向陣列關(guān)系，建立21個方程組，使之滿足如下函數(shù)優(yōu)化問題：式中R為地球半徑，表示飛行器到地球質(zhì)心的距離大于地球半徑，利用最小二乘法就可以求解出飛行器位置，測向陣列和這9個未知參數(shù)，此時，所得解處于最逼近滿足21個方程組的情況。4.3基于聚類分析的模型優(yōu)化 由于存在9個未知參數(shù)，因此需要9個方程組才能解出，觀察關(guān)系式可知，3顆衛(wèi)星反饋的的信號，利用最小二乘法，即能確定飛行器的最佳位置參數(shù)。因?yàn)閷?shí)際過程中存在測向陣列和地球同步衛(wèi)星夾角的測量誤差，故極大可能存在每次選取3顆衛(wèi)星所確定的飛行器位置參數(shù)是不同的，根據(jù)9顆衛(wèi)星可得個飛行器位置參數(shù)，對84個參數(shù)解按夾角余弦法聚類，將數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，因

12、此，可以得到標(biāo)準(zhǔn)化的飛行器位置參數(shù)（，），故其與地心遠(yuǎn)點(diǎn)構(gòu)成的向量，設(shè)論域，的觀測值為，即數(shù)據(jù)矩陣?yán)脢A角余弦法，可得相似系數(shù)：根據(jù)構(gòu)造模糊相似矩陣，在矩陣為基礎(chǔ)，用平方法求出的傳遞包閉，模糊等價矩陣即等于。然后，由大到小取一組，確定相應(yīng)的截矩陣，可以進(jìn)行分類。根據(jù)分類結(jié)果，剔除個體解或數(shù)目較少的類，將其視作測量、計算誤差較大引起偏離真實(shí)位置的情況，利用所剩下的參數(shù)解，求得其坐標(biāo)均值，該坐標(biāo)值與真實(shí)位置參數(shù)基本吻合，其中為剩余參數(shù)解個數(shù)。然后代入利用測向陣列和地球同步衛(wèi)星夾角得出的矢量關(guān)系式，求解出測向陣列和。4.4模型的求解利用基于最小二乘法建立的模型，將9顆衛(wèi)星數(shù)據(jù)同時帶進(jìn)方程組，利用MA

13、TLAB求解出飛行器的位置為 ,測向陣列和。4.5模型的結(jié)果評價基于最小二乘法的求解模型，求解出比較接近真實(shí)參數(shù)的數(shù)據(jù)，充分利用所有衛(wèi)星提供的參數(shù)，全面考慮測量、計算存在的誤差，方法簡明易懂，結(jié)果較為合理；再進(jìn)一步優(yōu)化的以夾角余弦法為基礎(chǔ)的聚類分析模型，篩選出明顯誤差的數(shù)據(jù)并剔除，在此基礎(chǔ)上求解，使結(jié)果更加合理，符合實(shí)際。5 問題二模型的建立與求解5.1問題二的分析由問題一分析可知，三個衛(wèi)星就可以確定飛行器的位置，但精度存在問題。所以當(dāng)同步衛(wèi)星數(shù)量較少時，可利用飛行器在勻速飛行過程中多次檢測的結(jié)果來提高定位的精度。利用問題二所給出的五組數(shù)據(jù)，可以分別確定五個時刻的位置坐標(biāo)，然后分別對五個時刻位

14、置坐標(biāo)分別進(jìn)行曲線擬合，得到飛行器的運(yùn)動軌跡曲線方程，則可得衛(wèi)星在時的位置參數(shù)。再利用此運(yùn)動軌跡對前五組數(shù)據(jù)進(jìn)行檢驗(yàn)，即可得到一組殘差值，將此值帶入灰色預(yù)測模型，得到其殘差預(yù)測值與擬合值相加即為所要預(yù)測飛行器更為精確的位置參數(shù)。通過計算前五個殘差值的方差與飛行器在相應(yīng)階段內(nèi)飛行距離的比值，即可算的飛行器位置參數(shù)的可靠性。5.2.1利用最小二乘法進(jìn)行擬合最小二乘法曲線擬合原理：給定平面上的點(diǎn)，求，使在處函數(shù)值與實(shí)驗(yàn)數(shù)據(jù)的偏差的平方和最小，即吻合度最高。Step1：利用問題二不同時刻檢測到的地球同步衛(wèi)星的相關(guān)數(shù)據(jù)，運(yùn)用問題一所建立的數(shù)學(xué)模型，分別解得時刻空間飛行器的位置坐標(biāo)（見表 2 ）：表 2

15、不同時刻飛行器的坐標(biāo)時刻x坐標(biāo)y坐標(biāo)z坐標(biāo)-4185.024606.536070.90-4154.034621.586083.15-4110.794667.036102.62-4088.174687.826126.47-4036.974711.756145.26Step2：通過以上表格數(shù)據(jù)，以時間為橫坐標(biāo)，分別以軸、軸、軸數(shù)據(jù)為縱坐標(biāo)，建立平面坐標(biāo)軸，運(yùn)用MATLAB軟件分別擬合出值、值、值與時間的函數(shù)關(guān)系，如下圖所示：圖 3 的坐標(biāo)和時間的擬合曲線圖其中的坐標(biāo)值擬合方程：圖 4 的坐標(biāo)和時間的擬合曲線圖其中的坐標(biāo)值擬合方程：圖 5 的坐標(biāo)和時間的擬合曲線圖其中的坐標(biāo)值擬合方程：Step3：皮

16、爾遜相關(guān)系數(shù)是趨勢線擬合程度的指標(biāo)，它的數(shù)值大小可以反映趨勢線的估計值與對應(yīng)的實(shí)際數(shù)據(jù)之間的擬合程度，擬合程度越高，趨勢線的可靠性就越高。經(jīng)計算：坐標(biāo)的皮爾遜相關(guān)系數(shù)為0.9848，坐標(biāo)的皮爾遜相關(guān)系數(shù)為0.9848，坐標(biāo)的皮爾遜相關(guān)系數(shù)為0.9848，各坐標(biāo)的皮爾遜相關(guān)系數(shù)都接近1，說明曲線擬合吻合度較高，可靠性強(qiáng)。Step4：根據(jù)以上函數(shù)關(guān)系，當(dāng)時，空中飛行器的位置坐標(biāo)為(-3934.03、4797.28、6201.73)5.2.2建立模型對距離進(jìn)行灰色預(yù)測灰色預(yù)測模型是通過少量的、不完全的信息，建立數(shù)學(xué)模型并做出預(yù)測的一種預(yù)測方法。對于灰色量的處理不是尋求它的統(tǒng)計規(guī)律和概率分布，而是將無

17、規(guī)律的原始數(shù)據(jù)，通過一定的方法處理，變成比較有規(guī)律的時間序列數(shù)據(jù)。即以數(shù)找數(shù)的規(guī)律，再建立動態(tài)模型。因?yàn)榭陀^系統(tǒng)無論多么復(fù)雜，它總是聯(lián)系的，有序的，有整體功能的。所以作為系統(tǒng)行為特征的數(shù)據(jù)，總是蘊(yùn)含著某種規(guī)律?；疑蓴?shù)列：對灰數(shù)的處理主要是利用數(shù)據(jù)處理方法去尋求數(shù)據(jù)間的內(nèi)在規(guī)律，通過對已知數(shù)據(jù)列中的數(shù)據(jù)進(jìn)行處理而產(chǎn)生新的數(shù)據(jù)列，以此來研究尋找數(shù)據(jù)的規(guī)律性，這種方法稱為數(shù)據(jù)的生成。建立模型：設(shè)原始非負(fù)時間序列為，為累加生成序列，即：的白化微分方程為：為發(fā)展灰數(shù)；為內(nèi)生控制灰數(shù)，亦稱為灰作用量，設(shè)待辨識向量，將上式寫成矩陣的形式：按最小二乘法得式中 于是可得到灰色預(yù)測的 離散時間相應(yīng)函數(shù)為：為所

18、得累加的預(yù)測值，將預(yù)測值還原即為：殘差序列相對誤差序列相應(yīng)的預(yù)測模型模擬根據(jù)灰色預(yù)測模型得到時距離殘差的預(yù)測值如表 3 所示：表 3 距離殘差預(yù)測值（單位：）時間（s）距離殘差03.711022.432025.953029.934034.705040.106046.457053.75由此可知，時，空中飛行器的可信位置范圍為以(-3934.03、4797.28、6201.73)為圓心，以距離殘差預(yù)測值為半徑的球體內(nèi)，位置可靠度即為。表 4 灰色預(yù)測模型的檢驗(yàn)：時間實(shí)際數(shù)據(jù)模擬數(shù)據(jù)殘差相對誤差22.6822.400.26 1.12%21.1420.640.482.27%37.8237.070.78

19、2.55%31.4431.750.321.07% 平均相對誤差：說明灰色預(yù)測模型相對誤差較小，可靠性強(qiáng)，對于空中飛行器有良好的定位性。5.3問題二的結(jié)果分析問題二，利用5個時刻的測量數(shù)據(jù)，運(yùn)用問題一中數(shù)學(xué)模型，解得5個時刻空間飛行器的空間位置坐標(biāo)，在利用最小二乘法進(jìn)行曲線擬合，得到其函數(shù)關(guān)系，進(jìn)而得到所求時刻空間飛行器的位置坐標(biāo)，最后再利用灰色模型求得距離誤差值，并將其作為可靠性指標(biāo)，說明灰色模型可靠性強(qiáng)，對飛行器有較好的定位性。6 問題三模型的建立與求解6.1問題三的分析 當(dāng)可用同步衛(wèi)星數(shù)量較多時，為了提高定位精度和定位效率，需要對可用的同步衛(wèi)星進(jìn)行一定的優(yōu)選，因此需要建立具體的優(yōu)選策略，也

20、就是如何選，選幾顆，選哪幾顆。而選擇幾何分布好的衛(wèi)星組合可以提高定位精度和定位效率。通常我們用幾何精度因子來表征用戶和可見衛(wèi)星在空間幾何分布的好壞。從顆衛(wèi)星中選擇顆衛(wèi)星，共有種選擇方案。從可見衛(wèi)星中選擇參加導(dǎo)航定位計算的衛(wèi)星數(shù)目不同, 的取值也不相同, 與衛(wèi)星數(shù)目之間有一定的變化規(guī)律，分析衛(wèi)星數(shù)目和之間的關(guān)系，從而找到一個合適的衛(wèi)星數(shù)目。同時，衛(wèi)星數(shù)目越多，意味著方程數(shù)量越多，計算量越大，定位效率也就越低。因此，為了兼顧定位的效率，應(yīng)該在保證定位精度的前提下，盡可能減小定位衛(wèi)星的數(shù)量。根據(jù)對定位精度和定位效率的要求，確定一個合適的衛(wèi)星數(shù)目。確定衛(wèi)星數(shù)目為之后，還需對衛(wèi)星的幾何布局進(jìn)行優(yōu)選，遍歷

21、從所有衛(wèi)星中選擇顆衛(wèi)星的所有方案，從中選擇最小的衛(wèi)星組合，即為優(yōu)選的結(jié)果。確定定位衛(wèi)星后，給測量角度加上最大為0.1°的高斯噪聲，進(jìn)行多次的仿真實(shí)驗(yàn)，統(tǒng)計定位精度。6.2基于幾何精度（）因子模型幾何精度GDOP：是衡量定位精度的很重要的一個系數(shù)，它代表GPS 測距誤差造成的接收機(jī)與空間衛(wèi)星間的距離矢量放大因子。實(shí)際表征參與定位解的從接收機(jī)至空間衛(wèi)星的單位矢量所勾勒的形體體積與成反比，故又稱為幾何精度因子。實(shí)際上，GDOP的數(shù)值越大，所代表的單位矢量形體體積越小，即接收機(jī)至空間衛(wèi)星的角度十分相似導(dǎo)致的結(jié)果，此時的GDOP會導(dǎo)致定位精度變差。好的GDOP, 是指其數(shù)值小，代表大的單位矢量

22、形體體積，導(dǎo)致高的定位精度。好的幾何因子實(shí)際上是指衛(wèi)星在空間分布不集中于一個區(qū)域，同時能在不同方位區(qū)域均勻分布。假設(shè)地球同步衛(wèi)星的數(shù)目為4個，定義下列向量和矩陣：，可將上式進(jìn)一步簡化，得到：由于陣階數(shù)為4*3，故在上式兩側(cè)左乘的轉(zhuǎn)置矩陣就能獲得最小二乘解。若將用分量形式表示，即：定義幾何精度因子如下：圖 6 GDOP計算流程圖6.3優(yōu)選方法的確定由于在最優(yōu)定位衛(wèi)星的選擇中，問題可抽象為一個組合問題，即已知顆地球同步衛(wèi)星，從中選出顆組合用于定位計算，以作為目標(biāo)函數(shù)，進(jìn)而計算出每一種組合對應(yīng)的值，即相同數(shù)量的衛(wèi)星條件下，有種組合。考慮衛(wèi)星空間分布對精度的影響和方便計算的要求，在地球同步衛(wèi)星數(shù)量相同

23、的條件下，選取兩組分布差別較大的組合，運(yùn)用問題（1）的方法聯(lián)立方程組即可求得飛行器的位置參數(shù)，并根據(jù)所構(gòu)建的模型。計算出不同數(shù)量下的衛(wèi)星不同組合得到的值，(見表 5 )：表 5 不同衛(wèi)星組合位置參數(shù)對照表衛(wèi)星組合位置坐標(biāo)(單位:km)（-5210.62，6608.68，3217.94）7.52(-5203.93， 6605.61， 3149.20)8.16(-5198.70， 6604.12，3122.95)4.69(-5202.86， 6605.96， 3166.78)4.13(-5198.26， 6605.17， 3124.96)0.76(-5201.33， 6605.27， 3157.3

24、6)0.72(-5200.26， 6605.29， 3139.00)0.36(-5202.51，6605.51，3147.92)0.31(-5202.40， 6605.72，3150.56)0.22(-5198.50， 6604.57，3121.25)0.19(-5188.74， 6605.39，3128.43)0.11根據(jù)上表數(shù)據(jù)計算出同一數(shù)量衛(wèi)星不同組合的平均值，見表6： 表 6 不同數(shù)量衛(wèi)星組合均值表衛(wèi)星數(shù)4567897.844.410.740.330.200.11圖 7 衛(wèi)星數(shù)量對GDOP的影響由上圖可知：值隨地球同步衛(wèi)星數(shù)目的增加而顯著減小，地球同步衛(wèi)星數(shù)目在4-6時，值顯著降低，下

25、降速率較大，而當(dāng)定位衛(wèi)星的數(shù)量達(dá)到6顆以后，即衛(wèi)星數(shù)量為6-9顆時，其值變化不大，即此時增加衛(wèi)星數(shù)量，對的改善已經(jīng)不明顯。因此綜合考慮定位精度和定位效率，6顆同步衛(wèi)星較為適宜。 從9顆衛(wèi)星選擇6顆衛(wèi)星，共有種選擇方案，根據(jù)問題一所構(gòu)建的數(shù)學(xué)模型計算出每一種方案對應(yīng)的值，我們選擇其中對應(yīng)最小的一種選擇方案。經(jīng)計算，當(dāng)選擇第1,2,3,5,7,9顆衛(wèi)星時，其值最小。經(jīng)度分別為，該組合對應(yīng)的為。6.4角度誤差影響的仿真正態(tài)分布：是一個在數(shù)學(xué)、物理及工程等領(lǐng)域都非常重要的概率分布。若隨機(jī)變量服從一格數(shù)學(xué)期望、方差為的高斯分布，記為。其期望值決定了其位置，其標(biāo)準(zhǔn)差決定了分布的幅度。在多次測量中，設(shè)備的測

26、量誤差是正態(tài)分布的，而且均值為。問題中給出的精度可以認(rèn)為是測量誤差的標(biāo)準(zhǔn)差。在這種理解下，用各自的誤差限對測量誤差進(jìn)行無量綱化(也可以看成是一種加權(quán)法)處理是合理的，即求解的無約束優(yōu)化問題更合理。由6.3的分析可知，當(dāng)六顆衛(wèi)星定位時，定位精度和效率最優(yōu)，故構(gòu)建數(shù)學(xué)模型如下：運(yùn)用計算機(jī)仿真技術(shù)進(jìn)行500次仿真實(shí)驗(yàn)得到實(shí)驗(yàn)定位位置與未加入誤差時定位位置的距離誤差，統(tǒng)計這些距離誤差，得到其分布圖(見 圖 8 )：圖 8 距離誤差分布圖由上圖可以看出，絕大部分距離誤差都分布在250km以內(nèi),因此，認(rèn)定定位精度為250km。6.4問題三的結(jié)果分析問題三，利用幾何精度因子來衡量衛(wèi)星幾何分布的好壞，遍歷所有

27、方案后，發(fā)現(xiàn)6星組合定位的平均精度和效率較好，其定位結(jié)果與9星定位的結(jié)果很接近，說明6星定位的組合是較為合理的。在給測量值加入最大為0.1°的高斯分布后，大量的仿真實(shí)驗(yàn)顯示，加入的誤差給定位帶來了250km的偏差。說明模型建立的定位方法對于角度測量有著很高的要求，否則就有可能造成較大的誤差，即說明本模型對測量誤差有很高的要求，測量誤差是影響定位準(zhǔn)確與否的重要因素。7 模型的優(yōu)缺點(diǎn)分析與改進(jìn)方向7.1優(yōu)點(diǎn)：（1） 問題一基于最小二乘法的求解模型，求解出比較接近真實(shí)參數(shù)的數(shù)據(jù)，充分利用所有衛(wèi)星提供的參數(shù)，全面考慮測量、計算存在的誤差，方法簡明易懂，結(jié)果較為合理；再進(jìn)一步優(yōu)化的以夾角余弦法

28、為基礎(chǔ)的聚類分析模型，篩選出明顯誤差的數(shù)據(jù)并剔除，在此基礎(chǔ)上求解，使結(jié)果更加合理，符合實(shí)際。（2） 問題二在問題一的基礎(chǔ)上利用同一時刻給出的三顆衛(wèi)星的參數(shù)算出各時刻對應(yīng)的飛行器位置參數(shù)，然后利用最小二乘法擬合成曲線，推算出運(yùn)動模型。并用模型預(yù)測擬合位置與實(shí)際位置距離，則其實(shí)際位置等于擬合值跟預(yù)測值的和，將誤差的變化考慮進(jìn)模型，提高了定位的精度。（3） 在解決問題三時，考慮衛(wèi)星優(yōu)選策略時使用了精度幾何因子GDOP的方法，兼顧定位的精度和效率，選擇了合適的衛(wèi)星個數(shù)和衛(wèi)星分布，這種方法精度比較高。當(dāng)測量角度發(fā)生偏差時，假設(shè)偏差符合正態(tài)分布，通過500次的仿真，得出距離誤差分布圖，更加形象直觀。 7.

29、2缺點(diǎn)：（1） 在問題一求解過程中沒有考慮到那些個別會引起測量偏差較大的衛(wèi)星對結(jié)果的影響，即沒有對衛(wèi)星進(jìn)行優(yōu)選。（2） 在問題二中僅僅利用三顆衛(wèi)星提供的參數(shù)得出飛行器在不同時刻的位置參數(shù)本身就存在較大誤差，求解過程中飛行器的位置參數(shù)在軸上的數(shù)據(jù)殘差過大。（3） 在問題三中，本文僅僅只考慮了衛(wèi)星個數(shù)帶來的運(yùn)算量對定位效率的影響 ，沒有沒有考慮其他因素的影響，這是不符合實(shí)際情況的。7.3改進(jìn)方向：問題二求解過程中，將飛行器各軸上的位置參數(shù)分別進(jìn)行數(shù)據(jù)擬合，最后得到時飛行器的位置在一個球體范圍內(nèi)，難以定位，可綜合各軸位置參數(shù)進(jìn)行擬合得出飛行器的空間軌跡隨時間得變化規(guī)律，可得預(yù)測到更為精確直觀的飛行器

30、位置參數(shù)。問題三中，考慮到選取衛(wèi)星組合的人為主觀因素，可盡量多的選取衛(wèi)星組合對飛行器進(jìn)行定位，然后根據(jù)需要賦予不同權(quán)重，從而得到最優(yōu)解。參考文獻(xiàn)1 孫仲康,郭福成,馮道旺.單站無源定位跟蹤技術(shù). 國防大學(xué)出版社. 2008.11第1版2 程銘東，劉利姣，黃光明.基于遺傳算法的多傳感器網(wǎng)絡(luò)中目標(biāo)定位算法J.中國科技論文在線3 陳玲，無源定位與跟蹤算法研究D. 北京：北京航空航天大學(xué)，2005。4 陳小平，騰云龍，康榮雷等.幾何精度因子改進(jìn)算法研究J.電子科技大學(xué)學(xué)報，2008，(6)5 卓金武.Matlab在數(shù)學(xué)建模中的應(yīng)用. 北京航天航空大學(xué)出版社. 2011.4第1版6 鄧聚龍.灰色預(yù)測與決

31、策,華中工學(xué)院出版社，1986.附錄：附錄1.第一問的求解的主函數(shù)function F=f1(x)R=6367;h=35800;theta1=142/360*2*pi;theta2=163.0/360*2*pi;theta3=172.0/360*2*pi;theta4=136.0/360*2*pi;theta5=130.0/360*2*pi;theta6=125.0/360*2*pi;theta7=110.0/360*2*pi;theta8=89.0/360*2*pi;theta9=76.0/360*2*pi; a1=100.32/360*2*pi;b1=11.58/360*2*pi;a2=1

32、25.58/360*2*pi;b2=36.04/360*2*pi;a3=135.93/360*2*pi;b3=46.34/360*2*pi;a4=92.91/360*2*pi;b4=5.99/360*2*pi;a5=85.46/360*2*pi;b5=6.95/360*2*pi;a6=79.24/360*2*pi;b6=11.99/360*2*pi;a7=60.73/360*2*pi;b7=29.81/360*2*pi;a8=35.82/360*2*pi;b8=54.63/360*2*pi;a9=21.33/360*2*pi;b9=69.25/360*2*pi; xm1=(R+h)*cos(t

33、heta1);ym1=(R+h)*sin(theta1);zm1=0;xm2=(R+h)*cos(theta2);ym2=(R+h)*sin(theta2);zm2=0;xm3=(R+h)*cos(theta3);ym3=(R+h)*sin(theta3);zm3=0;xm4=(R+h)*cos(theta4);ym4=(R+h)*sin(theta4);zm4=0;xm5=(R+h)*cos(theta5);ym5=(R+h)*sin(theta5);zm5=0;xm6=(R+h)*cos(theta6);ym6=(R+h)*sin(theta6);zm6=0;xm7=(R+h)*cos(t

34、heta7);ym7=(R+h)*sin(theta7);zm7=0;xm8=(R+h)*cos(theta8);ym8=(R+h)*sin(theta8);zm8=0;xm9=(R+h)*cos(theta9);ym9=(R+h)*sin(theta9);zm9=0; if(x(1)2+x(2)2+x(3)2<63672) F=inf;else F=(xm1-x(1)*x(4)+(ym1-x(2)*x(5)+(-x(3)*x(6)-(xm1-x(1)2+(ym1-x(2)2+(zm1-x(3)2)0.5*cos(a1)2+. (xm1-x(1)*x(7)+(ym1-x(2)*x(8)+

35、(-x(3)*x(9)-(xm1-x(1)2+(ym1-x(2)2+(zm1-x(3)2)0.5*cos(b1)2+. (xm2-x(1)*x(4)+(ym2-x(2)*x(5)+(-x(3)*x(6)-(xm2-x(1)2+(ym2-x(2)2+(zm2-x(3)2)0.5*cos(a2)2+. (xm2-x(1)*x(7)+(ym2-x(2)*x(8)+(-x(3)*x(9)-(xm2-x(1)2+(ym2-x(2)2+(zm2-x(3)2)0.5*cos(b2)2+. (xm3-x(1)*x(4)+(ym3-x(2)*x(5)+(-x(3)*x(6)-(xm3-x(1)2+(ym3-x(

36、2)2+(zm3-x(3)2)0.5*cos(a3)2+. (xm3-x(1)*x(7)+(ym3-x(2)*x(8)+(-x(3)*x(9)-(xm3-x(1)2+(ym3-x(2)2+(zm3-x(3)2)0.5*cos(b3)2+. (xm4-x(1)*x(4)+(ym4-x(2)*x(5)+(-x(3)*x(6)-(xm4-x(1)2+(ym4-x(2)2+(zm4-x(3)2)0.5*cos(a4)2+. (xm4-x(1)*x(7)+(ym4-x(2)*x(8)+(-x(3)*x(9)-(xm4-x(1)2+(ym4-x(2)2+(zm4-x(3)2)0.5*cos(b4)2+.

37、(xm5-x(1)*x(4)+(ym5-x(2)*x(5)+(-x(3)*x(6)-(xm5-x(1)2+(ym5-x(2)2+(zm5-x(3)2)0.5*cos(a5)2+. (xm5-x(1)*x(7)+(ym5-x(2)*x(8)+(-x(3)*x(9)-(xm5-x(1)2+(ym5-x(2)2+(zm5-x(3)2)0.5*cos(b5)2+. (xm6-x(1)*x(4)+(ym6-x(2)*x(5)+(-x(3)*x(6)-(xm6-x(1)2+(ym6-x(2)2+(zm6-x(3)2)0.5*cos(a6)2+. (xm6-x(1)*x(7)+(ym6-x(2)*x(8)+

38、(-x(3)*x(9)-(xm6-x(1)2+(ym6-x(2)2+(zm6-x(3)2)0.5*cos(b6)2+. (xm7-x(1)*x(4)+(ym7-x(2)*x(5)+(-x(3)*x(6)-(xm7-x(1)2+(ym7-x(2)2+(zm7-x(3)2)0.5*cos(a7)2+. (xm7-x(1)*x(7)+(ym7-x(2)*x(8)+(-x(3)*x(9)-(xm7-x(1)2+(ym7-x(2)2+(zm7-x(3)2)0.5*cos(b7)2+. (xm8-x(1)*x(4)+(ym8-x(2)*x(5)+(-x(3)*x(6)-(xm8-x(1)2+(ym8-x(

39、2)2+(zm8-x(3)2)0.5*cos(a8)2+. (xm8-x(1)*x(7)+(ym8-x(2)*x(8)+(-x(3)*x(9)-(xm8-x(1)2+(ym8-x(2)2+(zm8-x(3)2)0.5*cos(b8)2+. (xm9-x(1)*x(4)+(ym9-x(2)*x(5)+(-x(3)*x(6)-(xm9-x(1)2+(ym9-x(2)2+(zm9-x(3)2)0.5*cos(a9)2+. (xm9-x(1)*x(7)+(ym9-x(2)*x(8)+(-x(3)*x(9)-(xm9-x(1)2+(ym9-x(2)2+(zm9-x(3)2)0.5*cos(b9)2+.

40、(x(4)*x(7)+x(5)*x(8)+x(6)*x(9)2+. (x(4)2+x(5)2+x(6)2-1)2+. (x(7)2+x(8)2+x(9)2-1)2;endend附錄2.最小二乘法求解過程clc;clear; xx=-6000;-3000;-1000;yy=2000;4000;7000;zz=1000;3000;5000;vx=-30;0;30;vy=-30;0;30;vz=-30;0;30;xxx=;x0=-6000;2000;1000;1;0;1;2;1;6;30;30;30;mm=10;nn=x0;for i=1:3 for j=1:3 for k=1:3 for ii=1:3 for jj=1:3 for kk=1:3 aa1=xx(i);bb1=yy(j);cc1=zz(k);vx1=vx(ii);vy1=vy(jj);vz1=vz(kk); x0=aa1;bb1;cc1;1;0;1;2;1;6;vx1;vy1;vz1; options = optimset(optimset,'MaxFunEvals',100000); x,resnorm=lsqnonlin(f1,x0,options); xxx=xxx;x; if resnorm<mm nn=x; mm=resnorm; end end end en