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

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

下面簡(jiǎn)單介紹利用PCA模塊對(duì)時(shí)間序列的共模誤差進(jìn)行分析的過(guò)程（具體PCA的介紹請(qǐng)下載相關(guān)文獻(xiàn)）

QOCA各模塊調(diào)用最核心的部分在于驅(qū)動(dòng)文件，輸入輸出文件和對(duì)各參數(shù)的設(shè)置均在這個(gè)驅(qū)動(dòng)文件里面，本次以pca.drv文件為例，給出常規(guī)pca分析所用到的驅(qū)動(dòng)文件，并根據(jù)introduction文件給出每行參數(shù)的注釋。（以下文件刪除中文注釋部分再運(yùn)行）
注：本次pca運(yùn)算采用的輸入文件是經(jīng)過(guò)QOCA analyze_tseri模塊分析后，扣除了階躍、速度、周年項(xiàng)、半周年項(xiàng)等項(xiàng)的殘差序列。具體的分析過(guò)程將在另外一篇博文中介紹。

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

準(zhǔn)備上述驅(qū)動(dòng)文件中提到的先驗(yàn)文件coordinate.solutn、輸入文件（analyze_tseri模塊的殘差文件）、pca.site、pca_data.list文件，分別為：
(1). 先驗(yàn)文件coordinate.solutn，是在analyze_tseri模塊運(yùn)行完之后在輸出文件中提取的更新后的先驗(yàn)文件（在analyze_tseri模塊介紹中會(huì)介紹如何生成）

(2).輸入文件（analyze_tseri模塊的殘差文件）

(3).pca.site

(4).pca_data.list

設(shè)置好驅(qū)動(dòng)文件和前期準(zhǔn)備文件即可運(yùn)行pca pca.drv ，生成的文件有： pca.out、 pca.cpt、 pca.res、 pca.seign文件，分別為：
(1)pca.out

(2) pca.cpt

(3)pca.res

(4) pca.seign

為了方便后面解算完成提取各個(gè)測(cè)站的信息，在工程目錄下運(yùn)行完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提取前三階主成分對(duì)應(yīng)的標(biāo)準(zhǔn)化特征向量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

參考文獻(xiàn)：
[1]紀(jì)海源. 基于GAMIT/GLOBK與QOCA的汾渭斷陷帶地殼運(yùn)動(dòng)研究[D]. 西安科技大學(xué), 2014.
[2]QOCA官網(wǎng)online toturing
感謝昆明理工大學(xué)李建濤師兄提供的學(xué)習(xí)材料

總結(jié)

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

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