时间序列分析软件Hector用户手册（二）

5 模型

假设观测值是一个确定性模型和随机噪声的累加，其中确定性模型为：

$$\begin{array}{rl}\mathbf{x}=& \sum _{i=0}^{n_P}{\mathbf{p}}_i{\left(t-t_R\right)}^i+\sum _{i=1}^{n_J}{\mathbf{b}}_{j}H\left(t-t_j\right)\\ & +\sum _{i=1}^{n_F}{\mathbf{s}}_i\sin \left({\omega }_{i}t\right)+{\mathbf{c}}_i\cos \left({\omega }_{i}t\right)\\ & +\sum _{i=1}^{n_T}{\mathbf{e}}_i\left(1-\exp \left(-\left(t-t_i\right)/T_i\right)\right)\\ & +\sum _{j=1}^{n_L}{\mathbf{a}}_j\log \left(1+\left(t-t_j\right)/T_j\right)\\ & +\sum _{k=1}^{n_T}\frac{{\mathbf{u}}_k}{2}\left[\tanh \left(\left(t-t_k\right)/T_k\right)-1\right]\\ & +\sum _{l=1}^{n_l}{a}_l\times f_l(t)\end{array}\text{\hspace{1em}(6)}$$

其中${\mathbf{p}}_i$是$n_P$阶多项式的系数，默认$n_P=1$也就是线性趋势。按照惯例，速率是单位每年，其中一年定义为365.25天，加速度定义为二次项的两倍。$H\left(t\right)$是Heaviside阶梯函数，用来对振幅为${\mathbf{b}}_j$的偏移建模。周年和半年信号常常纳入模型中，所以有自己的关键字seasonal和halfseasonal，在公式中，它们的角速度用${\omega }_i$表示。其他周期信号需要用关键字PeriodicSignals+以天为单位的周期列表输入。为了计算周年半年的振幅，假定${\mathbf{s}}_i$和${\mathbf{c}}_i$的平均不确定度为$\sigma$，这样振幅就服从均值为$\sigma \sqrt{\pi /2}{L}_{1/2}\left(-{\nu }^2/\left(2{\sigma }^2\right)\right)$的莱斯分布，其中$\nu =\sqrt{c_i^2+s_i^2}$，$L_{1/2}(x)={}_1{F}_1(-1/2,1,x)$。

5.1 震后形变

从版本1.6之后Hector可以对震后形变建模，这需要在数据文件头添加形变开始和持续的时间，单位分别是MJD和天。一个例子是：

# log 51994.0 10.0 
# log 55044.0 10.0 
# exp 53544.0 100.0

可选指数函数或者对数函数，控制文件estimatetrend.ctl必须声明关键字estimatepostseismic为yes。注意用户无法选择东向、北向或者垂直分量，该选项对所有分量适用。

1.7版本之后支持用双曲正切函数对缓慢滑移事件建模，两个参数$t_k$和$T_k$分别代表缓慢滑移事件的中点和宽度，都必须在文件头输入为如下形式，并且在控制文件中设置关键字estimateslowslipevent为yes。

# tanh 51994.0 10.0

5.2 多趋势

一些GNSS时间序列在一场地震之后表现出了不同的线性趋势，为了估计单个时间序列中的不同线性趋势，可以将式6中的多项式趋势替换为：

$$p_0+\sum _{k=1}^M{p}_{1k}\left(t_i-b_k\right)\left(H\left(t_i-b_k\right)-H\left(t_i-b_{k+1}\right)\right)\text{\hspace{1em}(7)}$$

其中$p_{1k}$是第$k$段的线性趋势，假定有$M$段，第$k$段的初始时刻为$b_k$，$H\left(t\right)$仍然是Heaviside阶梯函数。

线性趋势的间断点在时间序列的文件头给出，假定前例中的两场地震伴随着北分量和东分量的线性趋势的变化，那么可以写成：

# break 51994.0 0 
# break 51994.0 1

当使用的文件只有一个分量时，对于分量的声明是非必要的。为了估计多趋势，控制文件还需要设置关键字estimatemultitrend为yes。

5.3 地球物理信号

