[电离层建模学习笔记]开源程序M_GIM学习记录
文章目录
- [电离层建模学习笔记]开源程序M_GIM学习记录
- 1. 程序相关信息
- 2. 程序学习记录
- 2.1 采用的数据说明
- 2.2 程序运行前
- 2.3 程序运行结果
- 3. 其他
1. 程序相关信息
开源程序M_GIM基于Matlab(Zhou et al., 2023)，用于实现全球和区域电离层建模，可支持多系统数据，建模方式也可以在多项式和球谐函数这种之间选择。该程序在M_DCB(Jin et al., 2012)的基础上发展而来，是便捷实用的电离层建模软件。
该程序发布在GPS Toolbox上，除了源程序还附带了对应论文中的测试数据。本文记录了该程序的学习过程，该程序在运行前还有一些小的错误需要改正，在此也做一个记录，以备日后查询。
相关的文献：
[1] Zhou C, Yang L, Li B, et al. M_GIM: a MATLAB-based software for multi-system global and regional ionospheric modeling[J]. GPS Solutions, 2023, 27(1): 1-7.
[2] Jin R, Jin S, Feng G. M_DCB: Matlab code for estimating GNSS satellite and receiver differential code biases[J]. GPS solutions, 2012, 16: 541-548.
2. 程序学习记录
2.1 采用的数据说明
本文采用的数据为2023年 DOY001 中国陆态网数据，共241个测站，包含GPS和GLONASS观测数据。本次解算仅采用GPS系统的数据，后发现所有测站GPS伪距观测值仅存在C1C、C1C和C2W，即C1、C2和P2数据，在计算中拟采用C1-P2来获取P4平滑值，需要在计算P4值之前加上码偏差改正，将C1归化到P1上，接下来将介绍具体的改动。
2.2 程序运行前
本次解算利用陆态网数据对中国区域电离层进行区域建模，采用6阶4次球谐函数模型，解算仅考虑GPS数据，故用到的主要程序为：Main_nonSH.m 和./Func_nonSH/Get_nonSH_G.m，其他子程序调用情况可以在这两个文件中看到。在程序运行前，首先准备好观测值文件、sp3轨道、DCB文件分别存放在对应文件夹下，运行前需要对上述程序以下部分做对应修改：
-
原程序里的小错误
(1) Get_nonSH_G.m （Line 42-50）
... if isempty(d_sat) G_n_s=32; else G_n_s=32-length(d_sat);%the number of satellites disp([doy ' PRN ',num2str(d_sat) ,' have no observations.']); for k=length(d_sat):-1:1 %% gps_d_sat -> d_sat gpsx(:,d_sat(k))=[];gpsy(:,d_sat(k))=[];gpsz(:,d_sat(k))=[]; % 没有数据的卫星坐标设为空 end end ...(2) Get_nonSH_G.m （Line 70-75）
... if ~isempty(d_sat) for k=length(d_sat):-1:1 %% GPSP4(:,k)=[] -> GPSP4(:,d_sat(k)) GPSP4(:,d_sat(k))=[]; end end ... -
加上C1归化到P1的部分
（1) 读取DCB文件中的C1C -> C1W：read_dcb
line 21-23 SDCB_REF.gps=zeros(n_d,32); SDCB_REF.gpsc1p1=zeros(n_d,32); % edit by He Rong 2023/06/08 SDCB_REF.glo=zeros(n_d,24); line 32 [GPS_DCB_rec,SDCB_REF.gps(i,:),SDCB_REF.gpsc1p1]=r_gps_dcb([i_ipath '/CAS0MGXRAP_20' num2str(sdoy) '0000_01D_01D_DCB.BSX'],Sites_Info.name(index2)); % edit by He Rong 2023/06/08 line 189-221 function [DCB_rec,DCB_sat,DCB_sat_C1P1]=r_gps_dcb(fpath,sites) % edit by He Rong 2023/06/08 fid=fopen(fpath,'r'); n_r=length(sites); DCB_rec=linspace(0,0,n_r); DCB_sat=linspace(0,0,32); DCB_sat_C1P1 = linspace(0,0,32); % edit by He Rong 2023/06/08 while 1 line=fgetl(fid); if ~ischar(line), break, end %--satellites' DCB if length(line)>100 && strcmpi(line(1:7),' DSB G') && strcmpi(line(26:33),'C1W C2W') && strcmpi(line(16:19),' ') prn=str2double(line(13:14)); DCB_sat(prn)=str2double(line(83:91)); continue; end % satelite DCB：C1C->C1W if length(line)>100 && strcmpi(line(1:7),' DSB G') && strcmpi(line(26:33),'C1C C1W') && strcmpi(line(16:19),' ') prn=str2double(line(13:14)); DCB_sat_C1P1(prn)=str2double(line(83:91)); continue; end % edit by He Rong 2023/06/08 %--receivers' DCB if length(line)>100 && strcmpi(line(1:12),' DSB G G') && strcmpi(line(26:33),'C1W C2W') index=find(strcmpi(line(16:19),sites), 1); if ~isempty(index) DCB_rec(index)=str2double(line(83:91)); end continue; end end fclose(fid); end(2) 在计算P4前进行改正：Get_P4
line 1 function [] = Get_P4(path_obs,path_p4,path_sp3,Sites_Info,lim,sate_mark, SDCB_REF) % 对应外部调用也需要修改 line 37-41 % C1W位置上的其实是C1C(C1)的数据，需要将其转到P1上(加上C1C->C1W) % 在参与计算前，将所有测站的所有 C1W 数据都进行改正 % SDCB_REF结构体里SDCB_REF.gpsc1p1记录所有卫星的C1P1 % edit by He Rong 2023/06/09 [obs] = correctC1C_C1W(SDCB_REF,obs); line 1277-1308 function [obs] = correctC1C_C1W(SDCB_REF,obs) % edit by He Rong 2023/06/08 %% correct the observations C1C to C1W use DCB C1C->C1W % INPUT: % SDCB_REF: satellites DCB for each system % obs: original observation structs % % % OUTPUT： % obs: updated observation structs % ________ correct the GPS C1C observation (in C1W position) ONLY _________ %%%% Constant defination %%%% V = 2.99792458E+08; NS2M = (V*1.0E-9); [row, col] = size(obs.GPSC1W); % 获取C1W数据的大小 for i = 1:1:col % 逐卫星循环 % satTmp = obs.GPSC1W(:,i); % 选取当前卫星的数据 prn = i; dcb_corr = SDCB_REF.gpsc1p1(prn); for j = 1:1:row % 当前卫星逐历元循环 if obs.GPSC1W(j,prn) == 0 continue end obs.GPSC1W(j,prn) = obs.GPSC1W(j,prn) - dcb_corr * NS2M; end end end -
提升代码运行速度（可选）：Main_nonSH.m
line 25-26 %--the cores number of computer processor Corenum=2; % 此处根据各人计算机而定
2.3 程序运行结果
程序运行需要较长时间，其间生成的观测值以mat格式保存在./OBS/regional/23001/目录下，精密轨道保存在./SP3，计算得到的P4平滑值保存在./P4/regional/GPS/23001/，最终最小二乘估计的结果保存在./M_Result/，接下来根据结果计算区域格网VTEC，并输出文件。由于本次解算仅采用GPS数据，故在原程序Plot_nonSH.m的基础上做了一些修改，命名为Plot_nonSH_ChinaRegion.m，可作为参考：
%% ==========Plot and write regional ionospheric maps=================
%% nonintergral SH model
clear all;
close all;
doy=23001; fig=24; K=6; M=4;
load('sate_mark.mat');
load(['M_Result/nonSH_G',num2str(doy),'.mat']);
load('sate_mark.mat');
warning off;
addpath('Tools/m_map','Tools/m_map/private');
lat2=0; lat1=60; lon1=70; lon2=140;
latlim=2.5; lonlim=5;
% latGrid = lat1:-latlim:lat2;
% lonGrid = lon1:lonlim:lon2;
VVTEC = Get_VTEC(fig, latlim, lonlim, IONC, NN, m0, K, M, lat2, lat1, lon1, lon2);
VTEC=VVTEC;VTEC(VTEC(:,4)<0,4)=0.05;
RMS=[VTEC(:,1:3),VTEC(:,5)];
% read CODE final GIMs (codg2750.19i) as reference
disp('--------> Read CODE final GIMs as reference !');
IGSData=read_ionex(fig,'TEC');
AreaTEC=Get_areaTEC(fig,lat2,lat1,lon1,lon2,IGSData);
DIFFTEC=[VTEC(:,1:3),VTEC(:,4)-AreaTEC(1:size(VTEC,1),4)];
% 保存变量
fname=['RIM_data_China\' num2str(doy) '.mat'];
save(fname,'VTEC','RMS','IGSData','AreaTEC','DIFFTEC','-mat');
% % 加载变量
% load(['RIM_data_China\',num2str(doy),'.mat'])
% %%%% 这里统计一下VTEC误差序列的结果 %%%%%
% sess = unique(DIFFTEC(:,1));
% staRes = zeros(length(sess), 3);
% for iis = 1:1:length(sess)
% diffTmp = [];
% for kk = 1:1:length(DIFFTEC(:,3))
% if DIFFTEC(kk,1) == sess(iis)
% diffTmp(end +1, :) = DIFFTEC(kk,:);
% end
% end
% % 统计
% biasVal = mean(diffTmp(:,4), 1);
% stdVal = std(diffTmp(:,4), 0, 1); % std(A，flag):flag用来区分std求法，如果是0，则代表除以N-1，如果是1代表的是除以N
% rmsVal = rms(diffTmp(:,4));
% fprintf('Session: %2d Bias: %6.3f STD: %6.3f RMS: %6.3f \n', sess(iis), biasVal, stdVal, rmsVal);
% staRes(iis,:) = [biasVal stdVal rmsVal];
%
% end
% biasVal = mean(staRes(:,1), 1);
% stdVal = mean(staRes(:,2), 1);
% rmsVal = mean(staRes(:,3), 1);
% fprintf('Average:Bias: %6.3f STD: %6.3f RMS: %6.3f \n', biasVal, stdVal, rmsVal);
Pname='2023-001-nonSH-';
Plot_TEC(fig,latlim,lonlim,Pname,VTEC,lat1,lat2,lon1,lon2,0,50);
Pname_D='2023-001-nonSH-D-';
Plot_TEC(fig,latlim,lonlim,Pname_D,DIFFTEC,lat1,lat2,lon1,lon2,-10,10);
Pname_R='2023-001-nonSH-RMS-';
Plot_TEC(fig,latlim,lonlim,Pname_R,RMS,lat1,lat2,lon1,lon2,0,5);
% write results to ionex files
% 由于只用了GPS数据，所以其他几种系统的sDCB和rDCB均赋值为0
C_R = zeros(23,1);
E_R = zeros(18,1);
EX_R = zeros(11,1);
R_R = zeros(28,1);
C_S = zeros(1,15);
E_S = zeros(1,24);
R_S = zeros(1,21);
Write_ionex(doy,fig,G_R,C_R,E_R,EX_R,R_R,G_S,C_S,E_S,R_S,Pname,VTEC,sate_mark,total_station,list_P4);
%% ++++++++++++++++PLOT OVER!!!+++++++++++++++++++
此处，对应的Write_ionex里面在输出每颗卫星的DCB之前，加上判断，若全为0则不输出。对应程序修改可自行完成。
以上，程序运行结束，下面贴出结果图：
以下四张图展示各时段的结果，依次为：估计得到的VTEC、与CODE作差得到的VTEC插值、估计VTEC的RMS、穿刺点分布（此为自己程序写的）
3. 其他
2.3 贴出的图是采用Python拼接而成的图，这里推荐大家使用，特别是需要将很多子图拼在一起的情况下，采用matlab直接画的话，需要微调很久（即使有循环画子图的脚本，colorbar、子图命名等均需要花费大量时间），采用adobe AI 手动拼接的话，在很多子图的情况下也是费时费力。此处参考了
给出的程序，对子图进行拼接，这里贴出程序，以备日后查用：
import os
import math
from PIL import Image
def merge_images(image_folder, output_file, n, m):
# 获取所有图像文件的列表
image_files = [os.path.join(image_folder, f) for f in os.listdir(image_folder) if f.endswith('.png')]
# 计算每个小图像的大小和大图像的大小
image_count = len(image_files)
if image_count == 0:
print('No image files found in the directory:', image_folder)
return
# 计算小图像的大小以及大图像的大小
img = Image.open(image_files[0])
img_size0 = img.size[0]
img_size1 = img.size[1]
new_img_size0 = img_size0 * n
new_img_size1 = img_size1 * m
# 创建一个新的大图像
new_img = Image.new('RGB', (new_img_size0, new_img_size1), 'white')
# 将所有小图像粘贴到新图像的正确位置
for i, f in enumerate(image_files):
row = int(i / n)
col = i % n
img = Image.open(f)
img = img.resize((img_size0, img_size1))
new_img.paste(img, (col * img_size0, row * img_size1))
# 保存大图像
new_img.save(output_file)
# 用法示例
# 图片目录
image_folder = r'x:/xxxxx/xxxxxxxx/M_GIM/M_PLOT/nonSH/'
# 输出目录
output = r'./nonSH/'
if not os.path.exists(output): # True/False
os.mkdir(output)
# 输出图片名（必须指定图片名.png）
output_file = r'./nonSH/comb-nonSH.png'
n = 4 # 每行显示的图像数
m = 6 # 每列显示的图像数
merge_images(image_folder, output_file, n, m)
mage_folder = r'x:/xxxxx/xxxxxxxx/M_GIM/M_PLOT/nonSH/'
# 输出目录
output = r'./nonSH/'
if not os.path.exists(output): # True/False
os.mkdir(output)
# 输出图片名（必须指定图片名.png）
output_file = r'./nonSH/comb-nonSH.png'
n = 4 # 每行显示的图像数
m = 6 # 每列显示的图像数
merge_images(image_folder, output_file, n, m)
