捷聯(lián)慣導(dǎo)算法與組合導(dǎo)航原理講義

-CAL-FENGHAI.-(YICAI)-CompanyOne1

-CAL-本頁僅作為文檔封面，使用請直接刪除-CAL-FENGHAI.-(YICAI)-CompanyOne1

-CAL-本頁僅作為文檔封面，使用請直接刪除捷聯(lián)慣導(dǎo)算法與組合導(dǎo)航原理講義(20170220)(總205頁)捷聯(lián)慣導(dǎo)算法與組合導(dǎo)航原理講義嚴(yán)恭敏，翁浚編著西北工業(yè)大學(xué)2016-9

前言近年來，慣性技術(shù)不論在軍事上、工業(yè)上，還是在民用上，特別是消費電子產(chǎn)品領(lǐng)域，都獲得了廣泛的應(yīng)用，大到潛艇、艦船、高鐵、客機、導(dǎo)彈和人造衛(wèi)星，小到醫(yī)療器械、電動獨輪車、小型四旋翼無人機、空中鼠標(biāo)和手機，都有慣性技術(shù)存在甚至大顯身手的身影。相應(yīng)地，慣性技術(shù)的研究和開發(fā)也獲得前所未有的蓬勃發(fā)展，越來越多的高校學(xué)生、愛好者和工程技術(shù)人員加入到慣性技術(shù)的研發(fā)隊伍中來。慣性技術(shù)涉及面廣，涵蓋元器件技術(shù)、測試設(shè)備和測試方法、系統(tǒng)集成技術(shù)和應(yīng)用開發(fā)技術(shù)等方面，囿于篇幅和作者知識面限制，本書主要討論捷聯(lián)慣導(dǎo)系統(tǒng)算法方面的有關(guān)問題，包括姿態(tài)算法基本理論、捷聯(lián)慣導(dǎo)更新算法與誤差分析、組合導(dǎo)航卡爾曼濾波原理、捷聯(lián)慣導(dǎo)系統(tǒng)的初始對準(zhǔn)技術(shù)、組合導(dǎo)航系統(tǒng)建模以及算法仿真等內(nèi)容。希望讀者參閱之后能夠?qū)萋?lián)慣導(dǎo)算法有個系統(tǒng)而深入的理解，并能快速而有效地將基本算法應(yīng)用于解決實際問題。本書在編寫和定稿過程中得到以下同行的熱心支持，指出了不少錯誤之處或提出了許多寶貴的修改建議，深表謝意：西北工業(yè)大學(xué)自動化學(xué)院：梅春波、趙彥明、劉洋、沈彥超、肖迅、牟夏、鄭江濤、劉士明、金竹、馮理成、趙雪華；航天科工第九總體設(shè)計部：王亞軍；遼寧工程技術(shù)大學(xué)：丁偉；北京騰盛科技有限公司：劉興華；東南大學(xué)：童金武；中國農(nóng)業(yè)大學(xué)：包建華；南京航空航天大學(xué)：趙宣懿；武漢大學(xué)：董翠軍；網(wǎng)友：Zoro；山東科技大學(xué)：王云鵬。書中缺點和錯誤在所難免，望讀者不吝批評指正。作者2016年9月

