基于權(quán)重變分模型的地基IDOAS條帶噪聲去除

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

1. 中國科學(xué)院合肥物質(zhì)科學(xué)研究院， 安徽光學(xué)精密機(jī)械研究所， 安徽 合肥 230031

2. 中國科學(xué)技術(shù)大學(xué)， 安徽 合肥 230026

引 言

成像差分吸收光譜技術(shù)(imaging differential optical absorption spectroscopy， IDOAS)是成像光譜技術(shù)和差分吸收光譜技術(shù)的結(jié)合， 作為一種測量大氣痕量氣體(SO2， NO2， O3等)的有效方法， 能夠獲取痕量氣體二維濃度分布信息， 近年來在地基、 機(jī)載和星載平臺(tái)均有較大的發(fā)展[1-3]。 受到傳感器不均勻性和外界環(huán)境變化等因素的影響， IDOAS儀器容易出現(xiàn)條帶噪聲。 作為一種在成像光譜技術(shù)中常見的數(shù)據(jù)降質(zhì)問題， 條帶噪聲會(huì)覆蓋目標(biāo)物的空間分布特征并產(chǎn)生相應(yīng)的偽結(jié)構(gòu)， 影響后續(xù)的反演分析和信息提取， 最終降低了數(shù)據(jù)的可用性。

為解決條帶噪聲問題， 國內(nèi)外已有學(xué)者提出了有效的處理算法。 在星載方面， Boersma等利用潔凈均勻區(qū)域的平均值作為真值， 計(jì)算不同探測像元和真值的差值來消除OMI NO2柱濃度中的條帶噪聲[4]； Cheng等為去除EMI NO2產(chǎn)品中的條帶噪聲， 使用了二維傅里葉變換加頻率濾波的方法[5]； 楊太平等將SICATRAN模型的模擬值作為參考， 實(shí)現(xiàn)了EMI O4產(chǎn)品穿軌方向的條帶校正[6]； Geffen等在TROPOMI NO2算法中， 將太平洋赤道地區(qū)的斜柱濃度與CTM/DA模型計(jì)算值的差異作為去條帶的校正系數(shù)[7]。 在機(jī)載方面， Popp等假定APEX觀測結(jié)果在穿軌方向上平滑變化， 對(duì)沿軌方向的均值進(jìn)行多項(xiàng)式擬合后， 最后從NO2柱濃度中減去殘差結(jié)構(gòu)以去除條帶[8]； Nowlan等在GeoTASO NO2柱濃度的去條帶算法中， 利用潔凈海灣區(qū)域的觀測值與CMAQ模擬值的偏差作為校正量[9]。

以上算法均有良好的去條帶效果， 然而在應(yīng)用到地基IDOAS中卻存在一些不適用的問題， 如均勻參考區(qū)域難以獲得， 地面遮擋物造成的異常區(qū)域較多， 條帶噪聲變化較大等。 本文將圖像處理中常用的變分模型應(yīng)用到地基IDOAS儀器的條帶去除中[10]， 在考慮了條帶噪聲的分布特性和地基儀器的的遮擋問題后， 提出了基于權(quán)重變分模型的去條帶算法。 經(jīng)過模擬實(shí)驗(yàn)和外場實(shí)驗(yàn)的測試， 該算法能夠有效去除地基IDOAS儀器中氣體濃度分布的條帶噪聲問題。

1 基本原理

1.1 地基IDOAS原理

如圖1所示， 地基IDOAS通過二維轉(zhuǎn)臺(tái)的水平旋轉(zhuǎn)實(shí)現(xiàn)擺掃成像， 其中箭頭表示掃描的方向，α表示水平方向上掃描的角度，β表示垂直方向上的視場角。 儀器采集天空中的太陽散射光， 并利用DOAS方法處理光譜得到痕量氣體沿光路的積分濃度。 DOAS技術(shù)的原理是Lambert-Beer定律[11]， 利用痕量氣體的窄帶吸收特性， 將其與Mie散射、 Rayleigh散射、 氣溶膠作用等寬帶信號(hào)相分離。 圖2為一組外場實(shí)驗(yàn)的觀測結(jié)果， 實(shí)驗(yàn)中IDOAS儀器水平掃描360°， 垂直方向覆蓋30°。 圖2(a)為360 nm處探測器接收到的數(shù)字值(digital number， DN)， 圖2(b)為地面障礙物分布， 圖2(c)和(d)為NO2和SO2的觀測結(jié)果， 圖2(e)和(g)為NO2觀測結(jié)果的垂直和水平梯度圖， 圖2(f)和(h)為SO2觀測結(jié)果的垂直和水平梯度圖。

