摘要:激光外差光谱探测由于其光谱分辨率高、体积小、重量轻等优点近年来得到了快速的发展, 可用于大气温室气体垂直廓线测量和碳卫星地面定标等. 本文报道了利用3.939 μm带间级联激光器作为本振光源的测量大气N₂O的激光外差系统, 自制高精度太阳跟踪仪收集太阳光作为激光外差的信号光源, 其跟踪精度达到7 arcsec, 激光外差系统的光谱分辨率达到0.004 cm^–1, 测量了合肥地区(31.902°N, 117.167°E)大气N₂O吸收光谱, 得到2838.336和2539.344 cm^–1两个强吸收峰, 并对吸收信号进行波长标定, 得到了N₂O分子的整层大气透过率谱, 信噪比为93. 将高分辨率光谱数据进行归一化处理和频率校正, 利用参考正向模型和最优估计算法得到N₂O大气整层浓度廓线, 标准偏差体积分数为0.000031 × 10^–6—0.0026 × 10^–6, 对应相对误差范围为0.009%—0.83%. 研究结果表明, 所搭建的激光外差系统能够实现对大气中N₂O的吸收光谱测量以及对N₂O的廓线反演, 为长期观测大气N₂O浓度提供保证.
关键词: 激光外差/
透过率光谱/
氧化亚氮(N₂O)/
反演计算

English Abstract

全文HTML

--> --> -->

氧化亚氮(N₂O)是大气中仅次于CO₂和CH₄的一种重要的温室气体. 全球大气中氧化亚氮的体积分数已从工业化前约270 × 10^–9, 增加到现在的340 × 10^–9[1]. 氧化亚氮总排放量中超过三分之一是来源于人为排放^[2], 主要来自于农业, 但是目前对N₂O的时空分布尚不明确. 因此, 持续准确地检测大气中N₂O的时空分布对气候变化的影响具有重要意义.
基于光谱技术的大气探测方法主要有傅里叶变换红外光谱技术^[3,4]、激光外差光谱技术^?[5]. 激光外差光谱技术因其光谱分辨率高、空间分辨率高、光谱采样时间短、体积小、灵敏度高等优点被广泛应用, 其利用窄线宽的本振激光与太阳光或其他光源得到携带整层大气分子吸收信息的吸收光谱. Palmer等^[6]、Weidmann等^[7]和Tsai等^[8]利用外腔量子级联激光器建立的中红外激光外差测量系统, 实现了对大气中CH₄, CO₂, H₂O, O₃和N₂O等多种温室气体的测量. 美国NASA的Wilson等^[9,10]利用小型的分布反馈式激光器建立了全光纤近红外激光外差测量系统, 获取了大气中CH₄与CO₂的吸收光谱, 同时对装置进行了集成化设计. Rodin等^[11]利用平衡探测器建立激光外差测量系统获取了大气中CH₄与CO₂的吸收光谱. 国内, 中国科学院安徽光学精密机械研究所也开展了以外腔量子级联激光器与带间级联激光器为本振光源的自由空间结构激光外差光谱技术的研究^[12,13].
本文基于激光外差探测技术结合反演算法, 设计了一套用于测量大气N₂O浓度的激光外差系统, 其中使用3.939 μm的带间级联激光器(interband cascade laser, ICL)作为光源, 采用N₂O位于2538.336 cm^–1的吸收谱线对大气中的N₂O浓度进行了探测. 实验中采用自制高精度太阳跟踪仪采集太阳光作为信号光源, 利用ICL作为本振光源, 两束空间光进行光学混频后, 差频信号经过光电探测器、射频电路及锁相放大器, 并最终由美国国家仪器(National Instruments, NI)采集卡采集得到高分辨率光谱数据. 将所采集的光谱数据进行归一化处理、频率校正之后, 结合参考正向模型和最优估计算法反演得到整层大气N₂O的浓度廓线.

