基于权重变分模型的地基IDOAS条带噪声去除

奚 亮， 司福祺， 江 宇， 周海金， 邱晓晗， 常 振

1. 中国科学院合肥物质科学研究院， 安徽光学精密机械研究所， 安徽 合肥 230031

2. 中国科学技术大学， 安徽 合肥 230026

引 言

成像差分吸收光谱技术(imaging differential optical absorption spectroscopy， IDOAS)是成像光谱技术和差分吸收光谱技术的结合， 作为一种测量大气痕量气体(SO2， NO2， O3等)的有效方法， 能够获取痕量气体二维浓度分布信息， 近年来在地基、 机载和星载平台均有较大的发展[1-3]。 受到传感器不均匀性和外界环境变化等因素的影响， IDOAS仪器容易出现条带噪声。 作为一种在成像光谱技术中常见的数据降质问题， 条带噪声会覆盖目标物的空间分布特征并产生相应的伪结构， 影响后续的反演分析和信息提取， 最终降低了数据的可用性。

为解决条带噪声问题， 国内外已有学者提出了有效的处理算法。 在星载方面， Boersma等利用洁净均匀区域的平均值作为真值， 计算不同探测像元和真值的差值来消除OMI NO2柱浓度中的条带噪声[4]； Cheng等为去除EMI NO2产品中的条带噪声， 使用了二维傅里叶变换加频率滤波的方法[5]； 杨太平等将SICATRAN模型的模拟值作为参考， 实现了EMI O4产品穿轨方向的条带校正[6]； Geffen等在TROPOMI NO2算法中， 将太平洋赤道地区的斜柱浓度与CTM/DA模型计算值的差异作为去条带的校正系数[7]。 在机载方面， Popp等假定APEX观测结果在穿轨方向上平滑变化， 对沿轨方向的均值进行多项式拟合后， 最后从NO2柱浓度中减去残差结构以去除条带[8]； Nowlan等在GeoTASO NO2柱浓度的去条带算法中， 利用洁净海湾区域的观测值与CMAQ模拟值的偏差作为校正量[9]。

以上算法均有良好的去条带效果， 然而在应用到地基IDOAS中却存在一些不适用的问题， 如均匀参考区域难以获得， 地面遮挡物造成的异常区域较多， 条带噪声变化较大等。 本文将图像处理中常用的变分模型应用到地基IDOAS仪器的条带去除中[10]， 在考虑了条带噪声的分布特性和地基仪器的的遮挡问题后， 提出了基于权重变分模型的去条带算法。 经过模拟实验和外场实验的测试， 该算法能够有效去除地基IDOAS仪器中气体浓度分布的条带噪声问题。

1 基本原理

1.1 地基IDOAS原理

如图1所示， 地基IDOAS通过二维转台的水平旋转实现摆扫成像， 其中箭头表示扫描的方向，α表示水平方向上扫描的角度，β表示垂直方向上的视场角。 仪器采集天空中的太阳散射光， 并利用DOAS方法处理光谱得到痕量气体沿光路的积分浓度。 DOAS技术的原理是Lambert-Beer定律[11]， 利用痕量气体的窄带吸收特性， 将其与Mie散射、 Rayleigh散射、 气溶胶作用等宽带信号相分离。 图2为一组外场实验的观测结果， 实验中IDOAS仪器水平扫描360°， 垂直方向覆盖30°。 图2(a)为360 nm处探测器接收到的数字值(digital number， DN)， 图2(b)为地面障碍物分布， 图2(c)和(d)为NO2和SO2的观测结果， 图2(e)和(g)为NO2观测结果的垂直和水平梯度图， 图2(f)和(h)为SO2观测结果的垂直和水平梯度图。

图1 地基IDOAS仪器示意图

1.2 条带噪声特性

从图2(c)和(d)可以看出， NO2和SO2观测结果有明显的条带噪声， 需要进行去条带处理。 与随机噪声相比， 条带噪声具有明显的空间分布特性， 如图2(e)—(h)所示， NO2和SO2观测结果在垂直方向的梯度远远大于水平方向的梯度。 同时条带噪声也具有稀疏性， 即条带噪声没有覆盖整幅图像， 且在很多区域接近于零。 此外， 对比图4中另外两组的观测结果， 可发现地基IDOAS中条带噪声不是固定的偏差值， 条带噪声的变化性较大。