公式6中的最后一项代表尺度因子$a_l$（或者称为回归系数）乘以一个函数$f_l$，这可以表示一个地球物理信号，示例2中显示了如何将海表面气压融入到分析中和如何获取接近于逆气压改正的尺度值。通过声明：

estimatemultivariate yes 
MultiVariateFile XXXX.mom

文件XXXX.mom为MOM格式，第一列是MJD，之后可以添加任意多列地球物理信号。

6 可接受格式

所有输入文件都应当是ASCII格式，时间单调且均匀增加，数据格式通过后缀声明，Hector也接受自由格式，这意味着数值的位数和空格数可变。

对于mom和enu格式，以天为单位的采样周期可以在文件头声明如下：

# sampling period 1.0

如果该信息缺失，hector将从头几个观测值中估计采样周期，其能自动探测的周期有：0.5小时、1小时、1天、7天。

6.1 mom格式

第一列是MJD格式的观测时间，第二列是观测值，第三列是估计的模型值，该列可选。允许缺失。

6.2 enu格式

类似于mom格式，但是有4列，第一列是MJD格式的观测时间，2-4列分别是东分量、北分量和垂直分量。允许缺失。

6.3 neu格式

该格式为SOPAC所用，并且为CATS软件接受。第一列是年为单位的观测时间，2-4列是米为单位的北分量、东分量和垂直分量。允许缺失。单位m很容易通过设置关键字ScaleFactor转换成mm。以年为单位的时间通过下式转换为MJD：

$$MJD=floor(365.25\ast (T-1970)+40587+0.1)-0.5\text{\hspace{1em}(8)}$$

6.4 rlrdata格式

Hector可以读取PSMSL的月均数据，为了产生间隔完全相等的数据，假定每个年份为精确的30.4375天即730.5小时。由于一些年份有365天而其他为366天，因而Hector内部使用的转换公式为：

$$MJD=30.4375(12(year-1859)+(month-1))+59\text{\hspace{1em}(9)}$$

7 噪声模型

Hector支持不同的噪声模型及其组合，GNSS时间序列中应用最广的是幂律噪声加白噪声，Williams (2008)将该特定组合的协方差矩阵$\mathbf{C}$写作：

$$\mathbf{C}={\sigma }^2\left({\cos }^2\varphi \mathbf{I}+{\sin }^2\varphi \mathbf{E}(d)\right)\text{\hspace{1em}(10)}$$

其中$\mathbf{I}$是单位阵（单位白噪声的协方差阵），$\mathbf{E}$是幂律噪声的协方差矩阵，取决于谱指数$d$，两个噪声振幅的分布由$\varphi$控制，取值在0和$\pi /2$之间，总的方差为${\sigma }^2$。该式可以推广到n维球坐标的形式：

$$\begin{array}{rl}{x}_1& =\cos \left({\varphi }_1\right)\\ x_2& =\sin \left({\varphi }_1\right)\cos \left({\varphi }_2\right)\\ x_3& =\sin \left({\varphi }_1\right)\sin \left({\varphi }_2\right)\cos \left({\varphi }_3\right)\\ ⋮& \\ x_{n-1}& =\sin \left({\varphi }_1\right)\dots \sin \left({\varphi }_{n-2}\right)\cos \left({\varphi }_{n-1}\right)\\ x_n& =\sin \left({\varphi }_1\right)\dots \sin \left({\varphi }_{n-2}\right)\sin \left({\varphi }_{n-1}\right)\end{array}\text{\hspace{1em}(11)}$$

对于$n=2$有：

$$\begin{array}{rl}{x}_1& =\cos \left({\varphi }_1\right)\\ x_2& =\sin \left({\varphi }_1\right)\end{array}\text{\hspace{1em}(12)}$$

对于$n=3$有：

$$\begin{array}{rl}{x}_1& =\cos \left({\varphi }_1\right)\\ x_2& =\sin \left({\varphi }_1\right)\cos \left({\varphi }_2\right)\\ x_3& =\sin \left({\varphi }_1\right)\sin \left({\varphi }_2\right)\end{array}\text{\hspace{1em}(13)}$$

