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