VTEC - 全球电离层垂直电子总量估计系统

项目简介

VTEC（Vertical Total Electron Content）是一个基于地基GNSS观测数据估计全球电离层垂直电子总量的系统。该系统利用多GNSS系统（GPS、GLONASS、Galileo、BeiDou等）的双频观测数据，计算电离层斜向电子总量（STEC），然后转换为垂直电子总量（VTEC），最终利用球谐函数建立全球电离层模型，并联合估计接收机和卫星差分码偏差（DCB）。

系统采用C++语言开发，结合LibTorch库实现高效的数学计算，具有高度模块化设计和丰富的配置选项，适用于科学研究、导航增强和空间天气监测等多种应用场景。

功能特性

• 多GNSS系统支持：

  • GPS（美国全球定位系统）
  • GLONASS（俄罗斯全球导航卫星系统）
  • Galileo（欧盟全球卫星导航系统）
  • BeiDou（中国北斗卫星导航系统）
  • QZSS（日本准天顶卫星系统）

• 文件格式支持：

  • RINEX 2.x/3.x格式观测数据文件
  • 广播星历和精密星历(SP3)文件
  • IONEX格式输出

• 完善的数据预处理功能：

  • 高度角筛选和信噪比筛选
  • 周跳检测与修复（多种算法联合检测）：
    • 几何无关组合(GF)检测
    • Melbourne-Wübbena组合(MW)检测
    • 时间差分检测
    • 电离层变化率检测
    • TurboEdit联合检测策略
  • 多路径效应检测与缓解
  • 多策略观测值加权：
    • 高度角加权
    • 信噪比加权
    • 多路径加权
    • 综合加权

• 高精度STEC/VTEC计算：

  • 电离层穿刺点精确定位
  • 多种映射函数支持：
    • 单层模型(SLM)
    • 修正单层模型(MSLM)
    • 薄壳映射函数(TSM)
    • 方位角梯度扩展映射函数
  • 相位对齐伪距水平的STEC精确计算：
    • 基于周跳探测的观测弧段自动划分
    • 滑动窗口平滑偏移量计算
    • 相位精度结合伪距绝对水平
  • DCB与球谐系数联合估计

• 全球电离层建模：

  • 基于球谐函数的全球电离层建模
  • 可配置的球谐函数阶数
  • 支持时间序列TEC地图生成
  • 正则化与约束条件支持

• 系统设计特点：

  • 高度模块化设计
  • 丰富的配置选项
  • 多线程并行计算
  • 优化的数值计算性能

系统要求

• C++17兼容的编译器（GCC 7+, Clang 5+, MSVC 2019+）
• CMake 3.12+
• LibTorch 2.0+
• Boost库 (可选)

安装步骤

1. 安装依赖

Ubuntu/Debian

sudo apt update
sudo apt install build-essential cmake git libboost-all-dev

LibTorch安装

从PyTorch官网下载适合您系统的LibTorch库：https://pytorch.org/get-started/locally/

解压到方便访问的位置，例如：

cd ~/Documents
tar -xzf libtorch-cxx11-abi-shared-with-deps-x.x.x+cpu.tar.gz

2. 编译项目

git clone https://github.com/your-username/vtec.git
cd vtec
mkdir build
cd build
cmake -DCMAKE_PREFIX_PATH=/path/to/libtorch ..
make -j4

使用方法

基本用法

cd build
./ionosphere_vtec --config ../config/default_config.ini

命令行选项

用法: ./ionosphere_vtec [选项]

