張豪偉，韓 舸，馬 昕，龔 威，王 琦

（1.武漢大學 電子信息學院，湖北 武漢 430079；2.武漢大學 遙感信息工程學院，湖北 武漢 430079；3.武漢大學測繪遙感信息工程國家重點實驗室，湖北 武漢 430079；4.國家國防科技工業(yè)局 重大專項工程中心，北京 100101）

0 引言

政府間氣候變化專門委員會（Intergovernmental Panel on Climate Change，IPCC）第五次報告確認人為碳排放是全球變暖的主要因素。近年來，世界多個主要碳排放國家相繼確立了碳中和目標，標志著人類社會再次在共同應對氣候變化的議題上達成一致。準確、有效的碳監(jiān)測是制定和審查減排措施的關鍵工具，也是國際間合作的互信基礎。利用大氣二氧化碳（CO2）濃度觀測數(shù)據(jù)推算CO2通量，正逐步成為標準方法，可以補充現(xiàn)有“自下而上”核算體系的不足［1-2］。目前全球CO2濃度監(jiān)測方式包括地基觀測和被動衛(wèi)星觀測。其中，CO2觀測地面站點可以提供長時間、高精度的觀測資料，但其數(shù)量有限、分布不均，無法為高分辨的CO2通量監(jiān)測提供充足的觀測資料。在過去10 年中，被動碳遙感衛(wèi)星的快速發(fā)展，使得人類初步具備了獲取大范圍CO2濃度觀測資料的能力［3-5］。該衛(wèi)星觀測資料促進了區(qū)域通量反演、植被碳匯監(jiān)測、人為強點源核查和城市排放核算等應用科學的快速發(fā)展。與此同時，被動遙感技術(shù)的內(nèi)在不足也引起人們的關注。由于依賴太陽光進行觀測，且反演是一個欠定的優(yōu)化過程，其CO2濃度產(chǎn)品性能易受氣溶膠和云層干擾。高緯度和夜間觀測資料的缺乏也限制了相關科學問題的研究［3-8］。

從2005 年起，美國、歐洲和日本等發(fā)達國家就開始了天基激光碳監(jiān)測激光雷達的研制工作，期望利用主動觀測機制的優(yōu)勢與被動遙感形成互補。ABSHIRE 等［9］在2011 年進行了1.572 μm 差分吸收激光雷達ASCENDS 科學飛行實驗，基于10 s 的CO2濃度平均值，在相對恒定的二氧化碳柱濃度（XCO2）區(qū)域獲得了2×10-6～3×10-6的誤差；REFAAT 等［10］在2014 年進行了2.0 μm 差分吸收激光雷達實驗，并基于獲得的海洋和陸地剖面進行建模，其結(jié)果降低了CO2脈沖的距離偏差。除此之外，SINGH 等［11］在2017 年進行模擬實驗，根據(jù)10 s的數(shù)據(jù)平均，項目的預期性能低于0.35×10-6。我國在該領域的起步較晚，但進展很快。朱亞丹等［12］在中國進行了1.572 μm 差分吸收激光雷達實驗，通過對整個探測系統(tǒng)硬件進行評估，驗證了機載實驗的有效性。相成志等基于我國1.572 μm 差分吸收飛行實驗對XCO2反演進行了方案設計。結(jié)果表明，對于即將進行的陸地［13］衛(wèi)星任務，XCO2的預測相對隨機誤差小于0.3%。2022 年4月，大氣環(huán)境監(jiān)測1 號星的成功發(fā)射，標志著我國成為世界上首個具有激光碳監(jiān)測能力的國家。

