[QOCA学习笔记]利用QOCA软件PCA模块进行共模误差分析

操作环境：win10+虚拟机ubuntu下QOCA（ver 1.33）
数据：通过解算得到的安徽省CORS参考站坐标时间序列（neu格式）

QOCA软件是优秀的GNSS坐标时间序列分析软件，其简介及PCA模块的使用方法在QOCA online tutor上已经有详细介绍：
附上链接：Principal component analysis.pdf

下面简单介绍利用PCA模块对时间序列的共模误差进行分析的过程（具体PCA的介绍请下载相关文献）

QOCA各模块调用最核心的部分在于驱动文件，输入输出文件和对各参数的设置均在这个驱动文件里面，本次以pca.drv文件为例，给出常规pca分析所用到的驱动文件，并根据introduction文件给出每行参数的注释。（以下文件删除中文注释部分再运行）
注：本次pca运算采用的输入文件是经过QOCA analyze_tseri模块分析后，扣除了阶跃、速度、周年项、半周年项等项的残差序列。具体的分析过程将在另外一篇博文中介绍。

=========================================================================
*          << key word controlled driving file format >>                *
* symbol ":" must exist in command lines as index pointor               *
* any non-blank character at first column means comment line            *
* 第一列存在任何非空的的字符，均认为该行为注释行，不执行                   *
* empty after ":" means comment line too, but line apears in out file   *
=========================================================================job-type:                     regional filtering
c  job-type:                     tectonic signal
*%本次pca分析的工作类型：“regional filtering”意味着所有的构造和季节性条件都被删除。我们的目标是获得共模误差（主要是高频信号）；“tectonic signal”意味着季节性信号被删除。时间序列包含构造信号（速度或跳跃或地震后）。*decomposition method:             PCA
*%选择分解算法：“PCA”：主成分分析；“KLE”:Karhunen-Loeve展开；“SVD”：奇异值分解*apriori value file:             coordinate.solutn
*%（站点位置和速度先验值文件）与之前QOCA各模块的先验文件格式一样*input site list file(site_list):         pca.site
*%（站点列表文件名）在站点列表文件中，第一行是站点总数。其余的行是站点名称。有时用户可能会意外地设置多余的站点名称。pca程序将检查站点名称列表并删除多余的站点名称。*input qob list file(in_list):          pca_data.list
*%（输入文件列表名）在In_list文件中，第一行是数据文件数。其余的行是文件名（包括路径）和文件类型。*
*%type=11:QOCA map文件。（analyze_tseri output residual file format）
%type=12:SOPAC neu文件。
%type=13:一个站点的单组件（n，e，u）文件（欧文的文件）
%type=14:GIPSY stacov format
%type=15:REASoN web service xyz tar file format
%type=16:high-rate coordinate file（map格式）
%默认数据类型为11。*output file:                         pca.out
*%（输出文件名）此文件记录每个站点的输入/输出信息、调整后的偏差、趋势和季节性术语、东、北和向上分量的本征值和功率百分比。如果“adj_6par”选项未激活，本征值和百分比功率仍列在输出文件中。*component file:                      pca.cpt   component
*%（输出主成分文件）第一个输入是文件名，第二个输入是输出类别。如果类别是“component”：输出缩放的主分量值。如果类别是“vector”：输出特征向量值。如果类别是“time series”：分别输出东、b北和天方向的时间序列。
%默认类别为“component”。其他类别仍在建设中。目前，PCA代码只输出前12个主分量。*regional filtered time series file:  pca.res   residual
*%（输出滤波结果文件）第一项是文件名，第二项是输出类别。如果类别是“residual”：输出未过滤、已过滤的值及其残差。*spatial eigenvector file:            pca.seign
*%（空间特征向量输出文件）该文件记录了前10个主分量每个站的标准化空间响应。*network range (sit_range, degree):    0 360 -90 90
*%（网络的地理范围）单位°*solution span (soln_span):           2013.0  2018.5
*%（分析的开始和结束时间点）*do horizontal jointly:              no
*%在“regional filtering”任务中，我们通常分别对东部和北部时间序列进行分析。在这种情况下，此选项应为“否”。然而，在“tectonic signal”任务中，许多用户希望联合分析东、北分量，得到水平主分量。在这种情况下，此选项应为“是”。*fill_gap option:                    yes
*%默认为“是”。如果答案是“否”，程序将禁用补齐缺失值功能。在大多数情况下，填补时间序列缺口是必要的。否则，间隙将被指定为零。它会产生人为的规律，造成误导。如果每个时间序列都没有间隙，则是理想的。然而，在实践中，许多时间序列都有大或小的缺失。如果我们在pca分析中丢弃有间隙的站点，我们将丢失许多站点，有时甚至是大多数站点。如何填补这一空白并不容易，仍然是一个悬而未决的问题。如果任务是“regional filtering”，即时间序列以高频共模误差为主，台站对cme的响应相似。在这种情况下，pca程序使用“叠加”值来填补空白。如果任务是“tectonic signal”，这意味着时间序列中的模式是由构造信号控制的，每个站点的构造信号可能会有一到两个数量级的差异。在这种情况下，pca程序首先将平滑的平均值分配给缺失值。然后进行pca分析。然后利用第一主成分时间序列和台站响应重新填补空白。在此过程之后，重新进行主成分分析，以获得更新的主成分和台站响应。利用更新后的第一主成分及其相应的台站响应来重新填补空白。*adj_6par option:                    no
*%如果答案是“是”，程序将在执行PCA分析之前为每个时间序列调整6个参数（偏差、趋势、周年和半周年）。一般来说，时间序列还有其他信号，比如跳跃。我们建议用户使用analyze_tseri来删除6个参数和其他不相关的信号（如跳跃）。然后使用残差时间序列（analyze_tseri程序的输出）作为pca程序的输入文件。在这种情况下，应关闭adj_6par选项。*outlier_sigma criterion (enu, mm):  30.0 30.0 50.0outlier_value criterion (enu, mm):  200.0 200.0 300.0minimum data percentage (cut_p):    30.0minimum sites for each epoch (cut_t,percentage):   15.0reference frame:                   WGS84reference coordinate, rtime:       geodetic 2013
c  iteration number:                  4end: exit:

准备上述驱动文件中提到的先验文件coordinate.solutn、输入文件（analyze_tseri模块的残差文件）、pca.site、pca_data.list文件，分别为：
(1). 先验文件coordinate.solutn，是在analyze_tseri模块运行完之后在输出文件中提取的更新后的先验文件（在analyze_tseri模块介绍中会介绍如何生成）

(2).输入文件（analyze_tseri模块的残差文件）

(3).pca.site

(4).pca_data.list

设置好驱动文件和前期准备文件即可运行pca pca.drv ，生成的文件有： pca.out、 pca.cpt、 pca.res、 pca.seign文件，分别为：
(1)pca.out

(2) pca.cpt

(3)pca.res

(4) pca.seign

为了方便后面解算完成提取各个测站的信息，在工程目录下运行完pca之后，再利用以下脚本分别新建文件夹，在利用脚本提取信息，现贴出各脚本，其功能在脚本中均有注释。
(1)

declare -a sites=(AQSS AQWJ AQYX CZLA CZMG CZTC FYFN FYJS HBSX HSQM LAHQ LAHS MASM SZDS SZDS SZLB WHWH XCJD XCLX XCNG)
for (( i=0 ; i<18 ; i++ ))
do  mkdir ${sites[$i]}chmod 777 ${sites[$i]}cp -r pca.cpt ./${sites[$i]}cp -r pca.out ./${sites[$i]}cp -r pca.res ./${sites[$i]}cp -r pca.seign ./${sites[$i]}cd ${sites[$i]}chmod 777 pca*cd ../
done

(2)

#从pca.seign提取前三阶主成分对应的标准化特征向量declare -a sites=(AQSS AQWJ AQYX CZLA CZMG CZTC FYFN FYJS HBSX HSQM LAHQ LAHS MASM SZDS SZDS SZLB WHWH XCJD XCLX XCNG)
for (( i=0 ; i<18 ; i++ ))
do cp ${sites[$i]}/pca.seign pca.seignsed -i '/*/d' pca.seign                  #删除含有“*”的那一行echo ${sites[$i]} >> 1st2nd3rd_eVector.txtcat pca.seign | awk '{print $3 " " $4 " " $5 " " $6}' >> 1st2nd3rd_eVector.txtrm pca.seign
done

(3)

#列出前三阶主成分特征值占总特征值的百分比echo '1st order' >> 1st2nd3rd_per.txt
grep '   1   ' pca.out >> 1st2nd3rd_per.txtecho '2nd order' >> 1st2nd3rd_per.txt
grep '   2   ' pca.out >> 1st2nd3rd_per.txtecho '3rd order' >> 1st2nd3rd_per.txt
grep '   3   ' pca.out >> 1st2nd3rd_per.txt

(4)