激光外差原理在文献[14]中已有详细描述, 在此只做简单描述. 微弱的输入信号与较强的本振(LO)信号混合在非线性电路元件中, 差频或中频信号(IF)由滤波器选择后, 在窄带固定频率放大器中放大并检测, 中频信号可以选择落在可用放大装置的频谱范围内. 图1给出了显示激光外差的原理示意图. 两束光束叠加在一个混频器上, 然后经过放大后检测中频信号.
图 1 激光外差原理图
Figure1. Schematic diagram of laser heterodyne.

大多数光学探测器是电场平方律器件, 响应与光照强度成正比, 假设两个振幅为A_s和A_LO、角频率为$ {w}_{\mathrm{s}} $和$ {w}_{\mathrm{L}\mathrm{O}} $且线性极化的完全对准平面波入射到探测器上. 探测器所产生的总光电流理论值为

$ i\left(t\right)=\alpha {e}^{2}\left(t\right)=\alpha ({A}_{\mathrm{s}}{\cos}{w}_{\mathrm{s}}t+{A}_{\mathrm{L}\mathrm{O}}{\cos}{w}_{\mathrm{L}\mathrm{O}}t{)}^{2} . $

经过三角函数和差化积公式简化为

$ \begin{split}i\left(t\right)=\;&\alpha \left({A}_{\mathrm{s}}^{2}\frac{1+ \cos 2{w}_{\mathrm{s}}t}{2}+{A}_{\mathrm{L}\mathrm{O}}^{2}\frac{1+ \cos 2{w}_{\mathrm{s}}t}{2}\right) \\&+\alpha {A}_{\mathrm{s}}{A}_{\mathrm{L}\mathrm{O}} \big[\cos ({w}_{\mathrm{s}}t + {w}_{\mathrm{L}\mathrm{O}}t) \\& + \cos({w}_{\mathrm{s}}t-{w}_{\mathrm{L}\mathrm{O}}t)\big] .\end{split} $

式中, ${\cos}2{w}_{\mathrm{s}}t, {\cos}2{w}_{\mathrm{L}\mathrm{O}}t, {\cos}({w}_{\mathrm{s}}t+{w}_{\mathrm{L}\mathrm{O}}t)$频率远远大于光电探测器电路带宽, 所以积分输出为0. 如果$ {w_{\text{s}}} $和$ {w_{{\text{LO}}}} $足够接近, $ {w_{\text{s}}} - {w_{{\text{LO}}}} $在光电探测器的低通截止频率之下. 那么(2)式可以进一步简化为

$ i\left(t\right)=\frac{\alpha \left({A}_{\mathrm{s}}^{2}+{A}_{\mathrm{L}\mathrm{O}}^{2}\right)}{2}+\alpha {A}_{\mathrm{s}}{A}_{\mathrm{L}\mathrm{O}}{\cos}({w}_{\mathrm{s}}t-{w}_{\mathrm{L}\mathrm{O}}t) \text{, } $

其中, 光电探测器输出信号包含两个部分: 直流信号和交流信号.
直流信号表达式:

$ {i}_{\mathrm{d}\mathrm{c}}= {\alpha ({A}_{\mathrm{s}}^{2}+{A}_{\mathrm{L}\mathrm{O}}^{2})}/{2} . $

交流信号表达式:

$ {i}_{\mathrm{a}\mathrm{c}}=\alpha {{A}}_{\mathrm{s}}{A}_{\mathrm{L}\mathrm{O}}\mathrm{cos}({w}_{\mathrm{s}}t-{w}_{\mathrm{L}\mathrm{O}}t) . $

直流信号反映了光电探测器对本振光和信号的直流响应, 交流信号即是激光外差信号, 如果本振光和信号光都是理想的单色光, 那么交流部分的信号是一个频率为两束光光频之差$ {w_{\text{s}}} - {w_{{\text{LO}}}} $的正弦波. 实际上, 外差的功率表达式为

