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