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