图2(b)为观测点附近地表障碍物的分布， 障碍物的遮挡会造成光谱结构的改变， 使图2(c)和(d)中DOAS反演结果出现异常， 并使条带噪声的结构遭到破坏。 地基IDOAS测量中受到障碍物遮挡和气体空间分布的影响， 常见去条带算法中所需的均匀参考区域可能难以获得。 本文采用基于权重变分模型的去条带算法， 该算法利用条带噪声本身的分布特性， 不需要均匀的参考区域， 并使用权重矩阵解决了障碍物遮挡的问题。

1.3 权重变分去条带算法

假设条带噪声为加性噪声， 去条带算法的各向异性变分模型为

λ2‖S‖0+λ3‖hS‖0}

(1)

式(1)中：S为条带噪声；Y为观测到的实际图像；h为水平梯度算子，v为竖直梯度算子； ‖·‖1为l1范数， 能够在最小化的过程中接受较大的突变， 适合刻画条带边界。 ‖·‖0为l0范数，l0范数能控制S的估算过程， 避免无条带区域信息的丢失； 式(1)中右边第一项表示条带噪声在水平方向的连续性， 第二项表示图像在垂直方向上的连续性， 第三和第四项表示条带噪声及其水平梯度的稀疏性；λ1，λ2和λ3为平衡方向性和稀疏性的参数。

式(1)利用图像的全局特性来估算条带噪声， 但在地基测量中由于异常区域的存在， 使得基于全局特性的算法不适用。 该异常区域由地表障碍物遮挡造成， 而在遮挡区域探测器的DN值通常要小于天空背景， 因此可通过阈值分割实现天空背景和地表障碍物的区分。 本文使用分块自适应阈值分割算法， 可得到表征遮挡区域的权重矩阵W。 如图2(b)所示，W在遮挡区域值为0， 在正常区域值为1。 在遮挡区域， 图像的方向特性遭到破坏， 因此可在式(1)中的前两项添加权重因子， 使其在遮挡区域的贡献为零。 结合以上分析， 权重变分模型可表示为

图2 地基IDOAS一组观测结果

λ2‖D‖0+λ3‖H‖0},

s.t.H=hS

(2)

V=v(Y-S)

D=S

式(2)中：W为权重矩阵； ∘为元素乘法；H，V和D为引入的变量， 使式(2)转化为有约束的优化问题。 为求解式(2)， 本文采用交替方向乘子算法(alternating direction method of multipliers， ADMM)[12]， 其对应的增广拉格朗日展开为

L(H,V,D,S,p1,p2,p3)=‖W∘H‖1+λ3‖H‖0+

(3)

式(3)中：p1，p2和p3为拉格朗日乘子；ρ1，ρ2和ρ3为惩罚参数。 原问题被分解为H，V，D和S四个相对简单的子问题， 并在迭代中交替求解

(4)

(5)

(6)

(7)

式中cshrink， softshrink和hardshrink分别为复合阈值收缩、 软阈值收缩和硬阈值收缩算子[13]， F和F-1分别表示二维快速傅里叶变换和逆变换； ∘为元素乘法； *为复共轭算子。

在每轮迭代的最后， 需要按以下公式对拉格朗日乘子进行更新

(8)

(9)

(10)

当迭代次数大于上限Nmax或S的变化量‖S(k+1)-S(k)‖/‖S(k)‖小于阈值ε时迭代停止。 从上述模型中估计出S后， 用实际图像Y减去S即可获得所需的去条带图像。

2 实验部分

为检验去条带算法的性能， 利用四种不同类型的模拟条带噪声进行了测试。 实验中理想图像为二维Savitzky-Golay滤波后模糊但无条带的实测数据， 图3展示了受不同模拟条带影响后的理想图像以及去条带结果， 其中图3(a)和图3(c)为位置周期、 强度随机的条带， 图3(a)中的条带较稀疏， 图3(c)中的条带较稠密； 图3(e)和图3(g)为位置周期、 强度随机的条带加上位置随机、 强度随机的条带， 图3(e)中的位置随机条带为部分条带， 即条带没有贯穿整行， 图3(g)中的位置随机条带为宽条跌， 即带宽度超过两行。 图3(a)—(h)中的底部灰色像元表示由遮挡所致的异常区域。