选项:
  --config <file>           指定配置文件路径
  --obs <file>              RINEX观测文件路径
  --nav <file>              RINEX导航文件路径
  --sp3 <file>              SP3精密星历文件路径
  --output <file>           输出IONEX文件路径
  --elev <degrees>          最小高度角阈值
  --degree <n>              球谐函数最高阶数
  --systems <sys>           要处理的GNSS系统 (如 "G,R,E,C")
  --interval <seconds>      时间间隔（秒）
  --mapping-function <type> 映射函数类型 (SLM, MSLM, TSM)
  --phase-aligned           启用相位对齐伪距水平STEC计算（默认）
  --no-phase-aligned        禁用相位对齐伪距水平STEC计算
  --smoothing-window <n>    相位对齐平滑窗口大小（默认30）
  --azimuthal-asymmetry     启用方位角非对称性（默认禁用）
  --lat-gradient <value>    纬度方向梯度系数
  --lon-gradient <value>    经度方向梯度系数
  --help                    显示帮助信息

批处理模式

系统支持通过模式匹配处理多个文件：

./ionosphere_vtec --config ../config/default_config.ini --obs-pattern "data/*.??o" --nav-pattern "data/*.??n"

配置文件说明

配置文件为INI格式，包含以下主要部分：

[General]

基本配置参数：

verbose = true          # 是否显示详细输出
num_threads = 4         # 使用线程数
output_dir = ./output   # 输出目录

[RINEX]

RINEX文件和GNSS系统选择：