在Hector中这些坐标的平方称为分数$f_i=x_i^2$，从而将式10推广为：

$$\mathbf{C}={\sigma }^2\left(f_1{\mathbf{E}}_1+f_2{\mathbf{E}}_2+f_3{\mathbf{E}}_3+\dots +f_n{\mathbf{E}}_n\right)\text{\hspace{1em}(14)}$$

对于$n$个噪声模型需要$n-1$个角度${\varphi }_1,\dots ,{\varphi }_{n-1}$，所有分数在0-1之间变化。如第一节提到的，只接受平稳噪声。该式将产生Toeplitz协方差阵，从而只需要存储第一列，用$\gamma$指代。

7.1 白噪声

对于白噪声，协方差阵就是单位阵，当$\sigma =1$，协方差阵$\mathbf{C}$的第一列就是：

$$\begin{array}{rlrlr}{\gamma }_i& =1& & \text{ for }& i=0\\ & =0& & & i\ne 0\end{array}\text{\hspace{1em}(15)}$$

单边功率谱密度为：

$$S(f)=2\frac{1}{f_s}\text{\hspace{1em}(16)}$$

其中$f_s$是采样频率（Hz）。如果从0频率积分到Nyquist频率将得到时间序列的方差。

7.2 幂律噪声

幂律噪声协方差阵的第一列为：

$${\gamma }_i=\frac{\Gamma (d+i)\Gamma (1-2d)}{\Gamma (d)\Gamma (1+i-d)\Gamma (1-d)}\text{\hspace{1em}(17)}$$

$\sigma =1$时单边功率谱密度为：

$$S(f)=2\frac{1}{f_s}\frac{1}{{\left(2\sin \left(\pi f/f_s\right)\right)}^{2d}}\text{\hspace{1em}(18)}$$

7.3 幂律噪声近似

Bos et al. (2013)介绍了幂律噪声的一种Toeplitz近似，可以在关键字Noisemodels后用标签PowerlawApprox选择该模型。此外，用户必须声明在关键字TimeNoiseStart后声明观测起始时间之前的天数，1000是一个较好的初值。该近似只对较弱的非平稳噪声有效，因而允许设置的最小谱指数$\kappa$为-1.6。现在我们倾向于使用$\varphi$接近于1的GGM噪声模型，该模型对幂律噪声的近似效果更好且对低至-2的谱指数有效。

7.4 ARFIMA和ARMA

ARFIMA噪声模型的定义为：

$$\Phi (L)(1-L{)}^d{z}_t=\Theta (L){ϵ}_t\text{\hspace{1em}(19)}$$

其中$L$是延迟算子，$z_t$是残差，$ϵ$是白噪声。上式意味着观测值的残差可以通过对白噪声做变换得到。算子$\Phi$和$\Theta$定义为：

$$\begin{array}{l}\Phi (L)=1-{\varphi }_{1}L-{\varphi }_2{L}^2-\dots -{\varphi }_p{L}^p\\ \Theta (L)=1+{\theta }_{1}L+{\theta }_2{L}^2+\dots +{\theta }_q{L}^q\end{array}\text{\hspace{1em}(20)}$$

整数$p$和$q$的值在控制文件estimatetrend.ctl中分别由关键字AR_p和MA_q设置，建议$p$小于5，以此保证MLE过程始终以产生平稳噪声的$\Phi (L)$的系数${\varphi }_1,\dots ,{\varphi }_p$开始。如果$p$和$q$均为0，将得到纯粹的幂律噪声过程。对于海平面的研究，AR(1)用的较多。对于ARMA噪声模型存在快速的最大似然方法。

7.5 GGM模型

一阶GGM模型参数为$\varphi$，修饰以额外的参数$d$可产生幂律噪声，该噪声的PSD斜率为$2d$，在极低频和极高频变平为白噪声，自协方差向量的分析表达式为：

$${\gamma }_i=\frac{\Gamma (d+i)(1-\varphi {)}^i}{\Gamma (d)\Gamma (1+i)}{}_2{F}_1\left(d,d+i;1+i;(1-\varphi {)}^2\right)\text{\hspace{1em}(21)}$$

