在本科阶段，很多同学可能会涉及到一些建模要求。笔者在这里公布自己的建模历程，希望可以帮到大家捏。
地磁场由不同场源的贡献组成，主要来自地球液态外核、地壳/上地幔、电离层、磁层等几个部分。若从地面台站数据出发，电离层和磁层可以归为外源，而从卫星数据出发则将电离层归为内源，磁层作为外源。地磁场是各种磁场成分的叠加，且各组成部分之间存在相互作用。其中，主磁场（即地核磁场）约占地球总磁场的95％，是地磁场中最主要的部分。考虑到地磁场的特点，即不随时间变化（在以年为时间尺度）；可以认为外部（地表附近）没有外源电流贡献等等。
大聪明高斯在1839年把球谐函数分析方法应用于地磁场，得出了地磁场的数学表达形式，奠定了地磁学的数理基础。高斯真就是大神，数学物理绕不开，地球物理行星科学也经常见到捏。
假定该区域没有磁源（净电流和变化的电场）。根据麦克斯韦方程组，磁感应强度𝑩满足：
等式右边等于0：
可以认为：
磁场是某个标量场的梯度。同时考虑磁场的散度等于0：
我们可以得到这个拉普拉斯方程的解：
式中，是地球平均半径，
是到地心的径向距离，
分别是余纬和经度，
是连带勒让德函数，
是内源场的高斯系数，直接反映了地球内部过程产生的磁场，
是外源场的高斯系数，直接反映了地球外部过程产生的磁场，外源场在研究区域内是有限的，项的存在表示它在距离地球非常远的地方对磁场的描述就不符合实际了。我们可以认为，每一个球谐函数都是一个磁场，而地球磁场是很多不同类型磁场的叠加。
在这里我们只关心地球本身产生的磁场：
磁场在球坐标系下三分量的表达式为
我们可以把上面的一坨写成下面的矩阵。矩阵是matlab的语言。来matlab就要说matrix话捏。
将上式简写为
通过构建核矩阵，我们就可以基于观测数据使用最小二乘法对模型高斯系数进行反演
从而由观测数据d和构造的矩阵A来解出高斯系数x.
下面可以开始代码环节了捏。
地磁台站的数据文件长这样子的捏。
把文件全部命名为zhsj（i）方便处理捏。处理数据的代码：
filename=strcat('zhsj (', num2str(i), ').min');
data{i}=importdata(filename);%按照文件的规律，循环提取文件
%-------------------------------------------------------------------------------------------
lat_1{i} = data{1,i}{5};
lon_1{i} = data{1,i}{6};
lat_2{i} = strsplit(lat_1{i});
lat_3{i} = lat_2{i}{4};
lat_4(i) = str2double(lat_3{i});%将文件中有用的部分分块，提取出来储存在变量中
lon_2{i} = strsplit(lon_1{i});
lon_3{i} = lon_2{i}{4};
lon_4(i) = str2double(lon_3{i});
adata{i} = data{1,i}(33:end);
end
%-------------------------------------------------------------------------------------------
for k=1:wenjianshu
for m=1:1420
undata{k,m}=adata{k}{m,1};
undatasp{k,m}=strsplit(undata{k,m});
xundata{k,m} = undatasp{k,m}{4};
rxundata(k,m) = str2double(xundata{k,m});
yundata{k,m} = undatasp{k,m}{5}; %将文件中有用的部分分块，提取出来储存在变量中
ryundata(k,m) = str2double(yundata{k,m});
zundata{k,m} = undatasp{k,m}{6};
rzundata(k,m) = str2double(zundata{k,m});
end
end
%-------------------------------------------------------------------------------------------
xm=rxundata;
xl=xm';
x=mean(xl);
ym=ryundata;
yl=ym';
y=mean(yl);
zm=rzundata;
zl=zm';
z=mean(zl);
for f=1:wenjianshu
last(1,3*f-2)=-z(f);
last(1,3*f-1)=-x(f);%对磁场数据进行求平均处理，并全部转化为一个向量last
last(1,3*f)=y(f);
end
r(1:wenjianshu,1)=1;
yuwei=pi*(90-lat_4)/180;
theta=yuwei';
jindu=pi*lon_4/180; %将经纬度数据导入计算球谐矩阵的函数中，计算出每个经纬度对应的A矩阵
phi=jindu';
[G_r, G_theta, G_phi] = cal_SHA_int(r, theta, phi, N);%计算A矩阵的函数
%-------------------------------------------------------------------------------------------
for mm=1:wenjianshu
A(3*mm-2,:)=G_r(mm,:);
A(3*mm-1,:)=G_theta(mm,:);%将计算得出来的每个位置的A矩阵打包，构成最后对应每个台站的A矩阵
A(3*mm,:)=G_phi(mm,:);
end
%-------------------------------------------------------------------------------------------
B=A'*A;
C=inv(B)*A';
%C=inv(A);
fin=C*last'; %使用最小二乘解方法解出高斯系数
%-------------------------------------------------------------------------------------------
关于球谐函数的部分代码如下：
function [A_r, A_theta, A_phi] = cal_SHA_int(r, theta, phi, N)
% Calculates design matrices A_i that connects the vector
% of (Schmidt-normalized) spherical harmonic expansion coefficients,
% x = (g_1^0; g_1^1; h_1^1; g_2^0; g_2^1; h_2^1; ... g_N^N; h_N^N)
% and the magnetic component B_i, where "i" is "r", "theta" or "phi":
% B_i = A_i*x
% Input: r(:) radius vector (in units of the reference radius a)
% theta(:) colatitude (in radians)
% phi(:) longitude (in radians)
% N maximum degree/order
% 'int' means internal sources (g_n^m and h_n^m)
% copied from design_SHA.m (January 2003, Nils Olsen, DSRI)
N_koeff=(N+1)^2-1;
cos_theta = cos(theta);
sin_theta = sin(theta);
N_data = length(theta);
A_r = zeros(N_data, N_koeff);
A_theta = zeros(N_data, N_koeff);
A_phi = zeros(N_data, N_koeff);
k=0;
for n = 1:N
r_n = r.^(-(n+2));
Pnm = legendre(n, cos_theta, 'sch')'; % P_n^m and derivatives vrt. theta
dPnm(:,n+1) = (sqrt(n/2).*Pnm(:,n)); % m=n
dPnm(:,1) = -sqrt(n*(n+1)/2.).*Pnm(:,2); % m=0
if n > 1; dPnm(:,2)=(sqrt(2*(n+1)*n).*Pnm(:,1)-sqrt((n+2)*(n-1)).*Pnm(:,3))/2; end % m=1
for m = 2:n-1 % m=2...n-1
dPnm(:,m+1)=(sqrt((n+m)*(n-m+1)).*Pnm(:,m)-sqrt((n+m+1)*(n-m)).*Pnm(:,m+2))/2;
end
if n == 1 dPnm(:,2) = sqrt(2)*dPnm(:,2); end
for m = 0:n
cos_phi = cos(m*phi);
sin_phi = sin(m*phi);
if m == 0
k = k+1; % terms corresponding to g_n^0
A_r(:,k) = (n+1).*r_n(:).*Pnm(:,1);
A_theta(:,k) = -r_n(:).*dPnm(:,1);
A_phi(:,k) = r_n(:)*0;
else
k = k+1; % terms corresponding to g_n^m
A_r(:,k) = (n+1).*r_n(:).*cos_phi.*Pnm(:,m+1);
A_theta(:,k) = -r_n(:).*cos_phi.*dPnm(:,m+1);
A_phi(:,k) = r_n(:).*m.*sin_phi.*Pnm(:,m+1)./sin_theta;
k = k+1; % terms corresponding to h_n^m
A_r(:,k) = (n+1).*r_n(:).*sin_phi.*Pnm(:,m+1);
A_theta(:,k) = -r_n(:).*sin_phi.*dPnm(:,m+1);
A_phi(:,k) = -r_n(:).*m.*cos_phi.*Pnm(:,m+1)./sin_theta;
end
end % m
end % n计算出来最后的高斯系数建模就完成了，然后就可以画图图哩。
我们先计算在经纬度网格里面的每一个网格的磁场：
dx = pi/120;
col = -pi/2:dx:pi/2;
az = -pi:2*dx:pi;
%az_r=az;
for aa=1:121
if az(aa)>=pi
az_r(aa)=2*pi+az(aa);
else az_r(aa)=az(aa);
end
end
col_r=pi/2-col;
az_r=az;
[Theta,Phi] = meshgrid(col,az);
lon=180.*Phi./pi;
lat=180.*Theta./pi;
r=1;
N=2;
for mm=1:121
for kk=1:121
phi=az_r(mm);
theta=col_r(kk);
[G_r, G_theta, G_phi] = cal_SHA_int(r, theta, phi, N);
B_R_M(mm,kk)=G_r*fin;
B_T_M(mm,kk)=G_theta*fin;
B_P_M(mm,kk)=G_phi*fin;
end
end然后就可以画图哩。这里使用的m_map画的，可以去学一下挺方便的。
m_proj('hammer','clongitude',-360);
MMM=B_r;
m_pcolor(lon-360,lat,MMM);
m_contourf(lon-360, lat, MMM, 'linewidth', 0.3);
shading interp;
colorbar;
colormap(flipud(get_cmap('PuOr')));
clim([100 800]);
%m_coast('linewidth',0.5,'color','r');
%m_coast('Color',[0 0 0],'LineWidth',0.5); % add coastline
m_grid('box','off','xticklabels',[],'xtick',12,'yticklabels',[],'ytick',12);
saveas(gcf, 'FSACBA.png')
最后的结果就是这样的，至于科学问题就自己去思考吧。
1074
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
