-CAL-FENGHAI.-(YICAI)-CompanyOne1

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

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

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

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