簡介這是一套面向電離層物理、GNSS信號處理及空間天氣研究領(lǐng)域的Fortran源碼工具專用于將標準RINEX格式的GNSS觀測數(shù)據(jù)如GPS/GLONASS高精度反演為電離層總電子含量TEC并輸出為GTEX專用格式服務(wù)于定位誤差校正、電離層建模與太陽活動影響分析等科研與工程場景。資源共42個文件含15個編譯目標文件.o、15個核心Fortran源碼.f涵蓋RINEX讀取、軌道插值、TEC計算、坐標轉(zhuǎn)換、周跳修正等模塊、3個Shell腳本含自動化處理流程、2個參數(shù)配置列表param.list、filename.list、1個可執(zhí)行程序rnx2gtex及Makefile、README、頭文件Define.h等完整構(gòu)建組件包體僅290KB結(jié)構(gòu)緊湊、依賴精簡。已有155人學習下載讀者可直接編譯運行獲得從原始RINEX觀測到GTEX格式TEC序列的端到端處理能力并通過源碼深入理解TEC反演中的幾何映射、衛(wèi)星鐘差修正、雙頻組合算法及Julian日歷轉(zhuǎn)換等關(guān)鍵技術(shù)實現(xiàn)細節(jié)。1. RNX2GTEX到底是什么一個老工具背后的硬需求做GNSS數(shù)據(jù)處理的人大概率繞不開一個詞TEC。搞電離層研究、做單頻定位誤差修正、評估空間天氣都需要從GNSS觀測數(shù)據(jù)里提取總電子含量。而RNX2GTEX這個Fortran源碼項目干的就是這件事——把RINEX格式的GNSS觀測文件轉(zhuǎn)換成電離層TEC值。先說清楚這個工具在整條數(shù)據(jù)處理鏈路里的位置。GNSS接收機采集到的原始數(shù)據(jù)通常是廠商私有格式比如Trimble的.dat、Leica的.m00經(jīng)過轉(zhuǎn)換后變成標準交換格式RINEX。RINEX文件里記錄的是偽距、載波相位、多普勒等觀測值但這些觀測值本身并不是TEC。TEC需要從L1和L2兩個頻率的觀測值差異中推算出來RNX2GTEX就是在這個環(huán)節(jié)工作的輸入RINEX觀測文件輸出每條信號路徑上的TEC時間序列供后續(xù)電離層建模、層析成像或精度評估使用。這個工具適合誰三類人最需要。第一類是電離層物理方向的研究生手里攢了一堆RINEX數(shù)據(jù)但還沒有系統(tǒng)的TEC提取工具鏈第二類是衛(wèi)星導航定位領(lǐng)域的工程師需要對觀測數(shù)據(jù)做預處理分析比如評估一條基線上的電離層活躍程度第三類是剛接觸Fortran科學計算的開發(fā)者想找一個結(jié)構(gòu)清晰、能直接編譯運行的GNSS處理源碼作為學習范本。無論哪類這個項目都能提供一個完整的、可運行的參考實現(xiàn)。我最初拿到這個源碼時有點意外因為現(xiàn)在Python做GNSS處理的生態(tài)已經(jīng)很成熟比如georinex解析RINEX、numpy做矩陣運算為什么還要用Fortran帶著這個疑問把源碼過了一遍之后我意識到這個選擇有它的道理。Fortran在科學計算領(lǐng)域積累了幾十年的代碼資產(chǎn)尤其是電離層、對流層這類地球物理模型中大量成熟的算法庫都是Fortran寫的RNX2GTEX作為其中一環(huán)用Fortran寫可以非常方便地嵌入到既有模型流程中。另外處理長時間序列的RINEX數(shù)據(jù)時Fortran的數(shù)組操作和I/O效率確實比腳本語言高出一截尤其是在沒有向量化優(yōu)化習慣的Python代碼面前差距很明顯。2. TEC計算原理代碼背后的物理與數(shù)學2.1 雙頻觀測值如何反演TEC要讀懂RNX2GTEX的源碼先得理解TEC反演的基本原理。GNSS信號從衛(wèi)星到接收機的傳播路徑上會穿過電離層而電離層中的自由電子會對無線電信號產(chǎn)生折射效應(yīng)。這個效應(yīng)的核心特征是折射量與信號頻率的平方成反比。所以用兩個不同頻率的信號同時傳播兩者之間的延遲差就直接反映了傳播路徑上的總電子含量。具體到數(shù)學表達偽距觀測方程可以寫成P1 ρ c(dt_r - dt_s) I1 T ε1 P2 ρ c(dt_r - dt_s) I2 T ε2其中電離層延遲I與頻率的關(guān)系是I1 40.3 × TEC / f12 I2 40.3 × TEC / f22兩式相減消去幾何距離ρ、鐘差和對流層延遲T得到P2 - P1 40.3 × TEC × (1/f22 - 1/f12)整理后就是代碼里最核心的公式TEC (P2 - P1) × f12 × f22 / (40.3 × (f12 - f22))對于GPS系統(tǒng)f1 1575.42 MHzf2 1227.60 MHz代入后系數(shù)大約是9.52也就是TEC ≈ 9.52 × (P2 - P1)。這里的TEC單位是TECU1 TECU 101? electrons/m2。RNX2GTEX源碼里必然有一段代碼在做這個換算但實際工程實現(xiàn)遠比公式復雜。因為偽距觀測值噪聲很大單歷元的P2 - P1可能包含幾TECU的噪聲直接使用效果很差。所以常見做法是用載波相位觀測值來平滑偽距或者直接用相位觀測值計算相對TEC變化再用偽距TEC來消除相位模糊度。2.2 載波相位TEC與模糊度處理載波相位觀測值的噪聲遠小于偽距大約只有毫米級對應(yīng)的TEC噪聲可以低到0.01 TECU以下。但相位觀測值有一個整數(shù)模糊度問題導致它只能給出TEC的相對變化無法給出絕對值。相位觀測方程類似L1 ρ c(dt_r - dt_s) - I1 T λ1 × N1 ε1 L2 ρ c(dt_r - dt_s) - I2 T λ2 × N2 ε2同樣做差得到L1λ1 - L2λ2 -(I1 - I2) λ1N1 - λ2N2整理后TEC (L1λ1 - L2λ2 λ2N2 - λ1N1) × f12 × f22 / (40.3 × (f12 - f22))這里的λ1N1 - λ2N2是一個常數(shù)偏差只要沒有周跳它在一段連續(xù)觀測弧段內(nèi)保持不變。所以RNX2GTEX這類工具的標準做法是用偽距TEC確定這個常數(shù)偏差的初值然后用相位TEC去跟蹤高精度的相對變化。這種組合方式兼顧了偽距的絕對性和相位的精密性。源碼里實現(xiàn)這個策略時最需要小心的就是周跳檢測。一旦發(fā)生周跳相位觀測值會跳變常數(shù)偏差也隨之改變?nèi)绻^續(xù)沿用之前的偏差值后續(xù)所有TEC估計都會整體偏移。我看到的實現(xiàn)里通常會使用電離層殘差組合或者M-W組合來檢測周跳碰到周跳就重新初始化偏差。2.3 硬件延遲偏差DCB和斜向TEC轉(zhuǎn)垂直TEC還有一個不可忽視的問題衛(wèi)星和接收機硬件對兩個頻率信號的處理延遲不同這種差分硬件延遲(DCBDifferential Code Bias)會直接疊加在TEC估計值上。如果不做修正TEC結(jié)果可能偏差幾TECU到十幾TECU在電離層平靜時期這個誤差已經(jīng)相當可觀了。處理DCB通常有兩條路線。一條是用外部產(chǎn)品修正比如CODE分析中心發(fā)布的衛(wèi)星和接收機DCB文件直接扣除另一條是在數(shù)據(jù)自身內(nèi)部估計假設(shè)一個區(qū)域內(nèi)所有接收機和衛(wèi)星的DCB之和在一定時間尺度內(nèi)恒定通過最小二乘平差同時解算TEC和DCB。RNX2GTEX這類單站處理工具如果源碼里沒有集成外部DCB文件的接口那它的輸出就是非修正TEC也就是包含DCB偏差的原始TEC值。這一點在用的時候要特別注意否則后續(xù)建模會引入系統(tǒng)性誤差。另外從觀測值中提取的TEC是沿衛(wèi)星到接收機視線方向的斜向TECSTEC而大多數(shù)應(yīng)用場景需要的是垂直TECVTEC。轉(zhuǎn)換需要知道電離層穿刺點IPP的位置和衛(wèi)星高度角常用的是單層映射函數(shù)VTEC STEC × cos(arcsin((Re×sin(90°-el)) / (Reh)))其中Re是地球半徑約為6371 kmh是電離層單層高度通常取350-450 km。RNX2GTEX如果輸出的是STEC那就需要在后處理中自己加這一層映射如果源碼里已經(jīng)包含了映射函數(shù)那它輸出的就是VTEC。我建議拿到源碼后第一件事就是確認這一點看一下主程序里是否調(diào)用了高度角計算和映射函數(shù)相關(guān)的子程序。3. 源碼結(jié)構(gòu)與核心模塊拆解3.1 RINEX文件解析模塊格式細節(jié)決定成敗RINEX格式是GNSS數(shù)據(jù)交換的通用語言但解析它并不輕松。RNX2GTEX的源碼首先要解決的就是讀取RINEX觀測文件通常是RINEX 2.11或2.12版本也可能支持3.x。這個模塊的代碼量往往占整個項目的三分之一以上屬于典型的看著簡單寫起來瑣碎的部分。RINEX觀測文件的結(jié)構(gòu)分兩部分頭部區(qū)Header和數(shù)據(jù)區(qū)Body。頭部區(qū)每一行都有固定的列位置比如第1-60列是描述信息第61-80列是標簽。程序需要從頭部讀取的關(guān)鍵信息包括接收機近似位置APPROX POSITION XYZ、觀測類型列表# / TYPES OF OBSERV、觀測歷元間隔INTERVAL、數(shù)據(jù)類型標志GGPSRGLONASSEGalileoCBDS等。這里有一個特別容易踩的坑不同版本的RINEX對觀測類型的命名規(guī)則不同。RINEX 2.x用的是兩位字符比如L1、L2、P1、P2、C1、D1、S1RINEX 3.x變成了三位字符比如L1C、L2W、C1C、C2W。RNX2GTEX源碼里解析模塊通常會維護一個映射表將不同版本的觀測碼統(tǒng)一映射到內(nèi)部定義的編號然后根據(jù)編號索引對應(yīng)的觀測值數(shù)組。如果你自己改源碼新增了某種觀測類型而沒有同步更新映射表程序會大概率在讀取數(shù)據(jù)時直接報錯。解析數(shù)據(jù)區(qū)時需要注意歷元行的格式。RINEX 2.x的歷元行格式為YY MM DD HH MM SS 標志 衛(wèi)星數(shù)然后是每顆衛(wèi)星的觀測值。觀測值按頭部聲明的類型順序排列缺失值用0.0或空白填充。源碼中常見的做法是逐行讀取后先判斷是否為歷元行通過正則或格式檢查再按衛(wèi)星數(shù)循環(huán)讀取觀測值每一顆衛(wèi)星的觀測值根據(jù)被跳過的空白字段對齊。I/O性能方面RNX2GTEX如果處理的是多天、多站的RINEX數(shù)據(jù)文件I/O的優(yōu)化就很重要。我看到一些實現(xiàn)里用了Fortran的直接訪問direct access或流式I/O來替代默認的順序讀實測在大文件上能快不少。3.2 觀測值預處理與質(zhì)量控制RINEX數(shù)據(jù)不等于干凈數(shù)據(jù)。由于多路徑效應(yīng)、接收機異常、衛(wèi)星故障等原因原始觀測值里混入了不少粗差RNX2GTEX的預處理模塊就是把這些異常值擋在TEC計算之前。預處理通常包含幾個環(huán)節(jié)。首先是衛(wèi)星高度角計算利用廣播星歷或精密星歷計算每顆衛(wèi)星在接收機處的仰角和方位角。高度角低于閾值比如10度的衛(wèi)星觀測值會被剔除因為低高度角信號穿過電離層的路徑長多路徑效應(yīng)嚴重TEC提取精度很難保證。其次是粗差剔除。常用的方法是考察相鄰歷元之間的觀測值變化率如果某個歷元的P2-P1值與前一個歷元的差值超過設(shè)定閾值比如50 TECU就標記為異常并剔除。更嚴格的做法是結(jié)合載波相位平滑偽距后的殘差來判斷殘差超過3倍標準差就剔除。再就是周跳檢測。我上面提到過相位觀測值的連續(xù)性直接決定了TEC跟蹤的精度所以RNX2GTEX里大概率會有一個專門檢測周跳的子程序。經(jīng)典的檢測手段是使用電離層殘差組合LI組合即L1 - L2的線性組合這個組合對周跳非常敏感——一個周跳在LI組合中會產(chǎn)生大約5.4 cm的跳變對GPS遠超噪聲水平。代碼中通常會給LI組合設(shè)置一個閾值比如10 cm超過就判定為周跳并觸發(fā)模糊度重新初始化。3.3 TEC計算與輸出格式預處理完成后TEC計算就是標準的公式推演了。RNX2GTEX主體循環(huán)偽代碼思路通常是打開RINEX觀測文件 讀取頭部信息 循環(huán)讀取每個歷元 循環(huán)讀取每顆衛(wèi)星 檢查觀測值有效性 計算衛(wèi)星高度角/方位角 如果高度角低于閾值跳過 計算偽距TECP2-P1 計算相位TECL1λ1-L2λ2 用偽距TEC初始化相位TEC常數(shù)偏差 檢測周跳若發(fā)生則重新初始化 輸出STEC或VTEC 關(guān)閉文件輸出部分一般有兩種設(shè)計。一種是直接輸出文本文件每行包含時間、衛(wèi)星PRN、高度角、方位角、STEC、VTEC等字段方便用繪圖工具展示另一種是直接輸出成IONEX格式或其他電離層地圖產(chǎn)品需要的中間格式方便與CODE、IGS等機構(gòu)的產(chǎn)品做對比。看RNX2GTEX這個名字GTEX大概率是GNSS TEC EXtraction的意思輸出的應(yīng)該是一種自定義的TEC交換格式后面可以做轉(zhuǎn)換器接到其他可視化工具上。這里我建議拿到代碼后先看輸出樣例確認各個字段的單位和坐標系。比如STEC是正的還是負的相位TEC的符號約定是什么時間是GPS時還是UTC這些細節(jié)如果不核對后面做時間序列對比時很容易出問題。4. 編譯與運行實戰(zhàn)4.1 Fortran環(huán)境準備gfortran安裝與驗證RNX2GTEX是Fortran源碼第一步是準備一個可用的Fortran編譯器。如果你在Linux或macOS上工作gfortran是最穩(wěn)妥的選擇它是GNU編譯器套件的一部分免費、跨平臺、兼容性好。大部分Linux發(fā)行版都可以通過包管理器直接安裝Ubuntu/Debian系用aptCentOS/RHEL系用yum或dnf。macOS上如果裝了Homebrewbrew install gcc就會一并安裝gfortran。Windows上可以安裝MinGW-w64或者使用WSLWindows Subsystem for Linux我個人的建議是直接用WSL跑Linux環(huán)境省去一堆環(huán)境變量和依賴的麻煩。安裝完成后用下面的命令驗證編譯器是否正常gfortran --version如果輸出了版本信息就說明環(huán)境就緒。接著寫一個hello world驗證編譯鏈路program hello print *, Fortran is ready end program hellogfortran -o hello hello.f90 ./hello這里有個小經(jīng)驗老版本的RNX2GTEX源碼可能用的是Fortran 77語法.f后綴而新版本可能是Fortran 90/95.f90后綴。gfortran對兩者都能編譯但處理.f文件時會按F77標準編譯如果你的代碼里用了F90的自由格式寫法記得把文件名改成.f90或者在編譯時加-ffree-form選項強制指定自由格式。4.2 源碼編譯流程從Makefile到手動編譯拿到RNX2GTEX源碼后先看看項目里有沒有Makefile。有Makefile的話直接執(zhí)行make就能出可執(zhí)行文件。但很多學術(shù)代碼的Makefile年代久遠編譯器選項和現(xiàn)代環(huán)境不完全兼容需要手動調(diào)整。我的建議是先把Makefile里的FC變量指到gfortran把FFLAGS里的優(yōu)化選項設(shè)置好FC gfortran FFLAGS -O2 -ffixed-line-length-none -fno-range-check LDFLAGS 這里-ffixed-line-length-none很重要因為老Fortran代碼默認固定格式下每行不能超過72個字符如果源碼里有長表達式不加這個選項會直接編譯報錯。-fno-range-check是應(yīng)對某些老代碼里整數(shù)溢出或類型轉(zhuǎn)換檢查過嚴的問題不過用的時候要謹慎它可能掩蓋真實的代碼bug。如果項目里沒有Makefile也可以手動編譯gfortran -O2 -ffixed-line-length-none -o rnx2gtex main.f90 tec_calc.f90 rinex_read.f90注意編譯順序依賴模塊一定要在主程序之前編譯。如果源碼用到了module.mod文件還需要保證module文件在編譯時能被找到必要時用-I參數(shù)指定module搜索路徑。4.3 運行參數(shù)與輸入輸出文件運行RNX2GTEX前需要準備好輸入文件。最基本的輸入是一份RINEX觀測文件文件名通常類似abcd0010.21o這種格式其中abcd是測站名001是年積日0表示該天內(nèi)第0個會話.21o表示2021年的觀測文件。如果你的源碼還依賴導航文件比如需要計算衛(wèi)星位置或高度角那還需要準備對應(yīng)的.21n或.21p文件。程序運行方式通常是命令行傳參./rnx2gtex abcd0010.21o或者交互式輸入文件名。運行后應(yīng)該會生成一個TEC輸出文件我建議第一次運行先用小數(shù)據(jù)量嘗試比如單站1小時的RINEX文件確認輸出符合預期后再跑完整數(shù)據(jù)。輸出文件的每一行大概長這樣2021 001 12 00 00 G01 45.3 182.6 12.34 11.87 2021 001 12 00 00 G02 38.7 75.2 18.92 17.44分別是時間、衛(wèi)星PRN、高度角、方位角、STEC、VTEC。這些字段名可能因源碼版本不同而有差異但核心信息是一致的。5. 常見問題與排查技巧實錄5.1 編譯期報錯的原因和解決RNX2GTEX這類老代碼編譯報錯幾乎是必然的不用慌大部分問題集中在幾類。最常見的是語法兼容性報錯。比如老代碼用了PAUSE語句這在現(xiàn)代Fortran標準里已經(jīng)被廢棄gfortran會報告Error: PAUSE statement is not standard可以直接注釋掉或者改成CONTINUE。還有OPEN語句里的STATUSNEW如果文件已經(jīng)存在會報錯改成STATUSREPLACE即可。另一個高頻問題是隱式類型沖突。老Fortran代碼默認I-N規(guī)則變量名以I、J、K、L、M、N開頭的是整型其他是實型但有些變量名在現(xiàn)代代碼里容易引起歧義。比如NMAX被當作整型如果代碼里給它賦了浮點值編譯時會報Type mismatch。解決辦法是在主程序和子程序里統(tǒng)一加上IMPLICIT NONE把所有變量顯式聲明雖然改起來麻煩但能一次性把這類問題清干凈。還有一種情況是整數(shù)溢出。處理RINEX數(shù)據(jù)時如果用到整型變量存儲時間戳或數(shù)據(jù)字節(jié)數(shù)在老代碼里可能是16位或32位整型數(shù)據(jù)量大了之后會溢出。gfortran默認整型是32位存儲文件大小沒問題但如果代碼里有從UNIX時間戳轉(zhuǎn)換的邏輯建議把相關(guān)變量改成INTEGER(8)。5.2 運行期數(shù)據(jù)異常的類型和應(yīng)對編譯過了不等于能拿到正確結(jié)果。我實際跑這類工具時經(jīng)常遇到幾類運行期異常。第一類是輸出TEC值全部為零。這通常是輸入的RINEX文件里沒有找到代碼期望的觀測類型。比如源碼默認讀取P1和P2觀測值但你的接收機輸出的是C1和L2源碼里與觀測類型匹配相關(guān)的邏輯就會失效導致所有觀測值被標記為無效。解決方法是查看源碼中解析觀測類型那段把實際存在的觀測碼加入匹配列表。第二類是TEC值出現(xiàn)明顯的系統(tǒng)性偏差或跳變。系統(tǒng)性偏差大概率是DCB未修正導致前面分析過這是預期內(nèi)的。跳變則要重點檢查周跳檢測是否生效可以對照該時段有沒有發(fā)生磁暴等空間天氣事件如果沒有就很有可能是周跳漏檢。第三類是程序運行到某個歷元直接崩潰或死循環(huán)。常見原因是RINEX文件數(shù)據(jù)區(qū)某行格式損壞導致解析邏輯進入錯誤分支。排查時可以打開調(diào)試選項重新編譯或者直接在讀取模塊里加打印語句輸出出問題歷元的時間再回到原始RINEX文件里檢查那一行數(shù)據(jù)。5.3 結(jié)果驗證與質(zhì)量評估拿到TEC結(jié)果后一定要做質(zhì)量檢查而不是直接進入后續(xù)建模。我的習慣做法是三個步驟。第一步畫出TEC時間序列圖看趨勢是否合理。正常情況下白天TEC值高中午前后達到峰值夜間TEC值低日出前達到谷值。一天內(nèi)STEC變化范圍通常在幾TECU到幾十TECU之間。如果曲線完全沒有日變化規(guī)律或者數(shù)值明顯超出合理范圍說明處理鏈路有問題。第二步與外部產(chǎn)品對比。IGS和CODE等機構(gòu)發(fā)布全球電離層地圖GIM可以從它們的網(wǎng)站下載對應(yīng)時刻的TEC數(shù)據(jù)在自己的測站位置插值再與RNX2GTEX的輸出對比。兩者差異在幾個TECU以內(nèi)是正常的如果差異巨大需要仔細檢查是否漏掉了DCB修正。第三步檢查多顆衛(wèi)星之間的一致性。同一測站同一時刻不同衛(wèi)星的STEC因為傳播路徑不同會有差異但VTEC經(jīng)過映射函數(shù)修正后應(yīng)該趨于一致。如果某顆衛(wèi)星的VTEC明顯偏離其他衛(wèi)星大概率是這顆衛(wèi)星數(shù)據(jù)本身有問題或者處理參數(shù)設(shè)置不對。以下整理了一份快速排查表方便復現(xiàn)時對照?，F(xiàn)象可能原因排查方法編譯報錯PAUSE/STOP老語法不適配注釋或替換為CONTINUE編譯報類型不匹配隱式類型沖突加IMPLICIT NONE顯式聲明輸出全零觀測類型不匹配檢查RINEX頭部的觀測類型列表TEC日變化無規(guī)律衛(wèi)星鐘差/軌道未處理確認是否輸入了導航文件TEC整體偏大DCB未修正引入外部DCB改正TEC出現(xiàn)階躍周跳漏檢調(diào)整周跳檢測閾值運行崩潰RINEX數(shù)據(jù)行損壞定位并修復損壞行6. 實際應(yīng)用場景與擴展思路6.1 電離層監(jiān)測與研究中的典型用法RNX2GTEX提取出的TEC時間序列在科研和工程中用途很廣。我自己最近的項目里用這類工具處理了一個中等緯度測站連續(xù)三天的數(shù)據(jù)研究磁暴期間TEC的擾動特征。方法是先提取STEC然后用滑動窗口計算TEC變化率指數(shù)ROTIRate of TEC Index能很清楚地看到磁暴期間電離層閃爍活動增強的信號。具體做法是對每個衛(wèi)星的連續(xù)弧段計算每30秒的TEC變化率再在5分鐘窗口內(nèi)統(tǒng)計標準差就得到了ROTI。用RNX2GTEX的輸出配合幾行Python腳本就能完成import numpy as np import pandas as pd data pd.read_csv(tec_output.txt, delim_whitespaceTrue, names[year,doy,hh,mm,ss,prn,el,az,stec,vtec]) data[time] pd.to_datetime(data[[year,doy,hh,mm,ss]].astype(str).agg(-.join, axis1), format%Y-%j-%H-%M-%S) data data.sort_values([prn,time]) for prn, grp in data.groupby(prn): grp grp.reset_index(dropTrue) dtec grp[stec].diff() / grp[time].diff().dt.total_seconds() * 30 roti dtec.rolling(10, min_periods1).std() data.loc[grp.index, roti] roti這類分析在空間天氣研究和導航增強系統(tǒng)性能評估中非常常見。6.2 從單站到網(wǎng)絡(luò)擴展方向RNX2GTEX的單站版本做出來后下一步很自然的擴展就是處理一個測站網(wǎng)絡(luò)的TEC數(shù)據(jù)做區(qū)域電離層地圖。常見做法是先把每個站的TEC輸出統(tǒng)一為VTEC然后選擇單層高度比如350 km計算每個觀測值的電離層穿刺點經(jīng)緯度再用克里金插值或球諧函數(shù)擬合生成區(qū)域TEC地圖。這個過程在Fortran源碼層面可以加但更便捷的做法是保留RNX2GTEX作為前端提取工具把輸出用Python的scipy和pykrige做插值。兩條路線的區(qū)別在于如果數(shù)據(jù)量巨大、要求實時處理就在Fortran側(cè)把插值也實現(xiàn)了減少I/O開銷如果只是離線研究腳本語言就夠用了沒必要改Fortran源碼。另外一個常見擴展是把TEC輸出與接收機DCB估計結(jié)合起來。如果你有同一個區(qū)域多個測站同時段的數(shù)據(jù)可以在TEC提取結(jié)果的基礎(chǔ)上用最小二乘平差同時估計每個測站的接收機DCB和每顆衛(wèi)星的衛(wèi)星DCB。這個內(nèi)容本身可以寫一篇單獨的博文這里提一句是因為RNX2GTEX輸出的STEC文件是這類平差的直接輸入做好格式對接能省不少事。6.3 源碼改造時的一些個人經(jīng)驗如果你打算在RNX2GTEX基礎(chǔ)上做二次開發(fā)我分享幾個實操中總結(jié)的經(jīng)驗。第一先跑通再改造。拿到源碼第一件事不是急著讀代碼而是先找一小段RINEX數(shù)據(jù)把原始代碼完整編譯運行一遍確認輸出正常。只有基線跑通了你之后做的任何改動才有對照。第二新增輸出字段時注意保持原有輸出格式穩(wěn)定。很多下游處理流程依賴固定列數(shù)的輸出文件如果你加了列下游腳本就會錯位。穩(wěn)妥的做法是新增一個輸出文件或者在文件末尾附加新列而不是改中間列的位置。第三RINEX 3.x的多系統(tǒng)支持值得優(yōu)先考慮?，F(xiàn)在越來越多的接收機輸出RINEX 3.x格式包含GPS、GLONASS、Galileo、BDS等多系統(tǒng)數(shù)據(jù)。如果RNX2GTEX只支持GPS可以照葫蘆畫瓢擴展支持其他系統(tǒng)。注意不同系統(tǒng)的頻率不同比如BDS的B1I是1561.098 MHzB2I是1207.14 MHzGalileo的E1是1575.42 MHzE5a是1176.45 MHz這些頻率常數(shù)必須對應(yīng)修改直接套GPS的公式會得到錯誤結(jié)果。第四單元測試斷不可少。TEC計算函數(shù)是整個工具的核心建議給它寫一組單元測試輸入已知的P1、P2差值驗證輸出TEC是否與理論值一致。比如P2-P1 1米時TEC應(yīng)該約為9.52 TECU幾個簡單用例就能把核心算法釘死后續(xù)改代碼不會改壞。7. 寫在最后的一點實操體會用RNX2GTEX處理TEC數(shù)據(jù)這條技術(shù)路線給我最大的感受是科學計算領(lǐng)域工具不論新舊能解決問題就是好工具。Fortran這門語言被很多人認為過時了但在GNSS電離層處理這種對數(shù)值精度和計算效率要求高的場景它依然有不可替代的價值。RNX2GTEX的價值不僅在于它本身能提取TEC更在于它是一個完整的、結(jié)構(gòu)清晰的樣例新手可以通過讀懂它理解TEC提取的全流程老手可以在它的基礎(chǔ)上做二次開發(fā)。我自己的做法是將它作為整個處理鏈路的源頭后面接Python做可視化和統(tǒng)計分析各取所長。最后再提醒一句跑任何GNSS處理工具之前先確認你對輸入數(shù)據(jù)的格式和數(shù)據(jù)質(zhì)量有充分的了解這是所有后續(xù)分析可信度的基礎(chǔ)。本文還有配套的精品資源點擊獲取