该噪声模型的使用就是在关键字Noisemodels后添加GGM，同时声明$1-\varphi$的参数：

GGM_1mphi XXXX

如果该值足够小，比如6.9e-6，那么GGM近似于一个幂律噪声模型。

7.6 闪烁噪声和随机游走噪声

闪烁噪声和随机游走噪声是两种简单的谱指数固定分别为0.5和1的幂律噪声，但是由于Hector只能处理平稳噪声（平稳噪声的协方差阵为Toeplitz形式从而允许快速求逆技术的应用），所以Hector实际上是通过给GGM赋予较小的$1-\varphi$参数值来实现，如：

NoiseModels FlickerGGM
GGM_1mphi  6.9e-06

8 AIC和BIC

为了实现最佳模型的选取，用户可以使用AIC或者BIC。两种准则都使用似然对数作为出发点，但是对参数的增加施加惩罚从而避免过拟合。

自然对数的定义为：

$$\ln (L)=-\frac{1}{2}\left[N\ln (2\pi )+\ln \det (\mathbf{C})+{\mathbf{r}}^T{\mathbf{C}}^{-1}\mathbf{r}\right]\text{\hspace{1em}(22)}$$

其中$N$是观测值的实际个数（缺失不计），协方差阵$\mathbf{C}$分解为：

$$\mathbf{C}={\sigma }^2\overline{\mathbf{C}}\text{\hspace{1em}(23)}$$

其中$\overline{\mathbf{C}}$是不同噪声模型之和，$\sigma$是驱动白噪声的标准差，从残差中估计：

$$\sigma =\sqrt{\frac{{\mathbf{r}}^T{\overline{\mathbf{C}}}^{-1}\mathbf{r}}{N}}\text{\hspace{1em}(24)}$$

利用上式及$\det c\mathbf{A}=c^N\det \mathbf{A}$，自然对数可以写成：

$$\ln (L)=-\frac{1}{2}\left[N\ln (2\pi )+\ln \det (\overline{\mathbf{C}})+2N\ln (\sigma )+N\right]\text{\hspace{1em}(25)}$$

参数个数$k$是设计矩阵$\mathbf{H}$、噪声模型和驱动白噪声过程的方差所有参数之和。例如：使用幂律噪声和白噪声估计线性趋势，包括5个参数：nominal bias、线性趋势、幂律噪声和白噪声之间方差的分布、幂律噪声的谱指数、驱动白噪声过程的方差。

AIC和BIC分别定义为：

$$\begin{array}{rl}\text{AIC}& =2k+2\ln (L)\\ \text{BIC}& =k\ln (N)+2\ln (L)\end{array}\text{\hspace{1em}(26)}$$

选取的模型应当具有最小的AIC/BIC值，但是注意这只是相对测度而不是绝对准则。

9 辅助python脚本

Hector的早期版本提供了一些TCL语言的脚本来实现时间序列的格式转换，但是由于该语言并没有被科学界接受所以使用python替代。下表列出了噪声模型的缩写，这在所有脚本中都一样。在某些情况下也允许噪声模型的组合，比如GNSS时间序列分析中PLWN是一个常用的组合。

缩写	描述
WN	White Noise
FN	Flicker Noise
PL	Power-law Noise
RW	Random Walk
GGM	Generalised Gauss-Markov noise model
AR1	ARMA(1,0)

9.1 analyse_timeseries.py

analyse_timeseries.py可以估计趋势，正确用法是：

./analyse_timeseries.py station_name GGMWN|PLWN|FNWN|RWFNWN|WN|AR1 [fraction]

其中station_name是存储在./obs_files/目录下没有后缀的文件名，第二个参数是噪声模型的组合，最后一个可选参数告诉脚本返回噪声分数而不是振幅，这在使用simulatenoise产生具有类似于真实数据的噪声性质的模拟时间序列的时候比较方便。

按照示例4.1，在./ex1/下运行命令将得到结果如下：

./analyse_timeseries.py TEST PLWN
16.417 0.705 1000 -1726.234 3478.469 3540.813 3578.018 4.103 -3.937 0.289 0.299 -0.002 -0.052 0.208 0.211 1.61287 1.19451 3.52647 0.3999

