據(jù)解析與顯示：從RadarYHM到PPI/CAPPI的完整流程)
簡介這是一份面向雷達數(shù)據(jù)處理與可視化開發(fā)者的C工程資源包聚焦基數(shù)據(jù)解析與雷達顯示兩大環(huán)節(jié)。壓縮包內(nèi)共80個文件以頭文件和源文件為主涵蓋VC工程配置、可執(zhí)行文件及圖標位圖等資源整體僅7.83MB。該資源已有344人瀏覽學習。其中包含PPI顯示模塊、CRadarBase基類封裝、圖像類型轉換與JPEG處理等核心代碼覆蓋從原始二進制數(shù)據(jù)讀取、解碼到圖像顯示的實現(xiàn)路徑可幫助理解脈沖壓縮、頻率分集等雷達信號處理技術在工程中的實際落地方式。PPI顯示模塊與圖像格式轉換部分對二維雷達成像和GUI更新有直接參考價值。同時完整的MFC工程架構也為讀者提供了雷達顯示界面與數(shù)據(jù)管理模塊的參考范式適合具備一定C與C編程和數(shù)字信號處理基礎、希望快速上手雷達數(shù)據(jù)解析或借鑒完整項目框架的開發(fā)者。1. RadarYHM拿到一個雷達基數(shù)據(jù)壓縮包之后你要做的第一件事做過氣象或者空管雷達數(shù)據(jù)處理的人大概率都有過這種經(jīng)歷甲方丟來一個名為“RadarYHM.rar”的壓縮包里面藏著一坨看起來像亂碼的二進制文件文件名寫著“基數(shù)據(jù)”。你說它是雷達基數(shù)據(jù)沒錯但你不知道它的仰角是怎么排的、方位分辨率是多少、強度場和速度場的存儲順序是什么更別說怎么把它變成一張能看的雷達圖。RadarYHM.rar 代表的就是這個場景雷達基數(shù)據(jù)解析程序包負責把廠家原始格式的基數(shù)據(jù)讀出來、換算成物理量再交給顯示模塊做 PPI 或 CAPPI 渲染。這篇筆記面向的是要用現(xiàn)成程序或自己寫腳本把基數(shù)據(jù)跑通的工作人員目標是讓你從拿到一個只有文件名的基數(shù)據(jù)包到看到第一張完整雷達圖全程不靠猜。2. 基數(shù)據(jù)到底長什么樣先搞懂極坐標存儲再聊解析很多人第一次打開基數(shù)據(jù)文件就懵了因為里面全是十六進制字節(jié)根本看不到“距離”“方位”“反射率”這些字眼。其實雷達基數(shù)據(jù)的存儲思路幾十年沒變過雷達繞垂直軸一圈一圈地掃描每掃一圈沿一條徑向線按距離庫range bin采樣把回波強度記錄下來。所以整包數(shù)據(jù)本質上是一個極坐標三維數(shù)組——方位角、距離庫、仰角三個維度定一個采樣點。2.1 徑向數(shù)據(jù)塊結構方位頭與距離庫的讀取方式不同廠家的基數(shù)據(jù)格式差異很大但只要拿到一份格式文檔解析邏輯幾乎都是同一個套路文件先是一個文件頭描述雷達站點經(jīng)緯度、雷達型號、體掃模式接著就是很多個“徑向數(shù)據(jù)塊”每一塊以一個徑向頭開頭徑向頭里寫著本徑仰角、方位角、徑向序號、距離庫數(shù)量、庫長等元信息徑向頭后面跟著三個等長的數(shù)據(jù)段分別對應反射率、速度、譜寬。常見做法是在解析前先讀一段診斷信息確認為大端還是小端字節(jié)序。國內(nèi)氣象雷達常見的基數(shù)據(jù)文件是 SA/SB 波段標準格式還有不少廠家在標準格式基礎上做了“私有擴展”所以程序里通常要留一個“格式探測開關”。用 Python 寫最直接的讀取邏輯時我一般先用struct.unpack按文檔把徑向頭拆出來。import struct # 假設已按文檔確認徑向頭 144 字節(jié)強度段在前速度段在中 def read_radial_head(buf, offset): # 這里以某常見格式為例各廠家的字段偏移以實際文檔為準 fields struct.unpack_from(HHI, buf, offset) # 仰角*64、方位角*64、徑向號 elev_raw, azim_raw, radial_no [x for x in fields] # 第 3 個位置是庫長單位米第 4 個位置是距離庫數(shù) range_len, bin_count struct.unpack_from(HH, buf, offset 6) return { elevation: elev_raw / 64.0, azimuth: azim_raw / 64.0, range_len_m: range_len, bin_count: bin_count } # 調用時開一個二進制文件按步長跳著讀 with open(radar_base_data.bin, rb) as f: raw f.read() offset 0 radials [] while offset len(raw): head read_radial_head(raw, offset) # 徑向頭長度不同數(shù)據(jù)段長度 庫數(shù) * 2short 存儲 data_len head[bin_count] * 2 offset 144 # 跳到強度段 # 強度、速度、譜寬依次排列這里一次性把三段讀走 ref struct.unpack_from(%dH % head[bin_count], raw, offset) vel struct.unpack_from(%dH % head[bin_count], raw, offset data_len) wid struct.unpack_from(%dH % head[bin_count], raw, offset data_len * 2) offset 144 data_len * 3先說明一個關鍵點elevation和azimuth經(jīng)常不是直接存浮點數(shù)而是存一個整數(shù)讀出來除以一個比例因子上面的 64 就是最常見的檔位有些雷達用 256。如果你讀了數(shù)據(jù)之后發(fā)現(xiàn)角度全部偏大或偏小幾十倍先懷疑這個比例因子而不是懷疑字節(jié)序。另一個常見情況是私有格式里每個距離庫不是一個 short 而是兩個 short前一個存原始計數(shù)后一個存質量標記這種情況可以把上面代碼的144改成64或其他實際頭長再把data_len乘以 2、4 或者 6。2.2 字節(jié)序、仰角掃描順序與雷達常數(shù)的影響字節(jié)序判斷不靠猜。一個很土但有效的辦法是讀一個小文件找到文件頭里寫的站點經(jīng)緯度如果解析出的緯度在 30 到 55 之間、經(jīng)度在 70 到 140 之間那字節(jié)序基本就對了如果讀了能算出緯度 11000 多這種數(shù)字那就是斷反了。注意一下這里說的基數(shù)據(jù)解析核心是“把字節(jié)翻譯成場”而不是“計算雷達氣象方程”。雷達常數(shù)比如發(fā)射功率、天線增益影響的是從原始計數(shù)值換算到 dBZ 的斜率反映在標定參數(shù)里解析程序通常只負責把已標定或者半標定的計數(shù)讀出來真正換算成物理量要再看格式文檔給的公式。仰角的組織方式也需要注意。雷達一次體掃一般包含 9 到 14 個仰角從低到高依次掃描。文件里的徑向數(shù)據(jù)不一定按仰角分組排列有的廠家按時間順序寫也就是先寫第 1 層仰角掃完的幾千個徑向再寫第 2 層也有的廠家先把第 1 層和第 2 層各自掃的中間結果交叉存放最后靠徑向頭里的“仰角編號”來做重排。我一般取出所有徑向頭后做一次排序排序的鍵不是一個 angle 值而是(elevation_index, radial_no)這兩個整數(shù)靠浮點數(shù)來排序容易因為 14.99 和 15.01 這種精度問題把順序打亂。提示不要在拿到新文件時直接把舊配置硬套進去。每個基數(shù)據(jù)的“數(shù)據(jù)段布局”都可能不同尤其是外廠雷達同一型號不同固件版本都能改偏移量。程序里留一個“radar_config”字典里面放頭長、庫長、字節(jié)序、比例因子、是否包含質量標記位換文件先改配置后解析能省掉大量反復打開 hexdump 的時間。3. 基數(shù)據(jù)解析與物理量換算把字節(jié)變成反射率和速度場解析字節(jié)只是第一步要把“原始計數(shù)”變成“雷達產(chǎn)品”還得做兩步換算一是把計數(shù)值轉換成物理量反射率因子、徑向速度、譜寬二是把極坐標網(wǎng)格貼到經(jīng)緯度或平面投影上。這兩個環(huán)節(jié)做不好出來的圖就是一塊花斑雷達顯示基本沒法用。3.1 反射率、速度與譜寬的換算公式反射率的下限和精度一般格式文檔里都會給一組閾值比如“小于 2 的計數(shù)值視為無效回波”“大于等于 10 的計數(shù)值剔除旁瓣干擾”。常見做法是程序啟動時讀一個標定文件里面包含cal_ref (count - threshold) * 0.5這種線性映射的斜率和偏置或者更常見的查表方式直接將計數(shù)值作為索引查一張 256 長度的查詢表表中存儲對應的 dBZ 值。查詢表的好處是雷達廠家可以針對硬件非線性單獨修正不至于動解析代碼。速度場換算更依賴硬件參數(shù)。多普勒速度的換算一般是velocity (count - velocity_nyquist_offset) * velocity_resolution但要注意速度的模糊Nyquist 速度問題如果解析出的速度場在某一徑向范圍內(nèi)出現(xiàn)整條從正到負的突變帶那基本就是速度模糊需要在顯示端做退模糊而數(shù)據(jù)解析端要保持原始值不被改動。# 反射率和速度換算的參考實現(xiàn) def convert_reflectivity(counts, lookup_table): # 256 長查詢表索引 0..255值單位 dBZ import numpy as np arr np.asarray(counts, dtypenp.uint16) # 無效閾值由文檔給出一般是 0 或 1 或 2 invalid_mask arr 255 dBZ lookup_table[np.clip(arr, 0, 255)] dBZ[invalid_mask] np.nan # 無效回波設為 NaN顯示時渲染色板跳過 return dBZ def convert_velocity(counts, invalid_threshold, nyquist_velocity): import numpy as np counts np.asarray(counts, dtypenp.float32) counts[counts invalid_threshold] np.nan # 典型公式速度 (count - 閾值) * 速度分辨率 - 奈奎斯特速度 vel (counts - invalid_threshold) * (2.0 * nyquist_velocity / 255.0) - nyquist_velocity return vel這段代碼里的lookup_table通常不是解析程序寫死的而是由雷達廠家提供的標定文件生成。搜 RadarYHM 相關資源時你會發(fā)現(xiàn)不同版本的基數(shù)據(jù)解析程序帶的標定文件有時是不一樣的老版本用ref.tbl固定文件新版本改成calibration.xml里面除了反射率還寫了仰角補償和方位增益修正系數(shù)。換標定文件后反射率整體偏移 1 到 2 個 dB 是常有的事所以雷達顯示程序在左下角會顯示標定文件名方便現(xiàn)場比對。3.2 體掃重組按仰角生成掃描層再做坐標變化徑向數(shù)據(jù)處理完之后下一步是“體掃重組”。你需要把整包徑向按仰角分組每個仰角組生成一個極坐標矩陣矩陣的形狀是(方位徑向數(shù), 距離庫數(shù))。方位徑向數(shù)取決于雷達掃描模式如果雷達按 0.5 度方位步進掃描一圈那一圈就是 720 條徑向如果是連續(xù)掃描并回傳角度頭徑向數(shù)可能不是整數(shù)最后一條徑向和第一條徑向的夾角覆蓋不到 360 度這時一般要做方位插值。體掃重組的代碼如下import numpy as np # radial_list: 每個元素是 {elevation: float, azimuth: float, ref: np.ndarray} def build_elevation_layer(radials, num_azimuths720, max_range_bins460): # 極坐標矩陣第 0 維度是方位第 1 維度是距離庫 layer np.full((num_azimuths, max_range_bins), np.nan, dtypenp.float32) for r in radials: azimuth_idx int(round(r[azimuth] / 0.5)) % num_azimuths # 防止某一徑向的庫數(shù)比矩陣寬 bins r[ref][:max_range_bins] layer[azimuth_idx, :len(bins)] bins # 對沒有徑向覆蓋的方位角做最近鄰填補 for az in range(num_azimuths): if np.isnan(layer[az]).all(): # 找左右兩邊最近的完整徑向拷貝過來 left az - 1 while left 0 and np.isnan(layer[left]).all(): left - 1 right az 1 while right num_azimuths and np.isnan(layer[right]).all(): right 1 if left 0: layer[az] layer[left] elif right num_azimuths: layer[az] layer[right] return layer這段代碼解決的是“掃描徑向數(shù)不夠一圈”的坑。注意這里int(round(r[azimuth] / 0.5))要求方位角頭是按 0.5 度為單位存儲的如果格式文檔里給的是 0.1 度步進那么 azimuth_idx 的計算要把除數(shù)改成 0.1。填補最近鄰這個操作看著粗暴但在雷達顯示里很實用因為不填補會在扇區(qū)邊界畫出一條黑色的“斷帶”影響對回波連續(xù)性的判斷。想要更平滑的效果可以改成線性插值但要注意方位角是環(huán)形的不能直接線性插值 359 度和 1 度之間的補數(shù)要么拆開角度要么把數(shù)據(jù)復制一份搬到 [0, 720) 再插。4. 從極坐標到屏幕坐標雷達顯示的基礎映射解析出極坐標場之后真正難的反而是顯示這步。雷達顯示不能把極坐標數(shù)據(jù)直接畫到屏幕上你必須為每個距離庫的經(jīng)緯度做坐標變換。這里有兩個選擇一是直接把數(shù)據(jù)網(wǎng)格化成經(jīng)緯度網(wǎng)格再做標準的等經(jīng)緯度渲染二是用 GPU 或圖形庫直接做極坐標扇區(qū)繪制。對從業(yè)人員來說前者實現(xiàn)門檻低更貼近大多數(shù)氣象雷達顯示的尺寸和精度。4.1 雷達站點為中心的極坐標到經(jīng)緯度網(wǎng)格插值把極坐標數(shù)據(jù)網(wǎng)格化到經(jīng)緯度核心是正變換已知雷達站點經(jīng)緯度、某距離庫中心的距離、某個徑向的方位角求該點的經(jīng)緯度。然后建立輸出網(wǎng)格比如 0.01 度分辨率在每個網(wǎng)格點上判斷它落在哪個距離庫和哪個方位角之間取對應值。逆插值比正變換更常用因為這相當于在查詢經(jīng)緯度點時直接算出源的極坐標位置。import numpy as np def polar_to_latlon(radar_lon, radar_lat, azimuth_deg, range_km): # 簡化的正變換假設地球為球體單位度、千米 R 6371.0 delta range_km / R * 180.0 / np.pi az_rad np.deg2rad(azimuth_deg) lat radar_lat delta * np.cos(az_rad) lon radar_lon delta * np.sin(az_rad) / np.cos(np.deg2rad(radar_lat)) return lon, lat # 逆插值從經(jīng)緯度網(wǎng)格點反查極坐標方位與距離 def latlon_to_polar(radar_lon, radar_lat, lon, lat): dx (lon - radar_lon) * np.cos(np.deg2rad(radar_lat)) dy lat - radar_lat az (np.degrees(np.arctan2(dx, dy)) 360.0) % 360.0 dist np.sqrt(dx * dx dy * dy) * 111.195 return az, dist這里的正變換公式?jīng)]有考慮地球曲率和雷達射線的彎曲。說實話對幾百公里內(nèi)的單雷達 PPI 顯示來說不做折射修正不會造成肉眼可辨的偏差但對遠距離超過 200 公里的強度值按地理位置填到圖上會有接近一個網(wǎng)格的偏移。真要做精確的話需要用標準大氣折射模型把雷達射線的高度算出來但這只影響高度不影響方位和水平距離。坐標轉換之后下一步是插值。距離庫方向用“就近取庫”方位方向用“線性插值”基本夠用。我自己一般先給每個經(jīng)緯度網(wǎng)格點算出方位角和距離然后直接算它在極坐標矩陣中的浮點坐標再做一次雙線性插值這比先正變換到經(jīng)緯度網(wǎng)格再重采樣少一個步驟速度也更快。注意橫跨 0 度方位線的處理方位角差要用(az1 - az2 180) % 360 - 180來算不然扇區(qū)邊界會出裂縫這個細節(jié)經(jīng)常讓圖上一圈一圈出現(xiàn)“拉絲”的翻車現(xiàn)場。4.2 顯示色標與底圖疊加的約定雷達顯示程序里反射率和速度的顏色是行業(yè)規(guī)范不是隨便挑的。反射率常用“藍綠黃紅紫”檔位0 dBZ 以下用藍黑漸變0 到 5 dBZ 用淺藍10 到 20 dBZ 用綠色到黃色40 以上直接用紅到紫速度場則用正負顏色區(qū)分朝向和遠離。解析程序只輸出物理量數(shù)組顯示模塊負責把物理量映射到顏色表。這里有一個容易做錯的地方不要讓顯示模塊把 NaN 值映射成“0”而要映射為透明色或背景色否則遠處沒回波的地區(qū)會顯示成一片綠色或紫色這會讓新手誤判為大范圍降水。底圖疊加建議使用標準矢量底圖分塊加載而不是現(xiàn)場去拼在線地圖服務。雷達顯示頁面一般只需要國界、省界、河流和城市點加載完經(jīng)緯度網(wǎng)格插值數(shù)據(jù)后直接用圖形庫按經(jīng)緯度疊加就行。最多見的定位錯誤是雷達站點經(jīng)緯度和底圖坐標系不一致基數(shù)據(jù)文件里如果是 CGCS2000 坐標系的經(jīng)緯度而底圖用的是 WGS84兩者在局部地區(qū)會差上百米雷達回波會被平移。解析層面要多一步坐標基準統(tǒng)一別把雷達中心點直接當作底圖原點來平移。5. 基數(shù)據(jù)解析和雷達顯示避坑五條血淚經(jīng)驗做雷達基數(shù)據(jù)解析的從業(yè)者十有八九都踩過下面這些坑。有些問題看著像數(shù)據(jù)壞了其實是解析配置錯了有些問題看似是顯示程序 bug實際上是插值算法的問題。5.1 反射率場全部偏大 10 dBZ標定文件沒加載對現(xiàn)象是雷達顯示的降水范圍看著特別夸張毛毛雨都顯示成中等強度降水強回波區(qū)幾乎全是紫色。你先懷疑標定文件而不是懷疑程序算法。原因很多基數(shù)據(jù)解析程序源碼里有標定文件的默認路徑換機器部署時路徑不對程序悄悄退回默認值導致查表用的還是出廠標定而雷達硬件近期已經(jīng)做過標定維護真實偏置變了。有些老格式里反射率計數(shù)和 dBZ 的換算斜率是 0.5新格式變成 1.0斜率不對直接把數(shù)值翻倍。解決解析程序加一行日志啟動時打印實際加載的標定文件路徑和查詢表第 128 個值對應的 dBZ。如果路徑不對去雷達控制軟件里讀當前標定參數(shù)重生成查詢表再啟動。線上最穩(wěn)妥的做法是每次雷達定標后導出的calibration.xml都執(zhí)行一次備份并且把當前正在用的路徑寫進配置文件的注釋里。5.2 雷達圖出現(xiàn)整圈稀疏的“齒梳”狀紋理徑向坐標映射錯誤現(xiàn)象是強回波區(qū)邊緣一圈一圈的鋸齒尤其是遠距離處一幀圖里同一個回波團的邊緣像被咬了一口每個齒間距完全一致。原因方位徑向數(shù)按 720 生成但實際基數(shù)據(jù)雷達是按 0.1 度步進掃描一圈有 3600 條徑向。做坐標映射時把 3600 條徑向強塞進 720 個方位槽里相當于每 5 條徑向數(shù)據(jù)只保留了 1 條漏掉了 80% 的方位精度。解決先掃描徑向頭統(tǒng)計方位角的最小差值和最大差值程序用統(tǒng)計結果自動決定方位槽數(shù)。不要寫死“720”。如果你拿到的是一個多普勒雷達同時做體掃和扇區(qū)掃描的數(shù)據(jù)還會出現(xiàn)某個仰角只有 120 度扇區(qū)覆蓋這時按 360 度槽位生成全圖扇區(qū)外全是補出來的假數(shù)據(jù)顯示為一條直線狀的假回波邊界。對策是要在坐標映射里加一個 coverage_mask記錄真實徑向覆蓋范圍顯示端只渲染有覆蓋的扇區(qū)其他置為背景。5.3 速度場出現(xiàn)“黑白棋盤”交錯折射修正和速度解模糊沖突現(xiàn)象是速度場圖上出現(xiàn)大面積的 0 值附近正負交錯看著像棋盤格。原因速度場比較敏感極坐標插值時若對速度場用了反射率場那套“最近鄰填補”會把雷達測量不到的遠距離庫填出一個接近 0 的速度值而這些地方實際沒有信號。特別是薄云或層狀云降水的邊緣插值出的 0 速度與真實的速度場混在一起看起來就是隨機正負交替。解決速度場插值時對無覆蓋區(qū)域一律保持 NaN不做最近鄰填充如果有業(yè)務需要風場連續(xù)性只做范圍不超過 2 個距離庫的小窗口平滑不要全圖插值。另一個坑是速度場先做了退模糊再插值退模糊后的場與原始場跳變特性不同插值引入的偽像更明顯。我現(xiàn)在的習慣是原始速度值只做掩碼在極坐標域完成退模糊最后才插值到經(jīng)緯度網(wǎng)格順序不能反。5.4 體掃數(shù)據(jù)各個仰角錯位仰角分組用錯字段現(xiàn)象是 PPI 高仰角圖上的回波和一個低仰角圖上的回波看起來毫無關聯(lián)甚至高低仰角的回波位置完全對不上。原因正常體掃中高層仰角看到的回波會相對低層仰角向雷達中心收縮因為雷達射線隨距離抬升遠距離處低層仰角能看到降水高層仰角已經(jīng)超出降水頂。如果仰角分組字段選錯比如把“俯仰角原始計數(shù)值”當作“仰角索引”使用會把不同掃描層的徑向混在一起顯示端畫出的是幾層數(shù)據(jù)的雜合體位置自然對不上。解決讀取徑向頭之后畫一張仰角索引-徑向序號的關系表比如打印(radial_no, elevation_index)的前 100 條正常情況應該是仰角索引不變、徑向序號遞增掃完一層后仰角索引跳變。如果前 100 條里仰角索引來回跳檢查一下解析的字段偏移量是不是偏了一位。經(jīng)驗是絕大多數(shù)解析錯誤都能通過打印前 10 條徑向頭元數(shù)據(jù)排查出來先看元數(shù)據(jù)再看數(shù)據(jù)。5.5 遠程顯示卡頓經(jīng)緯度重采樣太慢現(xiàn)象是雷達站本機顯示流暢但遠程查看端一加載體掃就卡死瀏覽器內(nèi)存直接拉滿。原因氣象雷達基數(shù)據(jù)一張體掃如果是 14 層、每層 3600 條徑向、每條 800 個庫那每個仰角矩陣是 3600 × 800 個 float14 層接近 1.6 億個數(shù)值。如果遠程顯示端在瀏覽器里逐個采樣點做經(jīng)緯度反算插值即使每層重采樣只算一次也要循環(huán)幾千萬次瀏覽器直接崩。解決常見的做法是服務器端先把體掃數(shù)據(jù)在極坐標域插值成經(jīng)緯度網(wǎng)格比如每層 600×700 網(wǎng)格然后對每個網(wǎng)格編碼成 PNG 瓦片或者 GeoTIFF 切片顯示端只拼圖不重算插值。這樣遠程顯示端只接收圖片內(nèi)存占用和 CPU 消耗都壓到很低。雷達顯示速度優(yōu)化越靠近數(shù)據(jù)源頭做的計算越多越靠近用戶做的計算越少這條原則不會錯。6. 進階用法用解析出的基數(shù)據(jù)做 CAPPI 和剖面圖驗證解析正確性當 PPI 能正常顯示之后你其實還可以用同一套解析結果繼續(xù)做 CAPPI等高平面位置顯示和垂直剖面。這兩個產(chǎn)品對解析正確性的驗證遠遠強于肉眼看 PPI因為 CAPPI 會把各個仰角的數(shù)據(jù)融合到同一個高度面上如果仰角元數(shù)據(jù)解析有誤CAPPI 會出現(xiàn)明顯的“接縫”或“斷崖”。CAPPI 的簡單做法是對每個經(jīng)緯度網(wǎng)格點計算該點在每個仰角上的高度用標準大氣折射近似公式h sqrt(R^2 r^2 2Rr) - R其中 R 是等效地球半徑 8500 公里r 是斜距然后找出距離目標高度最近上下兩個仰角的反射率值做高度線性插值。其實不需要對每個點都算三角函數(shù)可以先把每個仰角的距離庫先轉成高度層做成“高度-庫數(shù)”映射表經(jīng)緯度網(wǎng)格點只做一次徑向距離查找再查表取高度這樣性能能快五倍左右。剖面圖的實現(xiàn)更簡單但也更考驗解析質量選一條穿過雷達站的線段沿線段每個點用 4.1 里的逆變換找到它在每個仰角層中的位置把不同仰角的數(shù)值拼到一個二維剖面上。做剖面圖之前務必先驗證速度場和反射率場是否對齊——有的基數(shù)據(jù)格式里速度場距離庫的起點比反射率場晚一個庫不做對齊直接畫剖面會在近地面出現(xiàn)一條假的水平亮帶。最后想分享一個習慣不管從哪個渠道拿到 RadarYHM 這類基數(shù)據(jù)解析程序第一件事不是跑代碼而是先做一個最小驗證——生成一幀 PPI存為 GeoTIFF 帶坐標信息疊加在衛(wèi)星云圖或者上一幀雷達圖上對比一分鐘。這樣做的目的是確認坐標轉換和仰角分組沒問題而不是只是“看到圖了”。坐標差半度、仰角錯一層這類問題肉眼看實時動畫很難發(fā)現(xiàn)但疊加底圖或者疊加前一時次數(shù)據(jù)一眼就能看穿。希望這篇筆記幫你少走這一段彎路。本文還有配套的精品資源點擊獲取 雷達基數(shù)據(jù)解析與顯示：從RadarYHM到PPI/CAPPI的完整流程
簡介這是一份面向雷達數(shù)據(jù)處理與可視化開發(fā)者的C工程資源包聚焦基數(shù)據(jù)解析與雷達顯示兩大環(huán)節(jié)。壓縮包內(nèi)共80個文件以頭文件和源文件為主涵蓋VC工程配置、可執(zhí)行文件及圖標位圖等資源整體僅7.83MB。該資源已有344人瀏覽學習。其中包含PPI顯示模塊、CRadarBase基類封裝、圖像類型轉換與JPEG處理等核心代碼覆蓋從原始二進制數(shù)據(jù)讀取、解碼到圖像顯示的實現(xiàn)路徑可幫助理解脈沖壓縮、頻率分集等雷達信號處理技術在工程中的實際落地方式。PPI顯示模塊與圖像格式轉換部分對二維雷達成像和GUI更新有直接參考價值。同時完整的MFC工程架構也為讀者提供了雷達顯示界面與數(shù)據(jù)管理模塊的參考范式適合具備一定C與C編程和數(shù)字信號處理基礎、希望快速上手雷達數(shù)據(jù)解析或借鑒完整項目框架的開發(fā)者。1. RadarYHM拿到一個雷達基數(shù)據(jù)壓縮包之后你要做的第一件事做過氣象或者空管雷達數(shù)據(jù)處理的人大概率都有過這種經(jīng)歷甲方丟來一個名為“RadarYHM.rar”的壓縮包里面藏著一坨看起來像亂碼的二進制文件文件名寫著“基數(shù)據(jù)”。你說它是雷達基數(shù)據(jù)沒錯但你不知道它的仰角是怎么排的、方位分辨率是多少、強度場和速度場的存儲順序是什么更別說怎么把它變成一張能看的雷達圖。RadarYHM.rar 代表的就是這個場景雷達基數(shù)據(jù)解析程序包負責把廠家原始格式的基數(shù)據(jù)讀出來、換算成物理量再交給顯示模塊做 PPI 或 CAPPI 渲染。這篇筆記面向的是要用現(xiàn)成程序或自己寫腳本把基數(shù)據(jù)跑通的工作人員目標是讓你從拿到一個只有文件名的基數(shù)據(jù)包到看到第一張完整雷達圖全程不靠猜。2. 基數(shù)據(jù)到底長什么樣先搞懂極坐標存儲再聊解析很多人第一次打開基數(shù)據(jù)文件就懵了因為里面全是十六進制字節(jié)根本看不到“距離”“方位”“反射率”這些字眼。其實雷達基數(shù)據(jù)的存儲思路幾十年沒變過雷達繞垂直軸一圈一圈地掃描每掃一圈沿一條徑向線按距離庫range bin采樣把回波強度記錄下來。所以整包數(shù)據(jù)本質上是一個極坐標三維數(shù)組——方位角、距離庫、仰角三個維度定一個采樣點。2.1 徑向數(shù)據(jù)塊結構方位頭與距離庫的讀取方式不同廠家的基數(shù)據(jù)格式差異很大但只要拿到一份格式文檔解析邏輯幾乎都是同一個套路文件先是一個文件頭描述雷達站點經(jīng)緯度、雷達型號、體掃模式接著就是很多個“徑向數(shù)據(jù)塊”每一塊以一個徑向頭開頭徑向頭里寫著本徑仰角、方位角、徑向序號、距離庫數(shù)量、庫長等元信息徑向頭后面跟著三個等長的數(shù)據(jù)段分別對應反射率、速度、譜寬。常見做法是在解析前先讀一段診斷信息確認為大端還是小端字節(jié)序。國內(nèi)氣象雷達常見的基數(shù)據(jù)文件是 SA/SB 波段標準格式還有不少廠家在標準格式基礎上做了“私有擴展”所以程序里通常要留一個“格式探測開關”。用 Python 寫最直接的讀取邏輯時我一般先用struct.unpack按文檔把徑向頭拆出來。import struct # 假設已按文檔確認徑向頭 144 字節(jié)強度段在前速度段在中 def read_radial_head(buf, offset): # 這里以某常見格式為例各廠家的字段偏移以實際文檔為準 fields struct.unpack_from(HHI, buf, offset) # 仰角*64、方位角*64、徑向號 elev_raw, azim_raw, radial_no [x for x in fields] # 第 3 個位置是庫長單位米第 4 個位置是距離庫數(shù) range_len, bin_count struct.unpack_from(HH, buf, offset 6) return { elevation: elev_raw / 64.0, azimuth: azim_raw / 64.0, range_len_m: range_len, bin_count: bin_count } # 調用時開一個二進制文件按步長跳著讀 with open(radar_base_data.bin, rb) as f: raw f.read() offset 0 radials [] while offset len(raw): head read_radial_head(raw, offset) # 徑向頭長度不同數(shù)據(jù)段長度 庫數(shù) * 2short 存儲 data_len head[bin_count] * 2 offset 144 # 跳到強度段 # 強度、速度、譜寬依次排列這里一次性把三段讀走 ref struct.unpack_from(%dH % head[bin_count], raw, offset) vel struct.unpack_from(%dH % head[bin_count], raw, offset data_len) wid struct.unpack_from(%dH % head[bin_count], raw, offset data_len * 2) offset 144 data_len * 3先說明一個關鍵點elevation和azimuth經(jīng)常不是直接存浮點數(shù)而是存一個整數(shù)讀出來除以一個比例因子上面的 64 就是最常見的檔位有些雷達用 256。如果你讀了數(shù)據(jù)之后發(fā)現(xiàn)角度全部偏大或偏小幾十倍先懷疑這個比例因子而不是懷疑字節(jié)序。另一個常見情況是私有格式里每個距離庫不是一個 short 而是兩個 short前一個存原始計數(shù)后一個存質量標記這種情況可以把上面代碼的144改成64或其他實際頭長再把data_len乘以 2、4 或者 6。2.2 字節(jié)序、仰角掃描順序與雷達常數(shù)的影響字節(jié)序判斷不靠猜。一個很土但有效的辦法是讀一個小文件找到文件頭里寫的站點經(jīng)緯度如果解析出的緯度在 30 到 55 之間、經(jīng)度在 70 到 140 之間那字節(jié)序基本就對了如果讀了能算出緯度 11000 多這種數(shù)字那就是斷反了。注意一下這里說的基數(shù)據(jù)解析核心是“把字節(jié)翻譯成場”而不是“計算雷達氣象方程”。雷達常數(shù)比如發(fā)射功率、天線增益影響的是從原始計數(shù)值換算到 dBZ 的斜率反映在標定參數(shù)里解析程序通常只負責把已標定或者半標定的計數(shù)讀出來真正換算成物理量要再看格式文檔給的公式。仰角的組織方式也需要注意。雷達一次體掃一般包含 9 到 14 個仰角從低到高依次掃描。文件里的徑向數(shù)據(jù)不一定按仰角分組排列有的廠家按時間順序寫也就是先寫第 1 層仰角掃完的幾千個徑向再寫第 2 層也有的廠家先把第 1 層和第 2 層各自掃的中間結果交叉存放最后靠徑向頭里的“仰角編號”來做重排。我一般取出所有徑向頭后做一次排序排序的鍵不是一個 angle 值而是(elevation_index, radial_no)這兩個整數(shù)靠浮點數(shù)來排序容易因為 14.99 和 15.01 這種精度問題把順序打亂。提示不要在拿到新文件時直接把舊配置硬套進去。每個基數(shù)據(jù)的“數(shù)據(jù)段布局”都可能不同尤其是外廠雷達同一型號不同固件版本都能改偏移量。程序里留一個“radar_config”字典里面放頭長、庫長、字節(jié)序、比例因子、是否包含質量標記位換文件先改配置后解析能省掉大量反復打開 hexdump 的時間。3. 基數(shù)據(jù)解析與物理量換算把字節(jié)變成反射率和速度場解析字節(jié)只是第一步要把“原始計數(shù)”變成“雷達產(chǎn)品”還得做兩步換算一是把計數(shù)值轉換成物理量反射率因子、徑向速度、譜寬二是把極坐標網(wǎng)格貼到經(jīng)緯度或平面投影上。這兩個環(huán)節(jié)做不好出來的圖就是一塊花斑雷達顯示基本沒法用。3.1 反射率、速度與譜寬的換算公式反射率的下限和精度一般格式文檔里都會給一組閾值比如“小于 2 的計數(shù)值視為無效回波”“大于等于 10 的計數(shù)值剔除旁瓣干擾”。常見做法是程序啟動時讀一個標定文件里面包含cal_ref (count - threshold) * 0.5這種線性映射的斜率和偏置或者更常見的查表方式直接將計數(shù)值作為索引查一張 256 長度的查詢表表中存儲對應的 dBZ 值。查詢表的好處是雷達廠家可以針對硬件非線性單獨修正不至于動解析代碼。速度場換算更依賴硬件參數(shù)。多普勒速度的換算一般是velocity (count - velocity_nyquist_offset) * velocity_resolution但要注意速度的模糊Nyquist 速度問題如果解析出的速度場在某一徑向范圍內(nèi)出現(xiàn)整條從正到負的突變帶那基本就是速度模糊需要在顯示端做退模糊而數(shù)據(jù)解析端要保持原始值不被改動。# 反射率和速度換算的參考實現(xiàn) def convert_reflectivity(counts, lookup_table): # 256 長查詢表索引 0..255值單位 dBZ import numpy as np arr np.asarray(counts, dtypenp.uint16) # 無效閾值由文檔給出一般是 0 或 1 或 2 invalid_mask arr 255 dBZ lookup_table[np.clip(arr, 0, 255)] dBZ[invalid_mask] np.nan # 無效回波設為 NaN顯示時渲染色板跳過 return dBZ def convert_velocity(counts, invalid_threshold, nyquist_velocity): import numpy as np counts np.asarray(counts, dtypenp.float32) counts[counts invalid_threshold] np.nan # 典型公式速度 (count - 閾值) * 速度分辨率 - 奈奎斯特速度 vel (counts - invalid_threshold) * (2.0 * nyquist_velocity / 255.0) - nyquist_velocity return vel這段代碼里的lookup_table通常不是解析程序寫死的而是由雷達廠家提供的標定文件生成。搜 RadarYHM 相關資源時你會發(fā)現(xiàn)不同版本的基數(shù)據(jù)解析程序帶的標定文件有時是不一樣的老版本用ref.tbl固定文件新版本改成calibration.xml里面除了反射率還寫了仰角補償和方位增益修正系數(shù)。換標定文件后反射率整體偏移 1 到 2 個 dB 是常有的事所以雷達顯示程序在左下角會顯示標定文件名方便現(xiàn)場比對。3.2 體掃重組按仰角生成掃描層再做坐標變化徑向數(shù)據(jù)處理完之后下一步是“體掃重組”。你需要把整包徑向按仰角分組每個仰角組生成一個極坐標矩陣矩陣的形狀是(方位徑向數(shù), 距離庫數(shù))。方位徑向數(shù)取決于雷達掃描模式如果雷達按 0.5 度方位步進掃描一圈那一圈就是 720 條徑向如果是連續(xù)掃描并回傳角度頭徑向數(shù)可能不是整數(shù)最后一條徑向和第一條徑向的夾角覆蓋不到 360 度這時一般要做方位插值。體掃重組的代碼如下import numpy as np # radial_list: 每個元素是 {elevation: float, azimuth: float, ref: np.ndarray} def build_elevation_layer(radials, num_azimuths720, max_range_bins460): # 極坐標矩陣第 0 維度是方位第 1 維度是距離庫 layer np.full((num_azimuths, max_range_bins), np.nan, dtypenp.float32) for r in radials: azimuth_idx int(round(r[azimuth] / 0.5)) % num_azimuths # 防止某一徑向的庫數(shù)比矩陣寬 bins r[ref][:max_range_bins] layer[azimuth_idx, :len(bins)] bins # 對沒有徑向覆蓋的方位角做最近鄰填補 for az in range(num_azimuths): if np.isnan(layer[az]).all(): # 找左右兩邊最近的完整徑向拷貝過來 left az - 1 while left 0 and np.isnan(layer[left]).all(): left - 1 right az 1 while right num_azimuths and np.isnan(layer[right]).all(): right 1 if left 0: layer[az] layer[left] elif right num_azimuths: layer[az] layer[right] return layer這段代碼解決的是“掃描徑向數(shù)不夠一圈”的坑。注意這里int(round(r[azimuth] / 0.5))要求方位角頭是按 0.5 度為單位存儲的如果格式文檔里給的是 0.1 度步進那么 azimuth_idx 的計算要把除數(shù)改成 0.1。填補最近鄰這個操作看著粗暴但在雷達顯示里很實用因為不填補會在扇區(qū)邊界畫出一條黑色的“斷帶”影響對回波連續(xù)性的判斷。想要更平滑的效果可以改成線性插值但要注意方位角是環(huán)形的不能直接線性插值 359 度和 1 度之間的補數(shù)要么拆開角度要么把數(shù)據(jù)復制一份搬到 [0, 720) 再插。4. 從極坐標到屏幕坐標雷達顯示的基礎映射解析出極坐標場之后真正難的反而是顯示這步。雷達顯示不能把極坐標數(shù)據(jù)直接畫到屏幕上你必須為每個距離庫的經(jīng)緯度做坐標變換。這里有兩個選擇一是直接把數(shù)據(jù)網(wǎng)格化成經(jīng)緯度網(wǎng)格再做標準的等經(jīng)緯度渲染二是用 GPU 或圖形庫直接做極坐標扇區(qū)繪制。對從業(yè)人員來說前者實現(xiàn)門檻低更貼近大多數(shù)氣象雷達顯示的尺寸和精度。4.1 雷達站點為中心的極坐標到經(jīng)緯度網(wǎng)格插值把極坐標數(shù)據(jù)網(wǎng)格化到經(jīng)緯度核心是正變換已知雷達站點經(jīng)緯度、某距離庫中心的距離、某個徑向的方位角求該點的經(jīng)緯度。然后建立輸出網(wǎng)格比如 0.01 度分辨率在每個網(wǎng)格點上判斷它落在哪個距離庫和哪個方位角之間取對應值。逆插值比正變換更常用因為這相當于在查詢經(jīng)緯度點時直接算出源的極坐標位置。import numpy as np def polar_to_latlon(radar_lon, radar_lat, azimuth_deg, range_km): # 簡化的正變換假設地球為球體單位度、千米 R 6371.0 delta range_km / R * 180.0 / np.pi az_rad np.deg2rad(azimuth_deg) lat radar_lat delta * np.cos(az_rad) lon radar_lon delta * np.sin(az_rad) / np.cos(np.deg2rad(radar_lat)) return lon, lat # 逆插值從經(jīng)緯度網(wǎng)格點反查極坐標方位與距離 def latlon_to_polar(radar_lon, radar_lat, lon, lat): dx (lon - radar_lon) * np.cos(np.deg2rad(radar_lat)) dy lat - radar_lat az (np.degrees(np.arctan2(dx, dy)) 360.0) % 360.0 dist np.sqrt(dx * dx dy * dy) * 111.195 return az, dist這里的正變換公式?jīng)]有考慮地球曲率和雷達射線的彎曲。說實話對幾百公里內(nèi)的單雷達 PPI 顯示來說不做折射修正不會造成肉眼可辨的偏差但對遠距離超過 200 公里的強度值按地理位置填到圖上會有接近一個網(wǎng)格的偏移。真要做精確的話需要用標準大氣折射模型把雷達射線的高度算出來但這只影響高度不影響方位和水平距離。坐標轉換之后下一步是插值。距離庫方向用“就近取庫”方位方向用“線性插值”基本夠用。我自己一般先給每個經(jīng)緯度網(wǎng)格點算出方位角和距離然后直接算它在極坐標矩陣中的浮點坐標再做一次雙線性插值這比先正變換到經(jīng)緯度網(wǎng)格再重采樣少一個步驟速度也更快。注意橫跨 0 度方位線的處理方位角差要用(az1 - az2 180) % 360 - 180來算不然扇區(qū)邊界會出裂縫這個細節(jié)經(jīng)常讓圖上一圈一圈出現(xiàn)“拉絲”的翻車現(xiàn)場。4.2 顯示色標與底圖疊加的約定雷達顯示程序里反射率和速度的顏色是行業(yè)規(guī)范不是隨便挑的。反射率常用“藍綠黃紅紫”檔位0 dBZ 以下用藍黑漸變0 到 5 dBZ 用淺藍10 到 20 dBZ 用綠色到黃色40 以上直接用紅到紫速度場則用正負顏色區(qū)分朝向和遠離。解析程序只輸出物理量數(shù)組顯示模塊負責把物理量映射到顏色表。這里有一個容易做錯的地方不要讓顯示模塊把 NaN 值映射成“0”而要映射為透明色或背景色否則遠處沒回波的地區(qū)會顯示成一片綠色或紫色這會讓新手誤判為大范圍降水。底圖疊加建議使用標準矢量底圖分塊加載而不是現(xiàn)場去拼在線地圖服務。雷達顯示頁面一般只需要國界、省界、河流和城市點加載完經(jīng)緯度網(wǎng)格插值數(shù)據(jù)后直接用圖形庫按經(jīng)緯度疊加就行。最多見的定位錯誤是雷達站點經(jīng)緯度和底圖坐標系不一致基數(shù)據(jù)文件里如果是 CGCS2000 坐標系的經(jīng)緯度而底圖用的是 WGS84兩者在局部地區(qū)會差上百米雷達回波會被平移。解析層面要多一步坐標基準統(tǒng)一別把雷達中心點直接當作底圖原點來平移。5. 基數(shù)據(jù)解析和雷達顯示避坑五條血淚經(jīng)驗做雷達基數(shù)據(jù)解析的從業(yè)者十有八九都踩過下面這些坑。有些問題看著像數(shù)據(jù)壞了其實是解析配置錯了有些問題看似是顯示程序 bug實際上是插值算法的問題。5.1 反射率場全部偏大 10 dBZ標定文件沒加載對現(xiàn)象是雷達顯示的降水范圍看著特別夸張毛毛雨都顯示成中等強度降水強回波區(qū)幾乎全是紫色。你先懷疑標定文件而不是懷疑程序算法。原因很多基數(shù)據(jù)解析程序源碼里有標定文件的默認路徑換機器部署時路徑不對程序悄悄退回默認值導致查表用的還是出廠標定而雷達硬件近期已經(jīng)做過標定維護真實偏置變了。有些老格式里反射率計數(shù)和 dBZ 的換算斜率是 0.5新格式變成 1.0斜率不對直接把數(shù)值翻倍。解決解析程序加一行日志啟動時打印實際加載的標定文件路徑和查詢表第 128 個值對應的 dBZ。如果路徑不對去雷達控制軟件里讀當前標定參數(shù)重生成查詢表再啟動。線上最穩(wěn)妥的做法是每次雷達定標后導出的calibration.xml都執(zhí)行一次備份并且把當前正在用的路徑寫進配置文件的注釋里。5.2 雷達圖出現(xiàn)整圈稀疏的“齒梳”狀紋理徑向坐標映射錯誤現(xiàn)象是強回波區(qū)邊緣一圈一圈的鋸齒尤其是遠距離處一幀圖里同一個回波團的邊緣像被咬了一口每個齒間距完全一致。原因方位徑向數(shù)按 720 生成但實際基數(shù)據(jù)雷達是按 0.1 度步進掃描一圈有 3600 條徑向。做坐標映射時把 3600 條徑向強塞進 720 個方位槽里相當于每 5 條徑向數(shù)據(jù)只保留了 1 條漏掉了 80% 的方位精度。解決先掃描徑向頭統(tǒng)計方位角的最小差值和最大差值程序用統(tǒng)計結果自動決定方位槽數(shù)。不要寫死“720”。如果你拿到的是一個多普勒雷達同時做體掃和扇區(qū)掃描的數(shù)據(jù)還會出現(xiàn)某個仰角只有 120 度扇區(qū)覆蓋這時按 360 度槽位生成全圖扇區(qū)外全是補出來的假數(shù)據(jù)顯示為一條直線狀的假回波邊界。對策是要在坐標映射里加一個 coverage_mask記錄真實徑向覆蓋范圍顯示端只渲染有覆蓋的扇區(qū)其他置為背景。5.3 速度場出現(xiàn)“黑白棋盤”交錯折射修正和速度解模糊沖突現(xiàn)象是速度場圖上出現(xiàn)大面積的 0 值附近正負交錯看著像棋盤格。原因速度場比較敏感極坐標插值時若對速度場用了反射率場那套“最近鄰填補”會把雷達測量不到的遠距離庫填出一個接近 0 的速度值而這些地方實際沒有信號。特別是薄云或層狀云降水的邊緣插值出的 0 速度與真實的速度場混在一起看起來就是隨機正負交替。解決速度場插值時對無覆蓋區(qū)域一律保持 NaN不做最近鄰填充如果有業(yè)務需要風場連續(xù)性只做范圍不超過 2 個距離庫的小窗口平滑不要全圖插值。另一個坑是速度場先做了退模糊再插值退模糊后的場與原始場跳變特性不同插值引入的偽像更明顯。我現(xiàn)在的習慣是原始速度值只做掩碼在極坐標域完成退模糊最后才插值到經(jīng)緯度網(wǎng)格順序不能反。5.4 體掃數(shù)據(jù)各個仰角錯位仰角分組用錯字段現(xiàn)象是 PPI 高仰角圖上的回波和一個低仰角圖上的回波看起來毫無關聯(lián)甚至高低仰角的回波位置完全對不上。原因正常體掃中高層仰角看到的回波會相對低層仰角向雷達中心收縮因為雷達射線隨距離抬升遠距離處低層仰角能看到降水高層仰角已經(jīng)超出降水頂。如果仰角分組字段選錯比如把“俯仰角原始計數(shù)值”當作“仰角索引”使用會把不同掃描層的徑向混在一起顯示端畫出的是幾層數(shù)據(jù)的雜合體位置自然對不上。解決讀取徑向頭之后畫一張仰角索引-徑向序號的關系表比如打印(radial_no, elevation_index)的前 100 條正常情況應該是仰角索引不變、徑向序號遞增掃完一層后仰角索引跳變。如果前 100 條里仰角索引來回跳檢查一下解析的字段偏移量是不是偏了一位。經(jīng)驗是絕大多數(shù)解析錯誤都能通過打印前 10 條徑向頭元數(shù)據(jù)排查出來先看元數(shù)據(jù)再看數(shù)據(jù)。5.5 遠程顯示卡頓經(jīng)緯度重采樣太慢現(xiàn)象是雷達站本機顯示流暢但遠程查看端一加載體掃就卡死瀏覽器內(nèi)存直接拉滿。原因氣象雷達基數(shù)據(jù)一張體掃如果是 14 層、每層 3600 條徑向、每條 800 個庫那每個仰角矩陣是 3600 × 800 個 float14 層接近 1.6 億個數(shù)值。如果遠程顯示端在瀏覽器里逐個采樣點做經(jīng)緯度反算插值即使每層重采樣只算一次也要循環(huán)幾千萬次瀏覽器直接崩。解決常見的做法是服務器端先把體掃數(shù)據(jù)在極坐標域插值成經(jīng)緯度網(wǎng)格比如每層 600×700 網(wǎng)格然后對每個網(wǎng)格編碼成 PNG 瓦片或者 GeoTIFF 切片顯示端只拼圖不重算插值。這樣遠程顯示端只接收圖片內(nèi)存占用和 CPU 消耗都壓到很低。雷達顯示速度優(yōu)化越靠近數(shù)據(jù)源頭做的計算越多越靠近用戶做的計算越少這條原則不會錯。6. 進階用法用解析出的基數(shù)據(jù)做 CAPPI 和剖面圖驗證解析正確性當 PPI 能正常顯示之后你其實還可以用同一套解析結果繼續(xù)做 CAPPI等高平面位置顯示和垂直剖面。這兩個產(chǎn)品對解析正確性的驗證遠遠強于肉眼看 PPI因為 CAPPI 會把各個仰角的數(shù)據(jù)融合到同一個高度面上如果仰角元數(shù)據(jù)解析有誤CAPPI 會出現(xiàn)明顯的“接縫”或“斷崖”。CAPPI 的簡單做法是對每個經(jīng)緯度網(wǎng)格點計算該點在每個仰角上的高度用標準大氣折射近似公式h sqrt(R^2 r^2 2Rr) - R其中 R 是等效地球半徑 8500 公里r 是斜距然后找出距離目標高度最近上下兩個仰角的反射率值做高度線性插值。其實不需要對每個點都算三角函數(shù)可以先把每個仰角的距離庫先轉成高度層做成“高度-庫數(shù)”映射表經(jīng)緯度網(wǎng)格點只做一次徑向距離查找再查表取高度這樣性能能快五倍左右。剖面圖的實現(xiàn)更簡單但也更考驗解析質量選一條穿過雷達站的線段沿線段每個點用 4.1 里的逆變換找到它在每個仰角層中的位置把不同仰角的數(shù)值拼到一個二維剖面上。做剖面圖之前務必先驗證速度場和反射率場是否對齊——有的基數(shù)據(jù)格式里速度場距離庫的起點比反射率場晚一個庫不做對齊直接畫剖面會在近地面出現(xiàn)一條假的水平亮帶。最后想分享一個習慣不管從哪個渠道拿到 RadarYHM 這類基數(shù)據(jù)解析程序第一件事不是跑代碼而是先做一個最小驗證——生成一幀 PPI存為 GeoTIFF 帶坐標信息疊加在衛(wèi)星云圖或者上一幀雷達圖上對比一分鐘。這樣做的目的是確認坐標轉換和仰角分組沒問題而不是只是“看到圖了”。坐標差半度、仰角錯一層這類問題肉眼看實時動畫很難發(fā)現(xiàn)但疊加底圖或者疊加前一時次數(shù)據(jù)一眼就能看穿。希望這篇筆記幫你少走這一段彎路。本文還有配套的精品資源點擊獲取