天基遙感碳監(jiān)測的技術(shù)基礎是路徑積分差分吸收激光雷達（Integral Path Differential Absorption Lidar，IPDA Lidar）［14-23］。其通過差分吸收過程，抵消氣溶膠、地表反射等因素，精準的測量大氣CO2分子的吸收能力。該項技術(shù)是通過建立光學反演模型來反演大氣中的XCO2。本文利用2019 年在秦皇島地區(qū)獲得的機載激光雷達數(shù)據(jù)［24-25］，分析IPDA Lidar 反演XCO2需要解決的重要問題，并相應的提出改進方法。在IPDA Lidar 反演模型中積分權(quán)重函數(shù)上下界的精確確定，涉及氣象參數(shù)、飛機慣性導航信息、信號分解獲取測距距離這3 個部分的信息。在該3 個參數(shù)中，信號分解獲取到的CO2柱體長度信息在整個反演過程中起著非常重要的作用。因此，從準確測量CO2柱長的角度來降低CO2柱濃度反演誤差可行。除此之外，差分吸收光學厚度的獲取涉及對原始信號的處理。從參雜噪聲的回波信號中準確去除基線并估計出回波能量是測算差分吸收光學厚度的關鍵基礎。相對于常規(guī)的采用原始信號的離散數(shù)據(jù)積分，本文通過自適應的光譜能量進行建模，通過采用模型函數(shù)去積分獲取回波能量。提出一個基于光譜能量的模型，可以準確匹配非高斯信號波形。在實際數(shù)據(jù)采集中，特別是在山或云地表類型引起的多子波回波的信號，可以提高IPDA Lidar 系統(tǒng)的探測性能。

1 設計思路

為了獲得高精度的XCO2，提出的算法框架主要集中在2 個方面。首先，結(jié)合式（1）來優(yōu)化激光柱長的精度，從而從分母角度得到優(yōu)化積分權(quán)重函數(shù)；結(jié)合式（1）優(yōu)化波形能量積分，從分子的角度得到差分吸收光學厚度。算法框架流程如圖1所示。

圖1 算法流程Fig.1 Flow chart of the algorithm

1.1 模型函數(shù)構(gòu)建

Lidar 的回波信號被認為是幾個高斯信號的疊加，信號具有近似對稱的特點。但該特性通常會受到地形的影響而導致信號變寬，很難用一個方程來定量描述加寬的信號，且增大了差分吸收光學厚度和積分權(quán)重函數(shù)的計算誤差。針對1.572 μm IPDA激光雷達回波信號，因為常規(guī)的波形擬合不能用高斯模型準確描述，所以提出了一種基于頻譜能量模型的波形建模的新理論（式（3）），可以減少因信號展寬引起的DAOD 和IWF 的計算誤差。結(jié)合實驗中1.572 μm IPDA 激光雷達回波信號的實際情況，在式（2）的傅里葉變換形式的基礎上引入多個參數(shù)［26］，這樣修改后的模型效果可以減少波形的展寬效應。此外，信號的特征參數(shù)是通過整個回波信號公式的定量表達得到。頻譜能量公式和修改后的公式如下：

式中：a=1/4（1/τr-1/τf）；b=1/4（1/τr+1/τf），τr、τf分別為上升時間和下降時間。

式中：d為子波形參數(shù)中心，并且該參數(shù)也代表CO2柱體的長度；a、b、c和f需要從離散的信號數(shù)據(jù)中被擬合獲取基于波形方程。

此外，為了獲得多種土地利用類型回波信號的精確模型參數(shù)，采用了基于式（3）的列文伯格-馬夸爾特（LM）算法作為核心。LM 算法可以在以往模型參數(shù)的基礎上，結(jié)合相應的約束條件進行深度優(yōu)化，得到符合精度要求的非線性擬合結(jié)果。

1.2 模型函數(shù)對IWF 和DAOD 的優(yōu)化

對于積分權(quán)重函數(shù)的優(yōu)化，主要是對IWF 中的上限和下限進行修改。實驗結(jié)果表明，當CO2柱長度誤差小于3 m 時，在空間使用IPDA 激光雷達系統(tǒng)測得的XCO2柱誤差會小于1%。從精確測量CO2柱長的角度來看，通過結(jié)合上述模型的優(yōu)化參數(shù)d，改進IWF 反演是可行的。

