地磁场仿真与导航方案设计

地磁场仿真与导航方案设计

地磁场作为地球的固有资源，为航空、航天、航海提供了天然的坐标，可应用于航天器、航空飞行器和舰船等载体的定位、定向及姿态控制。本文依据地磁场的球谐模型，利用Matlab工具，设计了计算任意地点地磁场的程序，对仿真结果与实际数据进行对比分析；并选定一组卫星的数据，结合卫星的运动轨道，描绘90分钟内卫星处地磁总强度/矢量随轨道变化曲线，据此设计一种卫星地磁导航方案。

一、地磁场计算原理

1.1.球谐模型与球谐系数
地磁学的基本问题之一就是用数学表达式将地磁要素的地面分布表示成地理坐标的函数。1839 年 Gauss 提出了地磁场的球谐分析方法，1885 年 Schmidt 又发展了这一方法，形成了地磁学的高斯理论。国际地磁场随之建立起来，从而就有了全球统一的地磁场数学模型—球谐模型。
使用球谐模型计算地磁场，球谐系数必不可少。球谐系数是一个随时间变化的量。由于地球主磁场的源——高温液态铁镍环流处于不断的运动过程中，因此不同年份有不同的球谐系数。国际地磁和高空物理协会（IAGA）每隔五年提供一组球谐系数数据，计算特定年份特定时间的球谐系数，依据以下公式：

式中，$g_n^m(year)$和$h_n^m(year)$是year年的球谐系数（nT），$S{v}_{g_n^m}$和$S{v}_{h_n^m}$是球谐系数的年变率（nT/年），$g_n^m$和$h_n^m$是t年的球谐系数（nT）。
把所求年份t、对照年份year及该年份的球谐系数和年变率数据代入上式，就能得到任意年份的球谐系数。

1.2.地理坐标转地心坐标
由于IGRF的球谐系数是在半径为6371.2km的参考球中推导出来的，实际上，地球并非是理想的球体，而是椭球体。为了提高地磁计算的精度，需要考虑地球的扁率。假设地球某一测点的大地坐标为(ϕ,λ,h)，其中ϕ为大地(地理)纬度、λ为经度、h为 海拔高度；则该点的地心坐标为(φ,λ,r)，其中φ为地心纬度、λ为地心经度、r为测点离开地心的距离。则大地坐标和地心球坐标有如下关系：

上式中，a为地球长半轴(6378.137km)，b 为地球短半轴(6356.7523km)。

1.3.勒让德函数及其一阶导数求解
以矩阵$P_{(}n,m)(cos\theta )$表示勒让德函数，$\frac{d{P}_{n,m}(cos\theta )}{d\theta }$表示其一阶导数（（$cos\theta$）说明该函数以地心余纬θ的余弦值为自变量，$\theta =90°-\phi$），则

1.4.计算地磁场分量
使用球谐系数，并借助于高斯级数表达式，就可以计算出地心球坐标系下的地磁三分量 X ,Y 和Z 。IGRF三分量 X ,Y 和Z 的高斯级数表达式：

式中，R为地球半径(6371.2km)，r 为地心到计算点的径向距离(km)，θ为地心余纬(即θ=90°-φ)，φ为地心纬度，λ为地心或地理经度，$P_n^{m}cos(\theta )$为Schmidt-Legendre函数，$g_n^m$和$h_n^m$为球谐系数(nT)。

1.5.地心坐标转地理（大地）坐标
为了描述地磁向量场的空间分布特征，经常将地磁场F 在北━东━地地理坐标系中表示为7个地磁要素(F,H,X_大地，Y_大地，Z_大地，D,I)。下图给出了地理坐标系中各要素的定义和符号，其中XOY面为水平面，OZ为向下的铅锤方向，HOZ为当地磁子午面，XOZ为当地地理子午面，F为地磁场总强度，H为地磁场水平强度，X_大地为地磁场北向分量，Y_大地为地磁场东向分量，Z_大地为地磁场垂直分量，D为磁偏角，I为磁倾角。

