03. 数学反演：拟合方程、优化求解与残差评测

本章核心目标：构建 DOAS 实际光谱反演的完整数学模型；推导线性普通最小二乘法（OLS）正规方程与参数协方差矩阵；剖析波长漂移非线性影响与 Levenberg-Marquardt (LM) 阻尼迭代求解算法；建立斜柱浓度（SCD）、拟合残差（RMS）及物理检出限（Detection Limit）的严密量化评测体系。

---

1. DOAS 反演的数学建模与统一拟合方程

在实际测量中，光谱仪采用离散像元阵列（如包含 $N=1024$ 或 $2048$ 个像素的线阵/面阵 CCD 探测器）记录入射光子。波长轴被离散化为一系列采样点 ${\lambda }_i$（$i=1,2,\dots ,N$）。

1.1 对数差分拟合模型的离散表达

通常我们采集一条待分析的测量光谱 $I({\lambda }_i)$，以及一条选定的参考光谱 $I_{\text{ref}}({\lambda }_i)$（例如正午无污染天顶太阳散射光谱，或无吸收灯光光谱）。

根据上一章的差分吸收理论，测量谱与参考谱分别满足：

$$\ln I({\lambda }_i)=\ln I_{0,\text{meas}}^{\prime }({\lambda }_i)-\sum _{j=1}^M{\sigma }_j^{\prime }({\lambda }_i)\cdot {\text{SCD}}_{j,\text{meas}}$$$$\ln I_{\text{ref}}({\lambda }_i)=\ln I_{0,\text{ref}}^{\prime }({\lambda }_i)-\sum _{j=1}^M{\sigma }_j^{\prime }({\lambda }_i)\cdot {\text{SCD}}_{j,\text{ref}}$$

两式相减，定义差分斜柱浓度（Differential Slant Column Density, dSCD）：

$$\mathrm{\Delta }{\text{SCD}}_j\equiv {\text{SCD}}_{j,\text{meas}}-{\text{SCD}}_{j,\text{ref}}$$

两谱的慢变连续背景差值 $\ln I_{0,\text{meas}}^{\prime }-\ln I_{0,\text{ref}}^{\prime }$ 依然是一个慢变光滑函数，统一用低阶多项式 $P_m({\lambda }_i)$ 表征。同时引入仪器杂散光/未完全消除的暗电平强度偏移修正项 $\text{Offset}({\lambda }_i)$ 与随机测量噪声 ${\epsilon }_i$，得到 DOAS 统一拟合方程：

DOAS 离散对数差分统一拟合方程

$$\ln \left(\frac{I({\lambda }_i)}{I_{\text{ref}}({\lambda }_i)}\right)=-\sum _{j=1}^M{\sigma }_j^{\prime }({\lambda }_i)\cdot \mathrm{\Delta }{\text{SCD}}_j-\sum _{k=0}^m{a}_k\cdot ({\lambda }_i-{\lambda }_0{)}^k+\text{Offset}({\lambda }_i)+{\epsilon }_i$$

---

2. 纯线性条件下的普通最小二乘法 (OLS)

当假设仪器波长绝对准确、不存在波长微小温漂，且不考虑非线性强度偏移时，上述拟合模型关于所有待求解参数（各气体的 $\mathrm{\Delta }{\text{SCD}}_j$ 及多项式系数 $a_k$）均处于严格的线性关系。

2.1 矩阵与向量形式

定义观测向量 $\mathbf{y}\in {\mathbb{R}}^N$：

$$y_i=\ln \left(\frac{I({\lambda }_i)}{I_{\text{ref}}({\lambda }_i)}\right),i=1,2,\dots ,N$$

待求解状态参数向量 $\mathbf{x}\in {\mathbb{R}}^P$（参数总维度 $P=M+m+1$）：

$$\mathbf{x}={[\mathrm{\Delta }{\text{SCD}}_1,\mathrm{\Delta }{\text{SCD}}_2,\dots ,\mathrm{\Delta }{\text{SCD}}_M,a_0,a_1,\dots ,a_m]}^T$$

构建雅可比/设计矩阵（Design Matrix）$\mathbf{K}\in {\mathbb{R}}^{N\times P}$，其第 $i$ 行向量定义为各基底函数在波长 ${\lambda }_i$ 处的响应：