该脚本的输出就是长度最大可达20的数组，可以用于接下来脚本的输入，其含义列于表格3。这些输出对于用户并不友好，其主要目的是作为其他脚本的输入，

表格3 analyse_timeseries.py的输出说明

索引	含义	例1
1	趋势	16.417
2	趋势不确定度	0.705
3	观测值个数	1000
4	似然对数	-1726.234
5	AIC	3478.469
6	BIC	3540.813
7	BIC_c	3578.018
8	余弦年信号振幅	4.103
9	正弦年信号振幅	-3.937
10	余弦年信号振幅的不确定度	0.289
11	正弦年信号振幅的不确定度	0.299
12	半年信号余弦振幅	-0.002
13	半年信号正弦振幅	-0.052
14	半年信号余弦振幅的不确定度	0.208
15	半年信号正弦振幅的不确定度	0.211
16	驱动噪声的振幅	1.61287
	接下来的数值随不同的噪声模型变化
WN
17	白噪声的振幅/占比
AR1
17	AR1的振幅/占比
18	系数${\varphi }_1$
FLWN
17	白噪声的振幅/占比
18	闪烁噪声的的振幅/占比
RWFLWN
17	白噪声的振幅/占比
18	闪烁噪声的振幅/占比
19	随机游走噪声的振幅/占比
PLWN
17	白噪声的振幅/占比	1.19451
18	幂律噪声的振幅/占比	3.52647
19	谱指数$d(-\kappa /2)$	0.3999
GGMWN
17	白噪声的振幅/占比
18	幂律噪声的振幅/占比
19	谱指数$d(-\kappa /2)$
20	1mphi

使用fraction选项的输出为：

./analyse_timeseries.py TEST PLWN fraction
16.417 0.705 1000 -1726.234 3478.469 3540.813 3578.018 4.103 -3.937 0.289 0.299 -0.002 -0.052 0.208 0.211 1.61287 0.54850 0.45150 0.3999

数值1.61287 0.54850 0.45150可以作为simulatenoise的输入产生模拟时间序列。

9.2 analyse_and_plot.py

analyse_and_plot.py对./obs_files目录下的所有站点作如下处理：

1. 剔除异常值
2. 运行estimatetrend
3. 绘制时间序列及拟合模型值
4. 估计残差和拟合噪声模型的PSD
5. 绘制功率谱密度

运行命令如下：

./analyse_and_plot.py GGMWN|PLWN|FNWN|RWFNWN|WN|AR1 [station]

最后一个是可选项，如果给定了站点，将只对该站点分析并绘图，否则将对./obs_files下的所有文件分析；estimatetrend的输出存储到./mom_files，数据绘图存放到./data_figures，PSD绘图存放到./psd_figures。该脚本将产生3个csv文件：analysis_results.csv、seasonal.csv和noise.csv，这些文件包含每个站点analyse_timeseries.py的输出，索引1-7位于analysis_results.csv，8-15位于seasonal.csv，其他的位于noise.csv，异常值存放在outliers.csv。

9.3 find_offset.py

findoffset只寻找最有可能发生偏移的时刻，但是一般来说用户更想找到所有的偏移，因而需要多次调用findoffset，find_offset.py脚本即可实现该功能，其调用方式为：

./find_offset.py station_name PLWN|FNWN|RWFNWN|WN [3D]

和前面一样，假定时间序列以mom格式存放在./raw_files目录下，第3个可选参数表示在单个分析过程同时分析东分量、北分量和垂直分量，要让该选项生效，文件名应当符合惯例XXX_[012]，其中XXX是站点名称，0、1、2分别代表东分量、北分量和垂直分量。运行结果如下：

****** $ ./find_offset.py TEST PLWN
**************************************
    removeoutliers, version 1.7.2