$ {P}_{\mathrm{a}\mathrm{c}}={i}_{\mathrm{a}\mathrm{c}}^{2}{R}_{\mathrm{L}} \text{, } $

其中$ {R_{\text{L}}} $为光电探测器的负载电阻.

2

-->

3.1.太阳跟踪仪

激光外差技术可以在很小的仪器体积下获得很高的光谱分辨率, 光谱分辨率约为几十兆赫兹. 为了提高外差探测精度, 过去主要是在光谱仪本身进行改进, 例如仪器线性函数(ILS)表征^?[15], 在干涉图上应用直流信号(DC)校正、改进了采样精度和探测器非线性的表征. 为了达到足够的准确度, 必须考虑激光外差探测技术对信号光和本振光的空间模式有严格的要求, 而太阳跟踪仪为激光外差探测提供信号光, 必须保证其提供的信号光具有足够高的精度及稳定度, 才能满足测量的要求.
太阳跟踪仪跟踪方式有基于太阳运行轨迹方式和光电跟踪两种方式, 基于太阳运行轨迹方式在任何气象条件下都能正常跟踪, 但是没有形成环路反馈, 精度有所欠缺, 光电跟踪方式利用太阳光斑进行实时跟踪, 缺点是在阴天气候无法工作. 因此我们将两种跟踪方式结合起来提高跟踪精度, 首先利用已知天文数据算出太阳高度角h₀和方位角A.
太阳高度角计算公式为

$ {h}_{0}= \arcsin(\sin \theta \sin\beta +\cos \theta \cos \beta \cos \gamma ) . $

太阳方位角计算公式为

$ A= \arccos\left(\right(\sin{h}_{0}\sin\beta -\sin \theta )/\cos{h}_{0}\cos \beta ) \text{, } $

其中θ为赤纬角, β为纬度, γ为时角. 上位机Labview算法程序在计算获得太阳高度角和方位角之后, 得到与当前两角度差值, 转换成步进电机所需要的转动步数, 再由Labview上位机程序控制电机转动, 实现太阳运动轨迹跟踪. 上述跟踪完成后, 在图像传感器判断是否接收到太阳光斑信息, 如果没有接收到太阳光斑信息, 系统继续运行太阳运动轨迹跟踪程序, 一直到图像传感器接收到太阳光斑信息, Labview程序通过对太阳光斑进行二值化和中值滤波处理, 获得太阳光斑质心坐标, 计算太阳光斑质心坐标与视场中心坐标偏差, 并转换成步进电机转动步数, 用PID (proportion-integral-differential)控制算法实时不断的调整电机转动步数.
太阳跟踪仪部分的结构装置图如图2所示, 太阳光经过两面安装在电机上的高反镜反射, 白片反射的太阳光穿过巴德膜之后打到CCD光电探测器, 穿过氟化钙白片的光与本振激光拍频产生外差信号.
图 2 (a)太阳跟踪仪内部装置; (b)精跟踪模式下太阳光斑质心偏移量
Figure2. (a) Internal device of suntracker; (b) sunlight spot centroid offset in precision tracking mode.

如图3所示, 太阳跟踪仪在精跟踪模式下, X轴的偏移量在[–0.3, 0.3]像素之间, Y轴偏移量在[–0.5, 0.5]像素之间, 1 Pixel约等于14 arcsec, 所以整体跟踪精度约为7 arcsec.
图 3 系统结构原理图(LIA, 锁相放大器; PC, 电脑; DAQ, 数据采集卡; FL1, 2, 聚焦透镜1, 2; BS1, 2, 分束器1, 2)
Figure3. System structure schematic diagram. LIA, lock-in amplifier; PC, personal computer; DAQ, data acquisition card; FL1, 2, focus lens 1, 2; BS1, 2, beam splitter 1, 2.

2 -->

3.2.实验装置