$$\mathbf{K}=\left[\begin{array}{ccccccc}-{\sigma }_1^{\prime }({\lambda }_1)& \cdots & -{\sigma }_M^{\prime }({\lambda }_1)& -1& -({\lambda }_1-{\lambda }_0)& \cdots & -({\lambda }_1-{\lambda }_0{)}^m\\ -{\sigma }_1^{\prime }({\lambda }_2)& \cdots & -{\sigma }_M^{\prime }({\lambda }_2)& -1& -({\lambda }_2-{\lambda }_0)& \cdots & -({\lambda }_2-{\lambda }_0{)}^m\\ ⋮& \ddots & ⋮& ⋮& ⋮& \ddots & ⋮\\ -{\sigma }_1^{\prime }({\lambda }_N)& \cdots & -{\sigma }_M^{\prime }({\lambda }_N)& -1& -({\lambda }_N-{\lambda }_0)& \cdots & -({\lambda }_N-{\lambda }_0{)}^m\end{array}\right]$$

则整个离散拟合方程可写为标准线性矩阵形式：

$$\mathbf{y}=\mathbf{K}\mathbf{x}+\mathit{\epsilon }$$

2.2 正规方程推导与解析解

定义拟合目标函数为残差向量的二范数平方（无偏卡方统计量 ${\chi }^2$）：

$${\chi }^2(\mathbf{x})=\parallel \mathbf{y}-\mathbf{K}\mathbf{x}{\parallel }_2^2=(\mathbf{y}-\mathbf{K}\mathbf{x}{)}^T(\mathbf{y}-\mathbf{K}\mathbf{x})={\mathbf{y}}^T\mathbf{y}-2{\mathbf{x}}^T{\mathbf{K}}^T\mathbf{y}+{\mathbf{x}}^T{\mathbf{K}}^T\mathbf{K}\mathbf{x}$$

为了求使 ${\chi }^2$ 达到全局最小值的最优参数估计 $\hat{\mathbf{x}}$，对向量 $\mathbf{x}$ 计算一阶梯度并令其等于零向量：

$${\mathrm{\nabla }}_{\mathbf{x}}{\chi }^2=-2{\mathbf{K}}^T\mathbf{y}+2{\mathbf{K}}^T\mathbf{K}\mathbf{x}=\mathbf{0}$$

整理即得到经典的高斯正规方程（Normal Equations）：

$${\mathbf{K}}^T\mathbf{K}\mathbf{x}={\mathbf{K}}^T\mathbf{y}$$

当且仅当各分子的差分截面与多项式基底彼此线性无关时，信息矩阵 ${\mathbf{K}}^T\mathbf{K}\in {\mathbb{R}}^{P\times P}$ 为对称正定满秩矩阵，其唯一逆矩阵存在，得到最优参数的显式解析解：

线性最小二乘最优解析解

$$\hat{\mathbf{x}}=({\mathbf{K}}^T\mathbf{K}{)}^{-1}{\mathbf{K}}^T\mathbf{y}$$

2.3 参数误差估计与协方差矩阵

假设测量误差 $\mathit{\epsilon }$ 为零均值独立同分布高斯白噪声，其方差为 ${\sigma }_{\epsilon }^2$。根据线性统计误差传递定理，估计量 $\hat{\mathbf{x}}$ 的协方差矩阵 ${\mathbf{C}}_{\mathbf{x}}\in {\mathbb{R}}^{P\times P}$ 为：

$${\mathbf{C}}_{\mathbf{x}}=E[(\hat{\mathbf{x}}-\mathbf{x})(\hat{\mathbf{x}}-\mathbf{x}{)}^T]={\sigma }_{\epsilon }^2\cdot ({\mathbf{K}}^T\mathbf{K}{)}^{-1}$$

其中主对角线元素即为各参数拟合估计的方差：

• 第 $j$ 种气体的拟合标准差为：$\sigma (\mathrm{\Delta }{\text{SCD}}_j)=\sqrt{({\mathbf{C}}_{\mathbf{x}}{)}_{jj}}$；
• 非对角线元素 $({\mathbf{C}}_{\mathbf{x}}{)}_{jk}$ 反映了不同气体吸收截面之间的相关性干扰（Cross-Talk）。若两组分的差分峰位部分重叠，非对角元增大，会导致各自的独立反演不确定度急剧升高。

---

3. 非线性效应与 Levenberg-Marquardt (LM) 优化算法

