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