在模拟实验中，λ1，λ2和λ3分别设置为0.2， 0.001和0.2，ρ1，ρ2和ρ3设置为100倍的λ1， 最终去条带结果可见图3(b)， (d)， (f)和(h)。 从目视效果来看， 对四种不同类型的条带噪声， 权重变分算法都展现了良好的去条带性能。 采用四种常用的全参考评价指标对算法进行定量评价， 包括图像改善因子(improvement factor, IF)、 峰值信噪比(peak signal-to-noise ratio, PSNR)、 结构相似函数(structural similarity index, SSIM)和平均绝对误差(mean absolute error, MAE)[14]。 其中IF反映算法的降噪能力， 而PSNR， SSIM和MAE评价算法的保真能力。 一般来说， IF， PSNR的值越大， SSIM的值越接近于1， MAE的值越接近于0， 表明该算法的去条带效果越好。 表1列出了模拟实验中四种全参考指标计算结果， 不难发现权重变分算法都有良好的表现。 目视效果和评价指标均说明， 权重变分算法能够有效去除常见类型的条带噪声。

表1 模拟实验的评价结果

图3 不同模拟条带噪声和去条带结果

3 结果与讨论

地基IDOAS于2018年夏季在四川乐山进行了外场实验， 图2(b)中从左到右的障碍物分别为信号塔、 山丘、 工厂烟囱以及其他建筑物。 工业园区位于仪器的西南方位， 分布着多个污染气体的排放源。 外场实验中仪器的扫描间隔为1°， 最终每组观测图像的像素大小为360×48。 仪器的积分时间设置为500 ms， 每组全景扫描的总工作时间约为15 min。 在每组光谱采集完成后， 利用DOAS技术分别对NO2和SO2气体的斜柱浓度进行了反演， 具体的参数设置见表2。 图4展示了实验中两组观测的结果， 其中图4(a)和图4(c)为12:59开始采集的NO2和SO2浓度分布， 图4(e)和图4(g)为13:14开始采集的NO2和SO2浓度分布。 比较图4(a)， (c)， (e)和(g)可以发现， 条带噪声具有很大的变化性： 对于不同时间的观测结果， 条带噪声的空间分布和强度均不同； 对于不同反演气体， 条带噪声的空间分布和强度也不同。

表2 NO2和SO2 DOAS反演中的参数设置

图4 乐山实验中两组NO2和SO2的实际观测和去条带结果

在权重变分模型中，λ1，λ2和λ3用于平衡各个约束条件， 在实际应用中需根据图像的具体情况加以调整。 一般来说， 条带噪声的强度越大，λ1应设置越大； 条带噪声越稀疏， 需使用更大的λ2； 条带噪声越均匀， 适合用更大的λ3。 在外场实验中为达到最佳的去条带效果， 经多次测试后λ1，λ2和λ3分别设置为0.4， 0.001和0.2，ρ1，ρ2和ρ3设置为100倍的λ1。 去条带后的结果可见图4(b)， (d)， (f)和(h)， 总体上去条带算法取得了令人满意的效果。 如图4(a)方框所示的多个浓度高值区域， 在去条带之后保留了各自的空间分布信息； 图4(c)方框所示的观测区域有大量的随机噪声， 在去条带之后没有变得模糊。 结果表明， 权重变分算法能够有效去除条带噪声， 能够保留污染气体分布信息， 且未出现过度平滑的情况。

图5展示了去条带前后的权重列均值曲线， 其中虚线为去条带前的权重列均值， 实线为去条带后的权重列均值。 权重列均值表示在计算图像X(大小为m×n)的列均值c时， 需考虑权重因子W，c的第i行分量为

图5 去条带前后的NO2和SO2权重列均值对比

(11)

在低仰角区域会出现完全遮挡的情况， 此时定义分量ci为0。 从图5可以看出， 权重变分算法可以获得平滑的权重列均值曲线， 并且能保留原始曲线的变化趋势。

由于实测数据没有理想图像作为参考， 为评估外场实验中去条带算法的有效性， 采用Borsdorff等提出的方法对结果进行定量评价， 参量γ定义为[15]

(12)

表3 去条带前后的γ值

4 结 论

针对地基IDOAS仪器中存在的条带噪声问题， 提出了基于权重变分模型的去条带算法。 该算法基于条带噪声的方向性和稀疏性， 不需要均匀的参考区域， 并充分考虑了地基成像中障碍物遮挡的问题。 在测试算法性能的模拟实验中， 对四组不同类型条带的降噪效果证明了该算法具有良好的有效性和保真性。 在2018年夏季的外场实验中， 地基IDOAS反演结果为痕量气体斜柱浓度的二维分布。 运用该算法后， 多组NO2和SO2浓度分布中的条带噪声得到了有效去除。 结果表明， 该算法适用于多种气体以及变化条带的情况， 提高了地基IDOAS观测数据的可用性。