圖1 地基IDOAS儀器示意圖

1.2 條帶噪聲特性

從圖2(c)和(d)可以看出， NO2和SO2觀測結(jié)果有明顯的條帶噪聲， 需要進(jìn)行去條帶處理。 與隨機(jī)噪聲相比， 條帶噪聲具有明顯的空間分布特性， 如圖2(e)—(h)所示， NO2和SO2觀測結(jié)果在垂直方向的梯度遠(yuǎn)遠(yuǎn)大于水平方向的梯度。 同時(shí)條帶噪聲也具有稀疏性， 即條帶噪聲沒有覆蓋整幅圖像， 且在很多區(qū)域接近于零。 此外， 對(duì)比圖4中另外兩組的觀測結(jié)果， 可發(fā)現(xiàn)地基IDOAS中條帶噪聲不是固定的偏差值， 條帶噪聲的變化性較大。

圖2(b)為觀測點(diǎn)附近地表障礙物的分布， 障礙物的遮擋會(huì)造成光譜結(jié)構(gòu)的改變， 使圖2(c)和(d)中DOAS反演結(jié)果出現(xiàn)異常， 并使條帶噪聲的結(jié)構(gòu)遭到破壞。 地基IDOAS測量中受到障礙物遮擋和氣體空間分布的影響， 常見去條帶算法中所需的均勻參考區(qū)域可能難以獲得。 本文采用基于權(quán)重變分模型的去條帶算法， 該算法利用條帶噪聲本身的分布特性， 不需要均勻的參考區(qū)域， 并使用權(quán)重矩陣解決了障礙物遮擋的問題。

1.3 權(quán)重變分去條帶算法

假設(shè)條帶噪聲為加性噪聲， 去條帶算法的各向異性變分模型為

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

(1)

式(1)中：S為條帶噪聲；Y為觀測到的實(shí)際圖像；h為水平梯度算子，v為豎直梯度算子； ‖·‖1為l1范數(shù)， 能夠在最小化的過程中接受較大的突變， 適合刻畫條帶邊界。 ‖·‖0為l0范數(shù)，l0范數(shù)能控制S的估算過程， 避免無條帶區(qū)域信息的丟失； 式(1)中右邊第一項(xiàng)表示條帶噪聲在水平方向的連續(xù)性， 第二項(xiàng)表示圖像在垂直方向上的連續(xù)性， 第三和第四項(xiàng)表示條帶噪聲及其水平梯度的稀疏性；λ1，λ2和λ3為平衡方向性和稀疏性的參數(shù)。

式(1)利用圖像的全局特性來估算條帶噪聲， 但在地基測量中由于異常區(qū)域的存在， 使得基于全局特性的算法不適用。 該異常區(qū)域由地表障礙物遮擋造成， 而在遮擋區(qū)域探測器的DN值通常要小于天空背景， 因此可通過閾值分割實(shí)現(xiàn)天空背景和地表障礙物的區(qū)分。 本文使用分塊自適應(yīng)閾值分割算法， 可得到表征遮擋區(qū)域的權(quán)重矩陣W。 如圖2(b)所示，W在遮擋區(qū)域值為0， 在正常區(qū)域值為1。 在遮擋區(qū)域， 圖像的方向特性遭到破壞， 因此可在式(1)中的前兩項(xiàng)添加權(quán)重因子， 使其在遮擋區(qū)域的貢獻(xiàn)為零。 結(jié)合以上分析， 權(quán)重變分模型可表示為

圖2 地基IDOAS一組觀測結(jié)果

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

s.t.H=hS

(2)

V=v(Y-S)

D=S

式(2)中：W為權(quán)重矩陣； °為元素乘法；H，V和D為引入的變量， 使式(2)轉(zhuǎn)化為有約束的優(yōu)化問題。 為求解式(2)， 本文采用交替方向乘子算法(alternating direction method of multipliers， ADMM)[12]， 其對(duì)應(yīng)的增廣拉格朗日展開為

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

(3)