实验中基于大气中N₂O的吸收谱线, 选取了Nanoplus公司生产的3.9 μm ICL作为本振光源, 波数范围为2537.5—2540.5 cm^–1, 激光器输出约为1.8 mW, 采用LDC501激光控制器控制. 在该波段N₂O大气吸收占前七种大气吸收分子的95%以上, 其他气体吸收的干扰较小.
如图3所示, 使用上位机Labview程序控制LDC501进行逐点步进电流扫描, 同时上位机控制程序控制太阳跟踪仪实时采集太阳光, 经过斩波之后与激光在合束器上进行拍频, 通过滤光片和聚焦透镜后打到探测器上, 探测器输出的射频信号经过带通滤波器、检波器输入到锁相放大器, 锁相解调后, 通过采集卡采集光谱数据. 其中, 探测器的DC输出用来监测激光功率的稳定性, 透过半反镜的部分本振激光用波长计进行实时标定. 图4为激光外差实验整体结构的实物图.
图 4 系统结构实物图
Figure4. Physical picture of system structure.

2

-->

4.1.数据分析

实验过程中通过频谱密度图(图5所示)选用了10—70 MHz的带通滤波器, 避免了频谱分布中噪声较大的部分, 在保证光谱分辨率的同时为了提高系统的信噪比, 适当增大了积分时间和滤波带宽, 系统参数设置如表1所列.
图 5 (a) 射频信号功率谱; (b) 仪器函数
Figure5. (a) RF signal power spectrum; (b) instrument functions.

参数	数值
太阳光传输效率	0.305
滤波带宽/MHz	60
太阳温度/K	5773
测量波长/m	3.939
积分时间/s	0.3

表13.9 μm激光外差系统的参数设置
Table1.Parameter setting of 3.9 μm laser heterodyne system.

实验中, 本振激光使用逐点扫描方式控制激光器电流, 系统的仪器函数主要由电信号响应带宽决定. 太阳光作为一个宽带光源用来进行外差拍频时, 在探测器的整个响应带宽(0—200 MHz)都有均匀的信号分布. 实际测量时, 考虑电子噪声和光谱分辨率的要求, 需要选择添加特定波段的带通滤波器. 使用频谱分析仪(Agilent Technologies CXA Signal Analyzer N900A)采集并分析了系统的电子特性. 测量结果如图5(a)所示, 其中采集了本振激光3.939 μm的功率谱, 并分别选择加入45—98 MHz以及10—70 MHz带通滤波器进行采集, 在100 MHz左右有明显的系统电子噪声, 最后选择了噪声较低的区域(10—70 MHz).
根据激光外差原理以及频谱仪采集的功率谱归一化得到系统仪器函数, 如图5(b)所示, 系统实际双边带宽为120 MHz, 光谱分辨率为0.004 cm^–1.
信噪比是激光外差探测系统的1个重要参数, 可以反映系统的最小探测能力. 仪器函数则可体现外差系统可分辨的最窄谱线宽度, 所以信噪比和仪器函数是激光外差系统的两个重要参数, 提高信噪比和获得准确的仪器函数也是反演气体浓度廓线的必要工作.
对外差系统的信噪比(signal-to-noise ratio, SNR) S/N测量, 1966年Siegman^[16] 提出以黑体辐射光源为目标光源的的外差系统的理论信噪比:

$ \frac{S}{N}=\dfrac{2{T}^{*}\eta }{\mathrm{exp}\left[{hf}/({KT})\right]-1} $

外差探测的大部分目标光源都可以等效为一个特定温度的黑体, 所以这个信噪比公式可以简单的估计外差装置的信噪比等参数, 为激光外差系统的建立提供依据. 之后通过相关作者的研究发现通过锁相放大等探测器技术可大大提高系统的信噪比, 外差装置的信噪比可以写为

$ \frac{S}{N}=\dfrac{2{T}^{*}\eta} {\exp[{hf}/({KT})]-1}{\left(B\tau \right)}^{1/2} $