目錄第1章概述 61.1捷聯(lián)慣導(dǎo)算法簡介 61.2Kalman濾波與組合導(dǎo)航原理簡介 7第2章捷聯(lián)慣導(dǎo)姿態(tài)解算基礎(chǔ) 102.1反對稱陣及其矩陣指數(shù)函數(shù) 102.1.1反對稱陣 102.1.2反對稱陣的矩陣指數(shù)函數(shù) PAGEREF122.2方向余弦陣與等效旋轉(zhuǎn)矢量 132.2.1方向余弦陣 132.2.2等效旋轉(zhuǎn)矢量 \h142.3方向余弦陣微分方程及其求解 172.3.1方向余弦陣微分方程 172.3.2方向余弦陣微分方程的求解 172.4姿態(tài)更新的四元數(shù)表示 202.4.1四元數(shù)的基本概念 202.4.2四元數(shù)微分方程 232.4.3四元數(shù)微分方程的求解 252.5等效旋轉(zhuǎn)矢量微分方程及其泰勒級數(shù)解 262.5.1等效旋轉(zhuǎn)矢量微分方程 262.5.2等效旋轉(zhuǎn)矢量微分方程的泰勒級數(shù)解 292.6圓錐運動條件下的等效旋轉(zhuǎn)矢量算法 312.6.1圓錐運動的描述 312.6.2圓錐誤差補償算法 33第3章地球形狀與重力場基礎(chǔ) 403.1地球的形狀描述 403.2地球的正常重力場 463.3地球重力場的球諧函數(shù)模型 503.3.1球諧函數(shù)的基本概念 503.3.2地球引力位函數(shù) 583.3.3重力位及重力計算 63第4章捷聯(lián)慣導(dǎo)更新算法及誤差分析 694.1捷聯(lián)慣導(dǎo)數(shù)值更新算法 694.1.1姿態(tài)更新算法 694.1.2速度更新算法 704.1.3位置更新算法 PA764.2捷聯(lián)慣導(dǎo)誤差方程 764.2.1慣性傳感器測量誤差 764.2.2姿態(tài)誤差方程 784.2.3速度誤差方程 794.2.4位置誤差方程 794.2.5誤差方程的整理 804.3靜基座誤差特性分析 824.3.1靜基座誤差方程 824.3.2高度通道 83_Toc475862305"6.4車載慣性/里程儀組合導(dǎo)航 1596.4.1航位推算算法 1596.4.2航位推算誤差分析 1616.4.3慣性/里程儀組合 1646.5低成本姿態(tài)航向參考系統(tǒng)（AHRS） 1676.5.1簡化的慣導(dǎo)算法及誤差方程 1686.5.2地磁場測量及誤差方程 1696.5.3低成本組合導(dǎo)航系統(tǒng)模型 1706.5.4低成本慣導(dǎo)的姿態(tài)初始化 1716.5.5捷聯(lián)式地平儀的工作原理 173第7章捷聯(lián)慣導(dǎo)與組合導(dǎo)航仿真 1767.1飛行軌跡和慣性器件信息仿真 1767.1.1飛行軌跡設(shè)計 1767.1.2捷聯(lián)慣導(dǎo)反演算法 1777.1.3仿真 1787.2捷聯(lián)慣導(dǎo)仿真 PAGEREF_1807.2.1Matlab子函數(shù) 1807.2.2捷聯(lián)慣導(dǎo)仿真主程序 1857.3慣導(dǎo)/衛(wèi)星組合導(dǎo)航仿真 186子函數(shù) 1867.3.2組合導(dǎo)航仿真主程序 187附錄 190A一些重要的三維矢量運算關(guān)系 190B運載體姿態(tài)的歐拉角描述 192C姿態(tài)更新的畢卡算法、龍格—庫塔算法及精確數(shù)值解法 199D從非直角坐標(biāo)系到直角坐標(biāo)系的矩陣變換 207E線性系統(tǒng)基本理論 211F加權(quán)最小二乘估計 216G矩陣求逆引理 217H幾種矩陣分解方法（QR、Cholesky與UD） 219I二階濾波中的引理證明 223J方差陣上界的證明 225K三階非奇異方陣的奇異值分解 226LMatlab仿真程序 2311-3"\h\z\u第1章概述 11.1捷聯(lián)慣導(dǎo)算法簡介 11.2Kalman濾波與組合導(dǎo)航原理簡介 31.1捷聯(lián)慣導(dǎo)算法簡介在捷聯(lián)慣導(dǎo)系統(tǒng)（SINS）中慣性測量器件（陀螺和加速度計）直接與運載體固聯(lián)，通過導(dǎo)航計算機采集慣性器件的輸出信息并進行數(shù)值積分求解運載體的姿態(tài)、速度和位置等導(dǎo)航參數(shù)，這三組參數(shù)的求解過程即所謂的姿態(tài)更新算法、速度更新算法和位置更新算法。特別在惡劣的高動態(tài)環(huán)境下，高精度的SINS對慣性器件性能和導(dǎo)航算法精度的要求都非?？量?，由于高精度慣性器件往往價格昂貴并且進一步提升精度異常困難，所以在影響SINS精度的所有誤差源中要求因?qū)Ш剿惴ㄒ鸬恼`差比重必須很小，一般認(rèn)為應(yīng)小于5%。姿態(tài)更新算法是SINS算法的核心，對整個系統(tǒng)的解算精度影響最為突出，具有重要的研究和應(yīng)用價值。傳統(tǒng)的姿態(tài)更新算法有歐拉角法、方向余弦陣法和四元數(shù)法等方法，這些方法直接以陀螺采樣輸出作為輸入，使用泰勒級數(shù)展開或龍格—庫塔等方法求解姿態(tài)微分方程，未充分考慮轉(zhuǎn)動的不可交換性誤差問題。傳統(tǒng)姿態(tài)更新算法在理論上可以通過提高采樣和更新頻率來提高解算精度，但實際陀螺采樣頻率又受限于傳感器的帶寬和噪聲水平，因此傳統(tǒng)算法的精度提升空間相對有限，僅適用于對解算精度要求不太高的場合。早在1775年，歐拉就提出了等效旋轉(zhuǎn)矢量的概念，指出剛體的定點轉(zhuǎn)動（即繞固定點的任何有限角位移）均可用繞經(jīng)過該固定點的某軸的一次轉(zhuǎn)動來實現(xiàn)，建立了剛體上單位矢量在轉(zhuǎn)動前后的變換公式。1840年，羅德里格使用后人稱之為羅德里格參數(shù)的表示方法，推導(dǎo)了相繼兩次轉(zhuǎn)動的合成公式，它與W.R.Hamilton在1843年發(fā)明的四元數(shù)乘法表示是一致的。研究表明，相繼多次的定點轉(zhuǎn)動問題可用一系列的姿態(tài)變化量（變化四元數(shù)或變化矩陣）相乘來描述，每個姿態(tài)變化量與對應(yīng)轉(zhuǎn)動的等效旋轉(zhuǎn)矢量之間存在轉(zhuǎn)換公式，使用等效旋轉(zhuǎn)矢量計算姿態(tài)變化量不存在任何原理上的誤差。因此，現(xiàn)代的SINS姿態(tài)更新算法研究的關(guān)鍵就在于如何使用陀螺輸出構(gòu)造等效旋轉(zhuǎn)矢量，以盡量減小和避免不可交換性誤差，后續(xù)再使用等效旋轉(zhuǎn)矢量計算姿態(tài)變化量和進行姿態(tài)更新將變得非常簡單，而不像傳統(tǒng)方法那樣，直接使用陀螺輸出進行姿態(tài)更新容易引起不可交換性誤差。1949年，J.H.Laning在研究火控系統(tǒng)的過程中詳細(xì)地分析了空間轉(zhuǎn)動合成的性質(zhì)，推導(dǎo)了由等效旋轉(zhuǎn)矢量確定轉(zhuǎn)動角速度的公式，但是由于缺少更好的應(yīng)用背景驅(qū)動（比如后來SINS發(fā)展的迫切需求），未能獲得廣泛的研究重視。20世紀(jì)50年代是機械陀螺儀飛速發(fā)展的一個重要時期，也正是在那時發(fā)現(xiàn)了著名的圓錐運動現(xiàn)象，即當(dāng)陀螺儀在其旋轉(zhuǎn)軸和輸出軸出現(xiàn)同頻不同相的角振動時，盡管其輸入軸凈指向不變（在整體上沒有隨時間改變的趨勢項），但陀螺儀還是會敏感到并輸出常值角速率。1958年，為揭示圓錐運動現(xiàn)象產(chǎn)生的根源，L.E.Goodman建立了剛體轉(zhuǎn)動的等效旋轉(zhuǎn)矢量與角速度之間的關(guān)系式，后人稱之為Goodman-Robinson定理，該定理從幾何上將轉(zhuǎn)動不可交換性誤差的坐標(biāo)分量描述為單位球面上的一塊有向面積，其面積由對應(yīng)動坐標(biāo)軸在單位球面上掃過的曲線與連接該曲線端點的大圓圍成，Goodman借助二維Green積分理論獲得了不可交換性誤差的近似公式。1969年，基于Goodman近似公式，J.W.Jordan在假設(shè)陀螺角增量輸出為二次多項式條件下提出了等效旋轉(zhuǎn)矢量的“pre-processor”算法，它與后來發(fā)展的等效旋轉(zhuǎn)矢量二子樣算法完全一致。1969年，J.E.Bortz在其博士論文中詳細(xì)推導(dǎo)了等效旋轉(zhuǎn)矢量微分方程（1971年正式發(fā)表，后人稱之為Bortz方程），它是利用陀螺輸出求解等效旋轉(zhuǎn)矢量的基本公式，奠定了等效旋轉(zhuǎn)矢量多子樣算法的理論基礎(chǔ)。在實際應(yīng)用時一般需對較復(fù)雜的Bortz方程做近似處理，事實上，其簡化結(jié)果與Goodman公式完全一致，它也可以根據(jù)Laning公式簡化獲得。1983年，R.B.Miller采用在圓錐運動條件下使算法漂移誤差最小作為評價標(biāo)準(zhǔn)，推導(dǎo)了等效旋轉(zhuǎn)矢量三子樣優(yōu)化算法。1990年，J.E.Lee研究了四子樣優(yōu)化算法。1992年，Y.F.Jiang研究了利用本更新周期內(nèi)的三子樣及前更新周期內(nèi)的角增量計算旋轉(zhuǎn)矢量的優(yōu)化算法。1996年，M.B.Ignagni提出了由陀螺角增量構(gòu)造等效旋轉(zhuǎn)矢量的通式，并給出了多達(dá)10種類型的等效旋轉(zhuǎn)矢量算法。1999年，C.G.Park總結(jié)提出了各子樣下求解圓錐誤差補償系數(shù)和算法漂移誤差估計的通用公式。至此，從理論上看，在理想的圓錐運動條件下的不可交換性誤差補償問題得到了比較完美的解決。捷聯(lián)慣導(dǎo)的基本概念在20世紀(jì)50年代就已經(jīng)提出了，但是由于當(dāng)時計算機的運算能力極其有限，在算法發(fā)展的早期階段姿態(tài)更新通常采用雙速回路算法方案：高速回路（e.g.,400Hz-10kHz）使用簡單的一階算法補償由角振動引起的姿態(tài)不可交換性誤差；中速回路（e.g.,50Hz-200Hz）以高速回路的處理結(jié)果作為輸入再使用相對復(fù)雜的高階算法進行姿態(tài)矩陣或四元數(shù)更新。雙速回路算法的結(jié)構(gòu)設(shè)計和實現(xiàn)過程都稍顯繁瑣，它只是在計算機運算能力低下時期所采取的權(quán)宜之策，隨著通用計算機技術(shù)的飛速發(fā)展，尤其是80年代中后期之后，導(dǎo)航計算機的運算能力就不再是導(dǎo)航算法研究中需要著重關(guān)注的問題。雙速回路算法的結(jié)構(gòu)研究已經(jīng)成為歷史，目前的計算機完全能夠滿足高速高精度姿態(tài)更新解算的要求。1998年，P.G.Savage相繼發(fā)表的兩篇論文對整體捷聯(lián)慣導(dǎo)數(shù)值算法進行了比較全面的總結(jié)，但相對于普通技術(shù)人員而言，其算法描述過于繁雜，給具體實現(xiàn)帶來了很大的不便或困惑。1.2Kalman濾波與組合導(dǎo)航原理簡介如果信號受噪聲干擾，為了從量測中恢復(fù)出有用信號而又要盡量減少干擾的影響，常常采用濾波器進行信號處理。使用經(jīng)典濾波器時假定信號和干擾的頻率分布不同，通過設(shè)計特定的濾波器帶通和帶止頻段，實現(xiàn)有用信號和干擾的分離。但是，如果干擾的頻段很寬，比如白噪聲，在有用信號的頻段范圍內(nèi)也必然會存在干擾，這時經(jīng)典濾波器對濾除這部分干擾噪聲無能為力。若有用信號和干擾噪聲的頻帶相互重疊，信號處理時通常不再認(rèn)為有用信號是確定性的，而是帶有一定隨機性的。對于隨機信號不可能進行準(zhǔn)確無誤差的恢復(fù)，只能根據(jù)信號和噪聲的統(tǒng)計特性，利用數(shù)理統(tǒng)計方法進行估計，并且一般采取某種統(tǒng)計準(zhǔn)則使估計誤差盡可能小。借用經(jīng)典濾波器的術(shù)語，這種針對隨機信號的統(tǒng)計估計方法也常常稱為濾波器，或稱為現(xiàn)代濾波器以區(qū)別于經(jīng)典濾波器，但須注意經(jīng)典濾波器和現(xiàn)代濾波器之間是有本質(zhì)區(qū)別的。1Kalman濾波早在1632年，伽利略（GalileoGalilei）就嘗試用各種誤差函數(shù)最小化的方法提出了估計理論問題。1801年，數(shù)學(xué)家高斯（KarlGauss）將最小二乘估計法應(yīng)用于谷神星的軌道跟蹤和預(yù)測，取得了良好的效果。最小二乘估計以觀測殘差平方和最小作為估計準(zhǔn)則，它不需要關(guān)于量測的任何統(tǒng)計信息，算法簡單且實用性強，在參數(shù)估計領(lǐng)域獲得了廣泛的應(yīng)用。但是，通常情況下最小二乘估計只能應(yīng)用于靜態(tài)參數(shù)估計，而不適用于動態(tài)系統(tǒng)的狀態(tài)估計。20世紀(jì)40年代初期，維納（NorbertWiener）開始將統(tǒng)計方法應(yīng)用于通信系統(tǒng)和控制系統(tǒng)的研究中，提出了著名的維納濾波理論。同一時期，柯爾莫哥洛夫（AndreyKolmogorow）也進行了類似的研究。維納濾波是一種從頻域角度出發(fā)設(shè)計濾波器的方法，它根據(jù)有用信號和干擾信號的功率譜特性，通過構(gòu)造和求解維納—霍夫（Wiener-Hopf）方程得到最佳濾波器的傳遞函數(shù)，給出了最小均方誤差意義下的穩(wěn)態(tài)解。但是，在一般情況下求解維納—霍夫方程極為困難，甚至是不可能的。此外，維納濾波僅適用于低維平穩(wěn)隨機過程，人們試圖將它推廣到高維和非平穩(wěn)情況，但都因無法突破計算上的困難而難以實用，這嚴(yán)重限制了維納濾波的普及。維納濾波在歷史上有著非常重要的作用和獨特的地位，它首次將數(shù)理統(tǒng)計理論和線性系統(tǒng)理論有機結(jié)合起來，形成了對隨機信號進行估計的新理論，雖然維納濾波不適合用于狀態(tài)估計，但是它在信號處理和通信理論中依然十分有用。1960年，RudolfKalman將控制系統(tǒng)狀態(tài)空間的概念引入隨機估計理論中，建立了隨機狀態(tài)空間模型，利用了隨機狀態(tài)方程、量測方程以及激勵白噪聲的統(tǒng)計特性，構(gòu)造估計算法對隨機狀態(tài)進行濾波估計，后來被稱為Kalman（卡爾曼）濾波。在Kalman濾波中，所有利用的信息都是時域內(nèi)的參量，它不但可以應(yīng)用于一維平穩(wěn)的隨機過程，還可應(yīng)用于多維非平穩(wěn)過程，這就避免了Wiener濾波器設(shè)計的困境。Kalman濾波是一套由數(shù)字計算機實現(xiàn)的實時遞推算法，它以隨機系統(tǒng)的量測作為濾波器的輸入，濾波器的輸出是對系統(tǒng)狀態(tài)的最優(yōu)估計，這一特征與確定性控制系統(tǒng)中的狀態(tài)觀測器非常相似。在Kalman濾波器出現(xiàn)以后，估計理論的發(fā)展基本上都是以它為基礎(chǔ)的一些推廣和改進。20世紀(jì)60年代，Kalman濾波在美國的太空計劃中獲得了成功的應(yīng)用，但是由于當(dāng)時計算機字長較短，濾波器在實現(xiàn)過程中有時會出現(xiàn)一些問題，即計算機求解均方誤差陣容易出現(xiàn)無窮大情況，導(dǎo)致濾波發(fā)散。平方根濾波是一種在數(shù)學(xué)上增加Kalman濾波精度的方法，Potter為“阿波羅”太空計劃開發(fā)了第一個平方根濾波算法，它推動了后來一些其他平方根濾波方法的研究，比如Bierman提出的U-D分解濾波。平方根濾波精度性能的提升是以增加計算量為代價的，目前，隨著計算機硬件技術(shù)的發(fā)展，普遍采用雙精度浮點數(shù)進行計算和存儲，多數(shù)情況下不必再像過去那樣過于關(guān)注和擔(dān)心數(shù)值問題了。經(jīng)典Kalman濾波是基于線性系統(tǒng)的估計方法，一般只能適用于線性或者非常接近于線性的非線性問題，對于非線性比較明顯的問題，Kalman濾波往往不能給出滿意的結(jié)果，需要采用非線性估計方法。最為廣泛使用的非線性估計方法是EKF（擴展卡爾曼濾波），它通過泰勒級數(shù)展開，對非線性函數(shù)進行線性化近似。同樣，以泰勒級數(shù)展開為基礎(chǔ)，若保留二階項則稱為二階卡爾曼濾波方法，理論上二階濾波降低了EKF的線性化誤差，會得到比EKF稍好的估計性能，但這是以高復(fù)雜性和計算量為代價的。迭代濾波方法也是一種對EKF濾波的修正。隨著系統(tǒng)規(guī)模的不斷增大，如何有效處理多個傳感器測量信息的問題被提出并得到了廣泛的研究。傳統(tǒng)的方法是采用集中式Kalman濾波，將所有測量信息送到中心處理器進行集中處理，雖然它的處理結(jié)果是全局最優(yōu)的，但是這種處理方式存在通信負(fù)擔(dān)重、計算量大和容錯性能差等缺點。Speyer從分散控制的角度提出了多處理器結(jié)構(gòu)思想，每個局部傳感器都有自己的分處理器，處理包括自身在內(nèi)的所有傳感器的測量信息，得到的估計結(jié)果既是局部最優(yōu)的也是全局最優(yōu)的。Willsky對Speyer的方法進行了改進，提出了一個中心處理器（主）加多個局部處理器（子）的結(jié)構(gòu)方式，主處理器完成各個子處理器結(jié)果的合成，各子處理器間不要求通信聯(lián)系，因而是相互獨立的。Carlson對分散濾波算法做了重大改進，提出了聯(lián)邦濾波算法，采用信息分享原理，把全局狀態(tài)估計信息和系統(tǒng)噪聲信息分配給各個子濾波器，且不改變各子濾波器算法的形式，聯(lián)邦濾波具有實現(xiàn)簡單、信息分享方式靈活、容錯性能好的諸多優(yōu)點。2組合導(dǎo)航將運載體從起始點引導(dǎo)到目的地的技術(shù)或方法稱為導(dǎo)航，導(dǎo)航系統(tǒng)提供的信息主要有姿態(tài)、方位、速度和位置，甚至還包括加速度和角速率，這些信息可用于運載體的正確操縱和控制。隨著技術(shù)的發(fā)展，導(dǎo)航系統(tǒng)的種類越來越多，比如慣導(dǎo)系統(tǒng)、衛(wèi)星導(dǎo)航系統(tǒng)、磁羅盤、里程儀/多普勒測速儀/空速計、氣壓高度表/雷達(dá)高度表、地標(biāo)點/地圖匹配等。這些導(dǎo)航系統(tǒng)各有特色，優(yōu)缺點并存，比如慣導(dǎo)系統(tǒng)的優(yōu)點是自主性強、動態(tài)性能好、導(dǎo)航信息全面且輸出頻率高，但其缺點是誤差隨時間不斷累積，長期精度差；衛(wèi)星導(dǎo)航系統(tǒng)的優(yōu)點是精度高、誤差不隨時間增大，缺點是導(dǎo)航信息不夠全面、頻帶窄、信號容易受到干擾、在室內(nèi)等環(huán)境下接收不到衛(wèi)星信號而無法使用。在許多對導(dǎo)航性能要求苛刻的任務(wù)中，無論是精度要求高還是可靠性要求高，任何單一的導(dǎo)航系統(tǒng)可能都無法滿足要求，這就需要使用多種導(dǎo)航系統(tǒng)同時對運載體進行導(dǎo)航信息測量，再對所有測量信息作綜合處理（包括檢測、結(jié)合、相關(guān)和估計），從而得到更為準(zhǔn)確和可靠的導(dǎo)航結(jié)果。這種對多種導(dǎo)航信息作綜合處理的技術(shù)就稱為組合導(dǎo)航技術(shù)。從上述對慣導(dǎo)和衛(wèi)星導(dǎo)航的優(yōu)缺點描述中可以看出，兩者性能具有非常強的互補性，因而慣性/衛(wèi)星組合導(dǎo)航被公認(rèn)為是最佳的組合導(dǎo)航方案。組合導(dǎo)航系統(tǒng)的設(shè)計一般都采用Kalman濾波器，Kalman濾波器最早和最成功的應(yīng)用實例便是在導(dǎo)航領(lǐng)域。1960年卡爾曼在美國國家航空航天局埃姆斯研究中心（NASAAmesResearchCenter）訪問時，StanleySchmidt發(fā)現(xiàn)Kalman濾波方法對于解決阿波羅計劃的軌道預(yù)測很有用，后來阿波羅登月飛船的導(dǎo)航系統(tǒng)便使用了Kalman濾波器，通常認(rèn)為Schmidt首次實現(xiàn)了Kalman濾波器。此外，美國在航天飛機、潛艇和無人航空航天飛行器（比如巡航導(dǎo)彈）上均使用了Kalman濾波器。第2章捷聯(lián)慣導(dǎo)姿態(tài)解算基礎(chǔ)第2章捷聯(lián)慣導(dǎo)姿態(tài)解算基礎(chǔ) PAGER12.1反對稱陣及其矩陣指數(shù)函數(shù) 12.1.1反對稱陣 12.1.2反對稱陣的矩陣指數(shù)函數(shù) 32.2方向余弦陣與等效旋轉(zhuǎn)矢量 42.2.1方向余弦陣 42.2.2等效旋轉(zhuǎn)矢量 52.3方向余弦陣微分方程及其求解 82.3.1方向余弦陣微分方程 82.3.2方向余弦陣微分方程的求解 82.4姿態(tài)更新的四元數(shù)表示 112.4.1四元數(shù)的基本概念 112.4.2四元數(shù)微分方程 142.4.3四元數(shù)微分方程的求解 162.5等效旋轉(zhuǎn)矢量微分方程及其泰勒級數(shù)解 172.5.1等效旋轉(zhuǎn)矢量微分方程 172.5.2等效旋轉(zhuǎn)矢量微分方程的泰勒級數(shù)解 202.6圓錐運動條件下的等效旋轉(zhuǎn)矢量算法 222.6.1圓錐運動的描述 221-3"\h\z\u第章地球形狀與重力場基礎(chǔ) 13.1地球的形狀描述 13.2地球的正常重力場 73.3地球重力場的球諧函數(shù)展開 113.3.1球諧函數(shù) 113.3.2地球引力位函數(shù) 193.3.3重力位及重力計算 243.1地球的形狀描述實際的地球表面是一個凹凸不平、形狀十分復(fù)雜的物理面，難以準(zhǔn)確量化描述。為了研究方便，假想海洋表面靜止，并將其向陸地延伸，所得到的封閉曲面稱為大地水準(zhǔn)面，大地水準(zhǔn)面包圍的形體稱為大地水準(zhǔn)體。由于地球內(nèi)部密度分布不均勻和表面形狀起伏的影響，大地水準(zhǔn)體也是一個不規(guī)則的幾何體。實用中希望使用比較簡單的數(shù)學(xué)方程來擬合地球幾何形狀，按精度從低到高有如下三種近似：①近似成圓球體，中心選擇在地球質(zhì)心上，半徑約6371km，該描述比較粗略，適用于對精度要求不高的場合；②近似成旋轉(zhuǎn)橢球體，地球自轉(zhuǎn)軸（極軸）與一橢圓短半軸重合，橢圓的橢圓度約1/300（長短半軸相差約21km），橢圓繞其短半軸旋轉(zhuǎn)構(gòu)成旋轉(zhuǎn)橢球體的表面，該描述中地球赤道是圓形的；③近似成三軸橢球體，其中赤道是橢圓的，赤道橢圓度約1/100000（長短半軸相差約60m）。三軸橢球體描述雖然比旋轉(zhuǎn)橢球體更精確，但前者的相關(guān)計算比后者復(fù)雜許多，考慮到三軸橢球體的赤道橢圓度不大，可將其似成圓形的，因此旋轉(zhuǎn)橢球體應(yīng)用最為廣泛。下面先簡要介紹一下地球旋轉(zhuǎn)橢球上的一些基本概念。參見圖3.1-1，地球自轉(zhuǎn)軸的南端點和北端點分別稱為南極（）和北極（），包含南北極點的平面稱為子午面，子午面與旋轉(zhuǎn)橢球面的交線稱為子午圈（或經(jīng)圈）。通過英國格林尼治的經(jīng)線稱為本初子午線（或零度經(jīng)線）。任一經(jīng)線所在子午面與本初子午面之間的夾角，定義為經(jīng)度（記為），夾角方向與地球自轉(zhuǎn)軸同方向，取值范圍-180°～180°。包含旋轉(zhuǎn)橢球中心且垂直于自轉(zhuǎn)軸的平面稱為赤道面，赤道面與旋轉(zhuǎn)橢球面的交線稱為赤道，平行于赤道面的平面與橢球面的交線稱為緯圈。圖3.1-1旋轉(zhuǎn)橢球基本概念圖3.1-2子午圈橢圓對于地球旋轉(zhuǎn)橢球體而言，確定其三維形狀參數(shù)的關(guān)鍵在于確定二維子午圈橢圓。1子午圈橢圓參見圖3.1-1，建立地心右手直角坐標(biāo)系，常稱為地心地固坐標(biāo)系（EarthCenterEarthFixed，ECEF），坐標(biāo)原點選在地心，軸為自轉(zhuǎn)軸且指向北極，軸指向赤道與本初子午線的交點，軸在赤道平面內(nèi)且指向90°經(jīng)線，ECEF系與地球固聯(lián)，即跟隨地球自轉(zhuǎn)一起相對慣性空間轉(zhuǎn)動。對子午圈橢圓，不失一般性，選擇本初子午線橢圓作為研究對象，如圖3.1-2所示。橢圓上任一點與地心連線與軸的夾角稱為地心緯度，記為，取值范圍-90°～90°，南緯為負(fù)北緯為正。過點的橢圓法線與軸的夾角稱為地理緯度，簡稱緯度，記為，取值范圍-90°～90°。此外，與地心緯度對應(yīng)的方向稱為地心垂線，而與地理緯度對應(yīng)的方向稱為地理垂線。橢圓形狀完全由其長半軸和短半軸確定，但在涉及橢圓的計算中，為了方便常引入扁率和偏心率概念。橢圓方程為（3.1-1）其中，和分別為橢圓長半軸和短半軸。橢圓扁率（或稱橢圓度，flattening）定義為（3.1-2）橢圓偏心率（eccentricity）定義為（3.1-3）第二偏心率定義為且有（3.1-4）相對于第二偏心率而言，式（3.1-3）有時也稱為第一偏心率。由式（3.1-2）和（3.1-3）可分別得（3.1-5）（3.1-6）比較上述兩式，可得和（3.1-7）將橢圓方程（3.1-1）兩邊同時對求導(dǎo)，并考慮到式（3.1-6），得（3.1-8）上式移項整理得（3.1-9）式中表示橢圓在點的切線的斜率，顯然，切線與法線之間是相互垂直的，法線的斜率為，則有（3.1-10）從上式可解得（3.1-11）將式（3.1-6）和式（3.1-11）代入橢圓方程（3.1-1），可求得以地理緯度為參數(shù)的橢圓參數(shù)方程（3.1-12）參見圖3.1-2，記線段長度，則有（3.1-13）比較式（3.1-13）和式（3.1-12）中的第一式，可得（3.1-14）因而參數(shù)方程（3.1-12）可簡寫為（3.1-15）最后，比較一下地球表面上同一點的地理緯度與地心緯度之間的差別，或者說，地理垂線與地心垂線之間的偏差。對于地心緯度，注意到，根據(jù)式（3.1-11），有（3.1-16）記地理緯度與地心緯度之間的偏差量，則有（3.1-17）將和視為小量，上式可近似為（3.1-18）或者（3.1-19）其中，根據(jù)式（3.1-7）近似取。例如，若取橢球扁率，則在緯度=45°處可求得地理緯度與地心緯度的最大偏差值約為11.5′。2旋轉(zhuǎn)橢球表面的曲率半徑導(dǎo)航過程中，運載體在地球橢球表面附近移動，為了在合適的坐標(biāo)系（通常指地理坐標(biāo)系）下進行三維定位解算，旋轉(zhuǎn)橢球表面的幾何曲率半徑是一個非常重要的參數(shù)。首先給出法截線的概念。參見圖3.1-3，包含橢球面上某點法線的平面稱為法截面，法截面與子午面之間的夾角記為，法截面與橢球的交線稱為法截線。顯然，當(dāng)法截面包含南北極點時，法截線即為子午圈；當(dāng)法截面垂直于子午面時，法截線稱為卯酉圈。圖3.1-3法截線及其局部坐標(biāo)系在圖3.1-3中，不失一般性，假設(shè)點在本初子午線上，以為坐標(biāo)原點建立局部直角坐標(biāo)系，其中軸沿法線向外，軸沿法截線切線方向。不難看出，若坐標(biāo)系先繞軸轉(zhuǎn)動角度，然后繞軸轉(zhuǎn)動角度，再作相應(yīng)平移，則可得坐標(biāo)系（即地心直角坐標(biāo)系）。根據(jù)橢圓參數(shù)方程（3.1-12）知點在坐標(biāo)系下的坐標(biāo)為，此即前述平移的坐標(biāo)值，所以和兩坐標(biāo)系之間的坐標(biāo)變換關(guān)系為（3.1-20）將上式代入旋轉(zhuǎn)橢球方程，移項整理得（3.1-21）由于在局部坐標(biāo)系下表示的法截線方程必有，因而上式中不含項。將式（3.1-21）對坐標(biāo)求一次導(dǎo)和二次導(dǎo)，代入平面曲線的曲率計算公式，可得法截線在點的曲率（3.1-22）特別地，當(dāng)角度或時，分別有（3.1-23）（3.1-24）通常稱為子午圈主曲率半徑；而稱為卯酉圈主曲率半徑。在圖3.1-2中卯酉圈曲率半徑即對應(yīng)于線段的長度。除在地理緯度處有外，其它緯度處總有。3大地坐標(biāo)及其變化率在旋轉(zhuǎn)橢球表面上點處，緯圈切線與經(jīng)圈切線相互垂直，且兩切線同時垂直于橢球面的法線。在橢球表面上定義直角坐標(biāo)系：點為坐標(biāo)原點（重記為），緯圈切線指東、經(jīng)圈切線指北、橢球面法線指天分別為軸、軸和軸的正向。參見圖3.1-4，若某點在坐標(biāo)系中的直角坐標(biāo)為，顯然在橢球面點的法線上，稱為點的地理高度。以為原點建立坐標(biāo)系，其三軸分別平行于的同名坐標(biāo)軸，稱為當(dāng)?shù)氐乩碜鴺?biāo)系，簡記為系。點相對于地球橢球的空間位置可用所謂的大地坐標(biāo)（地理坐標(biāo)）表示，記為。（a）經(jīng)度變化率（b）緯度變化率圖3.1-4速度引起的經(jīng)緯度變化在圖3.1-4中，如果點對地球坐標(biāo)系的速度在系的投影記為。注意到，軸與緯圈相切，兩者在同一個平面內(nèi)，因而僅會引起經(jīng)度的變化，有（3.1-25）同理，軸與經(jīng)圈相切，兩者在同一平面內(nèi)，因而速度僅會引起緯度的變化，考慮到點所在子午圈的曲率半徑為，則有（3.1-26）對于地理高度為的點，假設(shè)其速度為，根據(jù)圖3.1-4中幾何關(guān)系，有（3.1-27）（3.1-28）上述兩式分別代入式（3.1-25）和（3.1-26），并考慮到天向速度僅引起地理高度變化，得速度與大地坐標(biāo)之間的關(guān)系，分別為（3.1-29）（3.1-30）（3.1-31）地理坐標(biāo)系與地球坐標(biāo)系之間的轉(zhuǎn)動關(guān)系可以用方向余弦陣表示，常稱之為位置矩陣，記為。參見圖3.1-5，系先繞軸轉(zhuǎn)動，接著繞軸轉(zhuǎn)動，再繞軸轉(zhuǎn)動，這時系三軸可與系相應(yīng)坐標(biāo)軸平行。據(jù)此，可計算得位置矩陣圖3.1-5系至系的三次轉(zhuǎn)動（3.1-32）對上式兩邊同時微分，得（3.1-33）上式與方向余弦陣微分公式對比，并將經(jīng)緯度公式（3.1-29）和（3.1-30）代入，分別記為，可得（3.1-34）這便是載體運動線速度引起導(dǎo)航系旋轉(zhuǎn)角速度的計算公式。4大地坐標(biāo)與地心直角坐標(biāo)之間轉(zhuǎn)換下面給出同一地點的地理坐標(biāo)與地心直角坐標(biāo)之間的相互轉(zhuǎn)換關(guān)系。（1）由求解根據(jù)式（3.1-15）和圖3.1-2，不難求得（3.1-35）其中，子午圈曲率半徑的計算見式（3.1-14）或重寫如下（3.1-36）（2）由求解首先，由式（3.1-35）中第二式除第一式，得（3.1-37）由上式可直接解得經(jīng)度（3.1-38）其次，對于緯度，不能求得顯式表示，通常采用迭代算法，推導(dǎo)過程如下：由式（3.1-35）中第一式和第二式實施如下運算（3.1-39）由式（3.1-35）中第三式移項整理，可得（3.1-40）式（3.1-40）除以式（3.1-39），得（3.1-41）將式（3.1-36）改寫成，再代入上式，得（3.1-42）如記，則由上式可構(gòu)造出求解緯度正切值的迭代公式，如下（3.1-43）令迭代初值，一般經(jīng)過5～6次迭代便可達(dá)到足夠的數(shù)值計算精度，再求反正切即可獲得緯度。最后，根據(jù)（3.1-39）求解高度，得（3.1-44）式（3.1-38）、（3.1-43）和（3.1-44）即為由求解的算法。3.2地球的正常重力場在地球的大地水準(zhǔn)體描述中，水準(zhǔn)體表面是地球?qū)嶋H重力場的一個等位面，每一點的重力方向均與該點所在等位面相垂直，實際的重力方向一般稱為天文垂線，或稱真垂線。由于實際地球內(nèi)部密度分布不均勻，并且表面凹凸不平，大地水準(zhǔn)面不規(guī)則，因而實際重力的大小和方向也不規(guī)則。與地球的幾何形狀描述類似，也希望使用一個簡單的數(shù)學(xué)函數(shù)來擬合地球重力場，這個簡單函數(shù)表示的重力場就稱為正常重力場。重力是地球萬有引力和離心力共同作用的結(jié)果，參見圖3.2-1，在點處，重力矢量是引力矢量和離心力矢量的合力。圖3.2-1重力矢量1．圓球假設(shè)下的地球重力若將地球視為圓球體并且認(rèn)為密度均勻分布，那么地球引力指向地心，根據(jù)牛頓萬有引力定律，地球?qū)ζ浔砻婊蛲獠繂挝毁|(zhì)點的引力大小為（3.2-1）其中，為萬有引力常數(shù)，為地球質(zhì)量，記為地心引力常數(shù)，為質(zhì)點至地心的距離。由于地球繞其極軸存在自轉(zhuǎn)角速率，使得與地球表面固連的單位質(zhì)點受到離心力的作用，其大小為（3.2-2）其中，為圓球半徑，是地理緯度（在圓球假設(shè)中即為地心緯度）。重力是引力與離心力的合力，引力與離心力之間的夾角為，根據(jù)余弦定理，在緯度為的地方上點的重力大小為（3.2-3）其中，為赤道上的重力大小，為赤道上的離心力與重力加速度的比值。實踐表明，基于圓球假設(shè)的重力公式（3.2-3）與實際橢球地球相比，在高緯度地區(qū)偏小將近2mg，部分原因歸結(jié)于實際橢球地球在高緯度地區(qū)半徑縮小，實際引力增大。為了更精確地計算正常重力值，需要在橢球條件下進行重力推導(dǎo)。2．旋轉(zhuǎn)橢球假設(shè)下的地球重力對于地球旋轉(zhuǎn)橢球體，假設(shè)在橢球表面上重力矢量處處垂直于表面，也就是說，旋轉(zhuǎn)橢球表面為重力的一個等位面，意大利人索密里安（Somigliana）于1929年經(jīng)過嚴(yán)密推導(dǎo)（過程比較復(fù)雜，從略），獲得了旋轉(zhuǎn)橢球體的正常重力公式，如下（3.2-4）其中，和分別為旋轉(zhuǎn)橢球的赤道長半軸和極軸短半軸，和分別為赤道重力和極點重力，為地理緯度處橢球表面的重力大小。 對于赤道重力和極點重力，近似有（3.2-5）其中，為旋轉(zhuǎn)橢球幾何形狀扁率，約等于1/298；并且（3.2-6）為赤道上的離心力與重力加速度的比值，約等于1/288。記地球重力扁率為（3.2-7）其實際值約等于1/189。將式（3.2-5）代入上式，并展開成關(guān)于和級數(shù)形式，可得（3.2-8）式（3.2-8）建立了重力扁率與橢球形狀扁率之間的關(guān)系，若忽略高階小量，近似有（3.2-9）上式是利用重力測量方法確定地球形狀參數(shù)的基本公式，它最早在1743年由法國數(shù)學(xué)家克萊羅（A.C.Clairaut）求得，通常稱為克萊羅公式?？巳R羅在推導(dǎo)公式（3.2-9）時作了如下前提假設(shè)：地球是由密度不同的均勻物質(zhì)層圈組成的橢球體，各橢球面都是重力等位面，且各層密度由地心向外有規(guī)律的減小。需要說明的是，克萊羅公式是近似成立的公式，而索密里安公式卻是理論上嚴(yán)格成立的，并且后者的前提條件更加寬松，對橢球體密度分布不做任何限制。若將和代入索密里安公式（3.2-4），展開為關(guān)于和的級數(shù)形式，可得（3.2-10）考慮到如下三角函數(shù)恒等式（3.2-11）將上式移項，有（3.2-12）再將式（3.2-12）代入式（3.2-10），忽略高階小量，可得（3.2-13）其中（3.2-14）若將克萊羅公式（3.2-9）代入上式，還可得（3.2-15）式（3.2-13）為索密里安公式的良好近似（最大誤差約為），實用中的正常重力模型常常是以該形式給出的，它比式（3.2-4）的計算量稍小些。歷史上曾給出過以下一些重要的正常重力模型。（1）1901年，德國人赫爾默（Helmert）根據(jù)當(dāng)時波斯坦系統(tǒng)的幾千個重力測量結(jié)果，給出正常重力公式為：（3.2-16）上式稱為赫爾默正常重力公式，其中重力扁率，利用克雷諾定理，可以計算出相應(yīng)的參考橢球的扁率。（2）1909年，美國人海福特（Hayford）根據(jù)美國當(dāng)時的大地測量結(jié)果給出了一個參考橢球，它的赤道半徑和幾何扁率；1928年，芬蘭人海斯卡寧（Heiskanen）根據(jù)當(dāng)時的重力測量結(jié)果計算出正常場地球模型在赤道上的重力為。若取地球自轉(zhuǎn)角速率，根據(jù)上述四個獨立參數(shù)，可計算得（3.2-17）1930年，國際大地測量與地球物理聯(lián)合會（IUGG）將上式定為國際正常重力公式。（3）利用現(xiàn)代衛(wèi)星測量技術(shù)，IUGG于1979年通過了1980大地參考系，與其對應(yīng)的正常重力公式為（3.2-18）（4）1987年，WGS-84（WorldGeodeticSystem1984）大地坐標(biāo)系給出的地球參數(shù)為：半長軸，扁率，地心引力常數(shù)（含大氣層），地球自轉(zhuǎn)角速率。如不考慮大氣層影響，可推導(dǎo)得正常重力公式（3.2-19）3．重力與海拔高度的關(guān)系在地球表面附近的重力場中，引力與離心力相比前者占主要成分，重力隨海拔高度增加而減小，其變化規(guī)律與引力隨距離增加而減小的規(guī)律近似相同。分析高度影響時，不妨將地球近似成圓球且質(zhì)量集中在地心，地球?qū)Ω叨葹榈膯挝毁|(zhì)點的引力為（3.2-20）對上式兩邊同時微分，得（3.2-21）其中（3.2-22）若設(shè)地心引力常數(shù)和地球平均半徑，則有，這表示在地球表面附近高度每升高1m，引力值（或重力值）將減小（約0.3）。綜合式（3.2-13）和式（3.2-21），可求得地球表面附近正常重力隨緯度和高度變化的實用計算公式，即在大地坐標(biāo)處的重力值，記為（3.2-23）值得注意的是，上式只給出了重力的大小，其方向應(yīng)該是點處的鉛垂向（真垂線），而不是該點的地理法向（地理垂線），如圖3.2-2所示。1967年，Heiskanen給出了真垂線和地理垂線之間夾角的近似公式（3.2-24）其中，。據(jù)上式計算，在緯度°處，高度每上升1km，約增加。顯然，除了赤道和極點外，只有在旋轉(zhuǎn)橢球表面上（）真垂線與地理垂線才能嚴(yán)格重合。圖3.2-2正常重力場的垂線偏差當(dāng)以地理坐標(biāo)系（系）作為慣性導(dǎo)航參考坐標(biāo)系時，為了扣除重力的影響，需將重力矢量投影到地理坐標(biāo)系，表示為（3.2-25）上式的三個分量依次表示重力矢量在東、北和天向的投影值，同樣在緯度°處，高度每上升1km，北向重力分量將變化，這一影響在高空高精度慣性導(dǎo)航中是需要考慮和補償?shù)?。更高精度的且與實際地球相關(guān)的重力場計算和補償詳見3.3節(jié)介紹。3.3地球重力場的球諧函數(shù)模型在3.2節(jié)中正常重力場描述的是規(guī)則地球產(chǎn)生的重力場，但實際地球并不規(guī)則，正常重力場只是實際重力場的一個較好的近似，為了更加細(xì)致地刻畫實際地球的重力場，需引入球諧函數(shù)和重力位等概念，這在高精度慣性導(dǎo)航系統(tǒng)的重力場建模和補償中十分有用。3.3.1球諧函數(shù)的基本概念1.拉普拉斯方程如果三元函數(shù)在空間區(qū)域上滿足如下偏微分方程（3.3-1）則稱是區(qū)域上的調(diào)和函數(shù)，或稱諧函數(shù)，上述方程稱為拉普拉斯方程（Laplaceequation）。方程（3.3-1）可簡記為（3.3-2）其中，算符（3.3-3）稱為拉普拉斯算子。參見圖3.3-1，球面上一點，其直角坐標(biāo)與極坐標(biāo)之間的坐標(biāo)轉(zhuǎn)換關(guān)系為圖3.3-1直角坐標(biāo)與極坐標(biāo)（3.3-4）其中，為點距坐標(biāo)中心的距離，為經(jīng)度，為極距或稱余緯。經(jīng)推導(dǎo)（過程繁瑣從略，附錄？），可求得拉普拉斯算子（3.3-3）在球坐標(biāo)下的表示（3.3-5）因而，用球坐標(biāo)表示的拉普拉斯方程為（3.3-6）2.球坐標(biāo)下拉普拉斯方程的求解求解拉普拉斯方程（3.3-6）一般采用分離變量法，假設(shè)它的解為（3.3-7）將式（3.3-7）代入式（3.3-6），得（3.3-8）上式兩邊同時乘以，移項可得（3.3-9）式（3.3-9）左邊只與函數(shù)有關(guān)，而右邊只與函數(shù)有關(guān)，欲使兩邊相等，必須使它們都等于同一常值，顯然該常值總可以表示為形式，所以有（3.3-10）（3.3-11）其中，、分別記為和，角標(biāo)表示與常值參數(shù)相關(guān)的解。注意到式（3.3-10）變成了常微分方程，稍微整理，可得（3.3-12）上式為歐拉常微分方程，其通解為（3.3-13）其中，和為任意常系數(shù)。對式（3.3-11）左邊的第一項進一步展開，該式可化為（3.3-14）上式兩邊同時乘以并移項整理，得（3.3-15）再次對式（3.3-15）采用分離變量法，令并代入，得（3.3-16）上式兩邊同時乘以并移項整理，得（3.3-17）式（3.3-17）的左邊只與函數(shù)有關(guān)，而右邊只與函數(shù)有關(guān)，同樣它們都應(yīng)等于同一常數(shù)，所以有（3.3-18）（3.3-19）在方程（3.3-18）中，考慮到經(jīng)度坐標(biāo)的自然周期條件，常數(shù)取值必須為，由此求得方程（3.3-18）的解為（3.3-20）其中，和為任意常系數(shù)。將式（3.3-19）整理，考慮到，可得（3.3-21）其中，的角標(biāo)表示與常值參數(shù)和相關(guān)的解。令，則有，并且有（3.3-22）（3.3-23）將式（3.3-22）和（3.3-23）代入式（3.3-21），并記為，可得（3.3-24）上式稱為連帶勒讓德方程（associatedLegendreequation），由可知自變量的取值范圍為。特別地，當(dāng)時稱其為勒讓德方程，重寫如下（3.3-25）至此，將拉普拉斯偏微分方程的求解問題轉(zhuǎn)化成了三個常微分方程的求解，即方程（3.3-12）、（3.3-18）和（3.3-25），其中前面兩個方程求解都比較容易，并已經(jīng)獲得解決，接下來主要介紹勒讓德方程（3.3-25）的求解及其解的特性。3.勒讓德多項式在勒讓德方程（3.3-25）中，實際使用時參數(shù)多為非負(fù)整數(shù)，求解該方程一般采用級數(shù)法，設(shè)級數(shù)解為（3.3-26）將上式及其一階和二階導(dǎo)數(shù)代入式（3.3-25），可求得待定系數(shù)（過程繁瑣從略），從而獲得通解（3.3-27）其中，和為任意常系數(shù)；稱為第一類勒讓德多項式，最高階次為，它在區(qū)間上為有限值；稱為第二類勒讓德多項式，為無窮級數(shù)，它在區(qū)間上是無界的。這里主要討論，它的級數(shù)表達(dá)式為（3.3-28）其中，表示取整運算。實際上，式（3.3-28）亦可表示成如下階導(dǎo)數(shù)的形式（3.3-29）在上式中，只需用二項式定理展開，然后逐項求階導(dǎo)，即可驗證其與式（3.3-28）相等。式（3.3-29）一般稱為勒讓德多項式的羅德里格（Rodrigues）表示。勒讓德多項式具有如下遞推公式（3.3-30a）（3.3-30b）（3.3-30c）上式為遞推求解高次的勒讓德多項式帶來了極大的方便。若將勒讓德方程（3.3-25）兩邊同時對微分次，并記，可得（3.3-31）引入新的變量，則有（3.3-32）顯然，式（3.3-32）與式（3.3-24）的含義完全相同，均表示連帶勒讓德方程。上式的推導(dǎo)過程說明，如果是勒讓德方程的解，則必定是連帶勒讓德方程的解，所以連帶勒讓德方程（3.3-32）的通解可以表示為（3.3-33）其中，和為任意常系數(shù)，并且有（3.3-34）（3.3-35）稱為次階第一類連帶勒讓德多項式，在區(qū)間上為有限值；稱為次階第二類連帶勒讓德多項式，在區(qū)間上是無界的。與式（3.3-27）一樣，這里主要討論。在式（3.3-34）中，當(dāng)時顯然有，因此可將勒讓德多項式記作。后面為了敘述簡便，在不引起混淆情況下，統(tǒng)一將勒讓德多項式和連帶勒讓德多項式稱為勒讓德多項式。為了便于直觀了解，表3.3-1給出一些低階的勒讓德多項式的顯式表達(dá)式，圖3.3-2給出了相應(yīng)的曲線。表3.3-1勒讓德多項式顯式0123450123--4-5注：表中簡記，若設(shè)則有。式（3.3-30）和（3.3-34）所示規(guī)律在表中得到了清楚的體現(xiàn)。圖3.3-2勒讓德多項式曲線連帶勒讓德多項式具有如下遞推公式（3.3-36a）（3.3-36b）（3.3-36c）對照表3.3-1，不難看出式（3.3-36a）給出了表中的對角線元素，式（3.3-36b）給出了次對角線元素，式（3.3-36c）給出剩余其它元素。顯然，式（3.3-30）可以看作式（3.3-36）中當(dāng)時的特例。此外，連帶勒讓德多項式還存在如下求導(dǎo)公式（3.3-37a）（3.3-37b）若固定某一非負(fù)整數(shù)，勒讓德多項式函數(shù)系在區(qū)間上是完備正交系，有（3.3-38）若記正規(guī)化勒讓德多項式，則是完備的正規(guī)正交系。如果一元函數(shù)在區(qū)間上單值連續(xù)，則可按勒讓德多項式展開且收斂，即有（3.3-39）其中（3.3-40）式（3.3-39）右端的級數(shù)稱為傅里葉—勒讓德級數(shù)，它是一種廣義傅里葉級數(shù)。最后，給出勒讓德多項式的一些重要特點：（1）當(dāng)時，在端點處取值為或，或者說，而當(dāng)時端點處均為0，即；（2）當(dāng)為偶數(shù)時，是偶函數(shù)，否則為奇函數(shù)；（3）在開區(qū)間上，有個零點，特別地，有個零點，而在上無零點。4.球諧函數(shù)（1）面球諧函數(shù)在勒讓德多項式中，做變量替換，考慮到的取值范圍為，則方程（3.3-21）的解為（3.3-41）因此，在微分方程（3.3-15）中，對于任一常數(shù)（），其解為（3.3-42）當(dāng)取遍時，所有解的線性組合記為方程（3.3-15）的通解，為（3.3-43）易知，在中共有個相互獨立的基本函數(shù)，二元自變量的定義域為單位球面，常稱這些函數(shù)為次面球諧函數(shù)（或球面調(diào)和函數(shù)，sphericalharmonicfunction），相應(yīng)地，式（3.3-15）稱為次球面函數(shù)方程。根據(jù)勒讓德多項式的特點，考慮到且在上是單調(diào)的，可知：當(dāng)為偶數(shù)時關(guān)于（赤道）是偶對稱的，而當(dāng)為奇數(shù)時關(guān)于赤道是奇對稱的；的零點個數(shù)與相同，即在上有個零點，特別地，無零點，而有個零點。還易知，在上均勻分布著個零點，特別地，當(dāng)時無零點。所以，對于余弦面球諧函數(shù)（）而言，其特點為：（1）當(dāng)時，有，這顯示面球諧函數(shù)與經(jīng)度無關(guān)，的個零點形成函數(shù)值等于0的條緯線，把球面分隔成個條帶，因而當(dāng)時面球諧函數(shù)也形象地稱為帶諧函數(shù)（或主諧函數(shù)）；（2）當(dāng)時，有，除極點和外無其它零點，而的個零點形成函數(shù)值等于0的條經(jīng)線，把球面分隔成個扇形，因而當(dāng)時面球諧函數(shù)也形象地稱為扇諧函數(shù)；（3）當(dāng)時，在中，一方面的個零點形成函數(shù)值等于0的條緯線，另一方面的個零點形成函數(shù)值等于0的條經(jīng)線，這些函數(shù)值等于0的緯線和經(jīng)線將球面分隔成個方塊（極點周圍為三角塊），因而當(dāng)時面球諧函數(shù)形象地稱為田諧函數(shù)。此外，對于任意球面諧函，由函數(shù)值為0的經(jīng)線（或緯線）分隔的相鄰區(qū)域的函數(shù)值正負(fù)符號正好相反。圖3.3-3給出了時的帶諧函數(shù)、扇諧函數(shù)和田諧函數(shù)示意圖。(a)帶諧函數(shù)(b)田諧函數(shù)（含中間三幅小圖）(c)扇諧函數(shù)圖3.3-3面球諧函數(shù)（，圖中紅色區(qū)域函數(shù)值為正而藍(lán)色為負(fù)）對于正弦面球諧函數(shù)，它與在經(jīng)度上相位相差正好為，即前者繞極軸轉(zhuǎn)動即可得后者，兩者其它特點完全一致，無需贅述。 在式（3.3-43）中，當(dāng)取遍時，所有的面球諧函數(shù)構(gòu)成面球諧函數(shù)系，記作。面球諧函數(shù)系是球面上的完備正交系，任意兩函數(shù)之間存在如下正交關(guān)系（3.3-44）其中在被積函數(shù)中，符號表示或任選其一，因而乘積表示四種可能結(jié)果之一。顯然，是球面上的完備正規(guī)正交系。在球面上的二元連續(xù)函數(shù)，總可按面球諧函數(shù)系展開成級數(shù)形式，即（3.3-45）其中系數(shù)在實際工作中，球面函數(shù)的具體表達(dá)式往往是未知的，只能在球面上測量獲得一些離散點處的函數(shù)值，一般使用有限階次的面球諧函數(shù)級數(shù)來擬合，便于后續(xù)插值等應(yīng)用。假設(shè)有個球面函數(shù)的測量數(shù)據(jù)，記為，代入式（3.3-45），略去高于次的面球諧函數(shù)，得（3.3-46）式中，簡記，表示測量誤差。若將個測量數(shù)據(jù)合并一起寫成矩陣形式，則有（3.3-47）其中在式（3.3-46）中，不大于次的面球諧函數(shù)待定系數(shù)共有個，欲使方程（3.3-47）有唯一解，測量數(shù)據(jù)不得少于待定系數(shù)個數(shù)，即必須滿足，通常為了提高待定系數(shù)的估計精度，選。采用最小二乘法求解式（3.3-47），得（3.3-48）最后指出，類似于一元函數(shù)的有限項傅里葉級數(shù)擬合，在次的面球諧函數(shù)擬合中，認(rèn)為球面函數(shù)以低頻成分（長波分量）為主，而忽略高頻成分（短波分量）的影響；帶諧函數(shù)和扇諧函數(shù)最多將半球分別分隔成和份，因此，次的面球諧函數(shù)擬合的最高角度分辨率近似為。（2）體球諧函數(shù)在方程（3.3-9）中，即球坐標(biāo)拉普拉斯方程（3.3-8）中，對于任一非負(fù)整數(shù)，其解為（3.3-49）當(dāng)遍歷，所有解的線性組合記為方程（3.3-8）的“形式通解”，如下（3.3-50）其中，、、和均為常系數(shù)，可見和兩者中的常系數(shù)相乘合并之后，不增加常系數(shù)的數(shù)目。為了簡寫方便，以后可以去掉這些常系數(shù)的上標(biāo)“”。在式（3.3-50）中，有（3.3-51）（3.3-52）式（3.3-51）一般適用于場合；而式（3.3-52）適用于場合，在地球引力場中主要應(yīng)用的是后者，為簡潔，后面討論時將去掉右上標(biāo)“”。在式（3.3-52）中，記（3.3-53）稱為次體球諧函數(shù)，它與面球諧函數(shù)之間僅相差了與球半徑有關(guān)的因子。式（3.3-52）說明，調(diào)諧函數(shù)可以展開成一系列體球諧函數(shù)之和，或者說，調(diào)諧函數(shù)可以展開成一系列與球半徑有關(guān)的面球諧函數(shù)之和。類似于球面函數(shù)情形，在實際工作中，諧函數(shù)的具體表達(dá)式也往往是未知的，只能在定義域上測量獲得一些離散點處的函數(shù)值，需使用有限階次的體球諧函數(shù)級數(shù)來擬合。同樣，假設(shè)測量數(shù)據(jù)為，代入式（3.3-52），略去高于次的面球諧函數(shù)，得（3.3-54）若將所有的測量數(shù)據(jù)合并一起，可得測量方程組，再求解出相應(yīng)系數(shù)向量，不再贅述。 如果單位球面是諧函數(shù)定義域的邊界，且測量數(shù)據(jù)均取樣在單位球面上，即為，則式（3.3-54）與（3.3-46）完全相同，這種情況下同樣可以求得系數(shù)向量，實際上這恰好反映了諧函數(shù)的性質(zhì)：諧函數(shù)可以由它在定義域邊界上的函數(shù)值唯一確定，這也正是高等數(shù)學(xué)中格林公式定理的反映。習(xí)慣上，常將面球諧函數(shù)和體球諧函數(shù)統(tǒng)稱為球諧函數(shù)。3.3.2地球引力位函數(shù)參見圖3.3-4，假設(shè)地球總質(zhì)量為，是其上的質(zhì)量微元。在地球外部空間上有一質(zhì)量為的質(zhì)點，和之間的距離為。圖3.3-4地球?qū)ζ渫獠恳毁|(zhì)點的引力根據(jù)萬有引力定律，和之間的引力大小為（3.3-55）式中，為萬有引力常數(shù)。若質(zhì)點沿矢徑方向移動，注意到和方向正好相反，可得引力做功為（3.3-56）若質(zhì)點從無窮遠(yuǎn)處移動至半徑處，則引力做功為（3.3-57）引力做功必然等于質(zhì)點位能（勢能）的減少量?？梢?，在質(zhì)量微元產(chǎn)生的引力場中，若將某一單位質(zhì)量質(zhì)點從無窮遠(yuǎn)處移動至半徑處，單位質(zhì)點的位能減少為，將該減少量定義為的引力位函數(shù)，簡稱位函數(shù)，記為（3.3-58）位函數(shù)是標(biāo)量，整個地球質(zhì)量產(chǎn)生的位函數(shù)等于其上各質(zhì)量微元的位函數(shù)之和，即（3.3-59）若建立直角坐標(biāo)系，和的坐標(biāo)分別記為和，則有兩點間距離公式，將對坐標(biāo)求一階和二階偏導(dǎo)，分別得（3.3-60）（3.3-61）再求位函數(shù)對坐標(biāo)的二階偏導(dǎo)數(shù)，得（3.3-62）同理，求位函數(shù)對坐標(biāo)和的二階偏導(dǎo)數(shù)，有（3.3-63）（3.3-64）上述三式（3.3-62）、（3.3-63）和（3.3-64）相加，得（3.3-65）可見，位函數(shù)在地球質(zhì)量的外部空間上是調(diào)和函數(shù)。參見圖3.3-5，假設(shè)球外質(zhì)點的球坐標(biāo)為，球內(nèi)質(zhì)量微元的球坐標(biāo)為，是質(zhì)點在所在球面上的投影，記向徑與之間的夾角為。圖3.3-5引力在極坐標(biāo)下的表示圖3.3-6球面三角在由三點構(gòu)成的三角形中，根據(jù)余弦定理，有（3.3-66）其中，記和。上式兩邊同時開方再倒數(shù)，得（3.3-67）當(dāng)時，可將展開成關(guān)于的級數(shù)形式，即（3.3-68）對比上式級數(shù)的系數(shù)與勒讓德多項式表3.3-1中的第一列數(shù)據(jù)，發(fā)現(xiàn)上式系數(shù)恰好等于各階勒讓德多項式，因而有（3.3-69）實際上，稱為勒讓德多項式的生成函數(shù)（或母函數(shù)）。將式（3.3-69）代入式（3.3-67），再代入式（3.3-59），有（3.3-70）式中，為地球上距坐標(biāo)原點最遠(yuǎn)的質(zhì)量微元的距離。參見圖3.3-6，根據(jù)球面三角形的余弦定理，有（3.3-71）可見，是四個球面角坐標(biāo)變量的函數(shù)，因而可以用這四個變量的球諧函數(shù)來表示，即（3.3-72）上式稱為球函數(shù)的加法公式（推導(dǎo)過程復(fù)雜，從略）。 在式（3.3-72）中，若暫且將視為固定坐標(biāo)，則可視為球面動坐標(biāo)的函數(shù)，將式（3.3-72）代入式（3.3-70），考慮到積分是針對空間動坐標(biāo)實施的，則有（3.3-73）其中（3.3-74）（3.3-75）記為地球引力常數(shù)。當(dāng)時，式（3.3-74）和（3.3-75）中各階系數(shù)具有明顯的物理意義。根據(jù)直角坐標(biāo)與球坐標(biāo)之間的轉(zhuǎn)換關(guān)系（3.3-4），重寫如下（3.3-76）將其代入式（3.3-74）和（3.3-75），容易得到以下結(jié)果：（3.3-77a）（3.3-77b）（3.3-77c）（3.3-77d）（3.3-77e）（3.3-77f）（3.3-77g）（3.3-77h）（3.3-77i）其中，式（3.3-77b）～（3.3-77d）中、和表示地球的質(zhì)心坐標(biāo)，若將坐標(biāo)系原點定義在地球質(zhì)心上，則有。另外，剛體的慣性張量定義為（3.3-78）其中，、和分別表示繞軸、軸和軸的轉(zhuǎn)動慣量；而、和為慣量積。比較式（3.3-77e）～（3.3-77i）和式（3.3-78），可求得（3.3-79a）（3.3-79b）（3.3-79c）（3.3-79d）（3.3-79e）若定義直角坐標(biāo)系時，使其坐標(biāo)軸與地球的慣性主軸重合，則有，即在式（3.3-79）中有。對于實際地球，選擇坐標(biāo)系一般優(yōu)先滿足坐標(biāo)原點與地球質(zhì)心重合，軸與地球自轉(zhuǎn)軸平行，軸在赤道面上且指向零度經(jīng)線，這時坐標(biāo)軸與地球慣性主軸一般不重合，因而實際上系數(shù)、和一般均不等于零。（3.3-79a）中反映了赤道與極軸轉(zhuǎn)動慣量的差別；（3.3-79c）中反映了赤道面上兩坐標(biāo)軸間轉(zhuǎn)動慣量的差別。利用三角函數(shù)恒等公式，式（3.3-73）可化為（3.3-80a）其中常稱為地球的動力扁率。在式（3.3-80a）中，是球形地球引起的引力位，約為，其它系數(shù)比小三個數(shù)量級，因而在有些應(yīng)用中，比如人造衛(wèi)星的運動，只需考慮和的影響，式（3.3-80a）可近似為（3.3-80b）根據(jù)牛頓第二運動定律，在地球引力場中衛(wèi)星的運動方程為（3.3-81）其中，為衛(wèi)星在地心慣性坐標(biāo)系中的位移矢量，且模值，為地球?qū)πl(wèi)星的引力。取式（3.3-81）的坐標(biāo)分量，并考慮到，可得（3.3-82a）同理，對于和坐標(biāo)分量有（3.3-82b）（3.3-82c）將式（3.3-82a）～（3.3-82c）整理成更為簡潔的矢量形式，為（3.3-83）其中，，表示上的單位矢量，表示地球自轉(zhuǎn)軸上的單位矢量，是與之間的夾角。上式便是衛(wèi)星在慣性坐標(biāo)系下的運動方程，它是矢量方程，與具體坐標(biāo)系選擇無關(guān)。值得說明的是，由于（3.3-80b）中只考慮了地球高階引力位中的主諧項，這時引力關(guān)于經(jīng)度是旋轉(zhuǎn)對稱的，因而無需考慮地球自轉(zhuǎn)的影響；如果出現(xiàn)引力的田諧項或扇諧項，就必須顧及地球自轉(zhuǎn)的影響了。3.3.3重力位及重力計算運載體在地球表面附近導(dǎo)航過程中，跟隨地球自轉(zhuǎn)一起轉(zhuǎn)動，通常選擇地固坐標(biāo)系作為參考坐標(biāo)系，這時運載體同時受地球引力和離心力作用，兩種力的合力即為重力。地心直角坐標(biāo)與球坐標(biāo)之間的關(guān)系（3.3-4），重寫如下（3.3-84）在球坐標(biāo)表示中，地球自轉(zhuǎn)僅會引起經(jīng)度相對慣性空間隨時間變化，記為（3.3-85）其中，表示初始經(jīng)度，其為常值，是地球自轉(zhuǎn)角速率。將式（3.3-85）代入式（3.3-84），并對時間求一階導(dǎo)，可得（3.3-86）繼續(xù)對上式求導(dǎo)，得（3.3-87）在式（3.3-87）左邊，坐標(biāo)對時間的二階導(dǎo)數(shù)，其物理含義是運動軌跡的向心加速度。 根據(jù)物理學(xué)知識，在引力場中單位質(zhì)點受力等于引力位函數(shù)的梯度，即（3.3-88）類似地，如果給定標(biāo)量函數(shù)（3.3-89）不難求得的梯度函數(shù)，為（3.3-90）比較式（3.3-90）與式（3.3-87），它們在數(shù)值大小上恰好相等但符號相反，可將式（3.3-90）視為單位質(zhì)點所受的慣性離心力，方向垂直于旋轉(zhuǎn)軸向外。因此，式（3.3-89）給出的標(biāo)量函數(shù)可稱為離心力位。本質(zhì)上，離心力不是物質(zhì)力（無外界施力物體），而是在旋轉(zhuǎn)坐標(biāo)系下引入的一種慣性力。容易驗證，離心力位的拉普拉斯運算不為零，即（3.3-91）可見，離心力位函數(shù)不是調(diào)和函數(shù)。 由于重力是引力和離心力之合力，因而重力位就等于引力位（3.3-73）和離心力位（3.3-89）之和，將重力位記為（3.3-92）重力位等于常值的閉合曲面稱為重力等位面，質(zhì)點沿等位面移動，重力不做功。當(dāng)重力位取不同的常值時，就得到一簇曲面，不同曲面之間既不平行也不相交。單位質(zhì)點受到的重力就等于重力位的梯度，即（3.3-93）在球坐標(biāo)系下，梯度的計算公式為（3.3-94）其中，、和分別為沿球的徑向、余緯和經(jīng)度方向的單位矢量，即分別對應(yīng)球面上的天向、南向和東向。將式（3.3-92）和（3.3-94）代入式（3.3-93），可求得地球重力在單位質(zhì)點處“東-北-天”坐標(biāo)系下的三個分量（特別注意，此處的天向是球心指向單位質(zhì)點的方向，而非重力鉛垂向或地理天向）：（3.3-95）在式（3.3-95）的第二式中，根據(jù)勒讓德求導(dǎo)公式（3.3-36c），可得參見圖3.3-7，記地理緯度與地心緯度的差值為，則由式（3.3-95）經(jīng)過旋轉(zhuǎn)變換，可獲得地理系下的重力值，如下（3.3-96）這便是由地球引力位函數(shù)系數(shù)及質(zhì)點的球坐標(biāo)求解地理系下重力的計算公式。圖3.3-7地理系下的重力計算值得注意的是，從式（3.3-38）中可以看出，當(dāng)和都比較大時（比如），勒讓德多項式的模值平方達(dá)到了量級，因此，直接使用勒讓德多項式在數(shù)值計算上存在問題，實際應(yīng)用中多采用正規(guī)化的勒讓德多項式。下面給出相關(guān)算法的正規(guī)化勒讓德多項式形式。首先，根據(jù)式（3.3-36a），不難求得正規(guī)化勒讓德多項式的遞推公式及求導(dǎo)公式，分別如下（3.3-97a）（3.3-97b）（3.3-97c）（3.3-98a）（3.3-98b）其次，將引力位式（3.3-73）改寫成如下形式（3.3-99）其中是正規(guī)化的勒讓德多項式，且有稱為正規(guī)化的球諧函數(shù)，其模方為，即單位半徑球面的面積。因此，重力位式（3.3-92）可表示為（3.3-100）類似于式（3.3-95），從上式可求得（3.3-101）其中至此，總結(jié)給出計算地理坐標(biāo)系下重力的四個步驟，如下：（1）由運載體地理位置計算球坐標(biāo)，這可通過ECEF直角坐標(biāo)作為過渡進行計算，并計算緯度差值，其中為地心緯度；（2）遞推計算勒讓德函數(shù)和；（3）讀入引力位球諧系數(shù)，根據(jù)式（3.3-101）計算；（4）通過式（3.3-96）計算。如果有或，則說明實際計算重力方向（若足夠精確的話可視為真垂線或天文垂線）與基于理想旋轉(zhuǎn)橢球的地理垂線不重合，兩者之間的角度偏差通常稱為垂線偏差，垂線偏差可分解為子午（南北）偏差分量和卯酉（東西）偏差分量，分別記為，（3.3-102）此外，天向重力與采用正常重力公式（3.2-23）計算結(jié)果之間的大小偏差一般稱為重力異常，可記為（3.3-103）在慣導(dǎo)系統(tǒng)中，高度通道本身就不穩(wěn)定，一般不會單獨使用，所以，多數(shù)情況下并不需要過多關(guān)注重力異常對慣導(dǎo)系統(tǒng)的影響；然而，在實際地球重力場中如果垂線偏差較大而又不能精確補償，將會帶來不利影響，它等效于加速度計偏值誤差或者水平失準(zhǔn)角誤差，因此，垂線偏差是影響慣導(dǎo)精度的一個重要因素。實際地球的垂線偏差約為數(shù)角秒，局部超過30角秒，最大甚至達(dá)100角秒，一般在地形起伏地區(qū)偏差較大，而在地形平坦的地方相對較小些。實際地球重力場比較復(fù)雜，據(jù)統(tǒng)計，如果僅僅采用正