简介本资源是一款面向GNSS基准站网坐标序列数据处理的MATLAB软件工具包专为计算机、电子信息工程及数学等专业本科生课程设计、毕业设计与科研实践开发解决高精度坐标时间序列建模、噪声分析与趋势提取等核心问题。压缩包共80个文件5.6MB含18个MATLAB源码.m、9个NEU格式坐标序列数据、33张结果可视化PNG图、3个可执行程序.exe及配套HTML说明页、PDF文档与日志文件覆盖数据读取、预处理、参数估计、残差分析与图形输出全流程。代码采用参数化设计关键参数集中配置、注释详尽逻辑清晰支持MATLAB 2014a/2019a/2024a多版本直接运行。附赠真实案例数据集开箱即用便于快速验证算法效果、理解GNSS坐标序列建模原理并为二次开发提供结构化脚本基础。 去年处理某个区域基准站网的坐标序列数据时我被一个情况磨得没脾气解算软件跑完之后N、E、U三个分量的时间序列里面既有离群得离谱的粗差又有断崖式台阶U分量上还叠着稳定的季节性起伏。那会儿如果直接拿最小二乘去拟合速度出来的结果和不确定度都不可信——粗差会把趋势线拉偏未建模的跳变会让残差失去意义而把噪声全当成白噪声来算速度误差更是会把不确定性低估到差不多5到10倍。所以我把整个流程整理成了一套带Matlab实现的GNSS基准站网坐标序列数据处理程序。它能完成从原始坐标序列文件读取、粗差剔除、跳变探测到趋势与周期项拟合、共模误差剔除、速度及不确定度估计的全流程处理。这套程序尤其适合做地壳形变研究、参考框架维持、滑坡与沉降监测这类需要长期坐标时间序列分析的工作也适合刚接触GNSS大地测量的研究生或工程师作为处理框架参考。1. 基准站网坐标序列处理的三个核心痛点先说清楚坐标序列是怎么来的。GNSS连续运行基准站每天观测经过GAMIT/GLOBK、Bernese或者PANDA这类高精度解算软件处理后会得到每个站在ITRF框架下的单日坐标解。把所有单日解按时间串起来就形成了每个站的坐标时间序列通常分开N北、E东、U垂直三个分量来分析。但解算软件给出的原始时间序列在研究尺度上几乎是没法直接用的。这里面的问题不是解算软件不行而是从卫星信号到坐标的整条链路里掺杂了太多与真实地壳运动无关的东西。1.1 粗差从哪里来粗差是指那些明显偏离正常值、没有物理规律可循的异常点。来源五花八门强电离层闪烁导致周跳未被完全修复、天线相位中心模型残差在低高度角时的反常表现、接收机短暂失锁后的伪距异常、甚至数据文件本身在传输过程中丢失了关键字节。这些粗差点在N、E方向通常表现为几毫米到几厘米的野值在U方向上可能更夸张。如果直接带着这些粗差做最小二乘拟合粗差会在目标函数里获得很大权重把趋势项、周期项的估计都带偏。而且粗差往往不是对称分布的正负误差用普通的标准差准则很难处理干净。1.2 跳变是时间序列里的台阶跳变指的是坐标序列在某一个时刻之后整体平移了一个量。常见原因包括天线被更换或发生位移、馈线或接收机硬件变更、解算策略调整。这种台阶在时间序列图上看起来非常典型就是突然往上或往下错开几毫米之后又平稳延续。跳变最麻烦的地方在于它和真实地壳运动的区分。地震引起的同震位移本质上也是一种跳变但是物理含义完全不同。处理时不能一股脑把所有跳变都删掉而是要把它们识别出来在建模时用阶跃函数去表示让跳变和趋势项、周期项一起被估计。1.3 有色噪声是速度不确定度被低估的元凶很多人在做坐标序列分析时有一个惯性思维残差就是随机误差用最小二乘残差的标准差来计算速度不确定度就万事大吉。但在GNSS坐标序列里残差根本不是白噪声而是典型的幂律噪声以闪烁噪声为主部分站还含随机游走噪声。这就意味着简单的白噪声假设会严重低估速度的不确定度。如果忽略时间相关性10年序列的速度不确定度可能被低估5到15倍。这一点在撰写学术论文或做形变解释时尤其致命因为“信号是否显著”完全取决于误差棒算得准不准。2. 第一步必须做对数据读取、时间轴整理与质量报告整套程序的第一步不是调用什么高级算法而是老老实实把数据读进来把时间轴理清楚搞清楚数据能支撑什么分析。这一步做不好后续所有环节都是在垃圾上盖楼。2.1 坐标序列文件的标准字段与Matlab读取方案目前各家机构提供的坐标序列文件格式不完全统一但核心字段差异不大。最常见的形式是类似于SOPAC的tseries文件每一行代表一个历元依次包含年、年积日DOY、月、日、N方向坐标偏离值、E方向坐标偏离值、U方向坐标偏离值有时还带各项标准差。一个典型的数据文件大概长这样2020 001 01 01 -0.42 1.13 3.27 2020 002 01 02 0.08 -0.67 -1.02 2020 003 01 03 -0.16 0.32 0.78 2020 005 01 05 0.25 -0.14 -2.14注意这里跳过了第004天说明数据存在缺失这在真实数据里非常常见。Matlab读取时我建议用textscan一次性把全部列读进来不要用load或csvread因为文件头部往往有若干行注释load处理不了。这里给出一个可用的读取函数骨架function data read_timeseries(filename) fid fopen(filename, r); C textscan(fid, %f %f %f %f %f %f %f %f %f, ... HeaderLines, 1, CommentStyle, #); fclose(fid); data.year C{1}; data.doy C{2}; data.N C{5}; % 单位 mm data.E C{6}; data.U C{7}; % 将年积日转为绝对时间用年为单位便于后续拟合 data.t data.year (data.doy - 1) / 365.25; end实际的读取函数还会处理不同列数、不同单位、NA字符串等情况但核心思路就是先把原文件里的所有信息完整读入再按列提取成结构体后续函数统一操作结构体不要一会儿传数组一会儿传cell容易乱。2.2 时间轴处理DOY转绝对时间、不等间隔与缺失处理坐标序列的时间轴处理里有个很容易被忽略的细节闰年。DOY转小数年时我用365.25作为平均年长精度足够日常分析但如果你需要精确到0.0001年的分辨率建议做一个完整的天数累积把每个历元换算成相对于某个起始时刻的小数年份。不等间隔本身不影响最小二乘拟合因为最小二乘并不要求数据等间隔。真正影响的是谱分析、噪声协方差矩阵构造这些对时间结构敏感的操作。所以程序里我会额外输出一个时间间隔分布向量看看这个站的数据是每天一个点还是存在大量缺数。此外数据缺失严重时比如连续三个月缺测这时候做跳变探测就要格外小心我在第6章会专门讲这个问题。2.3 质量报告程序跑批前先看三件事程序里我在主流程开始前会生成一个简单的质量报告输出三项指标数据完整率、最大连续缺失天数、初步的粗差疑似点数量。完整率用可用历元数除以理论历元数最大连续缺失天数直接决定了跳变探测等窗口类算法的最低窗口要求。这个报告我通常让程序打印到命令行同时写成一个txt方便批量处理多个站时快速筛出问题站。这一步有个很实际的价值几百个站的基准站网不可能每个站都人工目视检查一遍但机器自动报告标记出哪些站数据质量差、哪些站存在疑似跳变处理人员只需要重点查看这些站就够了。3. 核心算法链条粗差剔除、跳变探测与分段策略数据整理完就进入核心处理环节。我的处理顺序是先剔除粗差再识别跳变然后把跳变位置作为阶跃函数的断点连同趋势项、周期项一起放进最小二乘模型里求解。粗差剔除和跳变探测的顺序不能反因为跳变本身会让很多统计准则失效。3.1 用MAD代替3σ为什么标准差在这里不好使很多初学者直接对坐标序列做3σ剔除。3σ的缺陷在于标准差本身对粗差非常敏感——一个巨大的粗差会把均值拉向自己同时把标准差撑大结果就是小粗差混过去了而真正的极端粗差因为离均值太近反而没被识别出来。更稳健的做法是用中位数和MADMedian Absolute Deviation。MAD的定义是MAD median(|x_i - median(x)|)在正态分布假设下MAD乘以1.4826就可以作为标准差的稳健估计。判别准则改为|x_i - median(x)| k * 1.4826 * MADk一般取3到4。用中位数代替均值用MAD代替标准差粗差点无论多极端对中位数和MAD的影响都被大大抑制。实际剔除粗差时还有一个重要策略先做一次初步的趋势和周期拟合在残差序列上做MAD判别而不是直接对原始坐标序列做。因为坐标序列本身含有真实的趋势和季节性信号如果对原始值直接做MAD会把冬季和夏季之间的正常波动当成粗差剔除。3.2 跳变探测移动t检验与台站日志的组合跳变探测我采用了一个相对简单但效果很好的方法移动窗口t检验。思路是在时间序列上滑动一个窗口窗口分成前后两半检验两半的均值是否存在显著差异。如果p值很小就认为这个位置存在跳变。Matlab里可以直接用ttest2完整实现逻辑大概是这样function offsetIdx detect_offset(ts, halfWin) n length(ts); offsetIdx []; i halfWin 1; while i n - halfWin seg1 ts(i-halfWin : i-1); seg2 ts(i : ihalfWin-1); [h, p] ttest2(seg1, seg2); if h p 0.01 offsetIdx [offsetIdx; i]; i i halfWin; % 跳开已被标记的区域避免重复检测 else i i 1; end end end窗口长度一般取30天到90天太短容易对季节性信号敏感太长则可能抹掉真实跳变的检测能力。这里需要说明的是统计方法只能给出候选位置最终判定一定要结合测站日志antenna log、receiver log。如果一个候选跳变位点没有任何设备变更记录但统计上又非常显著那它可能就是真实的地壳运动信号需要结合地震目录或周边站网做联合判断。一旦确认了跳变位置处理方式不是把跳变之后的数据减掉一个常数而是在最小二乘模型中加入阶跃函数把跳变幅度作为未知参数估计出来。这样处理保留了跳变前后各自的真实信息比简单截断数据更科学。3.3 核心建模步骤趋势周期项的加权最小二乘拟合GNSS坐标序列的经典模型可以写成y(t) a bt Σ[A_isin(2πf_it) B_icos(2πf_it)] Σ[C_jH(t - t_j)] res(t)其中a是常数项b是线性速度项f_i通常取1 cpy年周期和2 cpy半年周期t_j是第j个跳变时刻H是Heaviside阶跃函数。在Matlab里构建设计矩阵的代码可以写成function A design_matrix(t, freq, offsetIdx) n length(t); A [ones(n,1), t(:)]; for k 1:length(freq) f freq(k) * 2 * pi; A [A, sin(f*t(:)), cos(f*t(:))]; end for i 1:length(offsetIdx) step zeros(n,1); step(offsetIdx(i):end) 1; A [A, step]; end end注意这里我把t归一化成以“年”为单位的绝对时间这样设计矩阵的条件数会比较合理不至于因为t的数值过大导致矩阵求逆时数值不稳定。拟合时如果数据文件里带了坐标解的标准差就用加权最小二乘如果没带就先用普通最小二乘残差估计出单位权中误差后再对速度和周期项的不确定度做尺度校准。拟合完成后趋势项b对应的是该站在N/E/U方向上的运动速度单位通常是mm/yr。到这一步粗差、跳变和趋势项已经被分离干净了剩下的残差才是后续共模误差和噪声分析的对象。4. 共模误差与噪声模型速度不确定度不能只靠白噪声粗差剔掉、趋势和周期项拟合完之后残差里仍然有很强的空间相干结构。同一区域内不同基准站的坐标残差往往表现出相似的波动形态这就是共模误差。如果不处理共模误差站间相关性将被带走速度不确定度估计和站间基线变化分析都会失真。4.1 PCA/EOF剔除共模误差的正确姿势共模误差剔除最常用的方法是PCA主成分分析/EOF分解。具体做法是把区域内所有站的残差序列排列成一个矩阵行是对应的时间历元列是各个站点。对N、E、U三个分量分别做PCA。以N分量为例写出核心处理流程function residue_new remove_cme(residue_matrix, ncm) % 输入residue_matrix行时间列站点已去均值、去趋势 [U, S, V] svd(residue_matrix, econ); % 取前ncm个主成分重构共模误差 cme_matrix U(:,1:ncm) * S(1:ncm,1:ncm) * V(:,1:ncm); residue_new residue_matrix - cme_matrix; endncm通常取1即第一主成分因为第一主成分在区域网数据完备时通常能解释50%以上的残差方差表现为整个区域同步抬升或水平错动。需要提醒的是做PCA之前必须先把每站的趋势项和周期项扣除干净否则PCA会把这些确定性信号当成主要方差来源吸进第一主成分结果把真实的地壳运动信号也一并当共模误差扣掉了。这个坑在工程应用里经常出现。区域站网空间跨度也是一个边界条件。如果区域范围太大比如超过500公里不同站点受到的水负荷、大气负荷不完全一致PCA提取的第一主成分可能混合了大尺度地球物理负载信号这时用PCA剔除共模误差就有争议。我的经验是用于CME分析的站网最好控制在200公里以内且站点分布均匀。4.2 为什么速度不确定度必须考虑闪烁噪声粗差、跳变、趋势、CME全都处理完之后残差保留了最纯粹的随机噪声部分。GNSS坐标序列的典型噪声特征按照Williams等人在2004年对全球IGS站的分析是白噪声闪烁噪声的组合部分站还混入随机游走成分。闪烁噪声是一种具有长程相关性的噪声它意味着今天残差的大小和几个月前残差的大小不是完全独立的。这种时间相关性会让有效样本量远小于数据点数。如果你用一个简单的白噪声模型来计算速度的标准差相当于假设每天的误差之间互不相关这显然不符合实际。一个粗略的经验数据同样一段10年序列用白噪声模型估计的速度标准差和用白噪声闪烁噪声模型估计的标准差差别往往是5到10倍。在做形变解释时用哪个误差直接影响结论——一个3mm/yr的站速度在纯白噪声下可能显示为显著运动但考虑闪烁噪声后可能并不显著。4.3 实际估计噪声参数的可行方案真正严格的噪声参数估计需要用最大似然估计MLE反演噪声模型。最完整的做法是同时估计白噪声、闪烁噪声、随机游走三种噪声的振幅以及谱指数。但自己从头实现闪烁噪声协方差矩阵并不是一件轻松的事尤其是当序列长度超过数千个历元时协方差矩阵的求逆计算量很大容易陷入数值不稳定。业界更常见的做法是使用Hector软件或CATS软件来完成噪声估计。Hector是开源的支持多分量同时估计内存管理也做得比较好。在Matlab程序里我的设计是把拟合与噪声估计解耦Matlab负责粗差剔除、跳变建模、趋势/周期拟合等前期工作把残差序列和协变量信息导出成Hector需要的输入格式然后调用Hector完成噪声模型反演最后把速度不确定度读回Matlab统一汇总。如果你的环境不具备安装外部软件的条件至少应该了解一个替代思路用bootstrap方法对残差序列进行重采样估计速度参数的采样分布。虽然理论上不如MLE严谨但比纯白噪声假设已经好很多而且实现难度低很多。5. Matlab代码框架从主函数到每个子模块前面讲了原理这部分说一下程序本身是怎么组织的。整个项目采用“主脚本功能函数”的结构主脚本负责流程编排功能函数负责独立模块实现。这样每个站只能跑一个脚本换数据时不用改代码批量处理时也方便并行。5.1 主处理脚本的工作流程主脚本main_processing.m的执行顺序非常直观% 1. 读取原始坐标序列 data read_timeseries(GNSS_station.txt); % 2. 粗差剔除默认MADk4 data remove_outliers(data, MAD, 4); % 3. 跳变探测窗口60天 offsetIdx detect_offset(data.U, 30); offsetIdx [offsetIdx; detect_offset(data.N, 30)]; % 4. 拟合趋势年/半年周期阶跃项 [params, A] fit_trend_seasonal(data, offsetIdx); % 5. 残差序列做PCA共模误差剔除批量多站时 residue_clean remove_cme(residue_matrix, 1); % 6. 导出结果与绘图 plot_timeseries(data, params, offsetIdx); export_velocity(velocity_output.txt, params);这里面每步的输入输出都是设计时想清楚的read_timeseries返回结构体data包含原始坐标和时间轴remove_outliers返回清洗后的data并附带一个log记录被剔除的历元fit_trend_seasonal返回参数向量包括速度、周期振幅、跳变幅度。后面无论加多少新功能都不需要改前面的函数接口。5.2 关键函数实现细节不仅仅是跑通design_matrix函数我前面已经给出了另外两个关键函数也说一下设计出发点。detect_offset函数里有一个很容易踩的细节当检测到跳变后代码要跳开一个窗口再继续否则一次跳变会在跳变点附近连续触发很多次检测。如果同时检测三个分量三个分量检测到的跳变位置可能不完全一致我的做法是做一次合并聚类把在10天以内的候选跳变归并为一个事件再结合测站日志做最终确认。粗差剔除函数remove_outliers的输入输出设计里我建议把被剔除点的信息返回出来而不是原地修改数据。原因是粗差剔除是一个有争议的决策保留log可以让后续分析阶段随时复查而不是数据一经处理就再也追溯不回来。PCA移除CME的函数里有一个细节不同方向N/E/U的残差应该各自独立做PCA而不是合在一起。因为这三个分量的物理噪声来源和幅度完全不同合在一起等于给U分量加了过大的权重。5.3 可视化输出从单站时序图到速度场程序里我也加入了可视化模块。单站输出是一张包含N、E、U三个子图的时序图原始数据画成灰点粗差被标记为红色圆圈跳变位置画一条竖直虚线拟合的趋势和周期信号叠加成一条实线残差画在最下面。这张图虽然简单但信息密度很高一个人几秒钟内就能判断一个站的处理质量是不是过关。批量处理完所有站后汇总速度结果并输出一个速度场文件格式大致是站名 经度 纬度 VN(mm/yr) VE(mm/yr) VU(mm/yr) sig_N sig_E sig_U AAA 116.13 39.52 1.23 -0.67 -2.11 0.15 0.18 0.32 BBB 117.14 40.11 -0.28 1.04 -1.33 0.14 0.16 0.28速度场文件可以直接导入GMT或Python的cartopy里绘制速度场箭头图这是做地壳形变研究最常用的成果图之一。6. 拿真实数据处理时踩过的坑处理真实数据的过程永远比原理复杂下面这几个坑是我在实际使用中反复踩过的都很有代表性。6.1 大地震的同震位移差点被当成粗差删掉有一年处理某区域站点数据时程序在某个日期附近标记出了大量粗差。我打开原始时序一看那个日期正是几百公里外一次强震发生的时间。同震位移在时间序列上确实表现为一个非常突兀的台阶从统计角度和粗差难以区分但物理含义完全不同——这是真实的构造信号必须保留并作为跳变建模拟合而不是当作粗差剔除。从那以后我在程序里加入了一个“人工复核”步骤被标记的粗差如果连续多站同时出现且集中在某个日期区间程序会弹出一个警告提示操作者检查地震目录或测站日志而不是静悄悄地剔除。6.2 连续缺失复测后的“假跳变”某个站因为设备供电问题连续停了40天恢复观测后数据整体似乎有微小偏移。移动t检验把复测点识别为跳变但这其实是错觉——因为缺失前后时间间隔太大而季节性信号在缺失区间内正常推进导致前后两段数据的平均值存在差异。这个问题的本质是任何基于均值的跳变探测方法都无法区分“真实跳变”和“数据缺口造成的相位差”。解决办法不复杂对每个候选跳变点先检查附近是否存在长于30天的数据缺失。如果存在需要谨慎对待最好结合log文件人工判断而不是盲目接受统计结果。6.3 短序列端点相位失真只有一年半载却要提取年周期和速度时端点处的拟合往往不稳定。因为正弦函数的拟合依赖于足够的完整周期来约束相位和振幅序列太短时端点的微小噪声会明显影响周期项的相位估计。如果你手里的序列不足两年我会建议只估计趋势项半年周期甚至只估计趋势项。否则拟合出的年周期振幅很可能是噪声驱动的还顺带污染速度项的估计。6.4 区域跨度太大时PCA吞掉真实信号之前把全省上百个站直接拿来跑PCA第一主成分出来之后N方向的残差被扣掉了一大块结果发现几个位于不同构造块体上的站似乎都不动了。后来意识到这些站跨越了不止一个地质单元地壳运动本来就存在区域差异PCA的“区域共模误差”在这个尺度上已经不再适用它把不同块体的真实相对运动也当成共模误差一起扣了。解决方式很简单要么按构造单元或地理距离把站网分组后再做PCA要么只对站间残差相关性高的区域做CME剔除。处理完对比速度和残差标准差效果明显不同。说一点我自己的体会。坐标序列处理真正花时间的地方从来不是实现某个统计过程而是把数据质量控制做扎实——粗差有没有剔除干净、跳变有没有正确建模、残差是否还存在明显的系统结构这些环节直接决定最后得到的速度值靠不靠谱。我每次拿到一批新站数据都会先花十来分钟随机挑几个站看原始时序图再让程序跑批用可视化和自动报告互相印证。这个习惯帮我避免了很多次“程序正常跑完但结果明显不对”的事情。做GNSS时间序列分析越是在前期把账算清楚后面做物理解释时才越踏实。本文还有配套的精品资源点击获取 GNSS基准站网坐标序列处理：Matlab实现粗差剔除、跳变探测与速度估计
简介本资源是一款面向GNSS基准站网坐标序列数据处理的MATLAB软件工具包专为计算机、电子信息工程及数学等专业本科生课程设计、毕业设计与科研实践开发解决高精度坐标时间序列建模、噪声分析与趋势提取等核心问题。压缩包共80个文件5.6MB含18个MATLAB源码.m、9个NEU格式坐标序列数据、33张结果可视化PNG图、3个可执行程序.exe及配套HTML说明页、PDF文档与日志文件覆盖数据读取、预处理、参数估计、残差分析与图形输出全流程。代码采用参数化设计关键参数集中配置、注释详尽逻辑清晰支持MATLAB 2014a/2019a/2024a多版本直接运行。附赠真实案例数据集开箱即用便于快速验证算法效果、理解GNSS坐标序列建模原理并为二次开发提供结构化脚本基础。 去年处理某个区域基准站网的坐标序列数据时我被一个情况磨得没脾气解算软件跑完之后N、E、U三个分量的时间序列里面既有离群得离谱的粗差又有断崖式台阶U分量上还叠着稳定的季节性起伏。那会儿如果直接拿最小二乘去拟合速度出来的结果和不确定度都不可信——粗差会把趋势线拉偏未建模的跳变会让残差失去意义而把噪声全当成白噪声来算速度误差更是会把不确定性低估到差不多5到10倍。所以我把整个流程整理成了一套带Matlab实现的GNSS基准站网坐标序列数据处理程序。它能完成从原始坐标序列文件读取、粗差剔除、跳变探测到趋势与周期项拟合、共模误差剔除、速度及不确定度估计的全流程处理。这套程序尤其适合做地壳形变研究、参考框架维持、滑坡与沉降监测这类需要长期坐标时间序列分析的工作也适合刚接触GNSS大地测量的研究生或工程师作为处理框架参考。1. 基准站网坐标序列处理的三个核心痛点先说清楚坐标序列是怎么来的。GNSS连续运行基准站每天观测经过GAMIT/GLOBK、Bernese或者PANDA这类高精度解算软件处理后会得到每个站在ITRF框架下的单日坐标解。把所有单日解按时间串起来就形成了每个站的坐标时间序列通常分开N北、E东、U垂直三个分量来分析。但解算软件给出的原始时间序列在研究尺度上几乎是没法直接用的。这里面的问题不是解算软件不行而是从卫星信号到坐标的整条链路里掺杂了太多与真实地壳运动无关的东西。1.1 粗差从哪里来粗差是指那些明显偏离正常值、没有物理规律可循的异常点。来源五花八门强电离层闪烁导致周跳未被完全修复、天线相位中心模型残差在低高度角时的反常表现、接收机短暂失锁后的伪距异常、甚至数据文件本身在传输过程中丢失了关键字节。这些粗差点在N、E方向通常表现为几毫米到几厘米的野值在U方向上可能更夸张。如果直接带着这些粗差做最小二乘拟合粗差会在目标函数里获得很大权重把趋势项、周期项的估计都带偏。而且粗差往往不是对称分布的正负误差用普通的标准差准则很难处理干净。1.2 跳变是时间序列里的台阶跳变指的是坐标序列在某一个时刻之后整体平移了一个量。常见原因包括天线被更换或发生位移、馈线或接收机硬件变更、解算策略调整。这种台阶在时间序列图上看起来非常典型就是突然往上或往下错开几毫米之后又平稳延续。跳变最麻烦的地方在于它和真实地壳运动的区分。地震引起的同震位移本质上也是一种跳变但是物理含义完全不同。处理时不能一股脑把所有跳变都删掉而是要把它们识别出来在建模时用阶跃函数去表示让跳变和趋势项、周期项一起被估计。1.3 有色噪声是速度不确定度被低估的元凶很多人在做坐标序列分析时有一个惯性思维残差就是随机误差用最小二乘残差的标准差来计算速度不确定度就万事大吉。但在GNSS坐标序列里残差根本不是白噪声而是典型的幂律噪声以闪烁噪声为主部分站还含随机游走噪声。这就意味着简单的白噪声假设会严重低估速度的不确定度。如果忽略时间相关性10年序列的速度不确定度可能被低估5到15倍。这一点在撰写学术论文或做形变解释时尤其致命因为“信号是否显著”完全取决于误差棒算得准不准。2. 第一步必须做对数据读取、时间轴整理与质量报告整套程序的第一步不是调用什么高级算法而是老老实实把数据读进来把时间轴理清楚搞清楚数据能支撑什么分析。这一步做不好后续所有环节都是在垃圾上盖楼。2.1 坐标序列文件的标准字段与Matlab读取方案目前各家机构提供的坐标序列文件格式不完全统一但核心字段差异不大。最常见的形式是类似于SOPAC的tseries文件每一行代表一个历元依次包含年、年积日DOY、月、日、N方向坐标偏离值、E方向坐标偏离值、U方向坐标偏离值有时还带各项标准差。一个典型的数据文件大概长这样2020 001 01 01 -0.42 1.13 3.27 2020 002 01 02 0.08 -0.67 -1.02 2020 003 01 03 -0.16 0.32 0.78 2020 005 01 05 0.25 -0.14 -2.14注意这里跳过了第004天说明数据存在缺失这在真实数据里非常常见。Matlab读取时我建议用textscan一次性把全部列读进来不要用load或csvread因为文件头部往往有若干行注释load处理不了。这里给出一个可用的读取函数骨架function data read_timeseries(filename) fid fopen(filename, r); C textscan(fid, %f %f %f %f %f %f %f %f %f, ... HeaderLines, 1, CommentStyle, #); fclose(fid); data.year C{1}; data.doy C{2}; data.N C{5}; % 单位 mm data.E C{6}; data.U C{7}; % 将年积日转为绝对时间用年为单位便于后续拟合 data.t data.year (data.doy - 1) / 365.25; end实际的读取函数还会处理不同列数、不同单位、NA字符串等情况但核心思路就是先把原文件里的所有信息完整读入再按列提取成结构体后续函数统一操作结构体不要一会儿传数组一会儿传cell容易乱。2.2 时间轴处理DOY转绝对时间、不等间隔与缺失处理坐标序列的时间轴处理里有个很容易被忽略的细节闰年。DOY转小数年时我用365.25作为平均年长精度足够日常分析但如果你需要精确到0.0001年的分辨率建议做一个完整的天数累积把每个历元换算成相对于某个起始时刻的小数年份。不等间隔本身不影响最小二乘拟合因为最小二乘并不要求数据等间隔。真正影响的是谱分析、噪声协方差矩阵构造这些对时间结构敏感的操作。所以程序里我会额外输出一个时间间隔分布向量看看这个站的数据是每天一个点还是存在大量缺数。此外数据缺失严重时比如连续三个月缺测这时候做跳变探测就要格外小心我在第6章会专门讲这个问题。2.3 质量报告程序跑批前先看三件事程序里我在主流程开始前会生成一个简单的质量报告输出三项指标数据完整率、最大连续缺失天数、初步的粗差疑似点数量。完整率用可用历元数除以理论历元数最大连续缺失天数直接决定了跳变探测等窗口类算法的最低窗口要求。这个报告我通常让程序打印到命令行同时写成一个txt方便批量处理多个站时快速筛出问题站。这一步有个很实际的价值几百个站的基准站网不可能每个站都人工目视检查一遍但机器自动报告标记出哪些站数据质量差、哪些站存在疑似跳变处理人员只需要重点查看这些站就够了。3. 核心算法链条粗差剔除、跳变探测与分段策略数据整理完就进入核心处理环节。我的处理顺序是先剔除粗差再识别跳变然后把跳变位置作为阶跃函数的断点连同趋势项、周期项一起放进最小二乘模型里求解。粗差剔除和跳变探测的顺序不能反因为跳变本身会让很多统计准则失效。3.1 用MAD代替3σ为什么标准差在这里不好使很多初学者直接对坐标序列做3σ剔除。3σ的缺陷在于标准差本身对粗差非常敏感——一个巨大的粗差会把均值拉向自己同时把标准差撑大结果就是小粗差混过去了而真正的极端粗差因为离均值太近反而没被识别出来。更稳健的做法是用中位数和MADMedian Absolute Deviation。MAD的定义是MAD median(|x_i - median(x)|)在正态分布假设下MAD乘以1.4826就可以作为标准差的稳健估计。判别准则改为|x_i - median(x)| k * 1.4826 * MADk一般取3到4。用中位数代替均值用MAD代替标准差粗差点无论多极端对中位数和MAD的影响都被大大抑制。实际剔除粗差时还有一个重要策略先做一次初步的趋势和周期拟合在残差序列上做MAD判别而不是直接对原始坐标序列做。因为坐标序列本身含有真实的趋势和季节性信号如果对原始值直接做MAD会把冬季和夏季之间的正常波动当成粗差剔除。3.2 跳变探测移动t检验与台站日志的组合跳变探测我采用了一个相对简单但效果很好的方法移动窗口t检验。思路是在时间序列上滑动一个窗口窗口分成前后两半检验两半的均值是否存在显著差异。如果p值很小就认为这个位置存在跳变。Matlab里可以直接用ttest2完整实现逻辑大概是这样function offsetIdx detect_offset(ts, halfWin) n length(ts); offsetIdx []; i halfWin 1; while i n - halfWin seg1 ts(i-halfWin : i-1); seg2 ts(i : ihalfWin-1); [h, p] ttest2(seg1, seg2); if h p 0.01 offsetIdx [offsetIdx; i]; i i halfWin; % 跳开已被标记的区域避免重复检测 else i i 1; end end end窗口长度一般取30天到90天太短容易对季节性信号敏感太长则可能抹掉真实跳变的检测能力。这里需要说明的是统计方法只能给出候选位置最终判定一定要结合测站日志antenna log、receiver log。如果一个候选跳变位点没有任何设备变更记录但统计上又非常显著那它可能就是真实的地壳运动信号需要结合地震目录或周边站网做联合判断。一旦确认了跳变位置处理方式不是把跳变之后的数据减掉一个常数而是在最小二乘模型中加入阶跃函数把跳变幅度作为未知参数估计出来。这样处理保留了跳变前后各自的真实信息比简单截断数据更科学。3.3 核心建模步骤趋势周期项的加权最小二乘拟合GNSS坐标序列的经典模型可以写成y(t) a bt Σ[A_isin(2πf_it) B_icos(2πf_it)] Σ[C_jH(t - t_j)] res(t)其中a是常数项b是线性速度项f_i通常取1 cpy年周期和2 cpy半年周期t_j是第j个跳变时刻H是Heaviside阶跃函数。在Matlab里构建设计矩阵的代码可以写成function A design_matrix(t, freq, offsetIdx) n length(t); A [ones(n,1), t(:)]; for k 1:length(freq) f freq(k) * 2 * pi; A [A, sin(f*t(:)), cos(f*t(:))]; end for i 1:length(offsetIdx) step zeros(n,1); step(offsetIdx(i):end) 1; A [A, step]; end end注意这里我把t归一化成以“年”为单位的绝对时间这样设计矩阵的条件数会比较合理不至于因为t的数值过大导致矩阵求逆时数值不稳定。拟合时如果数据文件里带了坐标解的标准差就用加权最小二乘如果没带就先用普通最小二乘残差估计出单位权中误差后再对速度和周期项的不确定度做尺度校准。拟合完成后趋势项b对应的是该站在N/E/U方向上的运动速度单位通常是mm/yr。到这一步粗差、跳变和趋势项已经被分离干净了剩下的残差才是后续共模误差和噪声分析的对象。4. 共模误差与噪声模型速度不确定度不能只靠白噪声粗差剔掉、趋势和周期项拟合完之后残差里仍然有很强的空间相干结构。同一区域内不同基准站的坐标残差往往表现出相似的波动形态这就是共模误差。如果不处理共模误差站间相关性将被带走速度不确定度估计和站间基线变化分析都会失真。4.1 PCA/EOF剔除共模误差的正确姿势共模误差剔除最常用的方法是PCA主成分分析/EOF分解。具体做法是把区域内所有站的残差序列排列成一个矩阵行是对应的时间历元列是各个站点。对N、E、U三个分量分别做PCA。以N分量为例写出核心处理流程function residue_new remove_cme(residue_matrix, ncm) % 输入residue_matrix行时间列站点已去均值、去趋势 [U, S, V] svd(residue_matrix, econ); % 取前ncm个主成分重构共模误差 cme_matrix U(:,1:ncm) * S(1:ncm,1:ncm) * V(:,1:ncm); residue_new residue_matrix - cme_matrix; endncm通常取1即第一主成分因为第一主成分在区域网数据完备时通常能解释50%以上的残差方差表现为整个区域同步抬升或水平错动。需要提醒的是做PCA之前必须先把每站的趋势项和周期项扣除干净否则PCA会把这些确定性信号当成主要方差来源吸进第一主成分结果把真实的地壳运动信号也一并当共模误差扣掉了。这个坑在工程应用里经常出现。区域站网空间跨度也是一个边界条件。如果区域范围太大比如超过500公里不同站点受到的水负荷、大气负荷不完全一致PCA提取的第一主成分可能混合了大尺度地球物理负载信号这时用PCA剔除共模误差就有争议。我的经验是用于CME分析的站网最好控制在200公里以内且站点分布均匀。4.2 为什么速度不确定度必须考虑闪烁噪声粗差、跳变、趋势、CME全都处理完之后残差保留了最纯粹的随机噪声部分。GNSS坐标序列的典型噪声特征按照Williams等人在2004年对全球IGS站的分析是白噪声闪烁噪声的组合部分站还混入随机游走成分。闪烁噪声是一种具有长程相关性的噪声它意味着今天残差的大小和几个月前残差的大小不是完全独立的。这种时间相关性会让有效样本量远小于数据点数。如果你用一个简单的白噪声模型来计算速度的标准差相当于假设每天的误差之间互不相关这显然不符合实际。一个粗略的经验数据同样一段10年序列用白噪声模型估计的速度标准差和用白噪声闪烁噪声模型估计的标准差差别往往是5到10倍。在做形变解释时用哪个误差直接影响结论——一个3mm/yr的站速度在纯白噪声下可能显示为显著运动但考虑闪烁噪声后可能并不显著。4.3 实际估计噪声参数的可行方案真正严格的噪声参数估计需要用最大似然估计MLE反演噪声模型。最完整的做法是同时估计白噪声、闪烁噪声、随机游走三种噪声的振幅以及谱指数。但自己从头实现闪烁噪声协方差矩阵并不是一件轻松的事尤其是当序列长度超过数千个历元时协方差矩阵的求逆计算量很大容易陷入数值不稳定。业界更常见的做法是使用Hector软件或CATS软件来完成噪声估计。Hector是开源的支持多分量同时估计内存管理也做得比较好。在Matlab程序里我的设计是把拟合与噪声估计解耦Matlab负责粗差剔除、跳变建模、趋势/周期拟合等前期工作把残差序列和协变量信息导出成Hector需要的输入格式然后调用Hector完成噪声模型反演最后把速度不确定度读回Matlab统一汇总。如果你的环境不具备安装外部软件的条件至少应该了解一个替代思路用bootstrap方法对残差序列进行重采样估计速度参数的采样分布。虽然理论上不如MLE严谨但比纯白噪声假设已经好很多而且实现难度低很多。5. Matlab代码框架从主函数到每个子模块前面讲了原理这部分说一下程序本身是怎么组织的。整个项目采用“主脚本功能函数”的结构主脚本负责流程编排功能函数负责独立模块实现。这样每个站只能跑一个脚本换数据时不用改代码批量处理时也方便并行。5.1 主处理脚本的工作流程主脚本main_processing.m的执行顺序非常直观% 1. 读取原始坐标序列 data read_timeseries(GNSS_station.txt); % 2. 粗差剔除默认MADk4 data remove_outliers(data, MAD, 4); % 3. 跳变探测窗口60天 offsetIdx detect_offset(data.U, 30); offsetIdx [offsetIdx; detect_offset(data.N, 30)]; % 4. 拟合趋势年/半年周期阶跃项 [params, A] fit_trend_seasonal(data, offsetIdx); % 5. 残差序列做PCA共模误差剔除批量多站时 residue_clean remove_cme(residue_matrix, 1); % 6. 导出结果与绘图 plot_timeseries(data, params, offsetIdx); export_velocity(velocity_output.txt, params);这里面每步的输入输出都是设计时想清楚的read_timeseries返回结构体data包含原始坐标和时间轴remove_outliers返回清洗后的data并附带一个log记录被剔除的历元fit_trend_seasonal返回参数向量包括速度、周期振幅、跳变幅度。后面无论加多少新功能都不需要改前面的函数接口。5.2 关键函数实现细节不仅仅是跑通design_matrix函数我前面已经给出了另外两个关键函数也说一下设计出发点。detect_offset函数里有一个很容易踩的细节当检测到跳变后代码要跳开一个窗口再继续否则一次跳变会在跳变点附近连续触发很多次检测。如果同时检测三个分量三个分量检测到的跳变位置可能不完全一致我的做法是做一次合并聚类把在10天以内的候选跳变归并为一个事件再结合测站日志做最终确认。粗差剔除函数remove_outliers的输入输出设计里我建议把被剔除点的信息返回出来而不是原地修改数据。原因是粗差剔除是一个有争议的决策保留log可以让后续分析阶段随时复查而不是数据一经处理就再也追溯不回来。PCA移除CME的函数里有一个细节不同方向N/E/U的残差应该各自独立做PCA而不是合在一起。因为这三个分量的物理噪声来源和幅度完全不同合在一起等于给U分量加了过大的权重。5.3 可视化输出从单站时序图到速度场程序里我也加入了可视化模块。单站输出是一张包含N、E、U三个子图的时序图原始数据画成灰点粗差被标记为红色圆圈跳变位置画一条竖直虚线拟合的趋势和周期信号叠加成一条实线残差画在最下面。这张图虽然简单但信息密度很高一个人几秒钟内就能判断一个站的处理质量是不是过关。批量处理完所有站后汇总速度结果并输出一个速度场文件格式大致是站名 经度 纬度 VN(mm/yr) VE(mm/yr) VU(mm/yr) sig_N sig_E sig_U AAA 116.13 39.52 1.23 -0.67 -2.11 0.15 0.18 0.32 BBB 117.14 40.11 -0.28 1.04 -1.33 0.14 0.16 0.28速度场文件可以直接导入GMT或Python的cartopy里绘制速度场箭头图这是做地壳形变研究最常用的成果图之一。6. 拿真实数据处理时踩过的坑处理真实数据的过程永远比原理复杂下面这几个坑是我在实际使用中反复踩过的都很有代表性。6.1 大地震的同震位移差点被当成粗差删掉有一年处理某区域站点数据时程序在某个日期附近标记出了大量粗差。我打开原始时序一看那个日期正是几百公里外一次强震发生的时间。同震位移在时间序列上确实表现为一个非常突兀的台阶从统计角度和粗差难以区分但物理含义完全不同——这是真实的构造信号必须保留并作为跳变建模拟合而不是当作粗差剔除。从那以后我在程序里加入了一个“人工复核”步骤被标记的粗差如果连续多站同时出现且集中在某个日期区间程序会弹出一个警告提示操作者检查地震目录或测站日志而不是静悄悄地剔除。6.2 连续缺失复测后的“假跳变”某个站因为设备供电问题连续停了40天恢复观测后数据整体似乎有微小偏移。移动t检验把复测点识别为跳变但这其实是错觉——因为缺失前后时间间隔太大而季节性信号在缺失区间内正常推进导致前后两段数据的平均值存在差异。这个问题的本质是任何基于均值的跳变探测方法都无法区分“真实跳变”和“数据缺口造成的相位差”。解决办法不复杂对每个候选跳变点先检查附近是否存在长于30天的数据缺失。如果存在需要谨慎对待最好结合log文件人工判断而不是盲目接受统计结果。6.3 短序列端点相位失真只有一年半载却要提取年周期和速度时端点处的拟合往往不稳定。因为正弦函数的拟合依赖于足够的完整周期来约束相位和振幅序列太短时端点的微小噪声会明显影响周期项的相位估计。如果你手里的序列不足两年我会建议只估计趋势项半年周期甚至只估计趋势项。否则拟合出的年周期振幅很可能是噪声驱动的还顺带污染速度项的估计。6.4 区域跨度太大时PCA吞掉真实信号之前把全省上百个站直接拿来跑PCA第一主成分出来之后N方向的残差被扣掉了一大块结果发现几个位于不同构造块体上的站似乎都不动了。后来意识到这些站跨越了不止一个地质单元地壳运动本来就存在区域差异PCA的“区域共模误差”在这个尺度上已经不再适用它把不同块体的真实相对运动也当成共模误差一起扣了。解决方式很简单要么按构造单元或地理距离把站网分组后再做PCA要么只对站间残差相关性高的区域做CME剔除。处理完对比速度和残差标准差效果明显不同。说一点我自己的体会。坐标序列处理真正花时间的地方从来不是实现某个统计过程而是把数据质量控制做扎实——粗差有没有剔除干净、跳变有没有正确建模、残差是否还存在明显的系统结构这些环节直接决定最后得到的速度值靠不靠谱。我每次拿到一批新站数据都会先花十来分钟随机挑几个站看原始时序图再让程序跑批用可视化和自动报告互相印证。这个习惯帮我避免了很多次“程序正常跑完但结果明显不对”的事情。做GNSS时间序列分析越是在前期把账算清楚后面做物理解释时才越踏实。本文还有配套的精品资源点击获取