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

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

引力位函数

$若外有引力做工(从无穷远移动到半径p处)为:A={\int }_{\infty }^p-\frac{Gdm\ast m}{p^2}=\frac{Gdm\ast m^{\prime }}{p}$
根据能量守恒，引力做工必然等于$m^{\prime }$减少的位能（势能），若单位质量的点(从无穷远移动到半径p处)，势能必然减少，将减少量定义为位函数:
$$dV=\frac{Gdm}{p}$$

地球外部空间的调和函数

地球产生的各质量微元函数之和为：
$$V={\int }_{M}dV=G{\int }_M\frac{dm}{p}（公式1）$$
又因为：
$$▽V=\frac{{\partial }^{2}V}{\partial x^2}+\frac{{\partial }^{2}V}{\partial y^2}+\frac{{\partial }^{2}V}{\partial z^2}=0(证略)$$
所以$V$是在地球质量M的外部空间上是调和函数

勒让德多项式的生成函数

参数说明：
$m^{\prime }：的球坐标为(r,\theta ,\lambda )$
$dm：的球坐标为(r,{\theta }^{\prime },{\lambda }^{\prime })$
$向量R与r夹角为\phi$
$m{0}^{\prime }是m^{\prime }在球面上的投影$
在,$o,dm,m^{\prime }$构成的三角形中，根据余弦定理,可以得到：
$$p^2=r^2+R^2-2Rrcos\phi =r^2(1-2ax+a^2)(a=\frac{R}{r},x=cos\phi )$$
再开平方，取倒数得到：
$$\frac{1}{p}=\frac{1}{r}(1-2ax+a^2{)}^{-\frac{1}{2}}（公式2）$$
又因为：
$$(1-2ax+a^2{)}^{-\frac{1}{2}}=\sum _{n=0}^{\infty }{P}_n(x)a^n(勒让德多项式的生成函数)（公式3）$$
所以$(1-2ax+a^2{)}^{-\frac{1}{2}}$为勒让德多项式的生成函数

球函数的加法公式

根据球面的三角余弦定理有:$cos\phi =cos\theta cos{\theta }^{\prime }+sin\theta sin{\theta }^{\prime }cos(\lambda -{\lambda }^{\prime })$
带入$$P_n(cos\phi )=\sum _{k=0}^n\frac{2(n-k)!}{(1+{\delta }_k)(n+k)!}{P}_n^k(cos\theta )P_n^k(cos{\theta }^{\prime })cosk{\lambda }^{\prime }+P_n^{k}cos\theta sin(k\lambda )P_n^k(cos{\theta }^{\prime })(sink{\lambda }^{\prime })(公式4)$$

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