#!/bin/bash
#从滤波后残差序列输出文件pca.res中提取各站的cme序列
mkdir cme
mkdir filtered
mkdir unfiltered
declare -a sites=(AQSS AQWJ AQYX CZLA CZMG CZTC FYFN FYJS HBSX HSQM LAHQ LAHS MASM SZDS SZDS SZLB WHWH XCJD XCLX XCNG)
for (( i=0 ; i<18 ; i++ ))
do  cd ${sites[$i]}grep ${sites[$i]}_GPS pca.res >> ${sites[$i]}_temp.txtgrep north ${sites[$i]}_temp.txt | awk '{print $1 " " $4}' >> ${sites[$i]}_n_cme.txtgrep east_ ${sites[$i]}_temp.txt | awk '{print $4}'        >> ${sites[$i]}_e_cme.txtgrep  up__ ${sites[$i]}_temp.txt | awk '{print $4}'        >> ${sites[$i]}_u_cme.txt          #cme序列paste -d" " ${sites[$i]}_n_cme.txt ${sites[$i]}_e_cme.txt ${sites[$i]}_u_cme.txt >> ${sites[$i]}_cme.txtmv ${sites[$i]}_cme.txt ../cmegrep north ${sites[$i]}_temp.txt | awk '{print $1 " " $3}' >> ${sites[$i]}_n_filtered.txtgrep east_ ${sites[$i]}_temp.txt | awk '{print $3}'        >> ${sites[$i]}_e_filtered.txtgrep  up__ ${sites[$i]}_temp.txt | awk '{print $3}'        >> ${sites[$i]}_u_filtered.txt     #滤波后的残差序列paste -d" " ${sites[$i]}_n_filtered.txt ${sites[$i]}_e_filtered.txt ${sites[$i]}_u_filtered.txt >> ${sites[$i]}_filtered.txtmv ${sites[$i]}_filtered.txt ../filteredgrep north ${sites[$i]}_temp.txt | awk '{print $1 " " $2}' >> ${sites[$i]}_n_unfiltered.txtgrep east_ ${sites[$i]}_temp.txt | awk '{print $2}'        >> ${sites[$i]}_e_unfiltered.txtgrep  up__ ${sites[$i]}_temp.txt | awk '{print $2}'        >> ${sites[$i]}_u_unfiltered.txt   #滤波前的序列paste -d" " ${sites[$i]}_n_unfiltered.txt ${sites[$i]}_e_unfiltered.txt ${sites[$i]}_u_unfiltered.txt >> ${sites[$i]}_unfiltered.txtmv ${sites[$i]}_unfiltered.txt ../unfilteredrm *.txtcd ..
done

参考文献：
[1]纪海源. 基于GAMIT/GLOBK与QOCA的汾渭断陷带地壳运动研究[D]. 西安科技大学, 2014.
[2]QOCA官网online toturing
感谢昆明理工大学李建涛师兄提供的学习材料

[QOCA学习笔记]利用QOCA软件PCA模块进行共模误差分析相关推荐

Python学习笔记：常用内建模块6 (urllib)
前言最近在学习深度学习,已经跑出了几个模型,但Pyhton的基础不够扎实,因此,开始补习Python了,大家都推荐廖雪峰的课程,因此,开始了学习,但光学有没有用,还要和大家讨论一下,因此,写下这些帖 ...
Python学习笔记：常用内建模块4：hmac
前言最近在学习深度学习,已经跑出了几个模型,但Pyhton的基础不够扎实,因此,开始补习Python了,大家都推荐廖雪峰的课程,因此,开始了学习,但光学有没有用,还要和大家讨论一下,因此,写下这些帖 ...
Python学习笔记：常用内建模块3：struct
前言最近在学习深度学习,已经跑出了几个模型,但Pyhton的基础不够扎实,因此,开始补习Python了,大家都推荐廖雪峰的课程,因此,开始了学习,但光学有没有用,还要和大家讨论一下,因此,写下这些帖 ...
【小猫爪】AUTOSAR学习笔记03-Communication Stack之CanIf模块
[小猫爪]AUTOSAR学习笔记03-Communication Stack之CanIf模块前言 1 CanIf简介 2 功能介绍 2.1 发送缓冲区 2.2 CAN Controller 模式控制 ...
【小猫爪】AUTOSAR学习笔记16-Memory Stack之Nvm模块
[小猫爪]AUTOSAR学习笔记16-Memory Stack之Nvm模块前言 1 Memory Stack结构 2 NvM模块 2.1 Block 2.2 数据结构 2.3 Redundant B ...
sizzle.js学习笔记利用闭包模拟实现数据结构：字典（Map）
sizzle.js学习笔记利用闭包模拟实现数据结构:字典(Map) 这几天学习和查看了jQuery和Property这两个很流行的前端库的御用选择器组件Sizzle.js的源代码,收获还是相对多的!之 ...
Arduino模块学习笔记（一）—GPS模块的使用
Arduino模块学习笔记(一)--GSP模块的使用文章目录 Arduino模块学习笔记(一)--GSP模块的使用所需组件一.模块使用介绍 1.GPS模块(在室内时,一般获取不到位置信息) 2. ...
DSP(TMSF280049C)学习笔记1：软件的安装与新工程的建立
DSP(TMSF280049C)学习笔记1:软件的安装与新工程的建立初次接触DSP(TMSF280049C)这款开发板,主要对学习的过程进行记录,因为是初次学习,可能有些地方存在错误,还请各位大神帮 ...
InSAR学习笔记之ISCE 软件安装
InSAR学习笔记之ISCE 软件安装 ISCE是一款常用的InSAR数据处理软件,2018年更新的版本基于ubuntu18.04系统安装过程简化了很多,本文分享一下安装过程.(之前在ubuntu16 ...

[QOCA学习笔记]利用QOCA软件PCA模块进行共模误差分析

[QOCA学习笔记]利用QOCA软件PCA模块进行共模误差分析相关推荐

最新文章

热门文章