IPDA 激光雷達發(fā)射器以200 μs 的間隔連續(xù)平行發(fā)射1.57 μm 的在線和離線光束。因此，該發(fā)射器可以被看作是對同一目標的兩個連續(xù)距離測量，從而產(chǎn)生了冗余的觀測測距信息。此外，盡管以式（3）為核心的LM 算法，但距離測量的觀測值相對于真實值仍有誤差。由于λon和λoff的波段距離是分別測量得到，所以兩個觀測值產(chǎn)生的距離測量值矛盾。為了平衡該矛盾，引入不等精度的直接調(diào)整理論。關于測量調(diào)整，一般認為權(quán)重是觀測結(jié)果可靠性的相對評價指標，其與中誤差的平方成反比。最有可能值公式用于獲得高精度范圍值，權(quán)重和最可能值公式如下：

式中：為λon或λoff對應的中值誤差；Pn為在λon或λoff對應的權(quán)重；L為上述數(shù)據(jù)處理后λon或λoff的測量值。

緊接著，結(jié)合積分權(quán)重函數(shù)式（6）獲取優(yōu)化后的參數(shù)IWF為

式中：IWF為積分權(quán)重函數(shù)計算的結(jié)果；PSFC和p=0為積分上、下限參數(shù)由優(yōu)化后的參數(shù)d、氣象參數(shù)和飛機慣性制導信息決定；為CO2積分權(quán)重函數(shù)。

對于差分吸收光學厚度的優(yōu)化，主要根據(jù)提出的模型方程（式（3））來擬合信號波形，以獲得對信號方程的準確描述。此外，將優(yōu)化后的參數(shù)輸入到DAOD 函數(shù)方程，然后對方程進行積分，以恢復信號的能量。因此，與傳統(tǒng)的數(shù)值離散積分和基于高斯核公式的積分相比，上述優(yōu)化方法可以更準確地獲得波形能量。其DAOD 函數(shù)方程如下：

式中：DAOD為差分吸收光學厚度函數(shù)；P、E的能量是由優(yōu)化后的公式參數(shù)組帶入到對應的式（3）計算所得。

在傳統(tǒng)的TSP問題中，解構(gòu)建圖為所有城市節(jié)點的連接網(wǎng)，邊的權(quán)值為兩城市節(jié)點間的距離。在本文中，首先需確定所有待編配車組與除晚高峰外的待編配車次的成本矩陣Cij，其中行代表車組號，列代表列車車次；車次按時間從早到晚依次排列，成本矩陣中的每一個值都作為一個節(jié)點vij(第i個車組擔任第j個車次)。在進行計劃編制時需按照車次時間順序依次進行編配，故將所有節(jié)點按從左到右的方向依次連接，當前節(jié)點僅可與右側(cè)相鄰列的所有節(jié)點連接，如圖1所示。

1.3 模型驗證

對于算法框架的驗證，主要是對模型參數(shù)d的驗證。如果模型參數(shù)得到準確，那么XCO2的精度也會有很大提高。在模型的參數(shù)組中，模型參數(shù)d也反映了CO2氣體的長度，而且參數(shù)d可以驗證。因此，模型參數(shù)d的驗證可以反映算法框架的準確性，也代表模型對XCO2反演的準確性。

模型驗證根據(jù)不同的地表類型，分別分為海洋和陸地。首先，在海面上用秦皇島的潮汐數(shù)據(jù)進行驗證；其次，用ICESat-2 的數(shù)據(jù)在平原和山區(qū)進行驗證。因此，需要重新評估經(jīng)緯度的偏移，以獲得激光雷達足跡的準確表面高程值。然后，根據(jù)圖2中D點的偏移激光雷達足跡坐標，將ICESat-2 的數(shù)據(jù)匹配為地表高程。受機載姿態(tài)角的影響情況，其導致真實的激光雷達足跡偏離了機下點如圖2 所示。其中，O點代表1.572 μm IPDA 激光雷達信號發(fā)射的位置；D點代表1.572 μm IPDA 激光雷達信號接觸到地面時的反射位置；而A點是機下點路徑（OA）上的一點，其高度與D點在同一平面上。C、A、B點在同一仰角平面上，A、O、C點形成的角度等于俯仰角。在圖2（c）中，N 和E 分別代表坐標軸上的北方和東方?？傊瑱C下點的位置（即GPS 值）并不代表激光雷達足跡的真實位置。