地理坐标系中地磁各要素示意图
地心坐标系中的地磁分量转换为大地坐标系中的地磁分量：
$X_{大地}=Xcos\sigma +Zsin\sigma$
$Y_{大地}=Y$
$Z_{大地}=-Xsin\sigma +Zcos\sigma$

1.6.卫星处地磁场矢量的描述
为了更方便地描述卫星处地磁总强度/矢量的变化，画出地磁随卫星轨道的三维矢量分布，需把地磁场从大地坐标系向地心地固坐标系做进一步转化。
可以借助之前的实验结果得到卫星在地心地固坐标系中的位置，记作$(X_D，Y_D，Z_D)$; 已知某时刻卫星运行点的地心经度λ，地心纬度φ，大地坐标系中的三地磁分量$X_{大地}、Y_{大地}、Z_{大地}$，设地心地固坐标系下卫星处地磁分量的矢量表示为$\overrightarrow{X_D}、\overrightarrow{Y_D}、\overrightarrow{Z_D}$，根据大地坐标系下各分量的定义：
$$\overrightarrow{X_D}=-Z_{大地}cos\phi cos\lambda -X_{大地}sin\lambda -Y_{大地}sin\phi cos\lambda$$
$$\overrightarrow{Y_D}=-Z_{大地}cos\phi sin\lambda +X_{大地}cos\lambda -Y_{大地}sin\phi sin\lambda$$
$$\overrightarrow{Z_D}=-Z_{大地}sin\phi +Y_{大地}cos\phi$$