式(3)中：p1，p2和p3為拉格朗日乘子；ρ1，ρ2和ρ3為懲罰參數(shù)。 原問題被分解為H，V，D和S四個(gè)相對(duì)簡單的子問題， 并在迭代中交替求解

(4)

(5)

(6)

(7)

式中cshrink， softshrink和hardshrink分別為復(fù)合閾值收縮、 軟閾值收縮和硬閾值收縮算子[13]， F和F-1分別表示二維快速傅里葉變換和逆變換； °為元素乘法； *為復(fù)共軛算子。

在每輪迭代的最后， 需要按以下公式對(duì)拉格朗日乘子進(jìn)行更新

(8)

(9)

(10)

當(dāng)?shù)螖?shù)大于上限Nmax或S的變化量‖S(k+1)-S(k)‖/‖S(k)‖小于閾值ε時(shí)迭代停止。 從上述模型中估計(jì)出S后， 用實(shí)際圖像Y減去S即可獲得所需的去條帶圖像。

2 實(shí)驗(yàn)部分

為檢驗(yàn)去條帶算法的性能， 利用四種不同類型的模擬條帶噪聲進(jìn)行了測試。 實(shí)驗(yàn)中理想圖像為二維Savitzky-Golay濾波后模糊但無條帶的實(shí)測數(shù)據(jù)， 圖3展示了受不同模擬條帶影響后的理想圖像以及去條帶結(jié)果， 其中圖3(a)和圖3(c)為位置周期、 強(qiáng)度隨機(jī)的條帶， 圖3(a)中的條帶較稀疏， 圖3(c)中的條帶較稠密； 圖3(e)和圖3(g)為位置周期、 強(qiáng)度隨機(jī)的條帶加上位置隨機(jī)、 強(qiáng)度隨機(jī)的條帶， 圖3(e)中的位置隨機(jī)條帶為部分條帶， 即條帶沒有貫穿整行， 圖3(g)中的位置隨機(jī)條帶為寬條跌， 即帶寬度超過兩行。 圖3(a)—(h)中的底部灰色像元表示由遮擋所致的異常區(qū)域。

在模擬實(shí)驗(yàn)中，λ1，λ2和λ3分別設(shè)置為0.2， 0.001和0.2，ρ1，ρ2和ρ3設(shè)置為100倍的λ1， 最終去條帶結(jié)果可見圖3(b)， (d)， (f)和(h)。 從目視效果來看， 對(duì)四種不同類型的條帶噪聲， 權(quán)重變分算法都展現(xiàn)了良好的去條帶性能。 采用四種常用的全參考評(píng)價(jià)指標(biāo)對(duì)算法進(jìn)行定量評(píng)價(jià)， 包括圖像改善因子(improvement factor, IF)、 峰值信噪比(peak signal-to-noise ratio, PSNR)、 結(jié)構(gòu)相似函數(shù)(structural similarity index, SSIM)和平均絕對(duì)誤差(mean absolute error, MAE)[14]。 其中IF反映算法的降噪能力， 而PSNR， SSIM和MAE評(píng)價(jià)算法的保真能力。 一般來說， IF， PSNR的值越大， SSIM的值越接近于1， MAE的值越接近于0， 表明該算法的去條帶效果越好。 表1列出了模擬實(shí)驗(yàn)中四種全參考指標(biāo)計(jì)算結(jié)果， 不難發(fā)現(xiàn)權(quán)重變分算法都有良好的表現(xiàn)。 目視效果和評(píng)價(jià)指標(biāo)均說明， 權(quán)重變分算法能夠有效去除常見類型的條帶噪聲。

表1 模擬實(shí)驗(yàn)的評(píng)價(jià)結(jié)果

圖3 不同模擬條帶噪聲和去條帶結(jié)果

3 結(jié)果與討論