在真实的户外或车载观测中，纯线性 OLS 假设往往会被打破。核心原因在于仪器光谱波长轴的微小物理漂移。

3.1 波长漂移（Shift）与展宽（Squeeze）

野外仪器内部温度波动（昼夜温差 ${15}^{\circ }\text{C}$）、外部震动会导致光栅基座或 CCD 靶面发生微米级热胀冷缩，引起像元对应的实际波长漂移：

$${\lambda }_i^{\prime }={\lambda }_i+\mathrm{\Delta }\lambda +s\cdot ({\lambda }_i-{\lambda }_0)$$

• $\mathrm{\Delta }\lambda$：平移参数（Shift），通常处于 $\pm 0.01\sim 0.1\text{ nm}$ 范围；
• $s$：拉伸/压缩参数（Squeeze），表征色散率的温度微变。

吸收截面关于波长漂移项满足：${\sigma }_j^{\prime }={\sigma }_j^{\prime }(\lambda +\mathrm{\Delta }\lambda )$。对其进行一阶泰勒展开：

$${\sigma }_j^{\prime }(\lambda +\mathrm{\Delta }\lambda )\approx {\sigma }_j^{\prime }(\lambda )+\frac{\partial {\sigma }_j^{\prime }}{\partial \lambda }\mathrm{\Delta }\lambda$$

由于参数 $\mathrm{\Delta }{\lambda }_j$ 与 ${\text{SCD}}_j$ 发生了乘积耦合（$\mathrm{\Delta }{\lambda }_j\cdot {\text{SCD}}_j$），模型蜕变为高度非线性拟合问题。

3.2 非线性最小二乘目标函数

设待拟合的总参数集合为 $\mathbf{p}\in {\mathbb{R}}^K$（既包含线性的 $\text{SCD}$，也包含非线性的 Shift 与 Squeeze 参数）。

模型在各波长处的预测函数记为 $f_i(\mathbf{p})=f({\lambda }_i,\mathbf{p})$。目标是求解最小化残差平方和的最优向量 $\hat{\mathbf{p}}$：

$$\underset{\mathbf{p}}{min}S(\mathbf{p})=\frac{1}{2}\sum _{i=1}^N{[y_i-f_i(\mathbf{p})]}^2=\frac{1}{2}\parallel \mathbf{r}(\mathbf{p}){\parallel }_2^2$$

其中残差向量 $\mathbf{r}(\mathbf{p})=\mathbf{y}-\mathbf{f}(\mathbf{p})$。定义雅可比矩阵 $\mathbf{J}(\mathbf{p})$：$J_{ij}=\frac{\partial f_i(\mathbf{p})}{\partial p_j}$。

3.3 Levenberg-Marquardt (LM) 算法的迭代推导与阻尼动力学

• 梯度下降法（Gradient Descent）：步长方向为负梯度 $-\mathbf{g}$。在远离极小值点时非常稳健，但在极值点附近收敛呈锯齿状，速度极慢；
• 高斯-牛顿法（Gauss-Newton）：迭代步长为 $\mathrm{\Delta }{\mathbf{p}}_{\text{GN}}=({\mathbf{J}}^T\mathbf{J}{)}^{-1}{\mathbf{J}}^T\mathbf{r}$。在极小值邻域享有二次收敛速度，但若初值较差导致 ${\mathbf{J}}^T\mathbf{J}$ 接近奇异时，会引发矩阵病态发散。

🧗‍♂️ 登山者的步履隐喻：LM 算法的阻尼艺术

想象一个在迷雾峡谷中寻找谷底最低点的登山者：

• 当阻尼因子 $\mu$ 极大时：迷雾极其浓重，地面情况完全不明，登山者只能极其小心翼翼地沿着脚下坡度最陡的负梯度方向摸索小步前进（梯度下降）——虽然速度缓慢，但绝不会失足坠下悬崖；
• 当阻尼因子 $\mu$ 极小时：迷雾散去，谷底清晰可见，登山者利用抛物线二次近似直接飞奔冲向碗底极小值点（高斯-牛顿）——几步之内超速收敛；
• LM 智能自适应调度：走完一步后，若海拔下降了（残差减小），就建立信心，大幅调小 $\mu$ 加速冲刺；若不小心踩空海拔上升了（残差变大），就撤回该步，果断调大 $\mu$ 切换回谨慎摸索！

