时间序列分析软件Hector用户手册（一）
时间序列分析软件Hector用户手册（一）
- 1 介绍
- 1.1 引用
- 1.2 主要特点
- 2 安装
- 3 Hector惯例
- 4 示例
- 4.1 示例1：带有异常值和跳变的模拟GNSS时间序列
- 4.1.1 剔除异常值
- 4.1.2 估计线性趋势
- 4.1.3 绘制功率谱密度图
- 4.2 示例2：Cascais的月均PSMSL验潮站数据
- 4.3 示例3：模拟有色噪声
- 4.4 示例4：震后回弹
- 4.5 示例5：自动分析+变化的季节信号
翻译自Hector用户手册，果然对我来说还是中文方便一些，有增删和修正，此外有一些不太会翻译就保留为英文。
1 介绍
Hector可以估计一个时间序列的趋势，尤其是在噪声时域相关的情况下。趋势估计在地球物理研究中应用广泛：温度、海面高、GNSS站点位置变化等等。Hector可以估计线性趋势或者高阶多项式趋势，还可以估计周期信号、offset（偏移？）、震后形变等等，使用的估计方法是极大似然估计。除了Hector，类似软件还有CATS。实际上，Hector就是快速版的CATS，之所以更快是因为Hector只接受平稳噪声，因而允许快速矩阵操作，至于非平稳噪声则用平稳噪声近似。另一个替代软件是est_noise，同样可以处理缺失数据，区别在于用了不同的方式建立协方差矩阵。
1.1 引用
Bos, M. S., Fernandes, R. M. S., Williams, S. D. P., and Bastos, L. (2013). Fast Error Analysis of Continuous GNSS Observations with Missing Data.J. Geod., Vol. 87(4), 351–360, doi:10.1007/s00190-012-0605-0.
1.2 主要特点
- 正确处理缺失数据。无需插值或者补零，也不用对协方差矩阵做近似（只需噪声平稳）
- 在估计线性趋势的过程中，允许包含周年、半周年和其他周期的信号
- 允许在给定时刻估计offset
- 包含幂律（power-law）噪声、ARFIMA、广义高斯-马尔科夫（Generalized Gauss-Markov，GGM）和白噪声模型等，以及这些噪声的任意组合
- 具有去除异常值，绘制功率谱密度图、模拟合成有色噪声的功能
- 可以自动探测offset
- 提供脚本以估计振幅变化的季节信号
2 安装
Hector可以运行在类Unix系统上，下载地址提供源码以及Ubuntu16.04和CentOS 5上编译好的版本，建议将二进制文件放在目录/usr/local/bin下，其他目录也行，只需注意将所在路径添加到环境变量中，该软件提供的程序如表1所示。
表1. Hector软件包提供的程序列表
| 名称 | 描述 |
|---|---|
| estimatetrend | 估计线性趋势 |
| estimatespectrum | 利用Welch周期图方法从数据或者残差中估计功率谱密度 |
| modelspectrum | 给定噪声模型和噪声参数，计算给定频率范围的功率谱密度 |
| findoffset | 寻找时间序列中可能存在的偏移 |
| simulatenoise | 模拟合成有色噪声 |
| date2mjd | 公历转换为简化儒略日（Modified Julian Date，MJD) |
| mjd2date | 简化儒略日转换为公历 |
如果是自己从源码编译安装的话需要GSL、LAPACK、BLAS和FFTW库。为了运行示例还需要安装python语言和gnuplot绘图软件。
3 Hector惯例
Hector包含一系列命令行程序：异常值剔除、趋势估计、功率谱密度估计、模拟有色噪声、寻找offset。每个程序有自己的控制文件，用以给定输入参数。比如程序removeoutliers的控制文件是removeoutliers.ctl，其他程序遵循同样的命名惯例。控制文件每一行都是关键字+参数值。
图1是地学时间序列分析步骤的流程图。观测值存放在./ori_files路径下，第一步是转换为Hector可识别的若干种格式，为了python脚本的最佳使用效果，推荐MOM格式（MJD、Observations、Model），这些文件存放在./raw_files，如果观测值分为东分量、北分量和垂直分量，那么可以在站点名称后添加_0、_1、_2分别代表各个分量，这有助于python脚本自动探测offset；对./raw_files中的文件探测可能的offset，然后拷贝到./obs_files，并在文件头附加可能的偏移所在日期的相关信息。
图1. 时间序列分析步骤。方框代表每个步骤采用的标准文件命名方式，红色部分是python脚本，_[012]后缀仅适用于含有三个分量（北向分量、东向分量、垂直分量）的时间序列
第二步是寻找异常值，预处理的结果将存放在./pre_files。实际上，异常值剔除和偏移探测可以迭代进行，一旦循环结束，就可以估计趋势，结果文件存放在./mom_files。用户可自行更改流程图中的所有文件路径，但是python脚本默认该结构，所以不推荐更改。所有目录惯例名称参考表2。
表2. 文件路径命名惯例
| 步骤 | 文件路径 | 描述 |
|---|---|---|
| 0 | ./ori_files | 原始格式的时间序列（非MOM格式） |
| 1 | ./raw_files | MOM格式，含有异常值，没有offset、break（偏移跳变？）等文件头信息 |
| 2 | ./obs_files | MOM格式，含有异常值，文件头附加偏移跳变等信息 |
| 3 | ./pre_files | 剔除完异常值的预处理观测值 |
| 4 | ./mom_files | 预处理文件的基础上增加了一列估计值（estimatetrend的输出） |
| 5 | ./sea_files | 估计得到的振幅变化的季节信号 |
| 6 | ./fil_files | 预处理观测值减去估计得振幅变化的季节信号 |
最后，如果需要结果可视化，可以绘制拟合模型的时间序列；如果需要查看噪声模型与实际数据的吻合程度，可以绘制残差和估计得噪声模型的功率谱密度图。
4 示例
所有代码和数据都在Hector网址上可以下载。
4.1 示例1：带有异常值和跳变的模拟GNSS时间序列
ASCII文件TEST.mom存放在目录ex1/obs_files，查看文件可以看到前几行如下：
# sampling period 1.0
# offset 50284.0
# offset 50334.0
# offset 50784.0
# offset 51034.0
第一行声明采样周期为天，其他几行声明需要估计的offset所在的时刻，格式为MJD，MJD对于绘图比较方便，用户可以使用程序date2mjd 和date2mjd在年/月/日/时/分/秒和MJD之间转换。文件头之后就是数据。
50084.0 -17.88951
50085.0 -16.88599
50086.0 -16.84916
...
利用gnuplot命令可产生图片：
set terminal png
set output 'myplot.png'
set xlabel 'years';
set ylabel 'mm';
plot './ex1/obs_files/TEST.mom' using (($1-51544)/365.25+2000):2 with lines t 'observed';
图2. 原始观测数据
4.1.1 剔除异常值
运行程序removeoutliers需要removeoutliers.ctl文件，其内容如下：
DataFile TEST.mom
DataDirectory ./obs_files
OutputFile ./pre_files/TEST.mom
interpolate no
seasonalsignal yes
halfseasonalsignal no
estimateoffsets yes
IQ_factor 3.0
PhysicalUnit mm
ScaleFactor 1.0
关键字ScaleFactor给定尺度变换因子，属于可选参数，默认1.0，其他大部分参数都很好懂，不予翻译，注意行与行之间的参数可以任意交换顺序。运行命令除了输出文件./pre_files/TEST.mom，还会产生outliers.out文件，包含剔除点的MJD。如果需要用其他控制文件可以运行命令“removeoutliers othercontrolfile.ctl”
直接运行程序“removeoutliers”结果如下。可以看到removeoutliers的原理是对原始数据进行普通最小二乘拟合得到线性趋势，然后从原始数据中减去该趋势得到残差，将这些残差排序得到四分位距（即75%处的值减去25%处的值），以中位数为参考，任何超出四分位距3倍的视为异常值，剔除异常值并继续迭代直到没有异常值为止。其中因子3由关键字IQ_factor指定，用户可以指定其他值。线性趋势的估计采用白噪声模型，周年信号也包括在估计过程中，如果需要声明其他周期信号可以使用关键字periodicsignals。
**************************************removeoutliers, version 1.7.2
**************************************
Data format: MJD, Observations, Model
Filename : ./obs_files/TEST.mom
Number of observations: 1000
Percentage of gaps : 10
after removal gaps, m=900
No Polynomial degree set, using offset + linear trend
No extra periodic signal is included.
Found 6 bad points.
after removal gaps, m=894
No Polynomial degree set, using offset + linear trend
No extra periodic signal is included.
Found 0 bad points.
--> ./pre_files/TEST.mom
将原始数据和预处理过的数据通过如下命令绘制为图3
set xlabel 'years';
set ylabel 'mm';
plot './ex1/obs_files/TEST.mom' using (($1-51544)/365.25+2000):2 with lines t 'observed', './ex1/pre_files/TEST.mom' using (($1-51544)/365.25+2000):2 with lines t 'preprocessed'
图3. 原始数据和预处理后的数据
4.1.2 估计线性趋势
剔除异常值之后即可估计线性趋势，输入参数在控制文件estimatetrend.ctl中，同上可以用别名，内容如下：
DataFile TEST.mom
DataDirectory ./pre_files
OutputFile ./mom_files/TEST.mom
interpolate no
seasonalsignal yes
halfseasonalsignal no
estimateoffsets yes
NoiseModels Powerlaw White
LikelihoodMethod AmmarGrag
PhysicalUnit mm
ScaleFactor 1
其中噪声模型还可以选择PowerlawApprox, Flicker, RandomWalk, ARMA, ARFIMA或者GGM 等等。注意Power-law只适用于谱指数大于1即平稳噪声的情况，对于非平稳幂律噪声，通常采用下面的模型，可以给出相近的结果：
NoiseModels GGM White
GGM_1mphi 6.9e-06
最大似然计算采用的是AmmarGrag，关键字LikelihoodMethod属于可选项，默认就是AmmarGrag，前提是缺失数据不超过50%，否则采用Full Covariance matrix (’FullCov’)方法。注意如果removeoutliers.ctl中的ScaleFactor不是1，那么最好在estimatetrend.ctl中设置为1，防止做两次尺度变换。偏移信息也可以通过使用关键字OffsetFile和component放到另一个文件中。estimatetrend的输出为：
************************************estimatetrend, version 1.7.2
************************************
0) Powerlaw
1) White
generator type: mt19937
seed = 0
first value = 2184496367
Data format: MJD, Observations, Model
Filename : ./pre_files/TEST.mom
Number of observations: 1000
Percentage of gaps : 10.6----------------AmmarGrag
----------------
No Polynomial degree set, using offset + linear trend
No extra periodic signal is included.Number of CPU's used (threads) = 41 0.40000 0.40000 f()= 1729.103280 size=0.2002 0.40000 0.40000 f()= 1729.103280 size=0.146
……50 0.37328 0.38351 f()= 1726.399399 size=0.00051 0.37328 0.38351 f()= 1726.399399 size=0.000
converged to minimum at52 0.37328 0.38351 f()= 1726.399399 size=0.000Likelihood value
--------------------
min log(L)=-1726.399
k =8 + 2 + 1 = 11
AIC =3474.799
BIC =3527.552
BIC_tp =3507.335
BIC_c =3558.758
ln_det_I =23.206Powerlaw:
fraction = 0.48921
sigma = 3.50785 mm/yr^0.19175
d = 0.3835 +/- 0.0867
kappa = -0.7670 +/- 0.1733White:
fraction = 0.51079
sigma = 1.15618 mm
No noise parameters to showSTD of the driving noise: 1.61773
bias : 0.838 +/- 1.001 mm (at 1997/5/15, 12:0:0.000)
trend: 16.400 +/- 0.689 mm/year
cos yearly : 4.092 +/- 0.281 mm
sin yearly : -3.937 +/- 0.286 mm
Amp yearly : 5.685 +/- 0.283 mm
Pha yearly : -43.896 degrees
offset at 50284.0000 : 24.49 +/- 0.69 mm
offset at 50334.0000 : -24.69 +/- 0.75 mm
offset at 50784.0000 : -39.76 +/- 0.70 mm
offset at 51034.0000 : 38.52 +/- 0.72 mm
--> ./mom_files/TEST.mom
Total computing time: 1.00000 sec
前几行的意义都很明确，显示CPU使用个数是为了确保多线程或者多核心处理器在正常运行，后面几行显示最小化似然值的自然对数的负值的过程（也即极大似然）。选定噪声模型的质量通过赤池信息准则（Akaike Information Criteria，AIC）和贝叶斯信息准则（Bayesian Information Criteria，BIC）评价，这些值在似然值的自然对数下面。
幂律噪声和白噪声的振幅估计过程为：首先将幂律噪声建模为分数布朗运动，这可以通过对白噪声滤波实现；第二步加入方差得到总方差
C=σwn2I+σpl2E(κ)(1)\mathbf{C}=\sigma_{w n}^{2} \mathbf{I}+\sigma_{p l}^{2} \mathbf{E}(\kappa) \tag{1} C=σwn2I+σpl2E(κ)(1)
其中σwn2\sigma_{w n}^{2}σwn2和σpl2\sigma_{p l}^{2}σpl2分别是待估的幂律噪声和白噪声的振幅，写成下式更为方便
C=σ2[ϕI+(1−ϕ)E(κ)](2)\mathbf{C}=\sigma^{2}[\phi \mathbf{I}+(1-\phi) \mathbf{E}(\kappa)] \tag{2} C=σ2[ϕI+(1−ϕ)E(κ)](2)
σ2\sigma^{2}σ2可以从残差中计算出来，所以只需要寻找ϕ\phiϕ使得似然值最大，需要搜寻参数空间的维数减少也有利于加快极大似然估计的速度。可以看到ϕ\phiϕ的值，以及σ2\sigma^{2}σ2的值为1.61773，则白噪声的标准差为
σwn=0.51079×1.61773=1.15618(3)\sigma_{w n}=\sqrt{0.51079} \times 1.61773=1.15618 \tag{3} σwn=0.51079×1.61773=1.15618(3)
Hector中，协方差矩阵E\mathbf{E}E没有用因子△T−κ/2{\bigtriangleup T}^{-\kappa/2}△T−κ/2尺度化，其中△T{\bigtriangleup T}△T是单位为年的采样周期。为了便于和文献里的幂律噪声振幅比较，将估计得到的振幅除以△T−κ/4{\bigtriangleup T}^{-\kappa/4}△T−κ/4
σpl=0.48921×1.61773(1/365.25)0.25×0.7670=3.50785(4)\sigma_{p l}=\frac{\sqrt{0.48921} \times 1.61773}{{(1/365.25)}^{0.25×0.7670}} =3.50785 \tag{4} σpl=(1/365.25)0.25×0.76700.48921×1.61773=3.50785(4)
幂律噪声的第二个参数是谱指数ddd也即其他文章里常用于GNSS时间序列的−κ/2{-\kappa/2}−κ/2，而白噪声没有额外的参数需要估计。剩余几行显示了模型的估计值，比如nominal bias（也即t0t_0t0时刻的截距）、线性趋势和季节信号。默认t0t_0t0为时间序列的时间中点。例如假定有5个观测值，设计矩阵H\mathbf{H}H为：
(1−21−1101112)\left( \begin{array}{cc} 1 & -2 \\ 1 & -1 \\1 & 0 \\1 & 1 \\1 & 2 \end{array} \right) ⎝⎜⎜⎜⎜⎛11111−2−1012⎠⎟⎟⎟⎟⎞
第一列估计nominal bias，第二列估计线性趋势，这两列是正交的因为HTH\mathbf{H}^T \mathbf{H}HTH产生一个对角阵，从而nominal bias的估计和线性趋势的估计互不影响。nominal bias对应第三行的时刻，如果需要声明另一个参考时刻，可以用关键字ReferenceEpoch+年月日，如：
ReferenceEpoch 2008 1 1
输出结果还包括偏移的估计值。estimatetrend的结果保存在./mom_files/TEST.mom文件中，输入gnuplot代码得到结果如图4所示。
set xlabel 'years';
set ylabel 'mm';
plot './obs_files/TEST.mom' using (($1-51544)/365.25+2000):2 with lines t 'observed' lw 1 lc rgb 'black', './pre_files/TEST.mom' using (($1-51544)/365.25+2000):2 with lines t 'preprocessed' lw 2 lc rgb 'red', './mom_files/TEST.mom' using (($1-51544)/365.25+2000):3 with lines t 'model' lw 2 lc rgb 'blue';
图4. 原始数据、预处理后数据和模型值
4.1.3 绘制功率谱密度图
为了测试模型是否正确，对残差绘制功率谱密度图，这可以通过程序estimatespectrum实现，该程序计算Welch周期图，结果存放在estimatespectrum.out，控制文件为estimatespectrum.ctl，其内容如下：
DataFile TEST.mom
DataDirectory ./mom_files
interpolate no
ScaleFactor 1.0
WindowFunction Parzen
Fraction 0.1
estimatespectrum默认将时间序列分为4份，由于50%的重叠，实际上一共有7块，本例中每一块的长度为250个数据点。用Parzen窗函数将每一块的前后各10%（关键字Fraction）平滑为0，窗函数的关键字是WindowFunction，还可以选择Hann。如果需要更多数据块，从而虽然损失频率范围但是获得更好的平滑效果，可以在命令行声明更大的分块个数如：
estimatespectrum 8
该命令将会把时间序列分为8份，由于重叠50%而产生15块。单边PSD下的面积应当等于时间序列的方差，这里面积的计算就是简单假定每个点代表一个宽度为1/△T1/{\bigtriangleup T}1/△T的矩形然后累加。程序运行如下：
************************************estimatespectrum, version 1.7.2
************************************
Data format: MJD, Observations, Model
Filename : ./mom_files/TEST.mom
Number of observations: 1000
Percentage of gaps : 10.6
window function : Parzen
window function fraction: 0.1
Number of data points n : 1000
Number of data used N : 1000
Number of segments K : 7
Length of segments L : 250
U : 0.9722
dt: 8.64e+04
scale for G to get Amplitude (mm): 0.3043
Total variance in signal (time domain): 3.295
Total variance in signal (spectrum) : 2.956
freq0: 4.6296e-08
freq1: 5.7870e-06
--> estimatespectrum.out
估计的噪声模型的PSD可以用程序modelspectrum计算，该程序从modelspectrum.ctl获取噪声模型的信息，并且用户必须手动输入噪声参数和频率范围，其输出保存在modelspectrum.out，运行如下：
************************************modelspectrum, version 1.7.2
************************************
---------------
ModelSpectrum
---------------Enter the standard deviation of the driving noise: 1.6190
Enter the sampling period in hours: 24
0) Powerlaw
1) White
generator type: mt19937
seed = 0
first value = 1667267897
Enter fraction for model Powerlaw: 0.490
Enter fraction for model White: 0.5103Powerlaw:
Enter value of fractional difference d:0.3831White:
1) Linear or 2) logarithmic scaling of frequency?: 2
Enter freq0 and freq1: 4.6296e-08 5.7870e-06
freq0 : 4.6296e-08
freq1 : 5.787e-06
--> modelspectrum.out
残差和噪声模型的PSD如图5所示：
set xlabel 'Frequency (cpy)' font '18';
set ylabel 'Power (mm^2/cpy)' offset -1,0 font '18'
set logscale xy;
set format y '10^{%T}';
set format x '10^{%T}';
s=31557600.0;
plot './ex1/estimatespectrum.out' using ($1*s):($2/s) linecolor rgb "red" t 'observed', './ex1/modelspectrum.out' using ($1*s):($2/s) with lines lc rgb "green" lw 2 t 'model'
图5. 残差和噪声模型的PSD
4.2 示例2：Cascais的月均PSMSL验潮站数据
ex2文件夹下存储了从PSMSL下载的Cascais站点数据，该时间序列没有异常值，所以可以直接估计线性趋势，控制文件estimatetrend.ctl如下：
DataFile 52.rlrdata
DataDirectory ./
OutputFile 52_out.mom
QuadraticTerm no
interpolate no
firstdifference no
seasonalsignal yes
halfseasonalsignal yes
estimateoffsets no
NoiseModels ARMA
PhysicalUnit mm
AR_p 1
MA_q 0
这里我们使用ARMA噪声模型，更精确地，设置了1个AR系数（关键字AR_p）和0个MA系数（关键字MA_q），从而可以将ARMA(1,0)简化为AR(1)。此时运行estimatetrend得到线性趋势1.270±0.075 mm/yr。如果将噪声模型分别改为AR(5)、ARFIMA（其中AR_p=1，MA_q=0）和GGM可以得到相应的趋势1.277±0.103、1.253±0.175和1.259±0.201 mm/yr，这些结果相比AR(1)都具有更小的BIC和AIC值。
利用modelspectrum和estimatespectrum可以得到PSD，给出的频率范围是1.1317e-09到1.9013e-07Hz，结果如图6所示，注意采样时间为730.5小时。另外，由于每次只用了一个噪声模型，因而分数始终为1。控制文件estimatespectrum.ctl如下：
DataFile 52_out.mom
DataDirectory ./
OutputFile estimatespectrum.out
interpolate no
set terminal postscript enhanced size 4,4 color solid "Helvetica"
set output './Cascais_psd.eps'
set xlabel 'Frequency (cpy)' font 'Helvetica, 18';
set ylabel 'Power (mm^2/cpy)' offset -1,0 font 'Helvetica, 18';
set logscale xy;
set key bottom left;
set format y '10^{%T}';
set format x '10^{%T}';
set xrange[*:6];
s=31557600.0
plot './ex2/estimatespectrum.out' using ($1*s):($2/s) with p pt 1 t 'observed',\'./ex2/modelspectrum_AR1.out' using ($1*s):($2/s) with l lw 3 t 'AR(1)',\'./ex2/modelspectrum_AR5.out' using ($1*s):($2/s) with l lw 3 t 'AR(5)',\'./ex2/modelspectrum_ARFI.out' using ($1*s):($2/s) with l lw 3 t 'ARFI(1,d)',\'./ex2/modelspectrum_GGM.out' using ($1*s):($2/s) with l lw 3 t 'GGM'
图6. Cascais验潮站数据的PSD
海平面变化研究有时候也关注加速度。这可以通过在控制文件estimatetrend.ctl中设置关键字DegreePolynomial为2实现（默认值为1），如果再将噪声模型设置为GGM，将得到二次项系数0.004±0.006 mm/yr2^22，该值为加速度的一半。
接下来的分析可以加入额外的地球物理信号。在控制文件中设置关键字estimatemultivatiate为yes并声明地球物理信号的文件，如下：
estimatemultivariate yes
MultiVariateFile hadslp2_cascais.mom
这里用到了月均哈德利中心海表气压数据集（Hadley Centre Sea Level Pressure dataset），再次设置噪声模型为GGM，此时运行estimatetrend得到
STD of the driving noise: 38.17881
bias : 19329.610 +/- 454.450 mm (at 1937/12/31, 12:45:0.000) trend: 1.292 +/- 0.186 mm/year
cos yearly : 20.003 +/- 1.991 mm
sin yearly : -27.878 +/- 1.941 mm
Amp yearly : 34.369 +/- 1.964 mm
Pha yearly : -54.340 degrees
cos hyearly : 0.650 +/- 1.625 mm
sin hyearly : -10.391 +/- 1.553 mm
Amp hyearly : 10.533 +/- 1.579 mm
Pha hyearly : -86.418 degrees
scale factor of channel 1 : -12.17 +/- 0.45
最后一行显示了multi-variate file第一列地球物理信号的回归系数（如尺度因子），该多变量文件遵循MOM格式，可以添加任意多列，每一列都包含不同的地球物理信号。在本例中回归系数是-12.17 mm/mbar，与标准逆气压改正值相近。
4.3 示例3：模拟有色噪声
为了进行Monte Carlo模拟，必须产生带有色噪声的时间序列，这可以通过程序simulatenoise实现，其思路是对每一个不同的噪声模型用不同的脉冲响应与白噪声时间序列做卷积，从而产生模拟噪声时间序列。上述过程可能产生spin-up effects，这是因为默认在第一次观测之前噪声为0。为了缓解该效应，可以给关键字TimeNoiseStart赋一个比较大的值来声明第一次观测前额外需要建模的点数，通常1000就够了。ex3目录下的控制文件simulatenoise.ctl为：
SimulationDir ./
SimulationLabel test_base
NumberOfSimulations 10
NumberOfPoints 5000
SamplingPeriod 1
TimeNoiseStart 1000
NoiseModels Flicker White
PhysicalUnit mm
关键字SimulationDir声明产生文件存放的位置，SimulationLabel声明文件的基本名称，NumberOfSimulations声明模拟次数，从而产生的文件名为基本名称+模拟次数，本例中为test_base0.mom、test_base1.mom、. . .、test_base9.mom。NumberOfPoints声明时间序列的点数。运行程序过程中，命令行会提示输入噪声参数，在本例中即为ϕ\phiϕ和ddd。
4.4 示例4：震后回弹
一场大型地震过后，地球表面会缓慢恢复到新的位置，这个震后回弹过程可以用指数函数或者对数函数建模。目录ex4下存放了一个带有多个回弹信号的模拟时间序列。已知地震发生的时刻，所以可以在./obs_files目录下的TEST.mom文件头添加偏移时刻，进一步地还可以给每个回弹事件添加信息。
# sampling period 1.0
# offset 51994.0
# offset 53544.0
# offset 55044.0
# log 51994.0 10.0
# log 55044.0 10.0
# exp 53544.0 100.0
为了让estimatetrend能够处理该文件，需要在控制文件中添加一行‘estimatepostseismic yes’
DataFile TEST.mom
DataDirectory ./obs_files
OutputFile ./mom_files/TEST.mom
interpolate no
seasonalsignal no
halfseasonalsignal no
estimateoffsets yes
estimatepostseismic yes
NoiseModels GGM White
GGM_1mphi 6.9e-06
PhysicalUnit mm
ScaleFactor 1.0
TimeNoiseStart 1000.0
观测值和模型值如图7所示。慢滑动事件的处理类似：在控制文件中添加一行’estimateslowslipevent yes’，在mom文件头添加’# tanh MMMM DD’，其中MMMM是简化儒略日，DD是回弹时间。
set xlabel 'years' font '18';
set ylabel 'vertical (mm)' offset -1,0 font '18';
plot './ex4/mom_files/TEST.mom' u (($1-51544)/365.25+2000):2 w l t 'observed',\'./ex4/mom_files/TEST.mom' u (($1-51544)/365.25+2000):3 w l lw 3 t 'fitted'
图7. 带有指数和对数形式震后回弹的模拟时间序列与拟合信号
4.5 示例5：自动分析+变化的季节信号
本例中使用站点MAS1、LHAZ、AUCK和ULAB的真实GNSS时间序列，这些站点都表现出振幅变化的季节信号。首先要做的就是探测偏移所在位置，这可以通过在ex5目录下运行python脚本文件find_all_offsets.py实现，该过程可能花费几小时，具体时间取决于电脑性能。另一个python脚本analyse_and_plot.py可用于查找异常值、估计线性趋势+季节信号+偏移、绘制时间序列及其功率谱密度，在需要分析大量数据的时候比较方便。
根据Bennett (2008)的工作，可以将周年和半年信号的振幅建模为Chebyshev多项式，控制文件estimatetrend_Bennett.ctl如下，结果如图8所示。
DataFile MAS1_2.mom
DataDirectory ./pre_files
OutputFile ./mom_files/MAS1_2.mom
interpolate no
PhysicalUnit mm
ScaleFactor 1.0
NoiseModels GGM White
GGM_1mphi 6.9e-06
seasonalsignal yes
halfseasonalsignal yes
estimatevaryingseasonal yes
varyingseasonal_N 3
estimateoffsets yes
ScaleFactor 1.0
PhysicalUnit mm
图8. MAS1的观测数据和拟合模型，红色是振幅固定的季节信号，黑色是变化的季节信号
对振幅变化的季节信号建模的另一个方式是用维纳滤波，前提是知道时间序列的噪声性质和变化的周年半周年信号的随机部分的值，前面一部分可以通过运行python脚本analyse_timeseries.py得到。用户也可以简单地使用apply_WF.py，输入
apply_WF.py MAS1_2 0.10 0.06 0.9999
其中0.1和0.06分别是周年和半周年信号振幅的AR(1)随机变化的标准差，0.9999是一阶AR系数（非常接近于1即随机游走）。该脚本将在sea_files和fil_files目录下创建文件，分别包括估计的变化季节信号和减去该变化信号的观测值，结果如图9所示，注意只要数据缺失少于10%左右就可以将数据缺失置为0。
图9. MAS1垂直分量的观测数据和拟合模型，红色是振幅固定的季节信号，黑色是维纳滤波得到的变化的季节信号
时间序列分析软件Hector用户手册（一）相关推荐
- 时间序列分析软件Hector用户手册（二）
时间序列分析软件Hector用户手册(二) 5 模型 5.1 震后形变 5.2 多趋势 5.3 地球物理信号 6 可接受格式 6.1 mom格式 6.2 enu格式 6.3 neu格式 6.4 rlr ...
- InSAR处理软件+时间序列分析软件
这里写目录标题 InSAR处理软件 gamma ISCE doris ROI_PAL GMTSAR SNAP PIE InSAR时间序列处理软件 MintPy PyRate StaMPS InSAR处 ...
- python 自相关函数_用Python的statsmodels库计算时间序列的自相关函数和画图
在时间序列分析课程中会需要用到自相关函数的计算,也就是当前期的值和滞后期的值之间的关系,这个指标的计算在计量软件中会比较容易实现,但是如果想要用python做怎么实现呢.代码如下:#导入库 impor ...
- 1.2 InSAR数据处理之软件介绍
本节内容主要分为两部分,分别是: 1.Windows操作系统下的软件介绍. 2.Linux系统下的软件介绍. 许多InSAR数据处理软件只能安装于Linux操作系统上,但有一部分软件既能够支持Linu ...
- [QOCA学习笔记]利用QOCA软件PCA模块进行共模误差分析
操作环境:win10+虚拟机ubuntu下QOCA(ver 1.33) 数据:通过解算得到的安徽省CORS参考站坐标时间序列(neu格式) QOCA软件是优秀的GNSS坐标时间序列分析软件,其简介及P ...
- 统计软件 matlab,统计软件列表 - MATLAB等数学软件专版 - 经管之家(原人大经济论坛)...
统计软件列表ActivStats 多媒体交互式学习软件包,统计学入门好帮手 Windows版本 ADE-4 一个多元数据分析软件. Windows版本 ALSCAL 多维等级分析(Multidimen ...
- [Hector学习笔记]GNSS时间序列处理软件Hector使用备忘（批处理脚本）
由Machiel Bos 和Rui Fernandes共同研发的Hector软件可以利用时间相关噪声来估计时间序列线性趋势,是与优秀的GNSS时间序列处理软件之一. 软件安装教程:Hector软件安装 ...
- abaqus分析用户手册单元卷_作用卷、分析卷、材料卷三件套，让你也能熟练应用Abaqus...
<Abaqus分析用户手册--指定条件.约束与相互作用卷>是"Abaqus用户手册大系"中的一册,包括指定条件.约束与相互作用三个部分.指定条件部分对各种物理过程中涉及 ...
- python 离散数据时间序列图_每个人都学的会的数据分析
数据分析已经成为数据时代各行各业突破各自行业发展瓶颈的最有效手段,无论是公司职员还是个体商户或大公司管理者,都需要有数据分析的能力.很多人认为数据分析能力就是对数据进行描述和做出漂亮的统计图形的能力, ...
- python 线性回归 统计检验 p值_PAST：最简便易用的统计学分析软件教程（一）软件基本信息介绍...
欢迎来到PAST的世界! PAST的全称为PAleontological STatistics,是1995年由P.D. Ryan.D.A.T. Harper和J.S. Whalley开发的软件PALS ...
最新文章
- 逆天了：Nature一篇论文57000位作者，更厉害的是，大多数作者都是游戏玩家
- Python基础教程（八）：日期和时间、文件I/O、异常处理
- 关于C++的extern关键字
- 一行SQL代码能做什么？
- 原来信用卡肉这么肥，怪不得银行天天给你发短信叫你办理
- b样条和三次样条_样条曲线
- python监控文件内容变化_Python监控文件内容变化
- Eclipse最全快捷键
- 控制kvm-qcow2增长空间-（一）
- 并发编程（十一）—— Java 线程池 实现原理与源码深度解析（一）
- 书籍-分布式系统常用技术及案例分析
- [量化学院]基于协整的配对交易
- 8、鼠标控制与32位模式切换
- 深度学习机器学习面试问题准备（必会）
- 【转载】专访罗升阳：老罗的Android之旅
- 时间复杂度的三种常见表示符号
- 如何实现移动端点击下拉箭头显示全部文字
- Web自动化测试面试
- （小白）Excel学习笔记
- Domino Web网页中更改密码比你想得简单得多