摘要:多普勒差分干涉仪通过测量干涉图相位变化反演大气风速. 干涉图信噪比是工程应用中影响测风干涉仪相位反演精度的核心指标之一, 基于观测数据分析多普勒差分干涉仪中相位不确定度与干涉图信噪比之间的定量关系建立解析模型, 对仪器设计、性能评价及风场数据应用具有重要意义. 本文基于傅里叶变换光谱学噪声传递理论和多普勒差分干涉仪相位反演模型, 建立干涉图信噪比与相位不确定度传递的理论模型. 利用多组不同信噪比仿真及实验实测数据, 将相位反演模型计算的统计结果与本文提出的理论模型分析结果进行比较验证, 结果表明, 多组实验数据平均残差小于0.07 mrad, 两者的一致性较好, 建立的模型为多普勒差分干涉仪效能评估、干涉仪设计优化提供理论依据.
关键词: 相位反演/
多普勒差分干涉仪/
大气风场探测/
光子噪声

English Abstract

全文HTML

--> --> -->

大气风场是描述空间环境特征的重要物理参数, 是空间天气与空间环境学研究领域的热点课题, 也是空间科学中最有应用价值的研究内容之一^[1-3]. 星载测风干涉仪获取中高层大气风场廓线的长期观测数据, 对空间天气预报、航天器发射、运行和返回轨道参数计算、无线电通信、空间科学实验等都具有重要意义, 为研究上下层大气之间的能量传输、大气成分变化以及大气波现象提供重要理论依据^[4-6].
多普勒差分干涉光谱技术是由Englert团队在2006年提出的一种新型的中高层大气风场探测技术^[6-9], 多普勒差分干涉仪是通过探测气辉谱线相位变化量反演多普勒频移, 进而计算大气风速, 相位反演精度决定风速测量精度^[10-15]. 在实际应用中, 仪器采集到的原始干涉图会存在光电转换、探测器暗电流、读出电路带来的各种随机噪声, 这些噪声会随相位反演模型传入相位, 导致反演相位随机误差. 在仪器工程设计中通常采用原始干涉图信噪比作为衡量仪器信号质量的重要指标, 是系统基本光学、电子学参数确定的依据^[16-20]. 因此, 研究多普勒差分干涉仪干涉图信噪比与反演相位不确定度之间的传递关系, 建立定量数理模型对指导仪器工程设计是非常重要的.
传统傅里叶变换光谱仪以反演光谱为目的, 已建立了从原始干涉图到复原光谱噪声传递关系, 但由于多普勒差分干涉仪目标是测量多普勒频移带来的相位漂移, 故已有的傅里叶变换型光谱仪数据反演的噪声传递理论模型不能直接用于衡量多普勒差分干涉仪相位不确定度与原始干涉图信噪比的关系. Luo^[17]分析了基于对称式的空间外差干涉光谱仪中以光子噪声为主要噪声源时的干涉域信噪比表达式, 提出空域与频域信噪比的数理模型, 但无法说明原始干涉图信噪比与相位不确定度之间的关系. Pritt^[19]针对传统的傅里叶变换光谱仪, 建立原始干涉图噪声与反演光谱噪声之间的传递模型, 但未进一步证明其与复干涉图噪声和相位噪声之间的传递关系.
本文基于傅里叶变换光谱学基础理论, 根据多普勒差分干涉仪的相位反演方法, 结合原始干涉图表达式、系统噪声模型及相位反演模型, 建立原始干涉图信噪比与相位不确定度传递的理论模型, 并利用实验室研制的多普勒差分干涉仪验证其正确性, 分析误差来源, 为仪器的工程设计与性能评价提供重要理论基础.

2

-->

2.1.多普勒差分干涉仪原理

多普勒差分干涉仪是将迈克尔逊干涉仪中两干涉臂末端的平面反射镜替换成闪耀光栅, 在一个干涉臂中加入非对称量$\Delta d$, 增大光程差, 提高干涉仪风速测量灵敏度. 如图1所示. 当任意单色光谱输入$B(\sigma )$进入仪器入瞳, 经过准直镜、滤光片和分束器后形成两束光, 并分别在两个干涉臂末端的闪耀光栅处发生衍射, 经过闪耀光栅后返回, 两束光的波阵面在分束器处相遇, 经过条纹成像系统L2、L3, 在探测器上形成的干涉条纹I满足^[21]:
图 1 多普勒差分干涉仪原理图
Figure1. Schematic of the Doppler asymmetric spatial heterodyne interferometer.