$$({\mathbf{J}}^T\mathbf{J}+\mu \cdot \mathbf{D})\mathrm{\Delta }\mathbf{p}={\mathbf{J}}^T\mathbf{r}(\mathbf{p})$$

LM 算法在 DOAS 光谱拟合软件（如 QDOAS、DOASIS）中作为标准反演内核，通常在 $5\sim 15$ 次迭代内即可高度稳健地收敛至全局最优解。

---

4. 斜柱浓度 (SCD) 物理单位与实用换算

反演算法直接输出的物理量是沿整个视线光程积分的粒子数密度总和，其标准 SI 衍生单位为：

$$\left[\text{SCD}\right]=\text{molecules}\cdot {\text{cm}}^{-2}\text{或}{\text{molec/cm}}^2$$

4.1 积分体积分数（$\text{ppm}\cdot \text{m}$）的工程换算

在主动式长光程 DOAS（LP-DOAS）或工业管道排放监测中，工程规范通常习惯使用体积分数与几何距离的乘积，如 $\text{ppm}\cdot \text{m}$（每百万分之一体积分数在 1 米光程内的积）或 $\text{mg}\cdot {\text{m}}^{-3}$。

根据理想气体状态方程：

$$p\cdot V=N_{\text{tot}}\cdot k_B\cdot T⟹n_{\text{air}}=\frac{N_{\text{tot}}}{V}=\frac{p}{k_B\cdot T}$$

其中玻尔兹曼常数 $k_B=1.380649\times {10}^{-23}\text{ J/K}$。在标准状态（STP: $T_0=273.15\text{ K}$, $p_0=101325\text{ Pa}$）下，标准空气分子数密度（洛施密特常数 $N_L$）为：

$$n_0=\frac{101325}{1.38065\times {10}^{-23}\times 273.15}\approx 2.6868\times {10}^{19}{\text{ molec/cm}}^3$$

工程实用换算等式 (标准状态 STP)

$$1\text{ ppm}\cdot \text{m}\approx 2.69\times {10}^{15}{\text{ molec/cm}}^2$$$$1\text{ ppb}\cdot \text{km}\approx 2.69\times {10}^{15}{\text{ molec/cm}}^2$$

若环境温度为室温 $T=298.15\text{ K}$（${25}^{\circ }\text{C}$），换算常数修正为：

$$1\text{ ppm}\cdot \text{m}\approx 2.46\times {10}^{15}{\text{ molec/cm}}^2(298.15\text{ K},1\text{ atm})$$

---

5. 反演质量评估体系：残差、RMS 与检出限

一次 DOAS 反演得到的一组浓度数值是否可信？必须依赖严谨的拟合统计学指标进行全自动校验。

5.1 残差谱（Residual Spectrum）与白噪声检验

在最优参数 $\hat{\mathbf{p}}$ 下，从实测对数光谱中减去所有拟合成功的成分，剩下的未被模型解释的部分构成残差向量（Residuals）：

$$\text{res}({\lambda }_i)=y_i-f({\lambda }_i,\hat{\mathbf{p}})$$

理想的反演残差必须表现为严格不相关的纯高斯白噪声：

• 均值检验：$\overline{\text{res}}\approx 0$；
• 自相关分析（Durbin-Watson 检验）：若残差谱中出现明显的周期性正弦波动或结构性尖峰，强烈表明：
  1. 拟合窗口内遗漏了某未知吸收气体（例如反演 $\text{HCHO}$ 窗口时未引入 $\text{BrO}$ 或 $\text{HONO}$ 截面）；
  2. 仪器谱线函数（Slit Function）存在未被校准的非对称彗形像差；
  3. 太阳光谱卷积存在微小波长错位或 Ring 效应未完全扣除。

5.2 拟合残差均方根（RMS）

残差均方根（Root Mean Square, RMS）是表征整个拟合窗口整体精度的单一关键指标：

$$\text{RMS}=\sqrt{\frac{1}{N-K}\sum _{i=1}^N{(\text{res}({\lambda }_i))}^2}$$

其中 $N$ 为有效像元采样点数，$K$ 为自由参数总数，$N-K$ 为系统的自由度。