将勒让德多项式的生成函数（公式(3)），将带入公式（2），再将公式2带入地球外部空间的调和函数公式（1）。
得到:
$R_e:旋转椭圆球的长半径$
$u=GM:地球引力常数$
${\delta }_k=\{\begin{array}{ll}1& k=0\\ 0& k\ne 1\end{array}$
$$V=\frac{G}{r}\sum _{n=0}^{\infty }(\frac{R_e}{r}{)}^n{\int }_M(\frac{R}{R_e}{)}^n{P}_n(cos\phi )dm$$
将公式（4）带入:
得到:
$$V=\frac{u}{r}\sum _{n=0}^{\infty }(\frac{R_e}{r}{)}^n\sum _{k=0}^n(C_n^{k}cosk\lambda +S_n^{k}sink\lambda )P_n^k(cos\theta )$$
$$C_n^k=\frac{2(n-k)!}{M(1+{\delta }_k)(n+k)!}{\int }_M(\frac{R}{R_e}{)}^n{P}_n^k(cos{\theta }^{\prime })cosk{\lambda }^{\prime }dm$$
$$S_n^k=\frac{2(n-k)!}{M(1+{\delta }_k)(n+k)!}{\int }_M(\frac{R}{R_e}{)}^n{P}_n^k(cos{\theta }^{\prime })sink{\lambda }^{\prime }dm$$
已知直角坐标系到球体坐标系的转换关系如下：
$$\{\begin{array}{l}x=Rsin{\theta }^{\prime }cos{\lambda }^{\prime }\\ y=Rsin{\theta }^{\prime }sin{\lambda }^{\prime }\\ z=Rcos{\theta }^{\prime }\end{array}$$
刚体的张量定义如下:刚体的惯性张量及其物理意义
$$I=\left[\begin{array}{ccc}{I}_x& -I_{xy}& -I_{xz}\\ -I_{xy}& I_y& -I_{yz}\\ -I_{xz}& -I_{yz}& I_{xz}\end{array}\right]={\int }_M\left[\begin{array}{ccc}{y}^2+z^2& -xy& -xz\\ -xy& x^2+z^2& -yz\\ -xz& -yz& x^2+y^2\end{array}\right]dm$$
可以得到：
$C_0^0=1$
$C_0^1=\frac{1}{R_e}(\frac{1}{M}{\int }_{M}zdm)$
$C_1^1=\frac{1}{R_e}(\frac{1}{M}{\int }_{M}xdm)$
$S_1^1=\frac{1}{R_e}(\frac{1}{M}{\int }_{M}ydm)$
$C_2^0=\frac{1}{M{R}_e^2}({\int }_M{z}^2-\frac{x^2+y^2}{2}dm)=\frac{1}{M{R}_e^2}(\frac{I_x+I_y}{2}-I_z)$
$C_2^1=\frac{1}{M{R}_e^2}{\int }_{M}xzdm=\frac{I_{xz}}{M{R}_e^2}$
$C_2^2=\frac{1}{4M{R}_e^2}({\int }_M{x}^2-y^{2}dm)=\frac{I_y-I_x}{4M{R}_e^2}$
$S_2^1=\frac{1}{M{R}_e^2}({\int }_{M}yzdm)=\frac{I_{yz}}{M{R}_e^2}$
$S_2^2=\frac{1}{4M{R}_e^2}({\int }_{M}xydm)=\frac{I_{xy}}{2M{R}_e^2}$
如果定义直角坐标与地球惯性惯性主轴重合，则：$I_{xy}=I_{yz}=I_{xz}=0$

对于实际应用时，地球坐标一边选择坐标原点与地球质心重合，oz与地球自转平行，ox轴在赤道面上且指向0度经线，这时坐标轴往往与惯性主轴不重合。
利用三角函数恒等式：
得到:
V = u r { 1 − ∑ n = 1 ∞ ( R e r ) n [ J n P n ( c o s θ ) + ∑ k = 1 n J n k P n k ( c o s θ ) c o s k ( λ + λ n k ) ] } 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[Jn​Pn​(cosθ)+k=1∑n​Jnk​Pnk​(cosθ)cosk(λ+λnk​)]}
$$J_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})/k$$
$J_2=C_2^0=\frac{1}{M{R}_e^2}(\frac{I_x+I_y}{2}-I_z)(动力扁率：反映赤道与极轴上转动惯量的差别)$
$\frac{u}{r}:球形地球引起的引力位$
在实际应用中因为其他系数比$J_2$小三个数量级，所以只需考虑$\frac{u}{r}和J_2$的影响，可以将$V$近似看为：
$$V=\frac{u}{r}[1-\frac{J_2{R}_e^2}{2r^2}(3co{s}^2\theta -1)]（只考虑主谐项）$$

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

已知
$位移矢量：r={\left[\begin{array}{ccc}x& y& z\end{array}\right]}^T$
$引力：f={\left[\begin{array}{ccc}{f}_x& f_y& f_z\end{array}\right]}^T$
$cos\theta =\frac{z}{r}$
对x,y,z求2阶片导，求得加速度方程为：
$$\ddot{r}=-\frac{u}{r^3}[I+\frac{3}{2}{J}_2(\frac{R_e}{r}{)}^2(D-5(u_r\ast u_p)I)]r$$
其中:$D=diag\left(\begin{array}{ccc}1& 1& 3\end{array}\right)$
$u_r表示r上的单位矢量$
$u_p表示自转轴上的的单位矢量$
$u_r\ast u_p表示u_p与u_r之间夹角的余弦值$