操作環(huán)境：win10+虛擬機ubuntu下QOCA（ver 1.33）
數(shù)據(jù)：通過解算得到的安徽省CORS參考站坐標時間序列（neu格式）

QOCA軟件是優(yōu)秀的GNSS坐標時間序列分析軟件，其簡介及PCA模塊的使用方法在QOCA online tutor上已經(jīng)有詳細介紹：
附上鏈接：Principal component analysis.pdf

下面簡單介紹利用PCA模塊對時間序列的共模誤差進行分析的過程（具體PCA的介紹請下載相關文獻）

QOCA各模塊調(diào)用最核心的部分在于驅(qū)動文件，輸入輸出文件和對各參數(shù)的設置均在這個驅(qū)動文件里面，本次以pca.drv文件為例，給出常規(guī)pca分析所用到的驅(qū)動文件，并根據(jù)introduction文件給出每行參數(shù)的注釋。（以下文件刪除中文注釋部分再運行）
注：本次pca運算采用的輸入文件是經(jīng)過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 * * 第一列存在任何非空的的字符，均認為該行為注釋行，不執(zhí)行 * * empty after ":" means comment line too, but line apears in out file * =========================================================================job-type: regional filtering c job-type: tectonic signal *%本次pca分析的工作類型：“regional filtering”意味著所有的構造和季節(jié)性條件都被刪除。我們的目標是獲得共模誤差（主要是高頻信號）；“tectonic signal”意味著季節(jié)性信號被刪除。時間序列包含構造信號（速度或跳躍或地震后）。*decomposition method: PCA *%選擇分解算法：“PCA”：主成分分析；“KLE”:Karhunen-Loeve展開；“SVD”：奇異值分解*apriori value file: coordinate.solutn *%（站點位置和速度先驗值文件）與之前QOCA各模塊的先驗文件格式一樣*input site list file(site_list): pca.site *%（站點列表文件名）在站點列表文件中，第一行是站點總數(shù)。其余的行是站點名稱。有時用戶可能會意外地設置多余的站點名稱。pca程序?qū)z查站點名稱列表并刪除多余的站點名稱。*input qob list file(in_list): pca_data.list *%（輸入文件列表名）在In_list文件中，第一行是數(shù)據(jù)文件數(shù)。其余的行是文件名（包括路徑）和文件類型。* *%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格式） %默認數(shù)據(jù)類型為11。*output file: pca.out *%（輸出文件名）此文件記錄每個站點的輸入/輸出信息、調(diào)整后的偏差、趨勢和季節(jié)性術語、東、北和向上分量的本征值和功率百分比。如果“adj_6par”選項未激活，本征值和百分比功率仍列在輸出文件中。*component file: pca.cpt component *%（輸出主成分文件）第一個輸入是文件名，第二個輸入是輸出類別。如果類別是“component”：輸出縮放的主分量值。如果類別是“vector”：輸出特征向量值。如果類別是“time series”：分別輸出東、b北和天方向的時間序列。 %默認類別為“component”。其他類別仍在建設中。目前，PCA代碼只輸出前12個主分量。*regional filtered time series file: pca.res residual *%（輸出濾波結(jié)果文件）第一項是文件名，第二項是輸出類別。如果類別是“residual”：輸出未過濾、已過濾的值及其殘差。*spatial eigenvector file: pca.seign *%（空間特征向量輸出文件）該文件記錄了前10個主分量每個站的標準化空間響應。*network range (sit_range, degree): 0 360 -90 90 *%（網(wǎng)絡的地理范圍）單位°*solution span (soln_span): 2013.0 2018.5 *%（分析的開始和結(jié)束時間點）*do horizontal jointly: no *%在“regional filtering”任務中，我們通常分別對東部和北部時間序列進行分析。在這種情況下，此選項應為“否”。然而，在“tectonic signal”任務中，許多用戶希望聯(lián)合分析東、北分量，得到水平主分量。在這種情況下，此選項應為“是”。*fill_gap option: yes *%默認為“是”。如果答案是“否”，程序?qū)⒔醚a齊缺失值功能。在大多數(shù)情況下，填補時間序列缺口是必要的。否則，間隙將被指定為零。它會產(chǎn)生人為的規(guī)律，造成誤導。如果每個時間序列都沒有間隙，則是理想的。然而，在實踐中，許多時間序列都有大或小的缺失。如果我們在pca分析中丟棄有間隙的站點，我們將丟失許多站點，有時甚至是大多數(shù)站點。如何填補這一空白并不容易，仍然是一個懸而未決的問題。如果任務是“regional filtering”，即時間序列以高頻共模誤差為主，臺站對cme的響應相似。在這種情況下，pca程序使用“疊加”值來填補空白。如果任務是“tectonic signal”，這意味著時間序列中的模式是由構造信號控制的，每個站點的構造信號可能會有一到兩個數(shù)量級的差異。在這種情況下，pca程序首先將平滑的平均值分配給缺失值。然后進行pca分析。然后利用第一主成分時間序列和臺站響應重新填補空白。在此過程之后，重新進行主成分分析，以獲得更新的主成分和臺站響應。利用更新后的第一主成分及其相應的臺站響應來重新填補空白。*adj_6par option: no *%如果答案是“是”，程序?qū)⒃趫?zhí)行PCA分析之前為每個時間序列調(diào)整6個參數(shù)（偏差、趨勢、周年和半周年）。一般來說，時間序列還有其他信號，比如跳躍。我們建議用戶使用analyze_tseri來刪除6個參數(shù)和其他不相關的信號（如跳躍）。然后使用殘差時間序列（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:

準備上述驅(qū)動文件中提到的先驗文件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

設置好驅(qū)動文件和前期準備文件即可運行pca pca.drv ，生成的文件有： pca.out、 pca.cpt、 pca.res、 pca.seign文件，分別為：
(1)pca.out

(2) pca.cpt

(3)pca.res

(4) pca.seign

為了方便后面解算完成提取各個測站的信息，在工程目錄下運行完pca之后，再利用以下腳本分別新建文件夾，在利用腳本提取信息，現(xiàn)貼出各腳本，其功能在腳本中均有注釋。
(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官網(wǎng)online toturing
感謝昆明理工大學李建濤師兄提供的學習材料

總結(jié)

以上是生活随笔為你收集整理的[QOCA学习笔记]利用QOCA软件PCA模块进行共模误差分析的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。