RMS 数量级	反演质量级别	典型应用与可信度
$<3\times {10}^{-4}$	顶尖科研级	完美消除光学杂散光与拉曼散射，可精确检出超痕量 $\text{IO},\text{BrO},\text{HONO}$
$3\times {10}^{-4}\sim 1\times {10}^{-3}$	优良标准级	典型常规 MAX-DOAS 观测水平，${\text{NO}}_2$ 与 ${\text{SO}}_2$ 结果高度精确可靠
$1\times {10}^{-3}\sim 3\times {10}^{-3}$	临界可用级	存在强气溶胶消光、浓云干扰或轻微光学元件温漂，弱吸收气体误差增大
$>5\times {10}^{-3}$	严重失真级	拟合发散或发生严重波长漂移，光谱数据应标记为无效并剔除

5.3 仪器物理检出限（Detection Limit, DL）评估

分子的光学检出限受制于光谱系统的综合噪声水平。根据 IUPAC 通用准则，能够被可信判别的最小差分光学厚度 ${\tau }_{min}^{\prime }$ 必须达到残差噪声 RMS 的 $k$ 倍（通常取 $k=2$ 对应 $95.4\mathrm{\%}$ 置信区间）：

$${\tau }_{min}^{\prime }=k\cdot \text{RMS}$$

对于吸收特征在选定窗口内的气体 $j$，定义其特征差分峰-谷幅值（Differential Peak-to-Peak Amplitude）：

$$\mathrm{\Delta }{\sigma }_j^{\prime }\equiv \underset{\lambda }{max}({\sigma }_j^{\prime }(\lambda ))-\underset{\lambda }{min}({\sigma }_j^{\prime }(\lambda ))$$

则该气体在当前光谱质量下的**最小可检出斜柱浓度（Detection Limit, ${\text{DL}}_j$）**严格由下式给出：

物理检出限（Detection Limit）解析式

$${\text{DL}}_j=\frac{k\cdot \text{RMS}}{\mathrm{\Delta }{\sigma }_j^{\prime }}$$

若考虑窗口内多个独立差分吸收峰的统计平均降噪效应：

$${\text{DL}}_j\approx \frac{k\cdot \text{RMS}}{\mathrm{\Delta }{\sigma }_j^{\prime }\cdot \sqrt{N_{\text{peaks}}}}$$

实例计算：${\text{NO}}_2$ 的典型检出限

• 在 $430-450\text{ nm}$ 拟合窗口中，${\text{NO}}_2$ 的差分截面峰谷差 $\mathrm{\Delta }{\sigma }_{{\text{NO}}_2}^{\prime }\approx 6.0\times {10}^{-19}{\text{ cm}}^2/\text{molecule}$；

• 设仪器优化良好，测量光谱的拟合残差 $\text{RMS}=5.0\times {10}^{-4}$，取置信系数 $k=2$：

  $${\text{DL}}_{{\text{NO}}_2}=\frac{2\times (5.0\times {10}^{-4})}{6.0\times {10}^{-19}{\text{ cm}}^2/\text{molec}}\approx 1.67\times {10}^{15}{\text{ molec/cm}}^2$$

• 若在光程 $L=5\text{ km}$ 的水平长光程观测中：

  $${\text{DL}}_{\text{conc}}=\frac{1.67\times {10}^{15}{\text{ molec/cm}}^2}{2.46\times {10}^{15}\text{ molec}\cdot {\text{cm}}^{-2}/(\text{ppm}\cdot \text{m})\times 5000\text{ m}}\approx 0.000136\text{ ppm}=0.14\text{ ppb}!$$

这以严谨的数学推导证明了：DOAS 凭借极小的残差控制与差分指纹特征，具备探测亚 $\text{ppb}$ 甚至数十 $\text{ppt}$ 级极限超痕量气体的超凡灵敏度。

---

6. 本章总结与下一章预告

本章推导了 DOAS 的矩阵拟合正规方程与 LM 优化算法，构建了从原始对数比值到斜柱浓度、拟合残差与检出限的完整数学反演闭环。

然而，所有精妙算法的前提是拥有一台能够采集高信噪比、低杂散光且波长稳定的物理光谱仪。地面如何利用凹面镜与衍射光栅构建紧凑型色散系统？探测器如何通过半导体热电深度制冷压制暗电流？MAX-DOAS 又如何利用空间多轴仰角几何将低空污染信号放大数十倍？

请进入下一章：04. 典型仪器光路与 MAX-DOAS 几何。