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