氣溶膠數(shù)據(jù)處理:從原始回波到可信廓線的關(guān)鍵一躍)
簡(jiǎn)介這份PDF文獻(xiàn)聚焦激光雷達(dá)探測(cè)大氣氣溶膠的數(shù)據(jù)處理研究面向大氣科學(xué)、環(huán)境監(jiān)測(cè)及遙感方向的學(xué)習(xí)者與科研人員幫助理解米散射激光雷達(dá)的系統(tǒng)構(gòu)成與反演算法原理。資源包內(nèi)含1個(gè)PDF文件大小約180KB內(nèi)容源自期刊論文涵蓋激光發(fā)射單元、接收望遠(yuǎn)鏡、Si:APD單光子計(jì)數(shù)器探測(cè)以及回波信號(hào)分析等關(guān)鍵環(huán)節(jié)。文中詳細(xì)介紹了云和氣溶膠消光系數(shù)及衰減后向散射系數(shù)的反演方法包括斜率法、Klett法與Fernald法等常用策略并給出24小時(shí)連續(xù)觀測(cè)的處理結(jié)果便于讀者對(duì)照激光雷達(dá)方程理解重疊因子修正與參數(shù)計(jì)算流程。目前已有208人學(xué)習(xí)適合需要夯實(shí)激光雷達(dá)數(shù)據(jù)處理基礎(chǔ)、撰寫相關(guān)論文或開展大氣探測(cè)實(shí)驗(yàn)的研究者參考借鑒。1. 激光雷達(dá)氣溶膠數(shù)據(jù)處理從原始回波到可信廓線的關(guān)鍵一躍拿到一臺(tái)激光雷達(dá)最讓人頭疼的往往不是硬件調(diào)試而是后面那一長(zhǎng)串原始數(shù)據(jù)。尤其是做大氣氣溶膠測(cè)量時(shí)回波信號(hào)里混著背景光、幾何重疊因子、距離平方衰減還有各種噪聲直接畫出來的曲線根本沒法看。激光雷達(dá)測(cè)量大氣氣溶膠的數(shù)據(jù)處理研究核心就是解決從原始光子計(jì)數(shù)或模擬回波到最終氣溶膠消光系數(shù)、后向散射系數(shù)廓線的全過程。這套流程決定了你能不能從一臺(tái)幾萬塊的激光雷達(dá)里榨出真正有物理意義的數(shù)據(jù)。適合誰看做大氣環(huán)境監(jiān)測(cè)的、搞氣象觀測(cè)的、以及剛接手激光雷達(dá)數(shù)據(jù)的研究生。如果你手里有數(shù)據(jù)但不知道怎么處理或者處理完發(fā)現(xiàn)負(fù)值滿天飛、邊界值亂跳那這篇筆記就是給你寫的。我會(huì)把整個(gè)數(shù)據(jù)處理框架拆開從信號(hào)預(yù)處理到反演算法再到踩過的坑一步步講清楚。2. 原始回波信號(hào)預(yù)處理把噪聲和背景光先摁住激光雷達(dá)原始信號(hào)通常有兩種形式模擬信號(hào)和光子計(jì)數(shù)信號(hào)。不管哪種第一步都是扣除背景噪聲。背景噪聲的來源很多太陽(yáng)光、探測(cè)器暗電流、大氣分子散射都會(huì)貢獻(xiàn)。常見做法是取遠(yuǎn)距離末端的一段信號(hào)做平均作為背景基線。但這里有個(gè)細(xì)節(jié)如果大氣邊界層很高遠(yuǎn)端信號(hào)可能還沒完全衰減到背景水平這時(shí)候直接扣背景就會(huì)把有效信號(hào)也扣掉。我一般會(huì)先畫一遍距離校正信號(hào)看看遠(yuǎn)端是否平坦再?zèng)Q定背景段的位置。2.1 背景扣除與距離平方校正背景扣除之后緊接著就是距離平方校正。激光雷達(dá)方程里信號(hào)強(qiáng)度隨距離平方衰減所以要把接收到的信號(hào)乘以距離的平方才能還原出大氣的后向散射特性。這一步看起來簡(jiǎn)單但距離的起點(diǎn)選錯(cuò)整個(gè)廓線形狀都會(huì)變。距離起點(diǎn)應(yīng)該是激光出射點(diǎn)到接收望遠(yuǎn)鏡光軸的幾何交點(diǎn)而不是簡(jiǎn)單的望遠(yuǎn)鏡位置。很多商用激光雷達(dá)會(huì)在元數(shù)據(jù)里給出這個(gè)值如果沒有就需要用重疊因子反推。import numpy as np def preprocess_lidar_signal(raw_signal, distance, bg_start_idx, bg_end_idx): 激光雷達(dá)原始信號(hào)預(yù)處理 raw_signal: 原始回波信號(hào)數(shù)組 distance: 對(duì)應(yīng)的距離數(shù)組單位米 bg_start_idx: 背景段起始索引 bg_end_idx: 背景段結(jié)束索引 # 計(jì)算背景均值 bg_mean np.mean(raw_signal[bg_start_idx:bg_end_idx]) # 扣除背景 signal_bg_corrected raw_signal - bg_mean # 避免負(fù)值影響后續(xù)對(duì)數(shù)運(yùn)算 signal_bg_corrected np.maximum(signal_bg_corrected, 0) # 距離平方校正 range_corrected signal_bg_corrected * distance ** 2 return range_corrected, bg_mean這段代碼里bg_start_idx和bg_end_idx的選擇直接決定背景扣除的效果。通常我會(huì)選最遠(yuǎn)端的 10% 到 20% 的數(shù)據(jù)點(diǎn)但前提是確認(rèn)這段信號(hào)已經(jīng)衰減到接近零。如果遠(yuǎn)端信號(hào)還有明顯起伏說明背景段選得太近需要往后挪。np.maximum那一步是為了防止扣除背景后出現(xiàn)負(fù)值雖然物理上不應(yīng)該有負(fù)值但實(shí)際數(shù)據(jù)里噪聲會(huì)導(dǎo)致負(fù)值出現(xiàn)直接取對(duì)數(shù)會(huì)報(bào)錯(cuò)。2.2 重疊因子校正與信號(hào)平滑幾何重疊因子是激光雷達(dá)近距離信號(hào)失真的主要原因。發(fā)射光束和接收視場(chǎng)在近距離沒有完全重合導(dǎo)致近場(chǎng)信號(hào)被低估。校正方法有實(shí)驗(yàn)法和理論計(jì)算法。實(shí)驗(yàn)法一般用水平均勻大氣假設(shè)通過比較不同距離的信號(hào)斜率來反推重疊因子。理論計(jì)算則需要知道激光發(fā)散角、望遠(yuǎn)鏡視場(chǎng)角、兩者間距等參數(shù)。我一般先用實(shí)驗(yàn)法快速評(píng)估如果重疊因子在幾百米內(nèi)就接近 1那后續(xù)反演可以忽略這段如果重疊區(qū)域延伸到 1 公里以上就必須做校正。信號(hào)平滑是另一個(gè)容易翻車的地方。平滑窗口太寬會(huì)把氣溶膠層的精細(xì)結(jié)構(gòu)抹掉窗口太窄噪聲又壓不住。常見做法是用滑動(dòng)平均或者小波變換?;瑒?dòng)平均簡(jiǎn)單但會(huì)引入相位偏移小波變換能保留突變特征但參數(shù)不好調(diào)。我的經(jīng)驗(yàn)是先用 5 到 9 點(diǎn)的滑動(dòng)平均試一下如果氣溶膠層邊界變模糊了就換小波。平滑之后一定要檢查信噪比如果平滑后信噪比還是低于 3那這段數(shù)據(jù)基本不可用。3. 氣溶膠消光系數(shù)反演Klett 法和 Fernald 法怎么選預(yù)處理完的信號(hào)下一步就是反演消光系數(shù)。激光雷達(dá)方程里有兩個(gè)未知數(shù)消光系數(shù)和后向散射系數(shù)直接求解是欠定的。所以需要假設(shè)兩者之間的關(guān)系也就是激光雷達(dá)比。Klett 法和 Fernald 法是兩種最常用的反演方法。Klett 法假設(shè)后向散射系數(shù)和消光系數(shù)成冪律關(guān)系適合氣溶膠為主的情況Fernald 法把分子散射和氣溶膠散射分開處理需要知道分子消光系數(shù)適合邊界層以上氣溶膠較少的場(chǎng)景。3.1 Klett 反演法的參數(shù)設(shè)置與邊界值選擇Klett 法的核心公式里邊界值的選擇至關(guān)重要。邊界值通常選在遠(yuǎn)端假設(shè)那里大氣均勻消光系數(shù)已知。如果邊界值選得太遠(yuǎn)信號(hào)噪聲會(huì)放大選得太近又可能把氣溶膠層截?cái)?。我一般?huì)先畫距離校正信號(hào)的對(duì)數(shù)曲線找一段斜率穩(wěn)定的區(qū)域作為邊界。邊界值的大小可以參考大氣能見度或者太陽(yáng)光度計(jì)數(shù)據(jù)如果沒有就用經(jīng)驗(yàn)值比如 1e-5 到 1e-4 每米。def klett_inversion(range_corrected, distance, lidar_ratio, boundary_value, boundary_idx): Klett 法反演氣溶膠消光系數(shù) range_corrected: 距離平方校正后的信號(hào) distance: 距離數(shù)組 lidar_ratio: 激光雷達(dá)比典型值 30-70 sr boundary_value: 邊界處的消光系數(shù) boundary_idx: 邊界點(diǎn)索引 # 取對(duì)數(shù) log_signal np.log(range_corrected) # 計(jì)算積分項(xiàng) integral np.cumsum(range_corrected) * (distance[1] - distance[0]) # 邊界處的積分值 integral_boundary integral[boundary_idx] # 反演消光系數(shù) extinction range_corrected / (range_corrected[boundary_idx] / boundary_value - 2 * lidar_ratio * (integral - integral_boundary)) return extinction這里lidar_ratio的取值直接影響反演結(jié)果。氣溶膠類型不同激光雷達(dá)比差別很大。城市氣溶膠通常在 40 到 60 之間沙塵氣溶膠可以到 50 以上海洋氣溶膠偏低。如果不知道具體類型先用 50 試然后根據(jù)反演出的消光系數(shù)廓線是否合理來調(diào)整。boundary_value和boundary_idx需要配合使用邊界點(diǎn)一般選在信號(hào)信噪比還不錯(cuò)的遠(yuǎn)端比如 3 到 5 公里處。3.2 Fernald 法分離分子與氣溶膠散射Fernald 法的思路是把分子散射和氣溶膠散射分開。分子消光系數(shù)可以通過標(biāo)準(zhǔn)大氣模型計(jì)算比如美國(guó)標(biāo)準(zhǔn)大氣。氣溶膠消光系數(shù)則通過迭代求解。Fernald 法對(duì)邊界值的要求比 Klett 法更敏感因?yàn)榉肿由⑸涞呢暙I(xiàn)在遠(yuǎn)端占比更大。如果邊界值選得不對(duì)反演出的氣溶膠消光系數(shù)會(huì)出現(xiàn)負(fù)值。def fernald_inversion(range_corrected, distance, molecular_extinction, molecular_backscatter, lidar_ratio_aerosol, boundary_value, boundary_idx): Fernald 法反演氣溶膠消光系數(shù) molecular_extinction: 分子消光系數(shù)廓線 molecular_backscatter: 分子后向散射系數(shù)廓線 lidar_ratio_aerosol: 氣溶膠激光雷達(dá)比 # 分子激光雷達(dá)比 lidar_ratio_molecular 8 * np.pi / 3 # 計(jì)算分子后向散射與消光比 ratio_molecular molecular_backscatter / molecular_extinction # 迭代求解 extinction_aerosol np.zeros_like(distance) extinction_aerosol[boundary_idx] boundary_value for i in range(boundary_idx - 1, -1, -1): # 簡(jiǎn)化迭代公式實(shí)際使用時(shí)需要根據(jù)具體文獻(xiàn)調(diào)整 numerator range_corrected[i] * np.exp(-2 * (lidar_ratio_aerosol - lidar_ratio_molecular) * molecular_extinction[i] * (distance[i1] - distance[i])) denominator range_corrected[boundary_idx] / boundary_value 2 * lidar_ratio_aerosol * numerator extinction_aerosol[i] numerator / denominator return extinction_aerosolFernald 法的迭代方向是從邊界點(diǎn)向近端推進(jìn)所以邊界點(diǎn)的選擇決定了整個(gè)廓線的基準(zhǔn)。如果邊界點(diǎn)選在氣溶膠層內(nèi)部反演結(jié)果會(huì)嚴(yán)重失真。我一般會(huì)選在氣溶膠層以上、分子散射為主的區(qū)域比如 5 到 8 公里。molecular_extinction和molecular_backscatter可以用標(biāo)準(zhǔn)大氣模型算也可以用地基微波輻射計(jì)或者探空數(shù)據(jù)。如果沒有實(shí)測(cè)數(shù)據(jù)用標(biāo)準(zhǔn)大氣也能湊合但精度會(huì)打折扣。4. 數(shù)據(jù)處理中的避坑指南那些讓你白干一整天的細(xì)節(jié)做激光雷達(dá)數(shù)據(jù)處理最怕的不是算法復(fù)雜而是細(xì)節(jié)沒注意結(jié)果全錯(cuò)。下面這幾條是我和同行們踩過的坑每條都按現(xiàn)象、原因、解決來寫。4.1 避坑一背景扣除后信號(hào)出現(xiàn)大面積負(fù)值現(xiàn)象扣除背景后遠(yuǎn)端信號(hào)變成負(fù)值距離平方校正后負(fù)值更明顯。原因背景段選得太近把還有效的信號(hào)當(dāng)成了背景?;蛘咛綔y(cè)器飽和導(dǎo)致遠(yuǎn)端信號(hào)被截?cái)啾尘熬邓愠鰜砥摺=鉀Q先畫原始信號(hào)確認(rèn)遠(yuǎn)端是否平坦。如果遠(yuǎn)端有起伏把背景段往后挪。如果探測(cè)器飽和需要換用低增益通道或者加衰減片重新測(cè)量。4.2 避坑二反演出的消光系數(shù)出現(xiàn)負(fù)值現(xiàn)象Klett 或 Fernald 反演后某些距離上的消光系數(shù)是負(fù)數(shù)。原因邊界值選得太小或者激光雷達(dá)比設(shè)得太大。也可能是信號(hào)預(yù)處理時(shí)平滑過度把真實(shí)信號(hào)抹掉了。解決先檢查邊界值用太陽(yáng)光度計(jì)或者能見度數(shù)據(jù)校準(zhǔn)。如果邊界值沒問題把激光雷達(dá)比調(diào)小 10% 到 20% 再試。平滑窗口不要超過 9 點(diǎn)否則會(huì)引入虛假的負(fù)值。4.3 避坑三重疊因子校正后近場(chǎng)信號(hào)反而更差現(xiàn)象做了重疊因子校正近場(chǎng)信號(hào)反而出現(xiàn)異常峰值或者凹陷。原因重疊因子曲線是用水平大氣假設(shè)反推的如果實(shí)際大氣不均勻反推的重疊因子就不準(zhǔn)?;蛘咝U龝r(shí)距離起點(diǎn)沒對(duì)齊導(dǎo)致校正曲線偏移。解決先用水平均勻天氣的數(shù)據(jù)做重疊因子比如清晨或者陰天。校正時(shí)確保距離起點(diǎn)和重疊因子曲線的起點(diǎn)一致。如果近場(chǎng)信號(hào)還是不對(duì)干脆把 500 米以內(nèi)的數(shù)據(jù)標(biāo)記為不可用不要強(qiáng)行校正。4.4 避坑四不同時(shí)間的數(shù)據(jù)拼在一起出現(xiàn)斷層現(xiàn)象把不同時(shí)間測(cè)的廓線拼成時(shí)間序列發(fā)現(xiàn)相鄰時(shí)刻消光系數(shù)跳變。原因每次測(cè)量的背景噪聲不一樣或者激光能量有波動(dòng)。也可能是反演時(shí)邊界值每次都在變。解決每次測(cè)量都單獨(dú)算背景不要用固定值。激光能量波動(dòng)可以用監(jiān)測(cè)通道歸一化。邊界值盡量固定如果必須變記錄變化原因。拼時(shí)間序列前先做一致性檢查把明顯異常的時(shí)刻剔除。4.5 避坑五信噪比低的數(shù)據(jù)強(qiáng)行反演現(xiàn)象反演出的廓線噪聲極大看不出任何氣溶膠層結(jié)構(gòu)。原因原始信號(hào)信噪比太低可能是天氣不好、激光能量下降或者探測(cè)器老化。解決先算信噪比如果低于 3這段數(shù)據(jù)直接放棄。不要試圖用強(qiáng)平滑來救平滑只會(huì)把噪聲變成虛假的結(jié)構(gòu)。如果必須用就做時(shí)間平均比如 10 分鐘平均但要注意氣溶膠層的變化尺度。5. 從廓線到應(yīng)用氣溶膠邊界層高度提取與驗(yàn)證技巧處理完消光系數(shù)廓線下一步通常是提取氣溶膠邊界層高度。這是激光雷達(dá)數(shù)據(jù)最常用的應(yīng)用之一。方法有很多梯度法、小波變換法、曲線擬合法。梯度法最簡(jiǎn)單找消光系數(shù)梯度最大的點(diǎn)。但梯度法對(duì)噪聲敏感容易把氣溶膠層內(nèi)部的波動(dòng)當(dāng)成邊界層頂。小波變換法抗噪能力強(qiáng)但尺度參數(shù)不好選。我一般先用梯度法快速看一下再用小波變換法驗(yàn)證。import pywt def boundary_layer_height(extinction, distance, waveletdb4, scale10): 用小波變換提取氣溶膠邊界層高度 extinction: 消光系數(shù)廓線 distance: 距離數(shù)組 wavelet: 小波基 scale: 尺度參數(shù) # 對(duì)消光系數(shù)做小波變換 coeffs pywt.cwt(extinction, scales[scale], waveletwavelet) # 取模極大值對(duì)應(yīng)的距離 modulus np.abs(coeffs[0]) blh_idx np.argmax(modulus) return distance[blh_idx]這里scale參數(shù)決定了小波變換的尺度。尺度太小會(huì)把氣溶膠層內(nèi)部的細(xì)節(jié)當(dāng)成邊界層頂尺度太大邊界層頂?shù)奈恢脮?huì)偏移。我一般會(huì)試 5 到 20 之間的幾個(gè)尺度看哪個(gè)尺度下提取的高度和探空數(shù)據(jù)最接近。如果沒有探空數(shù)據(jù)就看時(shí)間序列是否連續(xù)如果邊界層高度在一天內(nèi)變化平滑說明尺度選得合適。驗(yàn)證方法也很重要。我習(xí)慣用兩種方式交叉驗(yàn)證一是和微波輻射計(jì)或者探空數(shù)據(jù)對(duì)比二是看時(shí)間高度圖上的結(jié)構(gòu)是否合理。如果提取的邊界層高度在時(shí)間序列上出現(xiàn)頻繁跳變那多半是算法參數(shù)沒調(diào)好。另外氣溶膠邊界層高度和云底高度容易混淆如果消光系數(shù)廓線在邊界層以上還有明顯的峰值那可能是云層需要區(qū)分開。最后說一個(gè)我自己的習(xí)慣每次處理完數(shù)據(jù)我都會(huì)把原始信號(hào)、距離校正信號(hào)、消光系數(shù)廓線、邊界層高度畫在一張圖上從頭到尾看一遍。如果中間哪一步出現(xiàn)異常圖上會(huì)很明顯。這個(gè)習(xí)慣幫我省了很多后悔藥。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取