**************************************
Data format: MJD, Observations, Model
Filename              : ./dummy_raw.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 3 bad points.
after removal gaps, m=897
No Polynomial degree set, using offset + linear trend
No extra periodic signal is included.
Found 0 bad points.
--> dummy0_0.mom
MJD=50784.000000, trend=5.828000, BIC_c=5147.758000
For first round (no new offsets added) BIC_c: 5147.758000
--> mjd=50784.000000,  value=5039.837182
MJD=51034.000000, trend=26.520000, BIC_c=5031.352000
For offsets 1 BIC_c: 5031.352000
--> mjd=51034.000000,  value=4915.401033
MJD=50284.000000, trend=15.444000, BIC_c=4910.484000
For offsets 2 BIC_c: 4910.484000
--> mjd=50284.000000,  value=4849.095773
MJD=50335.000000, trend=7.058000, BIC_c=4840.925000
For offsets 3 BIC_c: 4840.925000
--> mjd=50335.000000,  value=4767.628136
MJD=51065.000000, trend=16.024000, BIC_c=4722.658000
For offsets 4 BIC_c: 4722.658000
--> mjd=51065.000000,  value=4732.498088
MJD=51066.000000, trend=15.868000, BIC_c=4731.529000
For offsets 5 BIC_c: 4731.529000
--> mjd=51066.000000,  value=4606.235149
0 MJD=       0.0 BIC_c=   5147.76
1 MJD=   50784.0 BIC_c=   5031.35
2 MJD=   51034.0 BIC_c=   4910.48
3 MJD=   50284.0 BIC_c=   4840.93
4 MJD=   50335.0 BIC_c=   4722.66
5 MJD=   51065.0 BIC_c=   4731.53
Computation time in seconds: 15.916143

最后几行列出了含有0-5个偏移时的MJD和BIC_c值，显然含有4个偏移的时候BIC_c最小，从而确定偏移个数和偏移时刻。偏移时刻添加到mom文件头中并保存在./obs_files目录下。

9.4 find_all_offsets.py

find_all_offsets.py用于处理./raw_files目录下的所有站点文件，对每一个站点调用find_offset.py，默认参数是PLWN噪声模型和3D选项。迄今为止，offset探测只实现了AmmarGrag方法而没有实现FullCov，如果控制文件中没有声明NumericalMethod，那么Hector将在缺失数据超过50%的情况下自动选择FullCov。因此，如果用户希望对含有大量缺失的数据探测offset的话就需要在estimatetrend的控制文件中添加一行NumericalMethod AmmarGrag。

9.5 apply_WF.py

维纳滤波方法提供了对于变化季节信号较好的一种估计，apply_WF.py实现了该方法。该脚本利用了核心脚本analyse_timeseries.py以估计常季节信号、去除趋势和估计PLWN的噪声参数，估计的时变季节信号存放在./sea_files目录下，减去时变季节信号的观测值（也即滤波后的信号）存放在./fil_files目录下，调用命令为：

apply_WF.py station_name sigma_a sigma_sa phi

其中站点名称是存放在./obs_files目录下的mom文件名，sigma_a和sigma_sa分别是描述周年和半年信号振幅变化的AR(1)过程的标准差，phi是AR(1)系数，该值通常接近于1也即等同于随机游走过程。注意该脚本中设置缺失为0，意味着季节信号没有随机变化。

10 控制文件

10.1 removeoutliers.ctl

关键字	参数值
DataFile	观测文件名
DataDirectory	观测文件所在目录
OffsetFile	offset信息文件（可选）
OutputFile	带有观测值和估计得模型值的mom格式的输出文件
component	只在neu和enu格式文件用到或者是使用了OffsetFile，可选East、North、Up
interpolate	yes/no
DegreePolynomial	多项式阶数0-6，默认1
estimatemultitrend	yes/no，
estimatepostseismic	yes/no（可选，默认no，如果yes则忽略关键字DegreePolynomial，只估计线性趋势）
estimateslowslipevent	yes/no，默认no
seasonalsignal	yes/no
halfseasonalsignal	yes/no
periodicsignals	一行数字代表以天为单位的周期（可选）
estimateoffsets	yes/no
ScaleFactor	尺度因子（可选，默认1）
PhysicalUnit	观测的物理单位
IQ_factor	四分位距因子

10.2 estimatetrend.ctl和findoffset.ctl