其中$ B $是探测带宽, $ \tau $是积分时间.
实验中首先将黑体温度设置为1273 K, 不经过吸收直接与本振激光(3.939 μm)拍频, 记录锁相输出的外差信号. 实验其他参数保持不变, 将黑体光完全遮挡, 再次记录输出信号. 两次实验中, 斩波频率为75 Hz, 滤波带宽为10—70 MHz, 锁相积分时间为1 s. 实验结果如图6所示, 此状态下的信噪比约为93.
图 6 外差装置的信噪比测量
Figure6. SNR measurement of heterodyne device.

2 -->

4.2.大气透过率谱

大气中氧化亚氮(N₂O)的吸收光谱由激光外差系统进行电流逐点步进扫描获得, 如图7所示, 激光逐点扫描范围内存在两处N₂O强吸收, 分别位于2838.336 和2539.344 cm^–1处. 实验中电流扫描间隔为0.1 mA, 对应波数变化小于0.004 cm^–1, 图7(a)为激光外差系统测量的N₂O的外差信号, 图7(b)为ICL波数与注入电流的关系曲线. 实验地点位于31.902°N, 117.167°E的安徽省合肥市科学岛.
图 7 (a) N₂O吸收的外差信号; (b) 激光器波数与注入电流的关系
Figure7. (a) Heterodyne signal absorbed by N₂O; (b) the relationship between laser wavenumber and injection current.

实验过程中, 太阳光波动会对外差信号产生极大影响, 因此实验过程中需要保持太阳光功率相对稳定. 对图7直接采集的外差信号进行一次多项式拟合谱线, 再做归一化处理后得到氧化亚氮(N₂O)的大气透过率谱. 由图8所示, 实验测得的N₂O吸收光谱与Hirtran数据库的吸收光谱吻合. 由测量结果(图7(a))得到此系统信噪比(SNR)可以达到97, 高于预期效果. 实验所得两个强吸收峰(2838.336和2539.344 cm^–1)可用于对本振激光进行波长校正, 为反演大气氧化亚氮廓线和柱总浓度提供准确的光谱数据.
图 8 实验所得N₂O吸收光谱与模拟吸收光谱的对比
Figure8. Comparison of the N₂O absorption spectrum obtained in the experiment and the simulation absorption spectrum.

2 -->

4.3.反演计算

图9为激光外差光谱(laser heterodyne spectroscopy, LHR)数据反演的流程, 由参考正向模型、实验数据处理和迭代反演等组成. 首先读取温度、压力、先验廓线及实验所得的光谱数据, 利用数据计算出初始透过率, 然后对实验数据进行预处理, 包括去除太阳光波动干扰和背景信号、波长校准^[17]等, 计算出初始Jacobian矩阵和增益系数, 用于迭代计算. 最后进行循环迭代, 这一步为反演算法的核心. 迭代反演采用Rodger^[18,19]提出的最佳估计方法(optimal estimation method, OEM). 对于OEM在基于LHR大气探测中的应用, Weidmann等^[20]已经做了详细的介绍. 在此仅简要说明, 测量数据与大气状态向量之间的关系被描述为
图 9 数据反演流程图
Figure9. Data inversion flow chart.

$ {{\boldsymbol{y}}_m} = F({\boldsymbol{x}},a,b,c) + {\boldsymbol{\varepsilon}} \text{, } $

其中y_m是m维的测量向量; x是n维参考正向模型中状态向量; ε是误差向量, 二阶多项式基线(系数为a, b和c)表示目标光谱区域中的非选择性吸收和未校正的功率.
OEM为最小化损耗函数χ²的迭代运算过程, 其中损耗函数表示为

