捷联惯导系统学习3.3(引力位函数)

参数说明：
M:地球总质量M:地球总质量M:地球总质量
dm:地球上的质量微元dm:地球上的质量微元dm:地球上的质量微元
m′:地球外部空间上的一个质点m':地球外部空间上的一个质点m′:地球外部空间上的一个质点
p:m与dm之间的距离p:m与dm之间的距离p:m与dm之间的距离
万有引力大小为：
f=Gdm∗m′p2f=\frac{Gdm*m'}{p^2}f=p2Gdm∗m′

引力位函数

若外有引力做工(从无穷远移动到半径p处)为:A=∫∞p−Gdm∗mp2=Gdm∗m′p若外有引力做工(从无穷远移动到半径p处)为:A=\int_\infty^p-\frac{Gdm*m}{p^2}=\frac{Gdm*m'}{p}若外有引力做工(从无穷远移动到半径p处)为:A=∫∞p−p2Gdm∗m=pGdm∗m′
根据能量守恒，引力做工必然等于m′m'm′减少的位能（势能），若单位质量的点(从无穷远移动到半径p处)，势能必然减少，将减少量定义为位函数:
dV=GdmpdV=\frac{Gdm}{p}dV=pGdm

地球外部空间的调和函数

地球产生的各质量微元函数之和为：
V=∫MdV=G∫Mdmp（公式1）V=\int_MdV=G\int_M\frac{dm}{p}（公式1）V=∫MdV=G∫Mpdm（公式1）
又因为：
▽V=∂2V∂x2+∂2V∂y2+∂2V∂z2=0(证略)\bigtriangledown V=\frac{\partial^2V}{\partial x^2}+\frac{\partial^2V}{\partial y^2}+\frac{\partial^2V}{\partial z^2}=0(证略)▽V=∂x2∂2V+∂y2∂2V+∂z2∂2V=0(证略)
所以VVV是在地球质量M的外部空间上是调和函数

勒让德多项式的生成函数

参数说明：
m′：的球坐标为(r,θ,λ)m'：的球坐标为(r,\theta,\lambda)m′：的球坐标为(r,θ,λ)
dm：的球坐标为(r,θ′,λ′)dm：的球坐标为(r,\theta',\lambda')dm：的球坐标为(r,θ′,λ′)
向量R与r夹角为φ向量R与r夹角为\varphi向量R与r夹角为φ
m0′是m′在球面上的投影m0'是m'在球面上的投影m0′是m′在球面上的投影
在,o,dm,m′o,dm,m'o,dm,m′构成的三角形中，根据余弦定理,可以得到：
p2=r2+R2−2Rrcosφ=r2(1−2ax+a2)(a=Rr,x=cosφ)p^2=r^2+R^2-2Rrcos\varphi=r^2(1-2ax+a^2)(a=\frac{R}{r},x=cos\varphi)p2=r2+R2−2Rrcosφ=r2(1−2ax+a2)(a=rR,x=cosφ)
再开平方，取倒数得到：
1p=1r(1−2ax+a2)−12（公式2）\frac{1}{p}=\frac{1}{r}(1-2ax+a^2)^{-\frac{1}{2}}（公式2）p1=r1(1−2ax+a2)−21（公式2）
又因为：
(1−2ax+a2)−12=∑n=0∞Pn(x)an(勒让德多项式的生成函数)（公式3）(1-2ax+a^2)^{-\frac{1}{2}}=\sum_{n=0}^\infty P_n(x)a^n(勒让德多项式的生成函数)（公式3）(1−2ax+a2)−21=n=0∑∞Pn(x)an(勒让德多项式的生成函数)（公式3）
所以(1−2ax+a2)−12(1-2ax+a^2)^{-\frac{1}{2}}(1−2ax+a2)−21为勒让德多项式的生成函数

球函数的加法公式

根据球面的三角余弦定理有:cosφ=cosθcosθ′+sinθsinθ′cos(λ−λ′)cos\varphi=cos\theta cos\theta'+sin\theta sin\theta'cos(\lambda-\lambda')cosφ=cosθcosθ′+sinθsinθ′cos(λ−λ′)
带入Pn(cosφ)=∑k=0n2(n−k)!(1+δk)(n+k)!Pnk(cosθ)Pnk(cosθ′)coskλ′+Pnkcosθsin(kλ)Pnk(cosθ′)(sinkλ′)(公式4)P_n(cos\varphi)=\sum_{k=0}^n\frac{2(n-k)!}{(1+\delta_k)(n+k)!}P_n^k(cos\theta)P_n^k(cos\theta')cosk\lambda'+P_n^kcos\theta sin(k\lambda) P_n^k(cos\theta')(sink\lambda') (公式4)Pn(cosφ)=k=0∑n(1+δk)(n+k)!2(n−k)!Pnk(cosθ)Pnk(cosθ′)coskλ′+Pnkcosθsin(kλ)Pnk(cosθ′)(sinkλ′)(公式4)

利用德让勒多项式求解V（调和函数）

将勒让德多项式的生成函数（公式(3)），将带入公式（2），再将公式2带入地球外部空间的调和函数公式（1）。
得到:
Re:旋转椭圆球的长半径R_e:旋转椭圆球的长半径Re:旋转椭圆球的长半径
u=GM:地球引力常数u=GM:地球引力常数u=GM:地球引力常数
δk={1k=00k≠1\delta_k= \begin{cases} 1&k=0\\ 0&k\neq1\\ \end{cases}δk={10k=0k=1
V=Gr∑n=0∞(Rer)n∫M(RRe)nPn(cosφ)dmV=\frac{G}{r}\sum_{n=0}^\infty(\frac{R_e}{r})^n\int_M(\frac{R}{R_e})^nP_n(cos\varphi)dmV=rGn=0∑∞(rRe)n∫M(ReR)nPn(cosφ)dm
将公式（4）带入:
得到:
V=ur∑n=0∞(Rer)n∑k=0n(Cnkcoskλ+Snksinkλ)Pnk(cosθ)V=\frac{u}{r}\sum_{n=0}^\infty(\frac{R_e}{r})^n\sum_{k=0}^n(C_n^kcosk\lambda+S_n^ksink\lambda)P_n^k(cos\theta)V=run=0∑∞(rRe)nk=0∑n(Cnkcoskλ+Snksinkλ)Pnk(cosθ)
Cnk=2(n−k)!M(1+δk)(n+k)!∫M(RRe)nPnk(cosθ′)coskλ′dmC_n^k=\frac{2(n-k)!}{M(1+\delta_k)(n+k)!}\int_M(\frac{R}{R_e})^nP_n^k(cos\theta')cosk\lambda'dmCnk=M(1+δk)(n+k)!2(n−k)!∫M(ReR)nPnk(cosθ′)coskλ′dm
Snk=2(n−k)!M(1+δk)(n+k)!∫M(RRe)nPnk(cosθ′)sinkλ′dmS_n^k=\frac{2(n-k)!}{M(1+\delta_k)(n+k)!}\int_M(\frac{R}{R_e})^nP_n^k(cos\theta')sink\lambda'dmSnk=M(1+δk)(n+k)!2(n−k)!∫M(ReR)nPnk(cosθ′)sinkλ′dm
已知直角坐标系到球体坐标系的转换关系如下：
{x=Rsinθ′cosλ′y=Rsinθ′sinλ′z=Rcosθ′\begin{cases} x=Rsin\theta'cos\lambda'\\ y=Rsin\theta'sin\lambda'\\ z=Rcos\theta'\\ \end{cases}⎩⎪⎨⎪⎧x=Rsinθ′cosλ′y=Rsinθ′sinλ′z=Rcosθ′
刚体的张量定义如下:刚体的惯性张量及其物理意义
I=[Ix−Ixy−Ixz−IxyIy−Iyz−Ixz−IyzIxz]=∫M[y2+z2−xy−xz−xyx2+z2−yz−xz−yzx2+y2]dmI=\left[\begin{matrix} I_{x}&-I_{xy}&-I_{xz}\\ -I_{xy}&I_{y}&-I_{yz}\\ -I_{xz}&-I_{yz}&I_{xz}\\ \end{matrix}\right]=\int_M\left[\begin{matrix} y^2+z^2&-xy&-xz\\ -xy&x^2+z^2&-yz\\ -xz&-yz&x^2+y^2\\ \end{matrix}\right]dmI=⎣⎡Ix−Ixy−Ixz−IxyIy−Iyz−Ixz−IyzIxz⎦⎤=∫M⎣⎡y2+z2−xy−xz−xyx2+z2−yz−xz−yzx2+y2⎦⎤dm
可以得到：
C00=1C_0^0=1C00=1
C01=1Re(1M∫Mzdm)C_0^1=\frac{1}{R_e}(\frac{1}{M}\int_Mzdm)C01=Re1(M1∫Mzdm)
C11=1Re(1M∫Mxdm)C_1^1=\frac{1}{R_e}(\frac{1}{M}\int_Mxdm)C11=Re1(M1∫Mxdm)
S11=1Re(1M∫Mydm)S_1^1=\frac{1}{R_e}(\frac{1}{M}\int_Mydm)S11=Re1(M1∫Mydm)
C20=1MRe2(∫Mz2−x2+y22dm)=1MRe2(Ix+Iy2−Iz)C_2^0=\frac{1}{MR_e^2}(\int_Mz^2-\frac{x^2+y^2}{2}dm)=\frac{1}{MR_e^2}(\frac{I_x+I_y}{2}-I_z)C20=MRe21(∫Mz2−2x2+y2dm)=MRe21(2Ix+Iy−Iz)
C21=1MRe2∫Mxzdm=IxzMRe2C_2^1=\frac{1}{MR_e^2}\int_Mxzdm=\frac{I_{xz}}{MR_e^2}C21=MRe21∫Mxzdm=MRe2Ixz
C22=14MRe2(∫Mx2−y2dm)=Iy−Ix4MRe2C_2^2=\frac{1}{4MR_e^2}(\int_Mx^2-y^2dm)=\frac{I_{y}-I_{x}}{4MR_e^2}C22=4MRe21(∫Mx2−y2dm)=4MRe2Iy−Ix
S21=1MRe2(∫Myzdm)=IyzMRe2S_2^1=\frac{1}{MR_e^2}(\int_Myzdm)=\frac{I_{yz}}{MR_e^2}S21=MRe21(∫Myzdm)=MRe2Iyz
S22=14MRe2(∫Mxydm)=Ixy2MRe2S_2^2=\frac{1}{4MR_e^2}(\int_Mxydm)=\frac{I_{xy}}{2MR_e^2}S22=4MRe21(∫Mxydm)=2MRe2Ixy
如果定义直角坐标与地球惯性惯性主轴重合，则：Ixy=Iyz=Ixz=0I_{xy}=I_{yz}=I_{xz}=0Ixy=Iyz=Ixz=0

对于实际应用时，地球坐标一边选择坐标原点与地球质心重合，oz与地球自转平行，ox轴在赤道面上且指向0度经线，这时坐标轴往往与惯性主轴不重合。
利用三角函数恒等式：
得到:
V=ur{1−∑n=1∞(Rer)n[JnPn(cosθ)+∑k=1nJnkPnk(cosθ)cosk(λ+λnk)]}V=\frac{u}{r}\{1-\sum_{n=1}^\infty(\frac{R_e}{r})^n[J_nP_n(cos\theta)+\sum_{k=1}^nJ_n^kP_n^k(cos\theta)cosk(\lambda+\lambda_n^k)]\}V=ru{1−n=1∑∞(rRe)n[JnPn(cosθ)+k=1∑nJnkPnk(cosθ)cosk(λ+λnk)]}
Jn=−Cn0,Jnk=(Cnk)2+(Snk)2,λnk=−arctan(SnkCnk)/kJ_n=-C_n^0,J_n^k=\sqrt{(C_n^k)^2+(S_n^k)^2},\lambda_n^k=-arctan(\frac{S_n^k}{C_n^k})/kJn=−Cn0,Jnk=(Cnk)2+(Snk)2,λnk=−arctan(CnkSnk)/k
J2=C20=1MRe2(Ix+Iy2−Iz)(动力扁率：反映赤道与极轴上转动惯量的差别)J_2=C_2^0=\frac{1}{MR_e^2}(\frac{I_x+I_y}{2}-I_z)(动力扁率：反映赤道与极轴上转动惯量的差别)J2=C20=MRe21(2Ix+Iy−Iz)(动力扁率：反映赤道与极轴上转动惯量的差别)
ur:球形地球引起的引力位\frac{u}{r}:球形地球引起的引力位ru:球形地球引起的引力位
在实际应用中因为其他系数比J2J_2J2小三个数量级，所以只需考虑ur和J2\frac{u}{r}和J_2ru和J2的影响，可以将VVV近似看为：
V=ur[1−J2Re22r2(3cos2θ−1)]（只考虑主谐项）V=\frac{u}{r}[1-\frac{J_2R_e^2}{2r^2}(3cos^2\theta-1)]（只考虑主谐项）V=ru[1−2r2J2Re2(3cos2θ−1)]（只考虑主谐项）

卫星在惯性坐标系下的运动方程（只考虑主谐项）

已知
位移矢量：r=[xyz]T位移矢量：r=\left[\begin{matrix} x&y&z\\ \end{matrix}\right]^T位移矢量：r=[xyz]T
引力：f=[fxfyfz]T引力：f=\left[\begin{matrix} f_x&f_y&f_z\\ \end{matrix}\right]^T引力：f=[fxfyfz]T
cosθ=zrcos\theta=\frac{z}{r}cosθ=rz
对x,y,z求2阶片导，求得加速度方程为：
r¨=−ur3[I+32J2(Rer)2(D−5(ur∗up)I)]r\ddot r=-\frac{u}{r^3}[I+\frac{3}{2}J_2(\frac{R_e}{r})^2(D-5(u_r*u_p)I)]rr¨=−r3u[I+23J2(rRe)2(D−5(ur∗up)I)]r
其中:D=diag(113)D=diag\left(\begin{matrix} 1&1&3\\ \end{matrix}\right)D=diag(113)
ur表示r上的单位矢量u_r表示r上的单位矢量ur表示r上的单位矢量
up表示自转轴上的的单位矢量u_p表示自转轴上的的单位矢量up表示自转轴上的的单位矢量
ur∗up表示up与ur之间夹角的余弦值u_r*u_p 表示u_p与u_r之间夹角的余弦值ur∗up表示up与ur之间夹角的余弦值