gnss_systems = G,R,E,C,J               # 启用的GNSS系统
observation_pattern = data/*.??o        # 观测文件模式
navigation_pattern = data/*.??n         # 导航文件模式
sp3_pattern = data/*.sp3                # 精密星历文件模式
min_elevation = 10.0                    # 最小高度角(度)

[Preprocessing]

数据预处理参数：

# 周跳检测方法
cycle_slip_method = All                 # 周跳检测方法: None, GeometryFree, MelbourneWubbena, TimeDifference, IonosphereRate, TurboEdit, All
gf_threshold = 0.25                     # 几何无关组合周跳检测阈值（周）
mw_threshold = 0.75                     # Melbourne-Wübbena组合周跳检测阈值（周）
td_threshold = 0.05                     # 时间差分周跳检测阈值（周/秒）
iono_rate_threshold = 0.02              # 电离层变化率周跳检测阈值（TECU/秒）
turbo_edit_strategy = Conservative      # TurboEdit策略: Conservative, Strict, Adaptive

# 周跳处理
fix_cycle_slips = true                  # 是否修复周跳
max_consecutive_slips = 5               # 最大连续周跳数（超过则拒绝该弧段）
max_time_gap = 300.0                    # 最大时间间隙（秒，超过则自动认为有周跳）

# 观测值加权
weighting_mode = Elevation              # 观测值加权模式: Uniform, Elevation, Snr, Multipath, Combined
elevation_power = 2.0                   # 高度角权重指数

[STEC]

斜向电子总量计算参数：

use_smoothed_code = true                # 是否使用平滑伪距
apply_phase_windup = true               # 是否应用相位缠绕修正
enable_phase_aligned = true             # 是否启用相位对齐伪距水平的STEC计算
smoothing_window_size = 30              # 相位对齐平滑窗口大小（观测历元数）
min_arc_length = 5.0                    # 最小弧段长度（分钟）
remove_outliers = true                  # 是否移除异常值
outlier_threshold = 3.0                 # 异常值识别阈值（标准差倍数）

[VTEC]

垂直电子总量计算参数：

iono_height = 450.0                     # 电离层等效高度（km）
mapping_function = SLM                  # 映射函数类型: SLM, MSLM, TSM
apply_azimuthal_asymmetry = false       # 是否应用方位角非对称性
latitude_gradient = 0.0                 # 纬度方向梯度
longitude_gradient = 0.0                # 经度方向梯度
use_azimuthal_weighting = false         # 是否使用方位角加权
north_weight = 1.0                      # 北方向权重
east_weight = 1.0                       # 东方向权重
south_weight = 1.0                      # 南方向权重
west_weight = 1.0                       # 西方向权重
latitude_cutoff = 80.0                  # 纬度截断值（度）

[SphericalHarmonic]

球谐函数建模参数：

max_degree = 8                          # 球谐函数最高阶数
max_order = 8                           # 球谐函数最高次数（通常等于阶数）
use_temporal = true                     # 是否使用时间变化项
temporal_degree = 2                     # 时间变化项的阶数
regularization = true                   # 是否使用正则化
regularization_factor = 0.001           # 正则化因子
constrain_nighttime = true              # 是否约束夜间电离层

[DCB]

差分码偏差估计参数：

estimate_satellite_dcb = true           # 是否估计卫星DCB
estimate_receiver_dcb = true            # 是否估计接收机DCB
constrain_zero_mean = true              # 是否约束DCB零均值
sat_dcb_std_threshold = 3.0             # 卫星DCB标准差阈值（ns）
rec_dcb_std_threshold = 5.0             # 接收机DCB标准差阈值（ns）

[IONEX]

IONEX输出参数：

grid_longitude_resolution = 5.0         # 输出格网经度分辨率（度）
grid_latitude_resolution = 2.5          # 输出格网纬度分辨率（度）
temporal_resolution = 2                 # 时间分辨率（小时）
output_satellite_dcb = true             # 是否输出卫星DCB
output_receiver_dcb = true              # 是否输出接收机DCB
output_rms = true                       # 是否输出RMS
ionex_version = 1                       # IONEX版本
output_file = output.ionex              # 输出文件名

详细配置示例请参考 config/default_config.ini。

技术实现详解

多种映射函数实现

本系统实现了多种电离层映射函数，以适应不同的应用场景和研究需求：

1. 单层模型(SLM)

最常用的映射函数，假设电离层集中在距地面一定高度的薄球壳上。

数学公式：

MF(z) = 1/cos(z')

其中z'是穿刺点处的天顶角，通过卫星高度角和地球曲率计算：

z' = arcsin((R_E/(R_E+h))*sin(z))

• R_E: 地球半径
• h: 电离层等效高度
• z: 观测站天顶角 = 90° - 高度角

实现代码片段:

double STECCalculator::calculateSLMMappingFunction(double elevation) {
    double z = PI / 2.0 - elevation;  // 天顶角
    double Re = EARTH_RADIUS;         // 地球半径
    double Ri = Re + ionoHeight_;     // 电离层高度
    
    double sinZp = (Re / Ri) * std::sin(z);
    double zp = std::asin(sinZp);     // 穿刺点天顶角
    
    return 1.0 / std::cos(zp);
}

适用场景：

• 一般电离层研究和建模
• 中高纬度地区

2. 修正单层模型(MSLM)

考虑了地球曲率的修正版本，提高了高度角较低时的映射准确性。

数学公式：

MF(z) = (Ri * cos(z')) / (cos(z) * h)

实现代码片段:

double STECCalculator::calculateMSLMMappingFunction(double elevation) {
    double z = PI / 2.0 - elevation;  // 天顶角
    double Re = EARTH_RADIUS;         // 地球半径
    double Ri = Re + ionoHeight_;     // 电离层高度
    
    // 计算地心角
    double psi = asin((Re / Ri) * sin(z));
    
    // 计算穿刺点天顶角
    double zp = z - psi;
    
    // 计算穿刺点到电离层路径长度与垂直路径的比率
    double mapFunc = (Ri * cos(zp)) / cos(z) / ionoHeight_;
    
    return mapFunc;
}

适用场景：

• 低高度角观测（<15°）
• 更精确的电离层VTEC建模
• 电离层层析成像

3. 薄壳映射函数(TSM)

通过多项式拟合的映射函数，在低高度角情况下表现更好。

数学公式：

MF(z) = (a0 + a1*z^1 + a2*z^2) / cos(z)

其中a0、a1、a2是经验参数，通常a0=1.0, a1=0.15, a2=0.0018。

实现代码片段:

double STECCalculator::calculateTSMMappingFunction(double elevation) {
    double z = PI / 2.0 - elevation;  // 天顶角
    
    // 添加纬度相关的修正项
    double a0 = 1.0;
    double a1 = 0.15;
    double a2 = 0.0018;
    
    return (a0 + a1 * pow(z, 1.0) + a2 * pow(z, 2.0)) / cos(z);
}

适用场景：

• 大范围区域的电离层建模
• 特殊区域（如赤道异常区）的映射

4. 方位角映射函数扩展

考虑电离层水平梯度，在标准映射函数基础上增加方位角依赖项。

数学公式：

MF(z,A) = MF(z) * (1 + G_NS*cos(A) + G_EW*sin(A))

其中：

• A是方位角
• G_NS是南北梯度系数
• G_EW是东西梯度系数

实现代码片段:

double STECCalculator::calculateMappingFunctionWithAzimuth(
    double elevation, double azimuth, double latGrad, double lonGrad) {
    
    double baseMapping = calculateMappingFunction(elevation);
    
    // 方位角转弧度
    double azimuthRad = azimuth * DEG2RAD;
    
    // 基于方位角添加水平梯度项
    double azimuthTerm = latGrad * cos(azimuthRad) + lonGrad * sin(azimuthRad);
    
    // 添加到基本映射函数
    return baseMapping * (1.0 + azimuthTerm);
}

适用场景：

• 存在强电离层扰动区域
• 赤道异常区域
• 太阳爆发期间的精细建模

相位对齐伪距水平STEC计算

系统实现了基于观测弧段的相位对齐伪距水平STEC计算方法，结合了载波相位的高精度和伪距的绝对水平优势：

1. 观测弧段划分

基于周跳标志和时间间隙将观测数据划分为连续的弧段：

std::vector<std::vector<GNSSObservation>> STECCalculator::divideObservationsIntoArcs(
    const std::vector<GNSSObservation>& observations) {
    
    std::vector<std::vector<GNSSObservation>> arcs;
    
    if (observations.empty()) {
        return arcs;
    }
    
    // 获取配置参数
    auto& config = Config::getInstance();
    auto& prepSettings = config.preprocessingSettings;
    double maxTimeGap = prepSettings.maxTimeGap;
    
    // 开始一个新的弧段
    std::vector<GNSSObservation> currentArc;
    currentArc.push_back(observations[0]);
    
    // 遍历观测数据，根据周跳标志或时间间隙划分弧段
    for (size_t i = 1; i < observations.size(); ++i) {
        const auto& obs = observations[i];
        const auto& prevObs = observations[i-1];
        
        // 检查是否有周跳或时间间隙过大
        double timeGap = obs.time - prevObs.time;
        
        if (obs.hasCycleSlip || timeGap > maxTimeGap) {
            // 结束当前弧段并开始新的弧段
            if (!currentArc.empty()) {
                arcs.push_back(currentArc);
                currentArc.clear();
            }
        }
        
        // 添加到当前弧段
        currentArc.push_back(obs);
    }
    
    // 添加最后一个弧段
    if (!currentArc.empty()) {
        arcs.push_back(currentArc);
    }
    
    return arcs;
}

2. 相位对齐伪距水平STEC计算

计算每个弧段内的STEC：

std::vector<STECObservation> STECCalculator::calculatePhaseAlignedSTEC(
    const std::vector<GNSSObservation>& arc) {
    
    // 计算伪距STEC和相位STEC
    std::vector<double> pseudorangeSTEC;
    std::vector<double> phaseSTEC;
    
    for (const auto& obs : arc) {
        double p_stec = calculatePseudorangeSTEC(obs);
        pseudorangeSTEC.push_back(p_stec);
        
        double ph_stec = calculatePhaseSTEC(obs);
        phaseSTEC.push_back(ph_stec);
    }
    
    // 计算偏移量并使用滑动窗口平滑
    std::vector<double> offsets;
    for (size_t i = 0; i < arc.size(); ++i) {
        offsets.push_back(pseudorangeSTEC[i] - phaseSTEC[i]);
    }
    
    // 使用滑动窗口平均计算稳定的偏移量
    std::vector<double> smoothedOffsets = calculateSmoothOffsets(offsets, arc.size());
    
    // 应用平滑偏移量，计算最终的STEC
    std::vector<STECObservation> stecObs;
    for (size_t i = 0; i < arc.size(); ++i) {
        // 相位STEC加上平均偏移量得到校准的STEC
        double calibratedSTEC = phaseSTEC[i] + smoothedOffsets[i];
        
        // 创建STEC观测结果
        STECObservation stecObservation;
        stecObservation.time = arc[i].time;
        stecObservation.satellite = arc[i].getSatellite();
        stecObservation.station = arc[i].getStation();
        stecObservation.elevation = arc[i].elevation;
        stecObservation.azimuth = arc[i].azimuth;
        
        // 计算电离层穿刺点
        auto ippCoords = calculateIonosphericPiercePoint(arc[i]);
        stecObservation.latitude = ippCoords[0];
        stecObservation.longitude = ippCoords[1];
        
        stecObservation.stec = calibratedSTEC;
        stecObservation.stecSigma = 0.1;  // 相位STEC通常比伪距STEC精度高
        
        stecObs.push_back(stecObservation);
    }
    
    return stecObs;
}

3. 计算滑动窗口平滑偏移量

std::vector<double> STECCalculator::calculateSmoothOffsets(
    const std::vector<double>& offsets, int arcSize) {
    
    // 获取配置参数
    auto& config = Config::getInstance();
    auto& stecSettings = config.stecSettings;
    int windowSize = stecSettings.smoothingWindowSize;
    
    std::vector<double> smoothedOffsets;
    
    for (int i = 0; i < arcSize; ++i) {
        // 确定窗口范围
        int startIdx = std::max(0, i - windowSize/2);
        int endIdx = std::min(arcSize - 1, i + windowSize/2);
        
        // 计算窗口内的平均偏移
        double sumOffset = 0.0;
        int count = 0;
        for (int j = startIdx; j <= endIdx; ++j) {
            sumOffset += offsets[j];
            count++;
        }
        
        double avgOffset = (count > 0) ? sumOffset / count : 0.0;
        smoothedOffsets.push_back(avgOffset);
    }
    
    return smoothedOffsets;
}

球谐函数与DCB联合估计

系统实现了球谐函数系数和DCB联合估计的方法：

void DCBEstimator::estimateVTECAndDCB(int maxDegree) {
    // 构建联合设计矩阵
    torch::Tensor designMatrix = buildJointDesignMatrix(maxDegree);
    
    // 构建观测向量
    torch::Tensor obsVector = torch::zeros({vtecObs_.size(), 1});
    for (size_t i = 0; i < vtecObs_.size(); ++i) {
        obsVector[i][0] = vtecObs_[i].stec;
    }
    
    // 构建权重矩阵
    torch::Tensor weights = torch::ones({vtecObs_.size()});
    for (size_t i = 0; i < vtecObs_.size(); ++i) {
        // 基于高度角加权
        weights[i] = std::pow(std::sin(vtecObs_[i].elevation * DEG2RAD), 2.0);
    }
    
    // 应用零均值约束
    applyZeroMeanConstraint(designMatrix, obsVector);
    
    // 求解联合参数
    torch::Tensor solution = Utils::svdSolve(designMatrix, obsVector);
    
    // 提取DCB值和球谐系数
    extractDCBAndSHValues(solution);
}

项目结构

vtec/
├── build/                    # 构建目录
├── config/                   # 配置文件
│   └── default_config.ini    # 默认配置文件
├── include/                  # 头文件
│   ├── Common.h              # 通用数据结构定义
│   ├── Config.h              # 配置管理
│   ├── DCBEstimator.h        # DCB估计
│   ├── GPSSatelliteOrbit.h   # 卫星轨道计算
│   ├── IONEXWriter.h         # IONEX文件输出
│   ├── Preprocessor.h        # 数据预处理
│   ├── RinexParser.h         # RINEX文件解析
│   ├── SP3Interpolator.h     # 精密星历插值
│   ├── SphericalHarmonic.h   # 球谐函数模型
│   ├── STECCalculator.h      # STEC计算
│   └── Utils.h               # 工具函数
├── src/                      # 源文件
│   ├── Common.cpp
│   ├── Config.cpp
│   ├── DCBEstimator.cpp
│   ├── GPSSatelliteOrbit.cpp
│   ├── IONEXWriter.cpp
│   ├── main.cpp              # 主程序
│   ├── Preprocessor.cpp
│   ├── RinexParser.cpp
│   ├── SP3Interpolator.cpp
│   ├── SphericalHarmonic.cpp
│   ├── STECCalculator.cpp
│   └── Utils.cpp
├── output/                   # 输出目录
├── doc/                      # 文档
│   ├── class.md              # 类设计文档
│   ├── report.md             # 项目报告
│   └── tech.md               # 技术报告
├── CMakeLists.txt            # CMake构建配置
└── README.md                 # 项目说明

数据可视化

VTEC系统生成的IONEX文件可以使用以下工具进行可视化：

1. GNSS-Toolbox：可读取IONEX文件，生成全球或区域电离层TEC地图。
2. MATLAB IONEX工具：提供了IONEX文件解析和可视化的MATLAB脚本。
3. Python绘图：通过Python的matplotlib和cartopy库可以绘制精美的全球TEC地图。

示例Python脚本：

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import numpy as np
from ionex_reader import read_ionex  # 自定义的IONEX文件读取模块

# 读取IONEX文件
tec_map, lats, lons, time_info = read_ionex('output.ionex')

# 创建全球地图
plt.figure(figsize=(12, 8))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()
ax.gridlines()

# 绘制TEC数据
im = ax.contourf(lons, lats, tec_map, 
                 levels=np.linspace(0, 100, 21),
                 transform=ccrs.PlateCarree(),
                 cmap='jet')

# 添加颜色条和标题
cbar = plt.colorbar(im, orientation='horizontal', pad=0.05)
cbar.set_label('VTEC (TECU)')
plt.title(f'全球VTEC地图 - {time_info}')

plt.savefig('vtec_map.png', dpi=300)
plt.show()

性能指标

VTEC系统达到了以下性能指标：

1. STEC计算精度：

   • 相位对齐STEC精度：优于0.5 TECU
   • 伪距STEC精度：2-5 TECU

2. VTEC建模精度：

   • 中低纬度区域：优于2 TECU
   • 高纬度区域：优于5 TECU
   • 赤道异常区：优于8 TECU

3. DCB估计精度：

   • 卫星DCB：优于0.5 ns
   • 接收机DCB：优于1.0 ns

4. 时空分辨率：

   • 空间分辨率：最高支持1°×1°全球格网
   • 时间分辨率：最高支持15分钟间隔

5. 计算性能：

   • 单次全球建模时间：<30秒(8阶球谐函数)
   • 24小时数据处理时间：<10分钟(单线程)

应用场景

VTEC系统可广泛应用于以下场景：

1. GNSS精密定位：提供高精度电离层延迟改正，应用于RTK、PPP等
2. 空间天气监测：监测电离层扰动、太阳耀斑影响、地磁暴响应
3. 电离层科学研究：研究电离层日变化、季节变化、太阳活动周期影响
4. 导航增强系统：支持SBAS、GBAS等导航增强系统的电离层改正
5. 通信系统：优化HF通信、卫星通信链路设计、雷达系统性能

贡献方式

欢迎通过Pull Request或Issue贡献代码、提出问题或建议。请确保您的代码遵循项目的编码规范，并通过适当的单元测试。

许可证

本项目采用 MIT 许可证。详情请参见 LICENSE 文件。