$\begin{split}I(x) =\;& \frac{1}{2} \int\nolimits_0^\infty {B(\sigma )} \left( {1 + \cos \left\{ {2{\text{π}} } \right.} \right.\\& \times\left. {\left. [4 (\sigma - {\sigma _{\rm{L}}}) \tan {\theta _{\rm{L}}}] x + 2 \Delta d \sigma\right\}} \right) {\rm{d}}\sigma \end{split}$

式中, $\sigma $为入射光波数, ${\sigma _{\rm{L}}}$为Littrow波数, ${\theta _{\rm{L}}}$为Littrow角, x为探测器上位置坐标, 定义探测器中心$x = 0$, $\Delta d$为两干涉臂非对称量.
2 -->

2.2.相位反演模型

多普勒差分干涉仪通过探测干涉图的相位变化反演气辉谱线的多普勒频移, 相位反演模型如图2所示, 按照参考文献[6], 将式(1)简化为
图 2 相位反演模型流程图
Figure2. Flow chart of phase retrieval method.

$\begin{split} {I_{\rm{D}}}(x) =\, & \sum\limits_j {{S_j}} [1 + {E_j}(x)\cos (2{\text{π}}{k_j}x + {\phi _j} + \text{δ} {\varphi _j})] \\ =& \!\sum\limits_j {S_j}\Big( 1+ \frac{1}{2}{E_j}(x) \big\{ \exp [{\rm i}(2{\text{π}}{k_j}x \!+\! {\phi _j} \!+ \!\text{δ} {\varphi _j})] \\& + \exp [ - {\rm i}(2{\text{π}}{k_j}x + {\phi _j} + \text{δ} {\varphi _j})]\big\}\!\Big),\\[-15pt]\end{split} $

式中, j为通过滤光片进入仪器的谱线, ${S_j}$为谱线强度系数, ${E_j}(x)$为第j根谱线的包络函数, ${k_j} \equiv $$ 4({\sigma _j} - {\sigma _L})\tan {\theta _L}$为差分频率, ${\phi _j} \equiv 4{\text{π}}({\sigma _j} - \sigma _{\rm{L}})\Delta d$为附加相位项, $\text{δ} {\varphi _j}$为多普勒频移引起的相移项. 对(2)式傅里叶变换可得空间频率在$ + {k_j}$和$ - {k_j}$有特征峰的反演光谱S, 用窗函数分离S中的一个特征峰($ + {k_j}$或者$ - {k_j}$), 其余的所有元素均赋0, 傅里叶逆变换后可得:

$\begin{split} I_{\rm{D}}^0(x) =\,& \frac{1}{2}{S_0}{E_0}(x)\exp [{\rm i} (2{\text{π}}{k_0}x + {\phi _0} + \text{δ} {\varphi _0}) \\=\, & \frac{1}{2}{S_0}{E_0}(x)[\cos (2{\text{π}}{k_0}x + {\phi _0} + \text{δ} {\varphi _0}) + \\& + {\rm i}\sin (2{\text{π}}{k_0}x + {\phi _0} + \text{δ} {\varphi _0})].\end{split} $

此时可计算任一采样点对应的相位值:

$2{\text{π}}{k_0}{x_0}+{\varphi _0}+ \text{δ} {\varphi _0} ={\rm{ arctan}}\left[ {\frac{{\Im (I_{\rm{D}}^0)}}{{\Re (I_{\rm{D}}^0)}}} \right]$

式中, $\Im (I_{\rm{D}}^0)$和$\Re (I_{\rm{D}}^{\rm{0}})$分别为${x_0}$对应的$I_{\rm{D}}^0$的虚部和实部.
2 -->

2.3.原始干涉图信噪比与相位不确定度理论推导分析

在多普勒差分干涉仪成像过程中, 噪声主要来源于电路噪声(读出噪声、暗电流噪声)和散粒噪声中的光子噪声(属于高斯白噪声). 因此叠加噪声的干涉图${I{'}}$可以表示为

${I{'}} = I + {\varepsilon _N}, $

式中${\varepsilon _N}$表示电路读出噪声、暗电流噪声和光子噪声引起的原始干涉图强度的随机变化, 由随机噪声特性可知, 数学期望为0, 且相邻空间位置不相关. 假设原始干涉图以光子噪声为主要噪声源, 光子噪声方差与探测器光电转换后生成的信号对应电子数相等^[16]. 分别对(5)式中${I{'}}$、I进行傅里叶变换、窗函数单频提取、傅里叶逆变换后得$I_{\rm{c}}^{\rm{\prime}}$、${I_{\rm{c}}}$, 结合随机信号分析方法, 用(4)式求得$I_{\rm{c}}^{\rm{\prime}}$、${I_{\rm{c}}}$对应的反演相位${\rm{Phas}}{{\rm{e}}^{\prime}}$、${\rm{Phase}}$, 两相位差$\Delta {\rm{Phase}}$即为原始干涉图中随机噪声引入的相位噪声:

$\begin{split} \, & \Delta {\rm{Phase}} = \arctan \left[ {\frac{{\Im (I_{\rm{c}}^{\prime})}}{{\Re (I_{\rm{c}}^{\prime})}}} \right] - \arctan \left[ {\frac{{\Im ({I_{\rm{c}}})}}{{\Re ({I_{\rm{c}}})}}} \right] \\ =& \arctan \left[ {\frac{{\Im (I_{\rm{c}}^{\rm{\prime}})*\Re ({I_{\rm{c}}}) - \Im ({I_{\rm{c}}})*\Re (I_{\rm{c}}^{\rm{\prime}})}}{{\Re (I_{\rm{c}}^{\rm{\prime}})*\Re ({I_{\rm{c}}}) + \Im (I_{\rm{c}}^{\prime})*\Im ({I_{\rm{c}}})}}} \right].\end{split}$

由于

$\frac{{\Im (I_{\rm{c}}^{\rm{\prime}})*\Re ({I_{\rm{c}}}) - \Im ({I_{\rm{c}}})*\Re (I_{\rm{c}}^{\rm{\prime}})}}{{\Re (I_{\rm{c}}^{\rm{\prime}})*\Re ({I_{\rm{c}}}) + \Im (I_{\rm{c}}^{'})*\Im ({I_{\rm{c}}})}} < 1, $

因此, 用反正切函数的一阶泰勒展开式求解(6)式, 即:

$\Delta {\rm{Phase}} = \frac{{\Im (I_{\rm{c}}^{\rm{\prime}})*\Re ({I_{\rm{c}}}) - \Im ({I_{\rm{c}}})*\Re (I_{\rm{c}}^{\rm{\prime}})}}{{\Re (I_{\rm{c}}^{\rm{\prime}})*\Re ({I_{\rm{c}}}) + \Im (I_{\rm{c}}^{\prime})*\Im ({I_{\rm{c}}})}}.$

求$\Delta {\rm{Phase}}$标准差得到由原始干涉图光子噪声引入的相位不确定度:

${\rm{Phase}}U = \frac{{\sqrt {E\left[ {\Im {{({\varepsilon _{\rm{c}}})}^2}} \right]} }}{{\sqrt {E[|I_{\rm{c}}^{\rm{'}}{|^2}]} }}, $

式中: $E\left[ {\Im {{({\varepsilon _{\rm{c}}})}^2}} \right]$, $E[|I_{\rm{c}}^{\prime}{|^2}]$分别表示${\varepsilon _{\rm{c}}}$虚部和$I_{\rm{c}}^{\rm{'}}$平方的均值; ${\varepsilon _{\rm{c}}}$为复干涉图中的噪声, 且实部与虚部同分布^[6,17], 则有:

$\sqrt {E\left[ {\Im {{({\varepsilon _{\rm{c}}})}^2}} \right]} = {\varepsilon _N} \cdot \sqrt {\frac{n}{N}},$

式中, n为窗函数尺寸, N为干涉图采样数, 将(10)式代入(9)式中可得:

${\rm{Phase}}U = \frac{{\sqrt n \cdot {\varepsilon _N}}}{{\sqrt {N \cdot E[|I_{\rm{c}}^{\rm{\prime}}{|^2}]} }}.$

原始干涉图任一采样点信噪比表达式为

${(S/N)_N} = \frac{{{I_N}}}{{{\varepsilon _N}}}, $

式中${\varepsilon _N}$和${I_N}$分别对应任一采样点的光子噪声和原始干涉图信号. 将(12)式代入(11)式可得原始干涉图信噪比与相位不确定度之间的解析表达式:

${\rm{Phase}}U = \frac{{\sqrt n \cdot {I_N}}}{{(S/N) \cdot \sqrt {N \cdot E[|I_c^{'}{|^2}]} }}. $

为充分验证推导的理论模型的正确性, 本节保证仿真干涉图与实验室实测干涉图有完全相同的边界条件, 使用计算机仿真生成信噪比与实测信噪比相同的干涉图, 每组数据1000帧. 首先, 采用第2节中的相位反演模型, 去除原始干涉图的低频基线, 傅里叶变换用长度为5的矩形窗函数对反演谱单频提取, 经傅里叶逆变换后, 反演每帧干涉图相位, 取第512个采样点相位统计标准差作为每组干涉图的相位不确定度; 其次, 再用本文推导的(13)式, 根据每帧干涉图信噪比计算相位不确定度, 统计平均后得到每组数据对应的结果. 两种方法计算结果如图3所示, 其中右图为两方法计算结果残差, 平均残差为0.0239 mrad. 造成残差的原因主要是在理论推导中为了方便计算, 用反正切函数的泰勒一阶展开式对相位方差近似计算. 从图中可看出, 随着信噪比的提高, 相位不确定度与残差均减小, 两种方法的计算结果一致性较好.
图 3 信噪比与相位不确定度仿真结果及残差　(a)为相位反演模型与理论模型计算的相位不确定度; (b)为两种计算结果残差
Figure3. Simulation results of retrieved phase uncertainty and residual with different SNR: (a) Phase uncertainty of two models; (b) the differences between two results.

为了进一步验证本文推导的原始干涉图信噪比与相位不确定度理论模型正确性, 结合实验室自研的多普勒差分干涉仪(主要参数参见表1)搭建实验平台, 实验器材包括多普勒差分干涉仪、波长为632.8 nm的稳频激光器、分光镜、衰减片、准直镜, 如图4所示.

干涉仪	基础光程差	50 mm
	光谱分辨率	0.78 cm–1
探测器	量化位数	16 bit
	满阱电荷	100000 e–1
	读出噪声(5 MHz)	18 e–1
	暗电流噪声(–75 ℃)	0.0003 e–1/pixel/sec

表1多普勒差分干涉仪主要参数
Table1.Principle components.

图 4 实验平台
Figure4. Photograph of laboratory setup.

为了提高系统的稳定性, 降低相位热漂移导致的系统性偏差, 对多普勒差分干涉仪进行被动隔热, 并对干涉仪组件主动温控, 精度为0.1 ℃, 探测器工作温度为–75 ℃, 并且在采集记录数据之前先使整个系统工作数小时到达稳定状态, 图5为积分时间为0.5 s时的原始干涉图.
图 5 原始干涉图
Figure5. Interferogram.

实验中用两种方法分别采集多组不同信噪比干涉图数据, 第一种方法使用同一衰减片, 调整积分时间控制原始干涉图强度, 获得不同信噪比干涉图; 第二种方法使用相同的积分时间, 调整衰减片获得不同信噪比干涉图. 这样设计主要有以下两种原因: 首先, 验证理论模型对不同噪声类型是否具有普适性; 其次, 为了使得实测干涉图信噪比覆盖更广. 由文献[16,17]可知, 实测干涉图噪声主要包含光子噪声、暗电流噪声、读出噪声, 因此当改变积分时间时, 三种噪声均发生改变, 当调整衰减片强度时, 仅光子噪声发生改变. 使用第一种方法选择两组不同透过率的衰减片, 其中第一组衰减片, 积分时间起始为0.2 s, 间隔为0.2 s; 第二组衰减片, 积分时间起始为0.5 s, 间隔为0.5 s. 每组衰减片实验中各采集9组数据, 每组100帧. 由于实验中没有精确测量入射光强, 因此选择多帧平均方法计算实测干涉图信噪比: 将每组100帧干涉图平均得到平均干涉图, 每帧干涉图与其偏差作为每帧干涉图噪声, 进而计算干涉图信噪比.
在用相位反演模型和本文推导的理论模型计算之前, 对原始干涉图预处理(去除冲击噪声和暗电平). 与仿真计算方法相同, 首先用第2节中的方法计算每帧干涉图相位, 得出每组干涉图的相位不确定度, 整个计算方法与前文计算仿真数据的数据处理方法相同. 其次, 用本文推导的理论模型计算每帧干涉图的理论相位不确定度, 再对每组结果统计平均得到每组干涉图的相位不确定度. 图6分别为两组衰减片对应的不同信噪比干涉图用相位反演模型与本文推导的理论模型计算结果及残差. 图7为相同积分时间, 不同衰减片对应的实测数据计算结果及残差.
图 6 变积分时间实测数据计算结果及残差　(a)为相位反演模型与理论模型计算的相位不确定度; (b)为两种计算结果残差
Figure6. Experimental results of retrieved phase uncertainty and residual with same attenuation coefficient and different integration time: (a) Phase uncertainty of two models; (b) the differences between two results.

图 7 变衰减片实测数据计算结果及残差　(a)为相位反演模型与理论模型计算的相位不确定度; (b)为两种计算结果残差
Figure7. Experimental results of retrieved phase uncertainty and residual with different attenuation coefficients: (a) Phase uncertainty of two models; (b) the differences between two results

图6中实线圈, 第3组数据为积分时间0.2 s的干涉图, 平均信噪比为73.488倍, 反演模型计算的不确定度为1.2213 mrad, 理论模型计算结果为1.0552 mrad, 仿真生成干涉图信噪比为72.8292倍, 相位反演模型计算结果为1.0396 mrad, 理论模型计算结果为1.0075 mrad. 残差偏大的原因是由于实测干涉图相位漂移较大, 图8中黑色为第3组实测数据反演相位差, 红色为仿真数据相同信噪比反演相位差, 从图中可以明显看出, 实测数据反演相位差的方差大于仿真数据.
图 8 相同信噪比实测数据与仿真数据对比分析
Figure8. Comparison and analysis of measured and simulated data with the same SNR.

图6中虚线圈, 第4组数据为积分时间1.5 s的干涉图, 平均信噪比为74.286倍, 反演模型计算的不确定度为1.1006 mrad, 理论模型计算结果为0.8663 mrad, 仿真生成的干涉图信噪比为72.8292倍时, 相位反演模型计算结果为1.0396 mrad, 理论模型计算结果为1.0075 mrad. 分析发现, 整组实测数据的信噪比波动较大, 如图9所示, 由于波动较大, 因此在平均计算时出现偏差, 造成计算结果残差较大.
图 9 积分时间1.5 s实测干涉图信噪比
Figure9. SNR of 1.5 s integration time interferograms

分析图6结果可得, 由于积分时间相对较短, 因此在使用不同方法采集干涉图对最终的计算结果影响甚小. 由图6、图7结果可得, 本文推导的理论模型与相位反演模型计算得到的相位不确定度一致性较好. 图6中两组不同衰减片计算结果平均残差分别为0.0634 mrad, 图7计算结果平均残差为0.0573 mrad, 对应风速误差分别为0.0373 m/s, 0.0337 m/s.(以波长为632.8 nm的激光作为输入光源, 实验室自研的多普勒差分干涉仪基础光程差为50 mm, 1 m/s的风速变化量对应的相位变化量为1.7 mrad). 两种方法的计算结果的一致性表明本文推导的理论模型正确.
对比图6、图7可发现, 在相同条件下, 两种计算方法的残差主要有以下几点: 首先, 在实际测量中, 随着积分时间的加长, 使得整组干涉图的采集时间加长, 信噪比波动范围变大, 如图6中当信噪比大于200倍时, 积分时间为4 s和4.5 s的两组实测数据, 两种方法的计算结果残差变大; 其次, 实验中对多普勒差分干涉仪中的成像系统被动隔热, 并未主动温控, 可能会对实测干涉图造成相位漂移. 后续研究中可将成像系统放大率带来的相位漂移修正后, 进一步提高反演相位精度.

基于多普勒差分干涉方程, 结合相位反演模型, 对干涉图进行傅里叶变换、窗函数单频提取、逆傅里叶变换、相位计算、反正切函数的一阶泰勒展开, 分步解析多普勒差分干涉仪中相位不确定度与干涉图信噪比之间的定量关系. 分别用多组不同信噪比的仿真与实测数据验证本文推导的理论模型, 与相位反演模型计算结果的高度一致性验证了本文提出模型的正确性. 本文成果在工程设计中结合目标源辐射特性分析、光学设计参数及电子学参数能够通过确定原始干涉图信噪比, 进而得出相位不确定度, 实现对多普勒差分干涉仪全链路观测能力分析的目的. 除此之外基于观测数据的信噪比分析, 还能够实现对实测观测数据的误差限分析, 为数据后续的科学应用提供不确定度输入.
感谢中国科学院西安光学精密机械研究所冯玉涛研究员和傅頔博士的讨论.