圖2 激光雷達足跡的偏離機下點的地理模型Fig.2 Geographical model of the deviation points of the lidar footprint

假設飛機的姿態(tài)角不發(fā)生變化，1.572 μm IPDA 激光雷達觀測到的機下點路徑是OA。然而，由于滾動角和俯仰角的綜合影響，實驗中1.572 μm IPDA 激光雷達觀測到的路徑是OD。所以AD是實際腳印和機下點的偏差距離（大約330 m）。因此，重新評估LiDAR 腳印的真實位置很重要，因為全球定位系統(tǒng)（GPS）數(shù)據(jù)并不代表激光雷達腳印的真實位置。對于該偏差距離的修正，使用以下公式來獲得經(jīng)緯度的偏差：

式中：LOD為被探測到的柱長距離，基于上述理論是圖2（a）中的OD；Pa為俯仰角；Ra為翻滾角度；VN為正北方向的速度；VE為正東方向的速度；LAD、LAC和LCD分別為圖2（a）中的AD、AC和CD距離；α為飛機的實際飛行方向與東面的偏差角度，并且α也代表飛機在實際飛行方向與東方之間的偏差角度；Δαβ為α和β角度差值；ΔLat和ΔLon為D相對于點A的偏移量。

1.4 模型評價指標

為了量化提出方法的合理性，選擇均方根誤差（RMSE），該指標用于量化模型的偏差。指標的計算公式為

式中：N為樣本總的數(shù)目；Pi為模型計算的值；Ri為真實觀測到的值。

2 實驗結(jié)果

2.1 基于提出的模型實驗結(jié)果

為了獲得高精度的XCO2，根據(jù)上述的算法框架對每組信號進行同樣的處理。在該過程中，建立模型方程，然后通過組合每個回波成分來獲得模型參數(shù)，以適應信號。單個信號中多個子組件的處理結(jié)果如圖3 所示。

圖3 模型處理結(jié)果Fig.3 Model processing results

2.2 模型驗證結(jié)果

為了進一步驗證算法框架，根據(jù)驗證思路和具體實驗條件，將驗證分為海面和陸地。對于海面的驗證思路，主要是將飛機提供的GPS 值與飛機平臺到海面的垂直CO2柱長和WGS 84 坐標系中的海面高度之差進行比較。而實時海面高度主要由2 部分組成，分別為驗潮站（秦皇島站）的高程和實時潮汐數(shù)據(jù)，如圖4 所示。驗潮站數(shù)據(jù)涉及1985 年高程系統(tǒng)數(shù)據(jù)與WGS 84 高程系統(tǒng)數(shù)據(jù)的轉(zhuǎn)換，轉(zhuǎn)換值是由EGM 2008 重力模型計算得到。

圖4 海面區(qū)域算法驗證結(jié)果Fig.4 Verification results of the sea area algorithm

在圖4 中，將GPS 高度的值與校正偏差后的值（該值為校正后的垂直CO2柱長與實時海面高度兩個值之和）進行比較。結(jié)果顯示，相對于約6 800 m 的測量高度，算法框架在海面上的精確度為0.94 m。此外，皮爾遜系數(shù)為1.00，表明數(shù)據(jù)之間在0.01 的水平上存在顯著的相關性，同時進一步表明算法框架在海面上的準確性具有很高的魯棒性。

為了驗證模型對陸地的影響，比較了GPS 數(shù)據(jù)與偏差修正值之間的差異如圖5 所示。

圖5 陸地區(qū)域算法驗證結(jié)果Fig.5 Verification results obtained by the land area algorithm