$\overrightarrow{X_D}、\overrightarrow{Y_D}、\overrightarrow{Z_D}$若取正值，说明它们的方向与地心地固坐标系三坐标轴正方向相同，反之则相反。
值得注意的是，由于$Z_{大地}$方向铅垂向下，故以上表达式中所有关于$Z_{大地}$的项都存在负号。
这样，卫星处地磁强度这一三维矢量，就能够由($(X_D，Y_D，Z_D，\overrightarrow{X_D}、\overrightarrow{Y_D}、\overrightarrow{Z_D})$完全确定。

二、地磁场Matlab仿真

2.1.地磁场参数的计算
根据我们的目的，需要实现的函数功能：

地磁场参数计算流程图
不妨记函数为main，它定义如下：

function [F,H,Xdadi,Ydadi,Zdadi,D,I]=main(longtitude,latitude,height,g,h)

输入参数longtitude,latitude,height分别表示地理经度、纬度、海拔高度，g,h是高斯球谐系数。由于该函数在后续实验中会被循环调用，最好把读取高斯球谐系数的语句写在函数之外，否则将产生很多不必要的运行时间。输出参数F,H,Xdadi,Ydadi,Zdadi,D,I分别表示大地坐标系下地磁场的总强度，水平强度，北向分量，东向分量，垂直分量，磁偏角，磁倾角。
2.1.1大地坐标和地心球坐标的坐标转换
根据1.2节中的公式，这一部分的程序可以表示为：

lamda=longtitude*pi/180;Fai=latitude*pi/180;h0=height;
a=6378.137;b=6356.7523;
N=a^2/(sqrt(a^2*cos(Fai)^2+b^2*sin(Fai)^2));
r=sqrt(h0*(h0+2*sqrt(a^2*cos(Fai)^2+b^2*sin(Fai)^2))+(a^4*cos(Fai)^2+b^4*sin(Fai)^2)/(a^2*cos(Fai)^2+b^2*sin(Fai)^2));
xigema=acos((h0+sqrt(a^2*cos(Fai)^2+b^2*sin(Fai)^2))/r);
fai=Fai-xigema;

细节：MATLAB中三角函数以弧度制进行运算，由于输入的经、纬度以角度值表示，故应先转成弧度（×π/180），再代入运算。
运算得到的变量，r为测点离开地心的距离，fai即φ地心纬度，lamda即λ地心纬度。
2.1.2勒让德函数计算
准备工作：

theta=pi/2-fai;
P=zeros(14,14);
dP=zeros(14,14);

theta即地心余纬θ=90°-φ；P为Schmidt-Legendre函数，dP为它的一阶导数。根据地磁分量的计算表达式，P、dP均为与球谐系数g_n^m、h_nm同大小的矩阵。由于给出的球谐系数数据中g、h是14×14大小的矩阵，因此令P、dP大小也为14×14。
根据3.5节的计算公式，P、dP的计算过程如下：

P(1,1)=1;P(1,2)=0;P(2,1)=cos(theta);P(2,2)=sin(theta);
for i=3:14
    for j=1:i
        n=i-1;m=j-1;
        if i~=j
            P(i,j)=(2*n-1)/sqrt(n^2-m^2)*cos(theta)*P(i-1,j)-sqrt(((n-1)^2-m^2)/(n^2-m^2))*P(i-2,j);
        else
            P(i,j)=(2*n-1)/sqrt((n+m)*(n+m-1))*sin(theta)*P(i-1,j-1);
        end
    end
end
dP(1,1)=0;
for i=2:14
    for j=1:i
        n=i-1;m=j-1;
        dP(i,j)=n*cos(theta)/sin(theta)*P(i,j)-sqrt(n^2-m^2)/sin(theta)*P(i-1,j);
    end
end

大体思路是先按照要求对矩阵赋初值，然后按照矩阵元素间的关系进行迭代。
须注意的是，P与dP都近似于下三角矩阵，按列进行迭代时应注意控制迭代的区间：列数≤行数，其余矩阵元素按0处理即可。
2.1.3 计算IGRF三分量
有了勒让德函数和球谐系数之后，就可以按照3.6节中的公式计算IGRF三分量：

X=0;Y=0;Z=0;
R=6371.2;
for i=2:14
    for j=1:i
        n=i-1;m=j-1;
      X=X+(R/r)^(n+2)*(g(i,j)*cos(m*lamda)+h(i,j)*sin(m*lamda))*dP(i,j);
        Y=Y+(R/r)^(n+2)*(g(i,j)*sin(m*lamda)-h(i,j)*cos(m*lamda))*m*P(i,j)/sin(theta);
        Z=Z-(n+1)*(R/r)^(n+2)*(g(i,j)*cos(m*lamda)+h(i,j)*sin(m*lamda))*P(i,j);
    end
end

X、Y、Z即为地心球坐标系下的地磁三分量，R为为地球平均半径(6371.2km)。由于是累加求和的形式，要在循环开始前对X、Y、Z赋初值0。
2.1.4 地心坐标转地理坐标
由于函数要求输出的是地理坐标下的参量，因此根据3.6节，进行这一步转换：

Xdadi=X*cos(xigema)+Z*sin(xigema);
Ydadi=Y;
Zdadi=-X*sin(xigema)+Z*cos(xigema);
F=sqrt(Xdadi^2+Ydadi^2+Zdadi^2);
H=sqrt(Xdadi^2+Ydadi^2);
D=acos(Xdadi/H)*180/pi;
I=acos(H/F)*180/pi;

磁偏角D与磁倾角I的返回值为角度，要进行弧度-角度转换（×180/π）。

2.2igrfmagm函数的使用
igrfmagm函数，可以认为是main函数功能的略缩版（输出参数较少，高度存在限制），在本实验中作为标准答案和对照组，同main的输出结果进行比较。它的使用方法：

XYZ=igrfmagm(height,latitude,longitude,decimalYear,generation)

Igrfmagm的输入参数：height=海拔高度，以米为单位；latitude=纬度；longitude=经度；decimalYear=回归年，输入待求年份；generation表示其使用的IGRF数据的版本，本实验中该版本为12，令它等于12即可。输出参数XYZ为3×1大小的向量，同main函数的输出Xdadi,Ydadi,Zdadi一致，分别为大地坐标系下地磁场的北向分量，东向分量，垂直分量。
Igrfmagm对输入参数height是有限制的，它的取值在-1000和600000 之间，意味着Igrfmagm不能给出地底1000m之下、海拔600000m之上的地磁场数据。

2.3绘制地磁总强度随海拔高度变化图
程序如下：

longtitude=108.8;latitude=34.11;
Hmax=36000;
F=zeros(Hmax:1);Real=zeros(Hmax:1);
[g,h]=Get_GM();
for height=1:1:Hmax
    [F(height),~,~,~,~,~,~]=main(longtitude,latitude,height,g,h);
end
height=1:1:Hmax;
plot(height,F,'r')

由上一节可知，Igrfmagm不能提供海拔较高处的地磁场信息，因此采取main函数。main函数的返回值有7个，这里我们只需选取地磁总强度F，并把不同高度的F存入向量中，就能完成图形绘制。其他返回参数由于不需要输出，用“~”代替即可。
在海拔较低处，我们还可以引入Igrfmagm的数据作为对比：

for height=1:1:Hmax
    XYZ=igrfmagm(1000*height,34.11,108.8,decyear(2015,1,1),12);
    Real(height)=sqrt(sum(XYZ.^2));
end
height=1:1:Hmax;
figure(2)
plot(height,Real,'g')

2.4全球地磁总强度二维图绘制
为了尽可能充分地表现全球地磁场分布，我们在如下区间对全球地磁场进行大致的绘图：纬度[-80°,80°]，经度[-180° ,180° ]，间隔0.5°，因此，存放各点地磁总强度数据的矩阵，其大小为321×721：

F=zeros(321,721);

以0.5为步长改变经纬度，利用main函数迭代：

for latitude=-80:0.5:80
    for longtitude=-180:0.5:180
    [F(latitude*2+161,longtitude*2+361),~,~,~,~,~,~]=main(longtitude,latitude,height,g,h);
    end
end

若能够进一步将地磁场分布画成三维曲面图，不仅更为直观地反映了地磁强度大小的相对值，而且将其旋转后，本身就包含了二维图的信息因此，在这里不妨使用meshgrid与mesh函数画出三维曲面图：

latitude=-80:0.5:80;longtitude=-180:0.5:180;
[LO,LA]=meshgrid(longtitude,latitude);
figure(1)
mesh(LO,LA,F);
xlabel('经度');
ylabel('纬度');
zlabel('地磁场强度');

对于等强度地磁图，使用contour函数，把经纬度和F矩阵作为参数输入即可，方法同画地理等高线：

figure(2)
[C,H]=contour(longtitude,latitude,F);
clabel(C,H);
xlabel('经度');
ylabel('纬度');

2.5地磁总强度/矢量的卫星轨道图绘制
根据卫星的TLE数据，我们容易求出卫星的地心地固坐标和地理坐标（经纬度和高度），以及它们随时间的变化。这一步在卫星坐标已求出的基础上进行。
90分钟内地磁总强度随轨道变化曲线不难绘制，只需把每一时刻卫星的经纬度、高度输入main函数：

[F(dt),~,~,~,~,~,~]=main(Longitude,Latitude,H/1000,g,h);

在循环结束后对向量F进行plot作图：

dt=1:1:5400;
plot(dt,F)
xlabel('运行时间');
ylabel('地磁总强度');

90分钟内地磁矢量随轨道变化曲线的绘制更为复杂。在1.6节中，我们已经提到要完全描述该矢量需要在地心地固坐标系下的六个变量：($(X_D，Y_D，Z_D，\overrightarrow{X_D}、\overrightarrow{Y_D}、\overrightarrow{Z_D})$已知，所以关键是利用1。6节中的坐标转换求解$\overrightarrow{X_D}、\overrightarrow{Y_D}、\overrightarrow{Z_D}$。首先得到大地坐标系下卫星处的磁场强度：

[~,~,X(dt/DT),Y(dt/DT),Z(dt/DT),~,~]=main(Longitude,Latitude,H/1000,g,h);

将其代入1.6节中的坐标转换式：

Longitude=Longitude*pi/180;Latitude=Latitude*pi/180;
XX(dt/DT)=-Z(dt/DT)*cos(Latitude)*cos(Longitude)-X(dt/DT)*sin(Longitude)-Y(dt/DT)*sin(Latitude)*cos(Longitude);
YY(dt/DT)=-Z(dt/DT)*cos(Latitude)*sin(Longitude)+X(dt/DT)*cos(Longitude)-Y(dt/DT)*sin(Latitude)*sin(Longitude);
ZZ(dt/DT)=-Z(dt/DT)*sin(Latitude)+Y(dt/DT)*cos(Latitude);

XX、YY、ZZ就表示矢量$\overrightarrow{X_D}、\overrightarrow{Y_D}、\overrightarrow{Z_D}$。迭代完成后，使用quiver3函数画出三维矢量场：

quiver3(XD,YD,ZD,XX,YY,ZZ,'b')

地磁矢量作图就此完成。

三、仿真结果、分析与地磁导航方案设计

3.1地磁场参数计算的分析与比较
依据main函数的计算结果与依据参考程序igrfmagm的计算结果对比如下表所示（M表示使用main函数计算出的结果，I表示使用igrfmagm函数计算出的结果）：

经度/°	纬度/°	海拔/km	地磁场北向分量		地磁场东向分量		地磁场垂直分量		地磁场总强度
			M	I	M	I	M	I	M	I
-171	-76	0	-3199	-2855	11138	11139	-60556	-60615	61654	61696
-153	-68	30	6418	6875	12698	12669	-56760	-56670	58516	58474
-135	-60	60	13572	14014	12312	12258	-48457	-48208	51806	51679
-117	-52	90	17810	18158	11004	10941	-37698	-37358	43121	42954
-99	-44	120	19571	19798	8561	8498	-25987	-25636	33640	33487
-81	-36	150	19477	19619	3967	3910	-15802	-15496	25393	25304
-63	-28	180	18042	18151	-2296	-2346	-10809	-10579	21157	21140
-45	-20	210	16478	16590	-6313	-6345	-11890	-11740	21278	21291
-27	-12	240	16815	16920	-6505	-6511	-14893	-14797	23385	23401
-9	-4	270	20980	21024	-3993	-3989	-15067	-15014	26136	26141
9	4	300	26557	26557	-947	-947	-8483	-8483	27895	27895
27	12	330	29608	29608	1038	1038	3247	3247	29803	29803
45	20	360	29417	29417	1048	1048	15920	15920	33465	33465
63	28	390	27481	27481	724	724	25888	25888	37761	37761
81	36	420	24476	24476	724	724	34386	34386	42214	42214
99	44	450	20209	20209	-401	-401	40960	40960	45676	45676
117	52	480	15899	15899	-2185	-2185	43774	43774	46623	46623
135	60	510	12603	12603	-2514	-2514	44112	44112	45946	45946
153	68	540	9486	9486	-1246	-1246	44449	44449	45468	45468
171	76	570	5630	5630	158	158	44940	44940	45292	45292

分析上表，发现在大部分地理条件下，两个函数计算得到的地磁场完全相同，仅有部分情况会有微小差别。差别在西经117°、南纬52°、海拔90km处达到最大，使用main函数计算出的磁场强度比实际对照值少167nT，且随着海拔高度增加，差别逐渐减小到0。从海拔高度的角度来看，该差别可能是随地面高度增加而衰减的地磁异常场与电离层扰动磁场叠加造成的结果。

3.2地磁总强度随海拔高度变化分析
在本次仿真中，选定地点陕西西安（东经108.8°，北纬34.11°）为分析点，画出的海拔在[0,36000km] 范围内地磁总强度曲线：

西安地磁总强度随海拔变化曲线1
上图曲线显示，随着海拔高度上升地磁总强度近似以反比例函数的形式减少。
将地磁总强度换做对数坐标后，曲线变为：

西安地磁总强度随海拔变化曲线2

西安地磁总强度随海拔变化曲线3
除了海拔较低的一段（可能受初值或地磁异常影响），该曲线基本线性。据此提出推测：海拔较高时，地磁总强度的对数值随√h呈线性关系。

3.3全球地磁总强度二维图绘制与分析
2015.1.1，海拔高度10km处的全球地磁总强度二维分布图：

全球地磁总强度二维分布图

全球等强度地磁图
若把二维分布图画成三维曲面，效果如下：

全球地磁总强度三维曲面图
观察分析以上图片，可以得到结论:

1. 地磁强度与纬度关系密切。总的来说，靠近两极的高纬度地区磁场较强，赤道附近地区磁场较弱
2. 南纬30°、西经50°附近的南大西洋和南美洲地区，相较于同纬度其他地区，磁场异常偏弱。
   南大西洋磁场异常区
   南大西洋异常区”是在1958年为科学家柯策（Pieter Kotze）所发现。科学家进一步研究南大西洋异常区在1590-2005年之间的变化，发现1750年是转折点，在1750年之前，现在的南大西洋异常区的地磁强度几乎没有变化。但1750年之后，地磁强度开始逐渐下降。据研究，全球磁场在1850-1860年开始明显减弱。换句说话，南大西洋异常区内磁力的减弱早于全球磁场减弱100年。许多科学家认为，异常区内的地磁减弱加速，带动了全球磁场开始弱化。
   南大西洋异常区不止磁场减弱，范围更逐渐扩大，而且其中心每年向西移动0.3经度。美国NASA的科学家甚至预测，如果南大西洋异常区以现在的速度继续扩大，到了2240年，将会涵盖半个南半球。
   地球磁场对地球有保护作用，它阻挡了许多来自外太空有害的宇宙射线和带电粒子，像X光之类的辐射线，还有带正电的质子或带负电的电子。南大西洋异常区上空的地球磁力保护变得薄弱之后，让更多来自外太空的辐射线得以穿透，到达大气层更低的地方，更接近地球表面，容易干扰到经过异常区的卫星、飞机和太空船的通信。
   例如，太空船经过”南大西洋异常区“时，太空人虽然闭起眼睛，眼前仍会感觉到闪光。而美国的哈勃太空望远镜绕行地球轨道飞过这里时，会出现仪器失灵、通讯错误，因此必须自动关闭感应器以防受损。这异常区的种种现象表示地球的磁北极和磁南极“即将”翻转，但就地球漫长历史而言，“即将”可能是明天，也可能是3,000年后。而“南大西洋异常区”的形成可能是异常区下方地球内部的变化所引起的。

3.4地磁总强度/矢量随轨道变化分析
选取amateur卫星，以2015年1月1日0点为计时起点，之后90分钟内卫星处地磁总强度随轨道变化曲线为：

amateur卫星地磁总强度随轨道变化曲线
90分钟内卫星处地磁矢量随轨道变化的三维示意图:

90分钟内amateur卫星地磁矢量随轨道变化的三维示意图
若延长卫星运行时间，三维矢量图如下：

图中箭头表示磁场矢量。大致可以看出，地磁场矢量从地理南极-地磁北极发出，在地理北极-地磁南极收回，与实际相符合。

3.5地磁导航方案的设计
地磁场是一个矢量场, 其强度大小和方向是位置的函数 ,同时地磁场具有丰富的例如总强度 、矢量强度、磁倾角 、磁偏角和强度梯度等特征, 为地磁匹配提供了充足的匹配信息 。因此 ,可以把地磁场当作一个天然的坐标系, 利用地磁场的测量信息来实现对卫星的导航定位 。
对本次仿真而言，地磁导航相当于做了一件与所编写的地磁场参数计算程序相反的事：已知一点地磁场参数，要求确定该点的经度、纬度、海拔高度（可以利用与地磁无关的方法得到），即进行一个逆向匹配的过程。由于整个地球磁场的全部信息数据库过于庞大，因此必须要对匹配方法进行优化。
我认为可以采取分级比对的方式：
1） 首先估计当前时刻太阳活动强度等产生宇宙射线的因素，利用随机信号处理的方法，滤除扰动磁场带来的噪声；
2） 测出地磁总强度，与全球等强度地磁图进行比对，估算卫星所处位置的大致区间；
3） 在选定的区间内，进一步比对磁场矢量强度，得到更精确的解。