关键字	参数值
DataFile	观测文件名
DataDirectory	观测文件所在目录
OffsetFile	offset信息文件（可选）
OutputFile	带有观测值和估计得模型值的mom格式的输出文件
component	只在neu和enu格式文件用到或者是使用了OffsetFile，可选East、North、Up
interpolate	yes/no
DegreePolynomial	多项式阶数0-6，默认1
estimatemultitrend	yes/no，
estimatepostseismic	yes/no（可选，默认no，如果yes则忽略关键字DegreePolynomial，只估计线性趋势）
estimateslowslipevent	yes/no，默认no
seasonalsignal	yes/no
halfseasonalsignal	yes/no
periodicsignals	一行数字代表以天为单位的周期（可选）
estimateoffsets	yes/no
ScaleFactor	尺度因子（可选，默认1）
PhysicalUnit	观测的物理单位
NoiseModels	从White、FlickerGGM、Ran- domWalkGGM、Powerlaw、PowerlawApprox、ARFIMA、ARMA、GGM中选取任意组合
LikelihoodMethod	可选项，从AmmarGrag、FullCov中选择，缺失少于50%情况下默认AmmarGrag，否则FullCov
AR_p	AR阶数（只对ARFIMA、ARMA有效）
MA_q	MA阶数（只对ARFIMA、ARMA有效）
GGM_1mphi	$1-\varphi$的参数值（只对GGM有效）
TimeNoiseStart	观测起始时间之前的天数，假定为噪声开始的时刻（只对PowerlawApprox有效）
RandomiseFirstGuess	yes/no（可选，默认no）
estimatevaryingseasonal	yes/no（可选，默认no）
varyingseasonal_N	多项式系数（1-8）
estimatemultivariate	yes/no（可选，默认no）
MultiVariateFile	待分析地球物理信号所在文件名称

10.3 estimatespectrum.ctl

关键字	参数值
DataFile	观测文件名
DataDirectory	观测文件所在目录
OffsetFile	offset信息文件（可选）
OutputFile	计算谱结果的输出文件（可选，默认输出到estimatespectrum.out）
interpolate	yes/no
ScaleFactor	尺度因子（可选，默认1）
WindowFunction	Parzen（默认）/Hann
Fraction	0-1（默认0.1）

10.4 modelspectrum.ctl

关键字	参数值
OutputFile	计算谱结果的输出文件（可选，默认输出到modelspectrum.out）
interpolate	yes/no
NoiseModels	从White、FlickerGGM、Ran- domWalkGGM、Powerlaw、PowerlawApprox、ARFIMA、ARMA、GGM中选取任意组合
AR_p	AR阶数（只对ARFIMA、ARMA有效）
MA_q	MA阶数（只对ARFIMA、ARMA有效）
GGM_1mphi	$1-\varphi$的参数值（只对GGM有效）
PhysicalUnit	观测的物理单位

10.5 simulatenoise.ctl

关键字	参数值
SimulationDir	产生的文件存放位置
SimulationLabel	产生文件的基准名称
NumberOfSimulations	模拟个数
NumberOfPoints	模拟长度
SamplingPeriod	以天为单位的采样周期
NoiseModels	从White、FlickerGGM、Ran- domWalkGGM、Powerlaw、PowerlawApprox、ARFIMA、ARMA、GGM中选取任意组合
AR_p	AR阶数（只对ARFIMA、ARMA有效）
MA_q	MA阶数（只对ARFIMA、ARMA有效）
GGM_1mphi	$1-\varphi$的参数值（只对GGM有效）
TimeNoiseStart	观测起始时间之前的天数，假定为噪声开始的时刻（只对PowerlawApprox有效）

11 高级功能

11.1 偏移文件

offset信息除了可以放在观测文件头还可以单独放在一个文件中，这么做的话就需要在控制文件中声明关键字OffsetFile和文件名，以及关键字component。由于大多数人喜欢简单的日-月-年格式，偏移信息文件的格式有所不同，示例如下：

20- 7-1996 NaN 0 0
8- 9-1996  NaN NaN 0
2-12-1997 NaN 0 0
9- 8-1998  NaN 0 0