偏差修正值主要包括2 個項目，即垂直CO2柱距（即圖5 中的算法結(jié)果）和相應激光雷達足跡的表面高程值。ICESat-2_ATL08 地形的高度誤差平均值在平原地區(qū)小于0.05 m，在山區(qū)小于0.5 m［27］?？紤]到IPDA 激光雷達數(shù)據(jù)和ICESat-2_ATL08 數(shù)據(jù)都是WGS 84 坐標系，使用ICESat-2_ATL08 數(shù)據(jù)作為表面高程輸入數(shù)據(jù)。在圖5 中，展示了算法模型在陸地上的驗證結(jié)果，精度為6.20 m。該算法框架在陸地上的低精度主要是由于驗證數(shù)據(jù)集的原因。為了增加高精度驗證數(shù)據(jù)的數(shù)量，采集了從2018 年12 月到2021 年12 月的ICESat-2_ATL08 數(shù)據(jù)。而且使用了8 m 的矩形框范圍，以匹配ICESat-2 和1.572 μm IPDA 激光雷達的足跡數(shù)據(jù)。但如果局部地形復雜時，會帶來一些誤差。此外，與2019年3 月采集的1.572 μm IPDA 激光雷達數(shù)據(jù)相比，驗證數(shù)據(jù)存在一定的誤差。因此，對于模型的驗證，在海面上的驗證更具有代表性。

2.3 二氧化碳反演結(jié)果和精度評估

不同積分時間下的XCO2反演結(jié)果如圖6 所示。在圖6（a）中XCO2反演在海洋上的波動明顯大于陸地。圖6（b）清晰地展示了XCO2在整個飛行軌跡上的變化趨勢。海洋區(qū)XCO2最低，平原區(qū)XCO2最高。山區(qū)的XCO2低于平原，但高于海洋。本文研究的平原地區(qū)是一個包括秦皇島市的城區(qū)。前人的研究表明，不同的估算方法對該地區(qū)大氣-海洋CO2通量的估算結(jié)果不同，但所有的方法都證實了渤海是冬季碳匯區(qū)，因此圖6 的結(jié)果合理。圖6（a）中黑點為正常的點數(shù)據(jù)，藍點表示云層數(shù)據(jù)；圖6（b）中藍線為每10 s 對黑點滑動平均的結(jié)果。

圖6 XCO2反演結(jié)果可視化Fig.6 Visualization of the XCO2 inversion results

檢索了研究區(qū)域2019 年3 月14 日前后3 天的OCO-2 產(chǎn)品數(shù)據(jù)集，發(fā)現(xiàn)2019 年3 月14 日和16 日的XCO2產(chǎn)品可用。如圖7 所示，OCO-2 在3 月16日的XCO2觀測值與研究區(qū)比較接近，且數(shù)量較高。3 月16 日OCO-2 的XCO2波動在401.4×10-6～420.6×10-6，平均值為412.83×10-6。3 月14日，OCO-2 的XCO2波動在410.8×10-6～420.0×10-6，平均值為414.65×10-6。OCO-2 的XCO2反演結(jié)果表明海洋和陸地之間存在CO2梯度。除此之外，機載差分吸收激光雷達的XCO2與OCO-2 的XCO2存在明顯的差異。機載差分吸收激光雷達的XCO2均值為416.26×10-6，而OCO-2 在3 月14 日和3 月16 日的XCO2均值分別為414.65×10-6、412.83×10-6。這些差異可以歸結(jié)為以下原因：1）OCO-2 的XCO2與機載IPDA 激光雷達路徑相比，OCO-2 的XCO2產(chǎn)物沒有覆蓋城市地區(qū)。人們普遍認為城市CO2高于郊區(qū)。2）相對于整個大氣柱的XCO2，本文中機載差分吸收激光雷達觀測的XCO2實際上只是部分XCO2。由于對流層低層的XCO2顯著高于高海拔地區(qū)，因此機載差分吸收激光雷達觀測的XCO2明顯高于全部OCO-2 觀測到的整個對流程的XCO2。另外，在基于差分吸收激光雷達觀測的XCO2和OCO-2 衛(wèi)星基于光譜儀觀測的XCO2的比較中，存在一個非常復雜的問題。目前OCO-2/3、GOSAT-1/2、TanSat 等類似傳感器的檢索算法都是使用平均核來計算最終的XCO2。然而，雙波長差分吸收激光雷達不產(chǎn)生平均核，而是直接計算柱加權(quán)的XCO2，因此，直接比較2 種產(chǎn)品存在問題。

