簡介這是一套面向電離層研究者、GNSS數(shù)據(jù)處理工程師及空間天氣分析人員的Fortran科學(xué)計算源碼專注于將標(biāo)準RINEX格式的GNSS觀測數(shù)據(jù)高精度反演為電離層總電子含量TEC并輸出為專用GTEX格式有效支撐定位誤差校正、電離層建模與太陽活動影響評估等關(guān)鍵任務(wù)。資源共42個文件含15個Fortran源文件.f實現(xiàn)核心算法如TEC反演、軌道讀取、周跳修正、15個編譯目標(biāo)文件.o、3個Shell腳本.sh用于自動化流程控制以及Makefile、頭文件、參數(shù)配置列表和詳細README文檔結(jié)構(gòu)完整、模塊清晰便于二次開發(fā)與教學(xué)實踐。壓縮包僅290KB輕量但功能完備。目前已有155人學(xué)習(xí)下載用戶可直接獲取可編譯運行的完整工程涵蓋從RINEX數(shù)據(jù)解析、衛(wèi)星幾何定位、偽距/相位組合TEC計算到GTEX格式寫入的全鏈路實現(xiàn)并附帶Julian日歷轉(zhuǎn)換、IGS周數(shù)計算等實用工具模塊。 從源碼讀懂GNSS電離層處理這事放到現(xiàn)在依然有吸引力。RNX2GTEX這個名字常做GNSS數(shù)據(jù)處理的人一看就明白它做的是把RINEX觀測文件轉(zhuǎn)成TEX格式輸出。這里的TEX不是LaTeX那套排版工具而是電離層TEC交換文件Total Electron Content Exchange的簡寫。整個程序的核心任務(wù)就是從雙頻GNSS觀測值中提取總電子含量TEC再把逐衛(wèi)星、逐歷元的TEC結(jié)果寫成標(biāo)準格式供后續(xù)電離層建模、單頻用戶電離層改正、空間天氣研究使用。我用這個程序處理過不少測站數(shù)據(jù)也把源碼從頭到尾翻過幾遍坦白說這不是一個代碼風(fēng)格很“現(xiàn)代”的項目但它把GNSS電離層測量的整條鏈路用最樸素的方式跑通了。這篇內(nèi)容我就從一個實際使用者的角度把RNX2GTEX涉及的物理原理、源碼結(jié)構(gòu)、編譯運行細節(jié)、輸出格式和常見坑都拆開講一遍適合剛接觸GNSS電離層數(shù)據(jù)處理、或者想通過老Fortran代碼理解觀測方程的同學(xué)參考。1. 先把名詞理清楚RINEX、TEX、TEC各自扮演什么角色很多初學(xué)者看到RNX2GTEX這個軟件名第一反應(yīng)是去查怎么編譯結(jié)果被RINEX文件名、TEX格式、TEC單位這些概念繞暈。我建議先花半小時把下面幾個名詞的關(guān)系理清后面看代碼會輕松很多。1.1 TEC是什么GNSS為什么能測到電離層TEC的全稱是Total Electron Content翻譯過來就是總電子含量指沿信號傳播路徑上單位截面積柱體內(nèi)的自由電子總數(shù)單位是電子數(shù)每平方米。平時我們更常用TECU做單位1 TECU等于10的16次方個電子每平方米。電離層里的自由電子會對GNSS信號產(chǎn)生折射延遲這個延遲的大小和信號頻率的平方成反比。對頻率為f的信號偽距上的電離層延遲近似為40.3乘STEC除以f的平方這里的STEC就是斜路徑上的TEC。因為兩個載波頻率不一樣同一顆衛(wèi)星發(fā)出的L1和L2信號穿過電離層時延遲量不同通過對比兩個頻率的觀測值就能把電離層影響單獨分離出來。我在實際給學(xué)生講的時候喜歡用一個類比電離層就像一塊有色玻璃不同顏色的光穿過時速度不一樣。GNSS接收機同時收到兩個頻率的信號相當(dāng)于拿兩束不同顏色的光同時穿過這塊玻璃通過比較它們的到達時間差就能反推玻璃的“厚度”。這里的“厚度”就是TEC。1.2 從RINEX到TEXRNX2GTEX在流水線中的位置RINEX是GNSS觀測數(shù)據(jù)的標(biāo)準交換格式接收機廠商輸出的原始數(shù)據(jù)經(jīng)過轉(zhuǎn)換后會以RINEX格式保存?zhèn)尉?、載波相位、多普勒等觀測值。RINEX文件里包含的信息非常豐富但直接拿RINEX文件做電離層研究并不方便因為觀測值里混著鐘差、軌道誤差、對流層延遲等一大堆無關(guān)量而且不同接收機輸出的文件格式細節(jié)也不一樣。TEX格式則是專門面向電離層TEC數(shù)據(jù)設(shè)計的交換格式它把處理好的TEC結(jié)果按站點、按歷元、按衛(wèi)星組織起來省去了用戶重復(fù)做預(yù)處理的工作。RNX2GTEX就是連接RINEX和TEX的橋梁輸入一個或多個RINEX觀測文件經(jīng)過質(zhì)量檢查和組合計算輸出TEX格式的TEC時間序列。整個GNSS電離層處理鏈路大致是RINEX觀測文件加精密星歷和鐘差經(jīng)過預(yù)處理、周跳探測、無幾何組合、相位平滑等步驟得到STEC再做硬件延遲校正和映射函數(shù)轉(zhuǎn)換輸出網(wǎng)格化的VTEC或直接輸出STEC序列。RNX2GTEX覆蓋的是從RINEX到STEC/TEC這一核心段落。1.3 單位轉(zhuǎn)換和數(shù)量級怎么判斷算出來的TEC是否合理用這個程序之前最好對TEC的數(shù)量級心里有數(shù)。中緯度地區(qū)平靜電離層情況下天頂方向VTEC通常在10到50 TECU之間低緯赤道異常區(qū)可以到80甚至100 TECU以上高緯和夜間會低很多夜間經(jīng)常只有幾個TECU。如果算出來的結(jié)果在一個測站、一整天內(nèi)變化幾百上千個TECU那基本可以判斷數(shù)據(jù)或處理方法出了問題。從幾何關(guān)系上也要有概念。斜路徑STEC一般是VTEC的好幾倍仰角越低路徑越長STEC越大。仰角30度左右時如果VTEC是30 TECUSTEC大概在50到60 TECU。程序輸出的如果明顯偏離這個范圍就要回頭檢查組合系數(shù)、單位換算是哪里出了問題。還有一個非常實用的換算關(guān)系必須掌握偽距組合P1減P2的1米差異大約對應(yīng)9.52 TECU具體推導(dǎo)依據(jù)是40.3乘以(1/f1的平方減1/f2的平方)的倒數(shù)f1是1575.42兆赫茲f2是1227.60兆赫茲。這個系數(shù)在驗證程序輸出時特別有用算完TEC以后可以反算一下P4殘差看看是否在合理范圍內(nèi)。2. 源碼核心邏輯STEC是怎么從觀測值里算出來的RNX2GTEX本質(zhì)上是把教科書里的雙頻電離層探測公式翻譯成了Fortran代碼。所以讀這份源碼數(shù)學(xué)上并不難難的是理解代碼里各種變量、常量和文件操作背后對應(yīng)的物理過程。我個人推薦在讀源碼之前先把下面幾個關(guān)鍵邏輯在紙上推導(dǎo)一遍。2.1 無幾何組合的推導(dǎo)和實現(xiàn)要點電離層探測的第一個關(guān)鍵組合叫無幾何組合geometry-free combination也叫電離層殘差組合。對于偽距組合形式是P4等于P1減P2對于載波相位組合形式是L4等于L1減L2。這個組合能消掉衛(wèi)星鐘差、接收機鐘差、對流層延遲、幾何距離等與頻率無關(guān)的項剩下的主要就是電離層延遲差異和硬件延遲偏差。在源碼里這個組合通常不是直接算兩個觀測值相減就完事還要考慮P1和P2的觀測值類型。有的接收機輸出C1和P2有的輸出P1和P2還有的只有C1和C2不同組合對應(yīng)的硬件延遲偏差不一樣代碼里一般會有對應(yīng)的分支判斷。這也是為什么源碼里出現(xiàn)一長串if條件判斷的原因。從組合值換算到STEC的時候符號問題特別容易搞暈。如果程序里寫的是P4等于P1減P2那么P4大約等于負的40.3乘STEC乘以(1除以f1平方減1除以f2平方)換算成TEC需要乘一個負系數(shù)如果程序里用的是P2減P1符號就反過來。我看到不少人在看源碼時卡在這里最后發(fā)現(xiàn)是符號理解反了。2.2 載波相位平滑偽距為什么需要代碼里怎么體現(xiàn)偽距觀測的噪聲比較大尤其C/A碼和P碼噪聲水平通常是分米級甚至米級直接用偽距組合算STEC結(jié)果會非常毛糙。載波相位觀測的噪聲小得多一個量級的差距但相位觀測值存在整周模糊度差分后依然有一個未知常數(shù)偏差無法直接給出絕對TEC。解決辦法是用相位組合的變化量來平滑偽距組合的絕對值這就是經(jīng)典的載波相位平滑偽距算法。具體實現(xiàn)思路是先對L4做周跳檢測把連續(xù)的、沒有周跳的弧段找出來在每個弧段內(nèi)計算(L4減去P4)的時間平均這個平均值包含了模糊度和硬件延遲的綜合常數(shù)然后用P4觀測值加上這個平均值得到平滑后的STEC。程序里通常會有一個循環(huán)逐歷元處理遇到周跳就重新初始化平滑窗口。我在源碼里看到平滑窗口長度設(shè)置時不同版本差別比較大。窗口太短平滑效果差窗口太長又容易把電離層的真實變化也抹平了。一般長弧段取20到30分鐘比較穩(wěn)妥但如果電離層活躍、TEC變化劇烈窗口要適當(dāng)縮短。這個參數(shù)值得根據(jù)你的數(shù)據(jù)和研究目的多試幾組。2.3 DCB偏差處理源碼的邊界在哪里這里要特別注意偽距組合P4里面除了電離層延遲還包含衛(wèi)星差分碼偏差和接收機差分碼偏差統(tǒng)稱DCB。即使做了相位平滑DCB依然保留在結(jié)果里。如果不修正STEC會出現(xiàn)一個系統(tǒng)性偏置中緯度地區(qū)這個偏置折算下來少則幾個TECU多則十幾二十個TECU對電離層建模影響很大。RNX2GTEX這個層級的程序通常是不做DCB估計的。它的定位是把原始觀測換算成不含幾何項的電離層組合量并輸出成TEXDCB修正一般留給后續(xù)處理鏈路比如用IGS或者CODE發(fā)布的DCB產(chǎn)品做后處理剔除。源碼里可能預(yù)留了DCB文件的讀取接口也可能完全沒有要看具體版本。你在把TEX數(shù)據(jù)用于定量研究之前一定要確認DCB處理在哪一步完成否則結(jié)果會整體偏移。這條邊界一定要搞清楚。如果你拿到一個RNX2GTEX輸出的TEX文件第一件事不是畫圖看趨勢而是要問這份數(shù)據(jù)是原始STEC還是已經(jīng)做了DCB修正做了衛(wèi)星端修正還是接收機端也修正了。我見過有人直接拿未修正DCB的數(shù)據(jù)做VTEC地圖出來的圖上整個測區(qū)都有一個固定偏置事后排查了半天才發(fā)現(xiàn)問題出在數(shù)據(jù)源頭。3. 把老Fortran代碼跑起來編譯與運行的全過程RNX2GTEX是Fortran寫的一般是Fortran 77風(fēng)格的固定格式代碼。這種老代碼在今天的Linux環(huán)境上編譯通常會遇到一些小問題但解決起來也不難關(guān)鍵是要知道幾個典型的坑。3.1 源碼文件組成與gfortran編譯命令從源碼庫拿到的RNX2GTEX可能是一個單獨的.f文件也可能是主程序加若干子程序文件的結(jié)構(gòu)。以常見版本為例主程序文件名一般是RNX2GTEX.f里面包含若干個subroutine比如讀取RINEX文件頭的子程序、讀取觀測記錄的子程序、計算TEC的子程序、寫出TEX文件的子程序每個子程序?qū)?yīng)一個獨立的處理階段。用gfortran編譯時最簡單的命令是gfortran -O2 -ffixed-line-length-132 -o rnx2gtex RNX2GTEX.f如果源碼拆成了多個文件就全部列在命令后面gfortran -O2 -ffixed-line-length-132 -o rnx2gtex RNX2GTEX.f SUBRTN.f CONST.f-ffixed-line-length-132這個選項容易忽略但經(jīng)常是編譯報錯的關(guān)鍵。很多老Fortran代碼的注釋和續(xù)行標(biāo)志在固定格式下默認只認72列超過的部分會被忽略如果源碼里某一行比較長不調(diào)整行長限制的話編譯時會報一堆莫名其妙的語法錯誤。如果是64位Linux系統(tǒng)一般不用加額外選項就能編譯但個別版本會用到非標(biāo)準的庫函數(shù)這時要根據(jù)報錯信息去源碼里查具體是哪個函數(shù)再決定是替換實現(xiàn)還是增加兼容代碼。我不建議一開始就改動源碼邏輯先試著原樣編譯遇到問題再對癥下藥。3.2 RINEX文件命名約定這個坑最容易栽RNX2GTEX對輸入文件的命名有約定基本上遵循標(biāo)準RINEX文件命名規(guī)則前四個字符是站名縮寫第五到第七個字符是年積日第八個字符是日內(nèi)文件序號第九第十個字符是年份點后面兩位是文件類型標(biāo)識。比如abmf0010.15o表示abmf站、年積日第001天、序號0、2015年、觀測文件。這個命名約定是程序正確運行的前提因為很多老程序不會做太智能的文件名解析而是直接按位置截取字符串來提取站名、年份和年積日。如果你把文件重命名成test_obs.15o之類的名字程序要么報錯要么給出完全錯誤的輸出。我建議在運行前把所有輸入文件統(tǒng)一改成標(biāo)準命名并放在同一個目錄下文件名全部用小寫或者全部用大寫不要混用。有的程序在文件系統(tǒng)大小寫敏感的環(huán)境里對文件名的判斷很嚴格稍微不一致就會找不到文件。3.3 運行、交互輸入與TEX輸出驗證編譯成功后運行方式通常有兩種取決于你拿到的版本。老版本一般是交互式提示輸入運行后程序會問你要RINEX文件名有的版本支持命令行參數(shù)直接指定文件名比如./rnx2gtex abmf0010.15o運行結(jié)束后目錄下會多出一個TEX文件。這時不要急著拿去用先打開文件看一下頭部信息對不對站名是否與輸入一致歷元數(shù)是否合理有沒有出現(xiàn)大量零值或負值。我一般習(xí)慣用head命令看前幾十行再用awk統(tǒng)計一下STEC列的數(shù)值范圍如果最小值是負幾十、最大值是正幾百說明數(shù)據(jù)預(yù)處理環(huán)節(jié)大概率有問題。如果輸出文件是空的或者程序中途崩潰最優(yōu)先檢查的永遠是RINEX文件本身可以先確認它能否被其他常用軟件正常讀取比如用teqc或者gfzrnx做一下質(zhì)量檢查排除RINEX文件損壞的可能再回頭查程序參數(shù)設(shè)置。4. 輸出文件長什么樣TEX格式解析與Python后處理TEX格式是RNX2GTEX的輸出也是后續(xù)處理的起點。雖然不同版本輸出格式存在差異但結(jié)構(gòu)上大體一致理解之后用腳本解析并不難。4.1 TEX文件頭部和正文結(jié)構(gòu)一個典型的TEX輸出文件頭部通常包含生成程序的標(biāo)識、站點名稱、站點坐標(biāo)、數(shù)據(jù)的時間范圍、觀測文件的相關(guān)信息等內(nèi)容。正文部分按歷元組織每個歷元下列出可見衛(wèi)星的TEC值可能帶有衛(wèi)星編號、仰角等信息。對于源碼里的寫語句建議逐個對照看。有的版本輸出STEC有的版本輸出VTEC有的版本會同時輸出仰角供你后續(xù)自己換算。源碼中每個格式描述符對應(yīng)的內(nèi)容最好對照RINEX文件和程序內(nèi)部變量Name來確認不要只看文件后綴就默認是VTEC。4.2 用Python快速解析并畫一條VTEC時間序列TEX文件雖然可以直接用文本編輯器打開看但要做時間序列分析或畫圖還是得寫腳本。下面給一個很基礎(chǔ)的Python解析示例具體列位置需要根據(jù)你那個版本的寫語句調(diào)整import matplotlib.pyplot as plt records [] with open(abmf0010.tex, r) as f: for line in f: if line.startswith(RNX2GTEX OUTPUT): parts line.split() station parts[2] elif line.strip() and line[0].isdigit(): doy int(line[0:3]) sec float(line[4:14]) prn int(line[15:17]) stec float(line[18:27]) records.append((doy, sec, prn, stec)) for prn in sorted(set(r[2] for r in records)): data [r for r in records if r[2] prn] times [r[1] / 3600.0 for r in data] values [r[3] for r in data] plt.plot(times, values, labelfPRN {prn}) plt.xlabel(Hour of Day) plt.ylabel(STEC (TECU)) plt.legend() plt.show()這段代碼把每個衛(wèi)星的STEC按小時畫成一條線可以快速看出各衛(wèi)星之間的系統(tǒng)偏差是否正常。如果某顆衛(wèi)星整體比別的衛(wèi)星高出一截大概率是衛(wèi)星DCB沒有修正如果出現(xiàn)鋸齒狀跳變說明周跳處理有問題。4.3 數(shù)據(jù)后處理的幾點建議我處理TEX數(shù)據(jù)時習(xí)慣做三步檢查。第一步看單顆衛(wèi)星連續(xù)弧段是否平滑有沒有突跳第二步把所有衛(wèi)星的STEC映射到天頂方向做VTEC按站點看日變化曲線是否合理正常情況中午高、夜間低第三步用同一天相鄰測站的數(shù)據(jù)做交叉驗證如果兩個測站距離很近VTEC應(yīng)該高度一致。映射STEC到VTEC時經(jīng)典做法是除以仰角的正弦值近似映射函數(shù)。但要注意低仰角時映射函數(shù)誤差很大一般會把仰角低于10度或15度的數(shù)據(jù)去掉再換算。RNX2GTEX如果本身不輸出仰角你可能還需要從RINEX文件或者星歷計算里補上這個信息這也是不少人在后續(xù)處理時額外寫模塊的原因。5. 常見問題與排查技巧實錄這個項目我用下來真正運行順利的情況其實不多大多數(shù)時間都在跟各種細節(jié)較勁。下面這些問題幾乎每個使用RNX2GTEX的人都會碰到我按階段整理成一張排查表方便你遇到問題時直接對照。5.1 編譯階段的問題編譯是第一個攔路虎也是最容易勸退新手的環(huán)節(jié)。我見過最多的報錯有兩類一類是固定格式行長問題用-ffixed-line-length-132基本能解決另一類是源碼里用了非標(biāo)準的擴展語法比如Tab開頭的代碼行、超過72列的續(xù)行、或者比較老的Fortran特性gfortran默認模式下不接受。遇到這類問題先把編譯器的警告信息完整看一遍重點找第一個報錯位置因為后續(xù)報錯往往是連鎖反應(yīng)。如果某個語法確實是老擴展最簡單的處理是把那幾行改寫成標(biāo)準Fortran 77語法。我不建議為了編譯通過而關(guān)掉所有警告容易埋下運行時隱患。5.2 運行階段的數(shù)據(jù)問題程序編譯通過只是開始運行結(jié)果異常的排查才更耗時間。以下幾種情況我都在實際數(shù)據(jù)里遇到過第一種輸出全部為0或者全部為負值。最常見的原因是觀測文件里沒有程序期望的觀測值類型比如程序默認讀P1和P2但現(xiàn)代接收機只輸出C1和C2或者雙頻數(shù)據(jù)里有大量L2觀測值缺失。這時需要檢查RINEX文件頭里的觀測類型列表確認實際包含哪些信號。第二種輸出的TEC在長時間段內(nèi)整體偏置。這是DCB沒修正的典型特征尤其是接收機端DCB在同一臺接收機的數(shù)據(jù)里是常數(shù)很容易被誤認為是真實電離層變化。如果相鄰兩天的數(shù)據(jù)在同一時刻都有固定差異大概率就是DCB問題。第三種單顆衛(wèi)星數(shù)據(jù)在某個歷元突然跳變之后又恢復(fù)正常這通常是觀測數(shù)據(jù)本身存在跳變或者周跳漏檢?？梢詸z查平滑窗口的周跳檢測閾值是否需要調(diào)整閾值設(shè)得太松會把小周跳放過去設(shè)得太緊又會頻繁重置平滑窗口導(dǎo)致結(jié)果噪聲增大。5.3 結(jié)果異常排查速查表現(xiàn)象可能原因處理建議輸出全部為零或負數(shù)觀測值類型不匹配、無P1/P2組合查看RINEX頭文件觀測類型確認輸入數(shù)據(jù)STEC整體偏置不同衛(wèi)星各自的基線不同未做衛(wèi)星端/接收機端DCB修正接入DCB產(chǎn)品做后處理修正個別衛(wèi)星弧段突跳周跳漏檢、平滑窗口參數(shù)不合理調(diào)低周跳檢測閾值縮短平滑窗口處理后VTEC夜間出現(xiàn)負值平滑噪聲、低仰角數(shù)據(jù)污染提高截止仰角檢查高度角輸入程序運行時提示文件無法打開文件名不符合RINEX命名規(guī)則按ssssdddf.yyt格式重命名輸入文件編譯時出現(xiàn)unclassifiable statement固定格式行長或非標(biāo)準擴展語法加編譯選項修改老式語法輸出文件只有頭部沒有正文RINEX觀測記錄讀取失敗檢查RINEX文件完整性先用teqc等工具質(zhì)檢表中列出的問題覆蓋了我看到的大多數(shù)求助。實際上每次排查這些問題都能加深對程序和數(shù)據(jù)格式的理解。你把這個表存下來遇到問題先按圖索驥能省不少時間。6. 源碼里值得多讀幾遍的幾個地方RNX2GTEX代碼量不大但里面有幾個片段非常值得反復(fù)讀。第一個是RINEX頭部解析部分這里能看到程序如何處理不同版本的RINEX格式老代碼往往用許多分支來兼容不同年代的文件格式讀這部分能學(xué)到不少處理歷史數(shù)據(jù)的經(jīng)驗。第二個是周跳檢測和平滑窗口的實現(xiàn)。這個模塊直接決定輸出STEC的質(zhì)量也是后來人改動最多的部分。有的版本用L4變化量超過固定閾值來判斷周跳有的版本用相鄰歷元差分的統(tǒng)計量動態(tài)設(shè)定閾值兩種方式各有優(yōu)劣。你完全可以在理解原邏輯后把平滑部分替換成更新的算法比如基于卡爾曼濾波的方式。第三個是TEX文件輸出的寫語句。這部分能直觀看出程序作者對輸出格式的考慮哪些信息被保留哪些信息被丟棄背后都有取舍。如果你要做更細致的分析可能需要在這里增加輸出內(nèi)容比如加上方位角、高度角、信噪比等。我在讀這份源碼過程中最大的體會是老程序雖然界面不友好、代碼風(fēng)格不現(xiàn)代但它的邏輯非常直接幾乎沒有多余的設(shè)計。這種直接性反而讓學(xué)習(xí)變得容易你能清楚地看到每一步算子對應(yīng)教科書里的哪個公式。對于想理解GNSS電離層處理全流程、或者想寫自己的TEC處理工具的人來說RNX2GTEX是一份非常合適的入門源碼。如果后續(xù)想擴展它的能力我建議優(yōu)先考慮這幾個方向一是增加對Galileo和北斗觀測值的支持老程序最初主要是針對GPS設(shè)計的二是把DCB修正模塊直接集成進去省去后處理環(huán)節(jié)三是增加輸出精度因子和高度角信息方便下游做質(zhì)量加權(quán)。改動過程中注意保持原有輸出格式的兼容性因為很多下游工具默認了TEX的舊版結(jié)構(gòu)。最后分享一個實際工作中的小習(xí)慣每次處理一個新的測站或新一天的數(shù)據(jù)時我會保留程序輸出的原始TEX文件用腳本自動生成一份包含最大值、最小值、均值、有效數(shù)據(jù)率的統(tǒng)計報告。這樣處理大量數(shù)據(jù)時能快速挑出異常那天。GNSS電離層數(shù)據(jù)處理麻煩往往不在程序本身而在數(shù)據(jù)質(zhì)量的把控上這個習(xí)慣幫我省了不少排查時間。本文還有配套的精品資源點擊獲取 從源碼解讀RNX2GTEX：GNSS電離層TEC提取與處理實踐
簡介這是一套面向電離層研究者、GNSS數(shù)據(jù)處理工程師及空間天氣分析人員的Fortran科學(xué)計算源碼專注于將標(biāo)準RINEX格式的GNSS觀測數(shù)據(jù)高精度反演為電離層總電子含量TEC并輸出為專用GTEX格式有效支撐定位誤差校正、電離層建模與太陽活動影響評估等關(guān)鍵任務(wù)。資源共42個文件含15個Fortran源文件.f實現(xiàn)核心算法如TEC反演、軌道讀取、周跳修正、15個編譯目標(biāo)文件.o、3個Shell腳本.sh用于自動化流程控制以及Makefile、頭文件、參數(shù)配置列表和詳細README文檔結(jié)構(gòu)完整、模塊清晰便于二次開發(fā)與教學(xué)實踐。壓縮包僅290KB輕量但功能完備。目前已有155人學(xué)習(xí)下載用戶可直接獲取可編譯運行的完整工程涵蓋從RINEX數(shù)據(jù)解析、衛(wèi)星幾何定位、偽距/相位組合TEC計算到GTEX格式寫入的全鏈路實現(xiàn)并附帶Julian日歷轉(zhuǎn)換、IGS周數(shù)計算等實用工具模塊。 從源碼讀懂GNSS電離層處理這事放到現(xiàn)在依然有吸引力。RNX2GTEX這個名字常做GNSS數(shù)據(jù)處理的人一看就明白它做的是把RINEX觀測文件轉(zhuǎn)成TEX格式輸出。這里的TEX不是LaTeX那套排版工具而是電離層TEC交換文件Total Electron Content Exchange的簡寫。整個程序的核心任務(wù)就是從雙頻GNSS觀測值中提取總電子含量TEC再把逐衛(wèi)星、逐歷元的TEC結(jié)果寫成標(biāo)準格式供后續(xù)電離層建模、單頻用戶電離層改正、空間天氣研究使用。我用這個程序處理過不少測站數(shù)據(jù)也把源碼從頭到尾翻過幾遍坦白說這不是一個代碼風(fēng)格很“現(xiàn)代”的項目但它把GNSS電離層測量的整條鏈路用最樸素的方式跑通了。這篇內(nèi)容我就從一個實際使用者的角度把RNX2GTEX涉及的物理原理、源碼結(jié)構(gòu)、編譯運行細節(jié)、輸出格式和常見坑都拆開講一遍適合剛接觸GNSS電離層數(shù)據(jù)處理、或者想通過老Fortran代碼理解觀測方程的同學(xué)參考。1. 先把名詞理清楚RINEX、TEX、TEC各自扮演什么角色很多初學(xué)者看到RNX2GTEX這個軟件名第一反應(yīng)是去查怎么編譯結(jié)果被RINEX文件名、TEX格式、TEC單位這些概念繞暈。我建議先花半小時把下面幾個名詞的關(guān)系理清后面看代碼會輕松很多。1.1 TEC是什么GNSS為什么能測到電離層TEC的全稱是Total Electron Content翻譯過來就是總電子含量指沿信號傳播路徑上單位截面積柱體內(nèi)的自由電子總數(shù)單位是電子數(shù)每平方米。平時我們更常用TECU做單位1 TECU等于10的16次方個電子每平方米。電離層里的自由電子會對GNSS信號產(chǎn)生折射延遲這個延遲的大小和信號頻率的平方成反比。對頻率為f的信號偽距上的電離層延遲近似為40.3乘STEC除以f的平方這里的STEC就是斜路徑上的TEC。因為兩個載波頻率不一樣同一顆衛(wèi)星發(fā)出的L1和L2信號穿過電離層時延遲量不同通過對比兩個頻率的觀測值就能把電離層影響單獨分離出來。我在實際給學(xué)生講的時候喜歡用一個類比電離層就像一塊有色玻璃不同顏色的光穿過時速度不一樣。GNSS接收機同時收到兩個頻率的信號相當(dāng)于拿兩束不同顏色的光同時穿過這塊玻璃通過比較它們的到達時間差就能反推玻璃的“厚度”。這里的“厚度”就是TEC。1.2 從RINEX到TEXRNX2GTEX在流水線中的位置RINEX是GNSS觀測數(shù)據(jù)的標(biāo)準交換格式接收機廠商輸出的原始數(shù)據(jù)經(jīng)過轉(zhuǎn)換后會以RINEX格式保存?zhèn)尉?、載波相位、多普勒等觀測值。RINEX文件里包含的信息非常豐富但直接拿RINEX文件做電離層研究并不方便因為觀測值里混著鐘差、軌道誤差、對流層延遲等一大堆無關(guān)量而且不同接收機輸出的文件格式細節(jié)也不一樣。TEX格式則是專門面向電離層TEC數(shù)據(jù)設(shè)計的交換格式它把處理好的TEC結(jié)果按站點、按歷元、按衛(wèi)星組織起來省去了用戶重復(fù)做預(yù)處理的工作。RNX2GTEX就是連接RINEX和TEX的橋梁輸入一個或多個RINEX觀測文件經(jīng)過質(zhì)量檢查和組合計算輸出TEX格式的TEC時間序列。整個GNSS電離層處理鏈路大致是RINEX觀測文件加精密星歷和鐘差經(jīng)過預(yù)處理、周跳探測、無幾何組合、相位平滑等步驟得到STEC再做硬件延遲校正和映射函數(shù)轉(zhuǎn)換輸出網(wǎng)格化的VTEC或直接輸出STEC序列。RNX2GTEX覆蓋的是從RINEX到STEC/TEC這一核心段落。1.3 單位轉(zhuǎn)換和數(shù)量級怎么判斷算出來的TEC是否合理用這個程序之前最好對TEC的數(shù)量級心里有數(shù)。中緯度地區(qū)平靜電離層情況下天頂方向VTEC通常在10到50 TECU之間低緯赤道異常區(qū)可以到80甚至100 TECU以上高緯和夜間會低很多夜間經(jīng)常只有幾個TECU。如果算出來的結(jié)果在一個測站、一整天內(nèi)變化幾百上千個TECU那基本可以判斷數(shù)據(jù)或處理方法出了問題。從幾何關(guān)系上也要有概念。斜路徑STEC一般是VTEC的好幾倍仰角越低路徑越長STEC越大。仰角30度左右時如果VTEC是30 TECUSTEC大概在50到60 TECU。程序輸出的如果明顯偏離這個范圍就要回頭檢查組合系數(shù)、單位換算是哪里出了問題。還有一個非常實用的換算關(guān)系必須掌握偽距組合P1減P2的1米差異大約對應(yīng)9.52 TECU具體推導(dǎo)依據(jù)是40.3乘以(1/f1的平方減1/f2的平方)的倒數(shù)f1是1575.42兆赫茲f2是1227.60兆赫茲。這個系數(shù)在驗證程序輸出時特別有用算完TEC以后可以反算一下P4殘差看看是否在合理范圍內(nèi)。2. 源碼核心邏輯STEC是怎么從觀測值里算出來的RNX2GTEX本質(zhì)上是把教科書里的雙頻電離層探測公式翻譯成了Fortran代碼。所以讀這份源碼數(shù)學(xué)上并不難難的是理解代碼里各種變量、常量和文件操作背后對應(yīng)的物理過程。我個人推薦在讀源碼之前先把下面幾個關(guān)鍵邏輯在紙上推導(dǎo)一遍。2.1 無幾何組合的推導(dǎo)和實現(xiàn)要點電離層探測的第一個關(guān)鍵組合叫無幾何組合geometry-free combination也叫電離層殘差組合。對于偽距組合形式是P4等于P1減P2對于載波相位組合形式是L4等于L1減L2。這個組合能消掉衛(wèi)星鐘差、接收機鐘差、對流層延遲、幾何距離等與頻率無關(guān)的項剩下的主要就是電離層延遲差異和硬件延遲偏差。在源碼里這個組合通常不是直接算兩個觀測值相減就完事還要考慮P1和P2的觀測值類型。有的接收機輸出C1和P2有的輸出P1和P2還有的只有C1和C2不同組合對應(yīng)的硬件延遲偏差不一樣代碼里一般會有對應(yīng)的分支判斷。這也是為什么源碼里出現(xiàn)一長串if條件判斷的原因。從組合值換算到STEC的時候符號問題特別容易搞暈。如果程序里寫的是P4等于P1減P2那么P4大約等于負的40.3乘STEC乘以(1除以f1平方減1除以f2平方)換算成TEC需要乘一個負系數(shù)如果程序里用的是P2減P1符號就反過來。我看到不少人在看源碼時卡在這里最后發(fā)現(xiàn)是符號理解反了。2.2 載波相位平滑偽距為什么需要代碼里怎么體現(xiàn)偽距觀測的噪聲比較大尤其C/A碼和P碼噪聲水平通常是分米級甚至米級直接用偽距組合算STEC結(jié)果會非常毛糙。載波相位觀測的噪聲小得多一個量級的差距但相位觀測值存在整周模糊度差分后依然有一個未知常數(shù)偏差無法直接給出絕對TEC。解決辦法是用相位組合的變化量來平滑偽距組合的絕對值這就是經(jīng)典的載波相位平滑偽距算法。具體實現(xiàn)思路是先對L4做周跳檢測把連續(xù)的、沒有周跳的弧段找出來在每個弧段內(nèi)計算(L4減去P4)的時間平均這個平均值包含了模糊度和硬件延遲的綜合常數(shù)然后用P4觀測值加上這個平均值得到平滑后的STEC。程序里通常會有一個循環(huán)逐歷元處理遇到周跳就重新初始化平滑窗口。我在源碼里看到平滑窗口長度設(shè)置時不同版本差別比較大。窗口太短平滑效果差窗口太長又容易把電離層的真實變化也抹平了。一般長弧段取20到30分鐘比較穩(wěn)妥但如果電離層活躍、TEC變化劇烈窗口要適當(dāng)縮短。這個參數(shù)值得根據(jù)你的數(shù)據(jù)和研究目的多試幾組。2.3 DCB偏差處理源碼的邊界在哪里這里要特別注意偽距組合P4里面除了電離層延遲還包含衛(wèi)星差分碼偏差和接收機差分碼偏差統(tǒng)稱DCB。即使做了相位平滑DCB依然保留在結(jié)果里。如果不修正STEC會出現(xiàn)一個系統(tǒng)性偏置中緯度地區(qū)這個偏置折算下來少則幾個TECU多則十幾二十個TECU對電離層建模影響很大。RNX2GTEX這個層級的程序通常是不做DCB估計的。它的定位是把原始觀測換算成不含幾何項的電離層組合量并輸出成TEXDCB修正一般留給后續(xù)處理鏈路比如用IGS或者CODE發(fā)布的DCB產(chǎn)品做后處理剔除。源碼里可能預(yù)留了DCB文件的讀取接口也可能完全沒有要看具體版本。你在把TEX數(shù)據(jù)用于定量研究之前一定要確認DCB處理在哪一步完成否則結(jié)果會整體偏移。這條邊界一定要搞清楚。如果你拿到一個RNX2GTEX輸出的TEX文件第一件事不是畫圖看趨勢而是要問這份數(shù)據(jù)是原始STEC還是已經(jīng)做了DCB修正做了衛(wèi)星端修正還是接收機端也修正了。我見過有人直接拿未修正DCB的數(shù)據(jù)做VTEC地圖出來的圖上整個測區(qū)都有一個固定偏置事后排查了半天才發(fā)現(xiàn)問題出在數(shù)據(jù)源頭。3. 把老Fortran代碼跑起來編譯與運行的全過程RNX2GTEX是Fortran寫的一般是Fortran 77風(fēng)格的固定格式代碼。這種老代碼在今天的Linux環(huán)境上編譯通常會遇到一些小問題但解決起來也不難關(guān)鍵是要知道幾個典型的坑。3.1 源碼文件組成與gfortran編譯命令從源碼庫拿到的RNX2GTEX可能是一個單獨的.f文件也可能是主程序加若干子程序文件的結(jié)構(gòu)。以常見版本為例主程序文件名一般是RNX2GTEX.f里面包含若干個subroutine比如讀取RINEX文件頭的子程序、讀取觀測記錄的子程序、計算TEC的子程序、寫出TEX文件的子程序每個子程序?qū)?yīng)一個獨立的處理階段。用gfortran編譯時最簡單的命令是gfortran -O2 -ffixed-line-length-132 -o rnx2gtex RNX2GTEX.f如果源碼拆成了多個文件就全部列在命令后面gfortran -O2 -ffixed-line-length-132 -o rnx2gtex RNX2GTEX.f SUBRTN.f CONST.f-ffixed-line-length-132這個選項容易忽略但經(jīng)常是編譯報錯的關(guān)鍵。很多老Fortran代碼的注釋和續(xù)行標(biāo)志在固定格式下默認只認72列超過的部分會被忽略如果源碼里某一行比較長不調(diào)整行長限制的話編譯時會報一堆莫名其妙的語法錯誤。如果是64位Linux系統(tǒng)一般不用加額外選項就能編譯但個別版本會用到非標(biāo)準的庫函數(shù)這時要根據(jù)報錯信息去源碼里查具體是哪個函數(shù)再決定是替換實現(xiàn)還是增加兼容代碼。我不建議一開始就改動源碼邏輯先試著原樣編譯遇到問題再對癥下藥。3.2 RINEX文件命名約定這個坑最容易栽RNX2GTEX對輸入文件的命名有約定基本上遵循標(biāo)準RINEX文件命名規(guī)則前四個字符是站名縮寫第五到第七個字符是年積日第八個字符是日內(nèi)文件序號第九第十個字符是年份點后面兩位是文件類型標(biāo)識。比如abmf0010.15o表示abmf站、年積日第001天、序號0、2015年、觀測文件。這個命名約定是程序正確運行的前提因為很多老程序不會做太智能的文件名解析而是直接按位置截取字符串來提取站名、年份和年積日。如果你把文件重命名成test_obs.15o之類的名字程序要么報錯要么給出完全錯誤的輸出。我建議在運行前把所有輸入文件統(tǒng)一改成標(biāo)準命名并放在同一個目錄下文件名全部用小寫或者全部用大寫不要混用。有的程序在文件系統(tǒng)大小寫敏感的環(huán)境里對文件名的判斷很嚴格稍微不一致就會找不到文件。3.3 運行、交互輸入與TEX輸出驗證編譯成功后運行方式通常有兩種取決于你拿到的版本。老版本一般是交互式提示輸入運行后程序會問你要RINEX文件名有的版本支持命令行參數(shù)直接指定文件名比如./rnx2gtex abmf0010.15o運行結(jié)束后目錄下會多出一個TEX文件。這時不要急著拿去用先打開文件看一下頭部信息對不對站名是否與輸入一致歷元數(shù)是否合理有沒有出現(xiàn)大量零值或負值。我一般習(xí)慣用head命令看前幾十行再用awk統(tǒng)計一下STEC列的數(shù)值范圍如果最小值是負幾十、最大值是正幾百說明數(shù)據(jù)預(yù)處理環(huán)節(jié)大概率有問題。如果輸出文件是空的或者程序中途崩潰最優(yōu)先檢查的永遠是RINEX文件本身可以先確認它能否被其他常用軟件正常讀取比如用teqc或者gfzrnx做一下質(zhì)量檢查排除RINEX文件損壞的可能再回頭查程序參數(shù)設(shè)置。4. 輸出文件長什么樣TEX格式解析與Python后處理TEX格式是RNX2GTEX的輸出也是后續(xù)處理的起點。雖然不同版本輸出格式存在差異但結(jié)構(gòu)上大體一致理解之后用腳本解析并不難。4.1 TEX文件頭部和正文結(jié)構(gòu)一個典型的TEX輸出文件頭部通常包含生成程序的標(biāo)識、站點名稱、站點坐標(biāo)、數(shù)據(jù)的時間范圍、觀測文件的相關(guān)信息等內(nèi)容。正文部分按歷元組織每個歷元下列出可見衛(wèi)星的TEC值可能帶有衛(wèi)星編號、仰角等信息。對于源碼里的寫語句建議逐個對照看。有的版本輸出STEC有的版本輸出VTEC有的版本會同時輸出仰角供你后續(xù)自己換算。源碼中每個格式描述符對應(yīng)的內(nèi)容最好對照RINEX文件和程序內(nèi)部變量Name來確認不要只看文件后綴就默認是VTEC。4.2 用Python快速解析并畫一條VTEC時間序列TEX文件雖然可以直接用文本編輯器打開看但要做時間序列分析或畫圖還是得寫腳本。下面給一個很基礎(chǔ)的Python解析示例具體列位置需要根據(jù)你那個版本的寫語句調(diào)整import matplotlib.pyplot as plt records [] with open(abmf0010.tex, r) as f: for line in f: if line.startswith(RNX2GTEX OUTPUT): parts line.split() station parts[2] elif line.strip() and line[0].isdigit(): doy int(line[0:3]) sec float(line[4:14]) prn int(line[15:17]) stec float(line[18:27]) records.append((doy, sec, prn, stec)) for prn in sorted(set(r[2] for r in records)): data [r for r in records if r[2] prn] times [r[1] / 3600.0 for r in data] values [r[3] for r in data] plt.plot(times, values, labelfPRN {prn}) plt.xlabel(Hour of Day) plt.ylabel(STEC (TECU)) plt.legend() plt.show()這段代碼把每個衛(wèi)星的STEC按小時畫成一條線可以快速看出各衛(wèi)星之間的系統(tǒng)偏差是否正常。如果某顆衛(wèi)星整體比別的衛(wèi)星高出一截大概率是衛(wèi)星DCB沒有修正如果出現(xiàn)鋸齒狀跳變說明周跳處理有問題。4.3 數(shù)據(jù)后處理的幾點建議我處理TEX數(shù)據(jù)時習(xí)慣做三步檢查。第一步看單顆衛(wèi)星連續(xù)弧段是否平滑有沒有突跳第二步把所有衛(wèi)星的STEC映射到天頂方向做VTEC按站點看日變化曲線是否合理正常情況中午高、夜間低第三步用同一天相鄰測站的數(shù)據(jù)做交叉驗證如果兩個測站距離很近VTEC應(yīng)該高度一致。映射STEC到VTEC時經(jīng)典做法是除以仰角的正弦值近似映射函數(shù)。但要注意低仰角時映射函數(shù)誤差很大一般會把仰角低于10度或15度的數(shù)據(jù)去掉再換算。RNX2GTEX如果本身不輸出仰角你可能還需要從RINEX文件或者星歷計算里補上這個信息這也是不少人在后續(xù)處理時額外寫模塊的原因。5. 常見問題與排查技巧實錄這個項目我用下來真正運行順利的情況其實不多大多數(shù)時間都在跟各種細節(jié)較勁。下面這些問題幾乎每個使用RNX2GTEX的人都會碰到我按階段整理成一張排查表方便你遇到問題時直接對照。5.1 編譯階段的問題編譯是第一個攔路虎也是最容易勸退新手的環(huán)節(jié)。我見過最多的報錯有兩類一類是固定格式行長問題用-ffixed-line-length-132基本能解決另一類是源碼里用了非標(biāo)準的擴展語法比如Tab開頭的代碼行、超過72列的續(xù)行、或者比較老的Fortran特性gfortran默認模式下不接受。遇到這類問題先把編譯器的警告信息完整看一遍重點找第一個報錯位置因為后續(xù)報錯往往是連鎖反應(yīng)。如果某個語法確實是老擴展最簡單的處理是把那幾行改寫成標(biāo)準Fortran 77語法。我不建議為了編譯通過而關(guān)掉所有警告容易埋下運行時隱患。5.2 運行階段的數(shù)據(jù)問題程序編譯通過只是開始運行結(jié)果異常的排查才更耗時間。以下幾種情況我都在實際數(shù)據(jù)里遇到過第一種輸出全部為0或者全部為負值。最常見的原因是觀測文件里沒有程序期望的觀測值類型比如程序默認讀P1和P2但現(xiàn)代接收機只輸出C1和C2或者雙頻數(shù)據(jù)里有大量L2觀測值缺失。這時需要檢查RINEX文件頭里的觀測類型列表確認實際包含哪些信號。第二種輸出的TEC在長時間段內(nèi)整體偏置。這是DCB沒修正的典型特征尤其是接收機端DCB在同一臺接收機的數(shù)據(jù)里是常數(shù)很容易被誤認為是真實電離層變化。如果相鄰兩天的數(shù)據(jù)在同一時刻都有固定差異大概率就是DCB問題。第三種單顆衛(wèi)星數(shù)據(jù)在某個歷元突然跳變之后又恢復(fù)正常這通常是觀測數(shù)據(jù)本身存在跳變或者周跳漏檢?？梢詸z查平滑窗口的周跳檢測閾值是否需要調(diào)整閾值設(shè)得太松會把小周跳放過去設(shè)得太緊又會頻繁重置平滑窗口導(dǎo)致結(jié)果噪聲增大。5.3 結(jié)果異常排查速查表現(xiàn)象可能原因處理建議輸出全部為零或負數(shù)觀測值類型不匹配、無P1/P2組合查看RINEX頭文件觀測類型確認輸入數(shù)據(jù)STEC整體偏置不同衛(wèi)星各自的基線不同未做衛(wèi)星端/接收機端DCB修正接入DCB產(chǎn)品做后處理修正個別衛(wèi)星弧段突跳周跳漏檢、平滑窗口參數(shù)不合理調(diào)低周跳檢測閾值縮短平滑窗口處理后VTEC夜間出現(xiàn)負值平滑噪聲、低仰角數(shù)據(jù)污染提高截止仰角檢查高度角輸入程序運行時提示文件無法打開文件名不符合RINEX命名規(guī)則按ssssdddf.yyt格式重命名輸入文件編譯時出現(xiàn)unclassifiable statement固定格式行長或非標(biāo)準擴展語法加編譯選項修改老式語法輸出文件只有頭部沒有正文RINEX觀測記錄讀取失敗檢查RINEX文件完整性先用teqc等工具質(zhì)檢表中列出的問題覆蓋了我看到的大多數(shù)求助。實際上每次排查這些問題都能加深對程序和數(shù)據(jù)格式的理解。你把這個表存下來遇到問題先按圖索驥能省不少時間。6. 源碼里值得多讀幾遍的幾個地方RNX2GTEX代碼量不大但里面有幾個片段非常值得反復(fù)讀。第一個是RINEX頭部解析部分這里能看到程序如何處理不同版本的RINEX格式老代碼往往用許多分支來兼容不同年代的文件格式讀這部分能學(xué)到不少處理歷史數(shù)據(jù)的經(jīng)驗。第二個是周跳檢測和平滑窗口的實現(xiàn)。這個模塊直接決定輸出STEC的質(zhì)量也是后來人改動最多的部分。有的版本用L4變化量超過固定閾值來判斷周跳有的版本用相鄰歷元差分的統(tǒng)計量動態(tài)設(shè)定閾值兩種方式各有優(yōu)劣。你完全可以在理解原邏輯后把平滑部分替換成更新的算法比如基于卡爾曼濾波的方式。第三個是TEX文件輸出的寫語句。這部分能直觀看出程序作者對輸出格式的考慮哪些信息被保留哪些信息被丟棄背后都有取舍。如果你要做更細致的分析可能需要在這里增加輸出內(nèi)容比如加上方位角、高度角、信噪比等。我在讀這份源碼過程中最大的體會是老程序雖然界面不友好、代碼風(fēng)格不現(xiàn)代但它的邏輯非常直接幾乎沒有多余的設(shè)計。這種直接性反而讓學(xué)習(xí)變得容易你能清楚地看到每一步算子對應(yīng)教科書里的哪個公式。對于想理解GNSS電離層處理全流程、或者想寫自己的TEC處理工具的人來說RNX2GTEX是一份非常合適的入門源碼。如果后續(xù)想擴展它的能力我建議優(yōu)先考慮這幾個方向一是增加對Galileo和北斗觀測值的支持老程序最初主要是針對GPS設(shè)計的二是把DCB修正模塊直接集成進去省去后處理環(huán)節(jié)三是增加輸出精度因子和高度角信息方便下游做質(zhì)量加權(quán)。改動過程中注意保持原有輸出格式的兼容性因為很多下游工具默認了TEX的舊版結(jié)構(gòu)。最后分享一個實際工作中的小習(xí)慣每次處理一個新的測站或新一天的數(shù)據(jù)時我會保留程序輸出的原始TEX文件用腳本自動生成一份包含最大值、最小值、均值、有效數(shù)據(jù)率的統(tǒng)計報告。這樣處理大量數(shù)據(jù)時能快速挑出異常那天。GNSS電離層數(shù)據(jù)處理麻煩往往不在程序本身而在數(shù)據(jù)質(zhì)量的把控上這個習(xí)慣幫我省了不少排查時間。本文還有配套的精品資源點擊獲取