$ {\chi ^2} = {({\boldsymbol{y}} - {\boldsymbol{F}})^T}S_\varepsilon ^{ - 1}({\boldsymbol{y}} - {\boldsymbol{F}}) + {({x_i} - {x_a})^T}{\boldsymbol{S}}_a^{ - 1}({x_i} - {x_a}) , $

其中, S_ε是测量协方差矩阵, S_a是先验协方差矩阵, x_a是先验N₂O廓线. 迭代过程采用Levenberg-Marquardt (LM)方法:

$\begin{split} \;&{x_{i + 1}} = {x_i} + {[(1 + \gamma ){\boldsymbol{S}}_a^{ - 1} + {\boldsymbol{K}}_i^{\text{T}}{\boldsymbol{S}}_\varepsilon ^{ - 1}{{\boldsymbol{K}}_i}]^{ - 1}}\\& \qquad \times [{\boldsymbol{K}}_i^T{\boldsymbol{S}}_\varepsilon ^{ - 1}({y_m} - {F_i}) - {\boldsymbol{S}}_a^{ - 1}({x_i} - {x_a})],\end{split}$

其中K为雅可比矩阵 (或加权函数), γ 为Levenberg-Marquardt (LM) 参数.
在反演过程中, 通过LM算法迭代计算调用参考正向模型以最小化损耗函数, 最终反演出N₂O的廓线. 图10(a)显示了经过预处理的激光外差实验数据和在2537.5—2539.5 cm^–1的光谱区域内的拟合光谱, 拟合结果之间的残差数据如图10(b)所示, 其随机分布在$ \pm $0.04左右表明拟合效果良好.
图 10 LHR数据反演结果　(a) 实验和拟合LHR谱图; (b) 残差
Figure10. Inversion results of LHR data: (a) Experimental and fitted LHR spectrogram; (b) residual.

如图11所示, 数据反演所采用的地表温度、压力数据来自于欧洲中期天气预报中心(ECMWF)对于N₂O的先验廓线, 使用了典型的中纬度冬天廓线数据. 大气分层采用40层, 从地表到70 km高空, 高度网格依次为1, 2.5和5 km.
图 11 温度和压力廓线
Figure11. Temperature and pressure profiles.

图12显示了2019年11月20日的连续测量反演结果, 实验对中午12:00左右的测量结果进行了分析, 这是因为欧洲中期天气预报中心(ECMWF)和中国气象数据服务中心CMCC提供的时间点为12:00, 该段时间内大气参数是稳定的. 图12显示了N₂O先验浓度廓线和反演结果所得浓度廓线, 其中先验浓度廓线结果选取在以观测点 (合肥, 31.9°N, 117.166°E ) 为中心, 6°半径范围内的数据. 在测量过程中, 体积分数标准偏差为0.000031 × 10^–6—0.0026 × 10^–6, 对应相对误差范围为0.009%—0.83%, 与ECWMF结果吻合较好, 可以对大气中的N₂O进行长期观测, 进而展开季节和年度的变化研究.
图 12 N₂O的先验和反演的垂直浓度分布图(VMR表示体积混合比)
Figure12. Prior and inversion vertical concentration profiles of N₂O. VMR, volume mixing ratio.

以3.939 μm激光器作为本振光源, 采用自制高精度太阳跟踪仪采集太阳光信号, 建立了一套测量整层大气氧化亚氮(N₂O)浓度廓线的中红外激光外差探测系统. 利用光谱分辨率为0.004 cm^–1的激光外差系统对实际大气中氧化亚氮的透过率谱测量, 得到信噪比为97的吸收光谱. 高分辨率大气吸收光谱经过信号预处理, 利用参考正向模型以及最优估计算法实现对整层N₂O浓度廓线的反演, 体积分数标准偏差为0.000031 × 10^–6—0.0026 × 10^–6, 对应相对误差范围为0.009%—0.83%. 后续将对测量系统的噪声、仪器函数进行优化处理, 从而提高系统的信噪比, 同时达到更高的分辨率, 从而进行长期观测.