地基IDOAS于2018年夏季在四川樂山進(jìn)行了外場實(shí)驗(yàn)， 圖2(b)中從左到右的障礙物分別為信號(hào)塔、 山丘、 工廠煙囪以及其他建筑物。 工業(yè)園區(qū)位于儀器的西南方位， 分布著多個(gè)污染氣體的排放源。 外場實(shí)驗(yàn)中儀器的掃描間隔為1°， 最終每組觀測圖像的像素大小為360×48。 儀器的積分時(shí)間設(shè)置為500 ms， 每組全景掃描的總工作時(shí)間約為15 min。 在每組光譜采集完成后， 利用DOAS技術(shù)分別對(duì)NO2和SO2氣體的斜柱濃度進(jìn)行了反演， 具體的參數(shù)設(shè)置見表2。 圖4展示了實(shí)驗(yàn)中兩組觀測的結(jié)果， 其中圖4(a)和圖4(c)為12:59開始采集的NO2和SO2濃度分布， 圖4(e)和圖4(g)為13:14開始采集的NO2和SO2濃度分布。 比較圖4(a)， (c)， (e)和(g)可以發(fā)現(xiàn)， 條帶噪聲具有很大的變化性： 對(duì)于不同時(shí)間的觀測結(jié)果， 條帶噪聲的空間分布和強(qiáng)度均不同； 對(duì)于不同反演氣體， 條帶噪聲的空間分布和強(qiáng)度也不同。

表2 NO2和SO2 DOAS反演中的參數(shù)設(shè)置

圖4 樂山實(shí)驗(yàn)中兩組NO2和SO2的實(shí)際觀測和去條帶結(jié)果

在權(quán)重變分模型中，λ1，λ2和λ3用于平衡各個(gè)約束條件， 在實(shí)際應(yīng)用中需根據(jù)圖像的具體情況加以調(diào)整。 一般來說， 條帶噪聲的強(qiáng)度越大，λ1應(yīng)設(shè)置越大； 條帶噪聲越稀疏， 需使用更大的λ2； 條帶噪聲越均勻， 適合用更大的λ3。 在外場實(shí)驗(yàn)中為達(dá)到最佳的去條帶效果， 經(jīng)多次測試后λ1，λ2和λ3分別設(shè)置為0.4， 0.001和0.2，ρ1，ρ2和ρ3設(shè)置為100倍的λ1。 去條帶后的結(jié)果可見圖4(b)， (d)， (f)和(h)， 總體上去條帶算法取得了令人滿意的效果。 如圖4(a)方框所示的多個(gè)濃度高值區(qū)域， 在去條帶之后保留了各自的空間分布信息； 圖4(c)方框所示的觀測區(qū)域有大量的隨機(jī)噪聲， 在去條帶之后沒有變得模糊。 結(jié)果表明， 權(quán)重變分算法能夠有效去除條帶噪聲， 能夠保留污染氣體分布信息， 且未出現(xiàn)過度平滑的情況。

圖5展示了去條帶前后的權(quán)重列均值曲線， 其中虛線為去條帶前的權(quán)重列均值， 實(shí)線為去條帶后的權(quán)重列均值。 權(quán)重列均值表示在計(jì)算圖像X(大小為m×n)的列均值c時(shí)， 需考慮權(quán)重因子W，c的第i行分量為

圖5 去條帶前后的NO2和SO2權(quán)重列均值對(duì)比

(11)

在低仰角區(qū)域會(huì)出現(xiàn)完全遮擋的情況， 此時(shí)定義分量ci為0。 從圖5可以看出， 權(quán)重變分算法可以獲得平滑的權(quán)重列均值曲線， 并且能保留原始曲線的變化趨勢(shì)。

由于實(shí)測數(shù)據(jù)沒有理想圖像作為參考， 為評(píng)估外場實(shí)驗(yàn)中去條帶算法的有效性， 采用Borsdorff等提出的方法對(duì)結(jié)果進(jìn)行定量評(píng)價(jià)， 參量γ定義為[15]

(12)

表3 去條帶前后的γ值

4 結(jié) 論

針對(duì)地基IDOAS儀器中存在的條帶噪聲問題， 提出了基于權(quán)重變分模型的去條帶算法。 該算法基于條帶噪聲的方向性和稀疏性， 不需要均勻的參考區(qū)域， 并充分考慮了地基成像中障礙物遮擋的問題。 在測試算法性能的模擬實(shí)驗(yàn)中， 對(duì)四組不同類型條帶的降噪效果證明了該算法具有良好的有效性和保真性。 在2018年夏季的外場實(shí)驗(yàn)中， 地基IDOAS反演結(jié)果為痕量氣體斜柱濃度的二維分布。 運(yùn)用該算法后， 多組NO2和SO2濃度分布中的條帶噪聲得到了有效去除。 結(jié)果表明， 該算法適用于多種氣體以及變化條帶的情況， 提高了地基IDOAS觀測數(shù)據(jù)的可用性。