[Hector学习笔记]GNSS时间序列处理软件Hector使用备忘（批处理脚本）

由Machiel Bos 和Rui Fernandes共同研发的Hector软件可以利用时间相关噪声来估计时间序列线性趋势，是与优秀的GNSS时间序列处理软件之一。
软件安装教程：Hector软件安装[附加软件环境变量设置]
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 BJSH BJYQ... )
for ((i=0;i<${#sites[*]};i++))
do
	site=${sites[$i]}
	n_raw=$site.neu
	c_raw=`grep "DataFile" removeoutliers.ctl | awk '{print $2}'`
#	终端显示此时正在处理的文件，以便修改
#	echo "$c_raw $n_raw"
#	将驱动控制文件中文件名修改成要处理的文件名
        sed -i "s/$c_raw/$n_raw/g" removeoutliers.ctl 
#	列出所有方向
	directions=(North East Up)
	for ((j=0;j<${#directions[*]};j++))
	do
		direction=${directions[$j]}
		c_dir=`grep "component" removeoutliers.ctl | awk '{print $2}'`
#		将驱动控制文件中坐标分量改成要处理的分量，并通过循环实现三个坐标分量的替换
		sed -i "s/$c_dir/$direction/g" removeoutliers.ctl
		out_file=$site.outliers.$direction.mom
#		修改控制文件中的输出文件，此时输出文件需要区分开站的各个方向
		c_out_file=`grep "OutputFile" removeoutliers.ctl | awk '{print $2}'`
		echo "$c_out_file"
		sed -i "s/$c_out_file/$out_file/g" removeoutliers.ctl
		removeoutliers removeoutliers.ctl
	done
done

【注】脚本中测站名的制作方法将单独介绍。如果要处理的测站有好几百个，列出测站名不能采用手动输入的方式，会有更好的办法进行制作。
更新：制作方法在这：GNSS时间序列分析笔记Ubuntu下制作站点列表
上述脚本可以放在与要处理的文件同一个文件夹下，当然也可以放在其他文件夹下，在removeoutliers.ctl文件中有对文件位置做详细的设置，比如我的是这样：（此时脚本与含有待处理文件的deformation_make在同一文件夹下，这样的逻辑关系可以自行体会）

经过处理后生成的文件与脚本在同一个文件夹下，生成的文件：

是Hector要求的mom格式。接下来进行趋势项估计。
2. 趋势项估计
输入文件格式：[site].outliers.[direction].mom
输出文件格式：[site].[direction].[noise].trend.mom
说明：在估计趋势项和周期项时，可以加入各种噪声以及噪声组合，此时输出文件需要对不同噪声进行区分。且输出问价后缀仍需是mom格式，以便后期的功率谱估计。
趋势项估计的脚本如下：

#!/bin/bash
#trend estimation by Hector
#列出所有测站(由于测站较多，省略了大部分)
sites=(AHAQ AHBB BJGB BJSH BJYQ... )
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.$m_name.trend.mom
			c_out_file=`grep "OutputFile" estimatetrend1.ctl | awk '{print $2}'`
			sed -i "s/$c_out_file/$out_file/g" estimatetrend.ctl
			echo "Outputting $out_file"
			c_model=`grep "NoiseModels" estimatetrend1.ctl | awk '{print $2}'`
			sed -i "s/$c_model/$model/g" estimatetrend.ctl
			re_file=$site.$direction.$m_name
			estimatetrend estimatetrend1.ctl > $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.outliers.Up.mom
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.Up.wh.trend.mom
Total computing time: 8.00000 sec

此时利用下面的脚本完成振幅相位的提取（我在解算时仅估计了周年项，若想估计半周年项，需要在ctl文件中进行相关设置）

#!/bin/bash
#grep amplitude and phase from log file 

sites=(AHAQ AHBB BJGB BJSH BJYQ CHUN CQCS CQWZ DLHA DXIN FJPT FJWY FJXP GDSG GDST GDZH GDZJ GSAX GSDH GSDX GSGL GSGT GSJN GSJT GSJY GSLX GSLZ GSMA GSML GSMQ GSMX GSPL GSQS GSTS GSWD GUAN GXBH GXBS GXGL GXHC GXNN GXWZ GZFG GZGY GZSC HAHB HAJY HAQS HBES HBJM HBXF HBZG HECC HECD HECX HELQ HELY HETS HEYY HEZJ HIHK HISY HLAR HLFY HLHG HLMH HLWD HNLY HNMY HRBN JIXN JLCB JLCL JLYJ JSLS JSLY JSNT JSYC JXHK JXJA KMIN LALB LALX LHAS LNDD LNHL LNJZ LNSY LNYK LUZH NMAG NMAL NMAY NMAZ NMBT NMDW NMEJ NMEL NMER NMTK NMWH NMWJ NMWL NMWT NMZL NXHY NXYC NXZW QHBM QHDL QHGC QHGE QHLH QHMD QHME QHMQ QHMY QHQL QHTR QHTT QHYS QION SCBZ SCDF SCGZ SCJL SCJU SCLH SCLT SCMB SCML SCMN SCMX SCNC SCNN SCPZ SCSM SCSN SCSP SCTQ SCXC SCXD SCXJ SCYX SCYY SDCY SDJX SDLY SDQD SDRC SDYT SDZB SHA2 SNAK SNHY SNMX SNTB SNXY SNYA SUIY SXCZ SXDT SXGX SXKL SXLF SXLQ SXTY SXXX SXYC TAIN TJBD TJBH TJWQ WUSH XIAG XIAM XJAL XJBC XJBE XJBL XJBY XJDS XJFY XJHT XJJJ XJKC XJKE XJKL XJML XJQH XJQM XJRQ XJSH XJSS XJTC XJTZ XJWL XJWQ XJWU XJXY XJYC XJYN XJYT XJZS XNIN XZAR XZBG XZCD XZCY XZDX XZGE XZGZ XZNM XZRK XZRT XZYD XZZB XZZF YANC YNCX YNDC YNGM YNHZ YNJD YNJP YNLA YNLC YNLJ YNMH YNMJ YNML YNMZ YNRL YNSD YNSM YNTC YNTH YNWS YNXP YNYA YNYL YNYM YNYS YNZD ZHNZ ZJJD ZJWZ ZJZS )
#for ((i=0;i<${#sites[*]};i++))
for ((i=0;i<2;i++))
do
	site=${sites[$i]}_GPS
	file_name=$site.Up.wh
	echo "new doing $file_name..."
	cos=`grep "cos yearly" $file_name | awk '{print $3 $4 $5}'`
	sin=`grep "sin yearly" $file_name | awk '{print $3 $4 $5}'`
	amp=`grep "Amp yearly" $file_name | awk '{print $3 $4 $5}'`
	pha=`grep "Pha yearly" $file_name | awk '{print $4}'`
	printf "%-1s %-1s %-1s %-1s\n" $file_name $cos $sin $amp $pha>> amp.dat
done

提取出来的结果文件如下：

可以对以上文件在excel中进行下一步分析。
3. 功率谱估计
按照之前的方法先用脚本运行estimatespectrum，脚本如下：（与上面脚本所用到的测站不一样）

#!/bin/bash

sites=(AQSS AQWJ AQYX CZLA CZMG CZTC FYFN FYJS HBSX HSQM LAHQ LAHS MASM SZDS SZLB WHWH XCJD XCLX XCNG)
for ((i=0;i<${#sites[*]};i++))
do
	site=${sites[$i]}
	directions=(North East Up)
	dirs=(N E U)
	for ((j=0;j<${#directions[*]};j++))
	do
		direction=${directions[$j]}
		dir=${dirs[$j]}
#		m_names=(wh wh_fn wh_rw wh_pl wh_fn_rw)
#		for((k=0;k<${#m_names[*]};k++))
#		do
#			m_name=${m_names[$k]}
			file_name=$site.$direction.wh_fn.trend.mom
			out_file=$site.$dir.wh_fn.spectrum
			log_file=$site.$direction.wh_fn.log
			c_file=`grep "DataFile" estimatespectrum.ctl | awk '{print $2}'`
			sed -i "s/$c_file/$file_name/g" estimatespectrum.ctl
			c_out_file=`grep "OutputFile" estimatespectrum.ctl | awk '{print $2}'`
			sed -i "s/$c_out_file/$out_file/g" estimatespectrum.ctl
			estimatespectrum estimatespectrum.ctl > $log_file
#		done
	done
done 

生成的文件有：