圖7 OCO-2 和路徑積分差分吸收激光雷達XCO2的可視化Fig.7 Visualization of the OCO-2 and path integral differential absorption lidar XCO2

3 結(jié)束語

隨著XCO2監(jiān)測精細化需求的提高，衛(wèi)星（如OCO-2/3、GOSAT-1/2 和Tan Sat）很難滿足基于被動探測理論高精度、全天和全天候的需求。目前全球迫切需要開發(fā)主動探測衛(wèi)星來支持全球碳源、碳匯和碳中和工作的開展。我國已經(jīng)在2022 年4 月16 日發(fā)射大氣環(huán)境監(jiān)測衛(wèi)星，該衛(wèi)星配備了CO2路徑積分差分吸激光雷達。星載差分吸收激光雷達在夜間以及高緯度地區(qū)和嚴重污染地區(qū)的有效觀測方面可以補充目前被動遙感技術(shù)。盡管夜間、白天和季節(jié)的CO2排放主動傳感、地球行星的高級空間碳和氣候觀測致力于填補該空白，但目前仍沒有基于激光雷達的衛(wèi)星在軌道上進行CO2探測。因此，完善機載（星載激光雷達的縮小原型）差分吸收激光雷達關鍵性技術(shù)開發(fā)，具有重要的科學意義。

基于該衛(wèi)星的縮比機載實驗，開展了相應關鍵算法技術(shù)的研究，其主要工作包括：1）提出了波形能量方程，針對差分吸收激光雷達信號非高斯性、復雜性和波動性，預先構(gòu)建了波形能量方程，其方程較好地解決了信號非高斯問題；然后基于波形能量方程式，分別提取出了后續(xù)計算CO2濃度所需要的模型參數(shù)組。解決了參數(shù)多方程傳遞的問題，也有效降低了參數(shù)間接攜帶的誤差。2）設計了差分吸收激光雷達信號模型參數(shù)優(yōu)化方法，針對差分吸收激光雷達信號非高斯性的特點，結(jié)合傳統(tǒng)的LM算法，通過替換LM 算法核心。有效地實現(xiàn)了波段on 和波段off 模型參數(shù)組的優(yōu)化，為后續(xù)計算光學厚度提供精確的參數(shù)輸入。3）設計了雙探測誤差消除方法，因為差分吸收激光雷達需對同一目標連續(xù)2 次距離測量，從而帶來測量信息的冗余觀測以及2 次測量結(jié)果相互不匹配的問題。將優(yōu)化后的模型參數(shù)組的參數(shù)don和doff進行不等精度的直接平差，從而獲取準確的參數(shù)d。為后續(xù)計算積分權(quán)重降低了很大的誤差，從而間接的提高了CO2濃度的計算。4）模型算法進行了驗證，通過在秦皇島的實驗，表明提出的模型在海上和陸地上的精度分別為0.74 和6.20 m。與1.572 μm 的差分吸收激光雷達空間分辨率1.2 m 相比，海面模型驗證的精度表明，算法精度已進入亞單位。而以OCO-2 衛(wèi)星的XCO2作為參考，本文算法可以得到趨勢一致的XCO2。特別是在10 s 滑動平均值下，海洋、城市和山區(qū)的平均XCO2值分別為411.07×10-6、425.71×10-6、417.87×10-6，標準差分別為1.93×10-6、0.85×10-6、0.96×10-6。綜上所述，該算法提高了XCO2反演的準確性和穩(wěn)健性，對我國近年來即將發(fā)射的衛(wèi)星搭載的差分吸收激光雷達具有重要的參考意義。