[Hector学习笔记]GNSS时间序列处理软件Hector使⽤备忘（批处理脚本） 由Machiel Bos 和Rui Fernandes共同研发的Hector软件可以利⽤时间相关噪声来估计时间序列线性趋势，是与优秀的GNSS时间序列处理软件之⼀。 软件安装教程：
Hector软件⽤户⼿册中包含了教程，⽰例⽂件等等，并且有⾮常详细的说明。在教程中，⼤多采⽤“设置驱动⽂件+在终端中逐⾏输⼊代码”进⾏计算的⽅式进⾏解算，这适⽤于新⼿，或者是需要处理的站点较少的情况。对于超过5个站点的时间序列的分析，如果还⽤逐个修改驱动⽂件和终端逐个输⼊的⽅式，未免太浪费时间。在学习该软件的过程中，逐步总结出批量处理站点时间序列的脚本，以备⽇后查询和学习。根据教程中的处理流程，分别从：剔除粗差、估计趋势项和估计功率谱密度三个⽅⾯进⾏总结，以上三部分具体的参数设置在⽤户⼿册中均有说明，此处不再赘述。
1. 剔除粗差
输⼊⽂件格式：[site].neu
输出⽂件格式：[site].outliers.[direction].mom（此处⽂件名称“outliers”可⾃定义，[direction]⽤来表⽰⽅向，mom是固定的）说明：由于Hector处理测站各个⽅向的时间序列是分开的，所以每处理⼀个⽅向所
⽣成的⽂件名对⽅向的标注是⾮常有必要的；另外上述输出格式中⽂件名称“outliers”部分可以⾃⼰修改成任何其他名称，只是为了后⾯更⽅便的辨认；最后，重要的是Hector后⾯的步骤的输⼊⽂件均是mom格式，所以最后⼀定要以mom格式结尾。
下⾯是脚本：
#!/bin/bash
#remove outliers use Hector
#列出所有测站(由于测站较多，省略了⼤部分)
sites=(AHAQ AHBB BJGB )
for((i=0;i<${#sites[*]};i++))
do
site=${sites[$i]}
n_raw=$u
c_raw=`grep"DataFile" l |awk'{print $2}'`
# 终端显⽰此时正在处理的⽂件，以便修改
# echo "$c_raw $n_raw"
# 将驱动控制⽂件中⽂件名修改成要处理的⽂件名
sed -i "s/$c_raw/$n_raw/g" l
# 列出所有⽅向
directions=(North East Up)
for((j=0;j<${#directions[*]};j++))
do
direction=${directions[$j]}
c_dir=`grep"component" l |awk'{print $2}'`
# 将驱动控制⽂件中坐标分量改成要处理的分量，并通过循环实现三个坐标分量的替换
sed -i "s/$c_dir/$direction/g" l
out_file=$site.outliers.$
scm文件# 修改控制⽂件中的输出⽂件，此时输出⽂件需要区分开站的各个⽅向
c_out_file=`grep"OutputFile" l |awk'{print $2}'`
echo"$c_out_file"
sed -i "s/$c_out_file/$out_file/g" l
l
done
done
【注】脚本中测站名的制作⽅法将单独介绍。如果要处理的测站有好⼏百个，列出测站名不能采⽤⼿动输⼊的⽅式，会有更好的办法进⾏制作。
更新：制作⽅法在这：
上述脚本可以放在与要处理的⽂件同⼀个⽂件夹下，当然也可以放在其他⽂件夹下，在l⽂件中有对⽂件位置做详细的设置，⽐如我的是这样：（此时脚本与含有待处理⽂件的deformation_make在同⼀⽂件夹下，这样的逻辑关系可以⾃⾏体会）
经过处理后⽣成的⽂件与脚本在同⼀个⽂件夹下，⽣成的⽂件：
是Hector要求的mom格式。接下来进⾏趋势项估计。
2. 趋势项估计
输⼊⽂件格式：[site].outliers.[direction].mom
输出⽂件格式：[site].[direction].[noise].
说明：在估计趋势项和周期项时，可以加⼊各种噪声以及噪声组合，此时输出⽂件需要对不同噪声进⾏区分。且输出问价后缀仍需是mom 格式，以便后期的功率谱估计。
趋势项估计的脚本如下：
#!/bin/bash
#trend estimation by Hector
#列出所有测站(由于测站较多，省略了⼤部分)
sites=(AHAQ AHBB BJGB )
for((i=0;i<${#sites[*]};i++))
do
site=${sites[$i]}
directions=(North East Up)
for((j=0;j<${#directions[*]};j++))
do
direction=${directions[$j]}
m_names=(wh wh_fn wh_rw wh_pl wh_fn_rw)
models=("White""White FlickerGGM""White RandomWalkGGM""White Powerlaw""White FlickerGGM RandomWalkGGM") for((k=0;k<${#m_names[*]};k++))
do
m_name=${m_names[$k]}
model=${models[$k]}
out_file=$site.$direction.$d.mom
c_out_file=`grep"OutputFile" l |awk'{print $2}'`
sed -i "s/$c_out_file/$out_file/g" l
echo"Outputting $out_file"
c_model=`grep"NoiseModels" l |awk'{print $2}'`
sed -i "s/$c_model/$model/g" l
re_file=$site.$direction.$m_name
l >$re_file
done
done
done
以上为完整的带坐标分量和噪声组合循环的脚本，我在解算是仅⽤到了测站的up⽅向，且仅考虑⽩噪声影响。如有类似这样的需要，可以将相应的循环部分取消，令其直接等于某个分量、噪声模型即可。此时解算得到的⽂件有mom⽂件和.wh⽂件（只考虑⽩噪声）。⽂件分别长这样：
这⾥的.wh⽂件是estimatetrend的log⽂件，⾥⾯记载了解算的结果，包括估计噪声的最⼤似然估计值、周年与半周年项等等信息，此时如果我们需要提取其中的信息，仍然可以采⽤脚本的⽅法。下⾯先以⼀个log⽂件为例：
************************************
estimatetrend, version 1.6
************************************
0) White
generator type: mt19937
seed = 0
first value = 3176401122
Data format: MJD, Observations, Model
Filename : ./AHAQ_GPS.
Number of observations: 5569
Percentage of gaps : 0.8799
----------------
AmmarGrag
----------------
No Polynomial degree set, using offset + linear trend
No extra periodic signal is included.
Number of CPU's used (threads) = 1
performing Ordinary Least-Squares
White:
fraction = 1.00000
sigma = 2.88355 mm
No noise parameters to show
Likelihood value
--------------------
min log(L)=-13678.350
k =4 + 0 + 1 = 5
AIC =27366.701
BIC =27399.781
STD of the driving noise: 2.884
bias : -0.085 +/- 0.039 mm (at 2009/11/13, 12:0:0.000)
trend: -0.343 +/- 0.009 mm/year
cos yearly : -2.982 +/- 0.055 mm
sin yearly : 1.843 +/- 0.055 mm
Amp yearly : 3.506 +/- 0.055 mm
Pha yearly : 148.276 degrees
--> AHAQ_GPS.d.mom
Total computing time: 8.00000 sec
此时利⽤下⾯的脚本完成振幅相位的提取（我在解算时仅估计了周年项，若想估计半周年项，需要在ctl⽂件中进⾏相关设置）
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议