[电离层建模学习笔记]开源程序M_GIM学习记录

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. 原程序里的小错误

   (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
   ...
   
   

2. 加上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
   

3. 提升代码运行速度（可选）：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 手动拼接的话，在很多子图的情况下也是费时费力。此处参考了使用Python批量拼接图片给出的程序，对子图进行拼接，这里贴出程序，以备日后查用：

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)

4. [2023.07.05补充]

解算之前数据准备工作详细流程

4.1 rinex文件

此程序只可以读取rinex3文件，如果自己的数据是rinex2版本，需要通过grzrnx或者rtkcov等转成rinex3，上述转换工具整理下放在这里：

链接: https://pan.baidu.com/s/1uUst-vbTDHKNEi3_M0Y_Xw 提取码: cpju 复制这段内容后打开百度网盘手机App，操作更方便哦

符合要求的rinex文件准备好后放在./Files_rinex/regional/（或者./Files_rinex/global/）目录下。

4.2 ionex文件

全球电离层TEC分布图文件，作为建模结果的外部参考，本实验中在运行Plot_nonSH_ChinaRegion.m时用到。不同机构发布的产品下载地址汇总包括但不限于下述地址：

1. https://cddis.nasa.gov/archive/gnss/products/ionosphere/（需注册账号）

2. ftp://gssc.esa.int/gnss/products/ionex/

3. ftp://gdc.cddis.eosdis.nasa.gov/pub/gps/products/ionex

4.3 dcb文件

该程序目前版本指定CAS发布的DCB产品，如需要使用其他机构的产品可以更改相应代码。CAS DCB产品可在以下地址下载：

1. ftp://cddis.gsfc.nasa.gov/pub/gps/products/mgex/dcb/

2. ftp://igs.ign.fr//pub/igs/products/mgex/dcb/

3. ftp://ftp.gipp.org.cn/product/dcb/

4.4 sp3文件

该程序需要多系统卫星的轨道产品，该产品可在以下地址下载：

ftp://igs.gnsswhu.cn/pub/gnss/products/mgex/