基于中頻域維納濾波的非視域成像算法研究*

唐佳瑤 羅一涵 謝宗良 夏詩燁 劉雅卿 徐少雄 馬浩統 曹雷

1) (中國科學院光束控制重點實驗室,成都 610209)

2) (中國科學院光電技術研究所,成都 610209)

3) (中國科學院大學,北京 100049)

4) (中國科學院大學電子電氣與通信工程學院,北京 100049)

非視域成像是對探測器視線外被遮擋的物體進行光學成像的新興技術,基于光錐變換反演法的非視域成像可以看作是一個反卷積的過程,傳統維納濾波反卷積方法是使用經驗值或者反復嘗試得到瞬態圖像的功率譜密度噪信比(power spectral density noise-to-signal ratio,PSDNSR)進行維納濾波反卷積,但非視域成像每個隱藏場景的PSDNSR 都不同,先驗估計難以適用.因此本文提出使用捕獲瞬態圖像的中頻域信息來估計PSDNSR 進行維納濾波從而實現非視域成像.實驗表明,基于中頻域維納濾波的非視域成像算法估計的PSDNSR 能夠落在一個重建效果較好的量級上.相比其他方法,本文算法能一步直接估計出PSDNSR,沒有迭代運算,計算復雜度低,能夠在保證重建效果的前提下提升重建效率.

1 引 言

非視域(non-line-of-sight,NLOS)成像是對視線外的隱藏物體進行光學探測和可視化的新興技術,類似于“視線拐彎”或“隔墻觀物”[1?7].在機器視覺、制造業、醫學成像、自動駕駛、軍事反恐等領域具有廣闊的應用前景.NLOS 基于飛行時間探測技術,主動向一個中介反射面發射激光,利用探測器捕獲從隱藏物體上發生散射返回的光子的空間時間信息來重建隱藏物體的形狀[8].

近年來,非視域成像技術成為研究熱點之一.2012 年,Velten 等[9]使用條紋相機和超快激光器搭建非視域成像系統,并提出了一種橢圓反投影算法對隱藏物體進行三維重建.2016 年,Klein 等[10]提出了基于合成瞬態渲染的NLOS,該研究將NLOS 轉化為了凸優化問題.2018 年,O’Toole 等[11]探索了一種適用于共焦NLOS 系統的光錐變換反演法(the light-cone transform,LCT),LCT 降低了重建復雜度.2019 年,Liu 等[5]提出了基于相量場虛擬波的NLOS,將NLOS 過程公式化為波成像問題,可以將經典光學中成熟的見解和技術應用在NLOS 領域,即將NLOS 成像系統模擬為一個虛擬的視線成像系統.在NLOS 定位和成像的研究中,為了提高重建質量和效率,除了對捕獲的飛行時間信息進行濾波處理外[12],主要從三維重建算法、成像裝置和中介反射面的選擇等方面進一步改進,如對反投影算法改進從而實現了多目標NLOS、遠距離NLOS 和快速反投影NLOS[2,13?15].也有研究引入深度學習來解決NLOS 問題[16?19].在之前的研究中,使用條紋相機、單光子雪崩二極管(single-photon avalanche diode,SPAD)、飛行時間(time of flight,TOF)相機等作為探測器做了大量研究[9,20,21].2021 年,中國科學技術大學研究團隊基于短激光脈沖泵浦技術構建了一個工作在近紅外波段的上轉換單光子探測器(upconversion single-photon detector,UCSPD),有效提高了探測器的時間分辨率[7].

目前主流的NLOS 是基于瞬態光傳輸進行研究的,主要包括橢圓反投影方法、光錐變換反演法和凸優化法等.橢圓反投影法對于存儲和處理的要求較高,并且接收到的回波信號弱[9,13,14,20],基于光錐變換的非視域成像通過使用共焦掃描系統解決了回波信號弱的問題.光錐變換反演法將探測器捕獲到的瞬態圖像表示為一個三維卷積,該卷積在變換域中模擬自由空間的光傳輸,光錐變換提供了快速且高效的方式計算反向光傳輸.可以將光錐變換反演法看作是反卷積的過程,維納濾波是經典的反卷積方法,其中功率譜密度噪信比(power spectral density noise-to-signal ratio,PSDNSR)決定了重建成像的質量.在非視域成像中每種隱藏場景的PSDNSR 不同,每次都需要調節,PSDNSR 通常采用經驗值或者反復嘗試來獲取.該方法進行反卷積很難一步就找到最佳的PSDNSR,需手動調節PSDNSR 進行多次實驗.而快速地進行反卷積對非視域成像的實時應用至關重要,因此本文重點研究改進的維納濾波算法使光錐變換反演法更高效.

2 基于光錐變換的非視域成像原理模型

本文同樣使用光錐變換反演法進行非視域成像,因此本文使用共焦光路,實驗設置如圖1 所示.通過使用分束鏡將激光發射光路與探測器接收光子的光路合成共焦非視域成像系統.激光器向中介反射面發射脈沖激光,經過三次散射后使用超導納米線單光子探測系統(superconducting nanowire single-photon detector,SNSPD)作為探測器捕獲返回的光子信息.時間相關單光子計數 (timecorrelated single-photon counting,TCSPC)模塊對SNSPD 捕獲到的每個像素的光子數進行直方圖統計,得到瞬態圖像.

圖1 實驗場景示意圖,激光通過振鏡對中介面掃描,探測器接收來自中介面反射的直接光和來自隱藏物體的間接光Fig.1.Experimental scene,the laser scans the intermediate surface through the galvanometer,and the detector receives direct light reflected from the intermediate surface and indirect light from a hidden object.

探測器捕獲到的瞬態圖像可以表示為[11]

式中,(x′,y′) 是激光照射在中介反射面上的點,(x,y,z) 是隱藏物體表面上的點,r表示這兩個點之間的距離.ρ是滿足z＞0 的三維半空間?中隱藏物體表面每個點的反照率,函數δ表示由x2+y2+z2=(tc/2)2給出的時空四維超圓錐的表面.

對瞬態成像模型做變量代換后可以寫成一個3D 卷積:Rt{g}=h ?Rz{ρ}.其中,h表示移位不變的卷積核,即Rt是對g沿t軸衰減及重采樣,Rz是對ρ沿z軸衰減及重采樣.接下來,對成像模型進行離散化,在空間域中經變量代換和離散后的瞬態圖像可以表示為

維納濾波是經典且有效的反卷積方法,文獻[11]正是使用維納濾波反卷積,最后得到隱藏物體的反照率為

一般來說,K值的選取決定著非視域成像的質量,通常采用經驗值或反復嘗試來手動調節,但該方法成像效率不高.并且常數項K值代替Π(u,v,w)進行維納濾波時,沒有用到足夠的先驗知識,導致難以快速找到最佳值.因此本文引入瞬態圖像中頻域的方法實現對K的快速計算.

3 基于中頻域維納濾波的非視域成像算法

瞬態圖像的頻譜G(u,v,w) 是一個M ×N ×P的復數矩陣,如圖2 所示為w=50 時瞬態圖像的頻譜圖 |G(u,v,w)|.可見瞬態圖像的大部分信息分布在頻譜原點G(u,v,w)=(0,0,50) 附近的“低頻”區域,卷積核的信息被完全淹沒在瞬態圖像的信息中,而頻譜的“高頻”區域幅度小且易被噪聲污染.由于“低頻”和“高頻”區域的頻譜幅度相差大,不會出現交疊,兩者之間必然存在一個中間過渡區域,即不含太多的瞬態圖像能量又沒有被噪聲淹沒,將之稱為“中頻域”,可以用于估計K值從而進行維納濾波反卷積[22].非視域成像系統中捕獲的瞬態圖像中大多數信號變化緩慢,只有少數信號變化大,頻譜往往為全局單減,因此瞬態圖像存在中頻域并且數值特征大都相似.

圖2 w =50 時瞬態圖像的頻譜圖|G(u,v,w)|Fig.2.Spectrum |G(u,v,w)|of the transient image at w=50.

如圖3 所示,取w=50 時瞬態圖像的幅度|G(u,v,w)|中過原點的一條線 |G(u,0,50)|,G(u,0,50) 是 共軛對稱的,所以只考慮前半個周期u ∈UT=[0,uT] 上的值,其中表示向零取整.那么低頻域位于u=0 附近,高頻域位于u=uT附近.通常該曲線中有兩個明顯的全局轉折點,即中頻域的邊界點.因為低頻域聚集了瞬態圖像的大部分信息,所以在曲線 |G(u,0,50)|中可以看到一個清晰的轉折點uLM,即低頻域和中頻域的分界點.而高頻域的值很小且被近似均勻分布的噪聲所污染,所以可以在曲線 l n|G(u,0,50)|中看到另一個清晰的轉折點uMH,即中頻域和高頻域的轉折點.但是,由于 l n|G(u,0,50)|通常是劇烈振蕩的,有很多局部轉折點,很難正確找到uMH,因此將曲線 l n|G(u,0,50)|平 滑之后得到Gr(u) 來計算轉折點uMH.如圖4 所示,根據曲線Gr(u) 的幾何特性,過(0,Gr(0)) 和 (uT,Gr(uT)) 兩點作直線,那么Gr(u) 上距離該直線最遠的距離就是uMH,因此可以計算出uMH為

圖3 幅度 |G(u,v,w)|中過原點的一條線|G(u,0,50)|Fig.3.A line |G(u,0,50)|through the origin in the spectral magnitude |G(u,v,w)|.

圖4 幅度 l n|G(u,0,50)|平滑之后的曲線Gr(u)Fig.4.Curve G r(u) after amplitude ln|G(u,0,50)|smoothing.

這時,對K值的估計問題就轉換為了獲取K值使得T1(u) 在抑制噪聲的基礎上最接近T0(u) .因為低頻域和中頻域中的噪聲很小可以忽略,有T0(u)≈T1(u)≈T2(u),u ∈ULM,其 中ULM為低頻域和中頻域.在高頻域中,噪聲近似為均勻分布,有T0(u)≈T0(uMH),u ∈UH.為了抑制噪聲的影響,使T1(u)≧T0(u) 在高頻區域成立,即:

其中u∈UT=[0,uT],上述引入的參數η稱為“噪聲抑制參數”,用來調節噪聲和圖像細節的平衡.

上文利用空間維度上瞬態圖像的頻譜圖|G(u,v,w)|的中頻域來估計K值.分析時間維度上瞬態圖像的頻譜圖,如圖5 所示為u=0 時瞬態圖像的頻譜圖 |G(u,v,w)|.取幅度 |G(u,v,w)|過原點的一條線 |G(0,0,w)|,如圖6 所示,也可計算出該曲線的中高頻轉折點wMH:

圖5 u =0 時瞬態圖像的頻譜圖|G(u,v,w)|Fig.5.Spectrum |G(u,v,w)|of the transient image at u=0.

圖6 幅度 |G(u,v,w)|中過原點的一條線|G(0,0,w)|Fig.6.A line |G(0,0,w)|through the origin in the spectral magnitude |G(u,v,w)|.

其中w∈WT=[0,wT] 是 曲線 |G(0,0,w)|的前半周期,Gr(w) 是曲線 l n|G(0,0,w)|平滑之后的曲線.

同理對K值進行估計:

分別使用空間維度上的頻譜圖和時間維度上的頻譜圖估計的K值沒有很大的差異,本文選擇使用時間維度上的頻譜圖 |G(0,0,w)|對K值進行估計.

針對基于光錐變換的非視域成像問題,使用基于中頻域的維納濾波算法對K值進行估計并完成重建的步驟如下:

1)在共焦非視域成像系統中,使用探測器捕獲瞬態圖像g(x′,y′,t) ;

3)計算出曲線 l n|G(0,0,w)|,再對該曲線平滑得到Gr(w) ;

4)根據(11)式計算出wMH;

5)根據(12)式計算出K值;

6)使用維納濾波算法對瞬態圖像反卷積并進行傅里葉反變換得到;

7)最終使用(4)式得到隱藏物體表面的反照率ρ?.

4 實驗及結果分析

搭建的實驗場景如圖7 所示,該共焦非視域成像系統包含的光源為波長1530 nm 的激光脈沖,脈沖寬度為70 ps,重復頻率為40 MHz,平均功率為750 mW.采用的探測器為超導納米線單光子探測系統(SNSPD),其探測效率約為70%.使用分束器(Thorlabs CCM1-BS015/M)將探測光路和激光發射光路共軸.使用振鏡(Thorlabs GVS012)來實現對中介面的掃描.時間相關單光子計數模塊以1 ps 的時間分辨率對探測事件進行時間標記.

圖7 實驗場景 (a)為在實驗室搭建的共焦非視域成像場景;(b)為共焦非視域成像系統Fig.7.Experimental scene: (a) The confocal non-horizontal imaging scene built in the laboratory;(b) demonstrates the confocal non-horizontal imaging system.

本文通過上述實驗裝置對如圖8 所示的5 個隱藏場景捕獲光子飛行時間信息,5 個隱藏場景分別為T 形紙板、人偶模型手臂向下、房屋形紙板、C 形紙板和人偶模型手臂向上.

圖8 隱藏場景 (a) T 形紙板;(b)人偶模型手臂向下;(c)房屋形紙板;(d)C 形紙板;(e)人偶模型手臂向上Fig.8.Hidden scene: (a) T-shaped cardboard;(b) puppet model with arms down;(c) house-shaped cardboard;(d) Cshaped cardboard;(e) puppet model with arms up.

(12)式中的參數η用來調節噪聲和圖像的細節,噪聲和圖像的清晰度會隨η減小而變弱.在評估了噪聲大小的基礎上,本文選擇η=1.1 來對K值進行估計,此時噪聲在可接受范圍內且細節更加豐富.基于中頻域的維納濾波算法對T 形紙板、人偶模型手臂向下、房屋形紙板、C 形紙板和人偶模型手臂向上5 個隱藏場景K值的估計值分別為2.8678,1.2353,3.6711,0.8096 和1.5939.

為了說明本文算法對非視域成像的有效性,分別取K為0.01,0.1,1,10,100,1000 對隱藏場景進行重建,將傳統維納濾波的重建結果與基于中頻域的維納濾波得到的結果進行對比.為方便起見,只對比重建結果的正視圖,見圖9—圖18.

圖9 隱藏物體T 形紙板和基于中頻域的維納濾波的重建結果Fig.9.T-shaped cardboard and reconstruction results of NLOS imaging based on mid-frequency Wiener filtering.

圖10 傳統維納濾波的重建結果與基于中頻域的維納濾波的重建結果對比 (a)?(f) T 形紙板的隱藏場景分別取K 為0.01,0.1,1,10,100,1000 進行維納濾波復原的結果;(h) T 形紙板基于中頻域的維納濾波重建的結果,估出K 為2.8678Fig.10.Comparison between the reconstruction results of traditional Wiener filtering and the reconstruction results based on Wiener filtering in the mid-frequency domain: (a)?(f) The results of Wiener filtering reconstruction with K as 0.01,0.1,1,10,100,and 1000 for T-shaped cardboard respectively;(h) the reconstruction result of T-shaped cardboard using the NLOS imaging algorithm based on Wiener filtering of mid-frequency domain,and the estimated K is 2.8678.

圖11 隱藏場景人偶模型手臂向下和基于中頻域的維納濾波的重建結果Fig.11.Puppet model with arms down and reconstruction results of NLOS imaging based on mid-frequency Wiener filtering.

圖12 傳統維納濾波的重建結果與基于中頻域的維納濾波的重建結果對比 (a)?(f)人偶模型隱藏場景分別取K 為0.01,0.1,1,10,100,1000 進行維納濾波復原的結果;(h)人偶模型隱藏場景基于中頻域的維納濾波重建的結果,估計出的K 為1.2353Fig.12.Comparison between the reconstruction results of traditional Wiener filtering and the reconstruction results based on Wiener filtering in the mid-frequency domain: (a)?(f) The results of Wiener filtering reconstruction with K as 0.01,0.1,1,10,100,and 1000 for puppet model with arms down respectively;(h) the reconstruction result of puppet model with arms down using the NLOS imaging algorithm based on Wiener filtering of mid-frequency domain,and the estimated K is 1.2353.

圖13 隱藏場景房屋形紙板和基于中頻域的維納濾波的重建結果Fig.13.House-shaped cardboard and reconstruction results of NLOS imaging based on mid-frequency Wiener filtering.

圖14 傳統維納濾波的重建結果與基于中頻域的維納濾波的重建結果對比 (a)?(f)為房屋形紙板的隱藏場景分別取K 為0.01,0.1,1,10,100,1000 進行維納濾波復原的結果;(h)房屋形紙板隱藏場景基于中頻域的維納濾波重建的結果,估計出的K 為3.6711Fig.14.Comparison between the reconstruction results of traditional Wiener filtering and the reconstruction results based on Wiener filtering in the mid-frequency domain: (a)?(f) The results of Wiener filtering reconstruction with K as 0.01,0.1,1,10,100,and 1000 for House-shaped cardboard respectively;(h) the reconstruction result of house-shaped cardboard using the NLOS imaging algorithm based on Wiener filtering of mid-frequency domain,and the estimated K is 3.6711.

圖15 隱藏場景C 形紙板和基于中頻域的維納濾波的重建結果Fig.15.C-shaped cardboard and reconstruction results of NLOS imaging based on mid-frequency Wiener filtering.

圖16 傳統維納濾波的重建結果與基于中頻域的維納濾波的重建結果對比 (a)?(f) C 形紙板的隱藏場景分別取K 為0.01,0.1,1,10,100,1000 進行維納濾波復原的結果;(h) C 形紙板隱藏場景基于中頻域的維納濾波重建的結果,估出K 值為0.8096Fig.16.Comparison between the reconstruction results of traditional Wiener filtering and the reconstruction results based on Wiener filtering in the mid-frequency domain: (a)?(f) The results of Wiener filtering reconstruction with K as 0.01,0.1,1,10,100,and 1000 for C-shaped cardboard respectively;(h) the reconstruction result of C-shaped cardboard using the NLOS imaging algorithm based on Wiener filtering of mid-frequency domain,and the estimated K is 0.8096.

圖17 隱藏場景人偶模型手臂向上和基于中頻域的維納濾波的重建結果Fig.17.Puppet model with arms up and reconstruction results of NLOS imaging based on mid-frequency Wiener filtering.

圖18 傳統維納濾波的重建結果與基于中頻域的維納濾波的重建結果對比 (a)?(f)人偶模型手臂向上的隱藏場景分別取K 為0.01,0.1,1,10,100,1000 進行維納濾波復原的結果;(h)人偶模型手臂向上隱藏場景基于中頻域的維納濾波重建的結果,估出的K 為1.5939Fig.18.Comparison between the reconstruction results of traditional Wiener filtering and the reconstruction results based on Wiener filtering in the mid-frequency domain: (a)?(f) The results of Wiener filtering reconstruction with K as 0.01,0.1,1,10,100,and 1000 for puppet model with arms up respectively;(h) the reconstruction result of puppet model with arms up using the NLOS imaging algorithm based on Wiener filtering of mid-frequency domain,and the estimated K is 1.5939.

考慮到人眼對于圖像的主觀評價不是從單一角度出發,因此本文選擇綜合圖像Tenengrad 梯度和結構相似性(structural similarity,SSIM)兩個指標來進行像質評價[23],從清晰度和空間結構性兩方面來評價重建圖像.從圖像像素的角度來看,圖像Tenengrad 梯度GRAD 越大,圖像內的邊緣越清晰,圖像質量越好;原始圖像和重建圖像的結構相似性越大,重建效果越好,圖像質量越好.考慮到圖像Tenengrad 梯度和SSIM 兩個指標的一致性和量級后,通過線性加權,得到最終圖像評價指標為

式中,α和β分別為圖像Tenengrad 梯度和SSIM 的權值,α+β=1 .本文根據對重建圖像的分析,分別取α=0.1 ,β=0.9 .Eval越大,圖像質量越好,Eval越小,圖像質量越差.圖16 所示為對不同K值維納濾波的重建結果和基于中頻域的維納濾波重建結果的客觀評價.

如圖19 所示,本文方法計算出5 個隱藏場景的K值分別在 1 00,1 00,1 00,1 0?1,1 00量級上.基于中頻域的維納濾波算法取得的K值總能落在一個最佳的量級上,K值在該量級時,圖像綜合質量評價值最高,圖像重建效果最好.可以看出,基于中頻域維納濾波的非視域成像算法重建圖像的Eval值在一個較高的范圍,說明本文方法估計出的K值進行非視域成像得到的重建圖像效果好,并且速度更快.實驗結果表明,使用本文方法能一步估出K值,并且接近最佳值.

圖19 重建圖像綜合質量評價,其中藍色線為設置不同K 得到的重建圖像的綜合質量評價值,紅色圓圈為基于中頻域的維納濾波重建圖像的綜合質量評價值 (a) T 形紙板;(b)人偶模型手臂向下;(c)房屋形紙板;(d) C 形紙板;(e)人偶模型手臂向上Fig.19.Comprehensive quality evaluation of reconstructed images,the blue line is the comprehensive quality evaluation value of the reconstructed image obtained by setting different K,and the red circle is the comprehensive quality evaluation value of the reconstructed image based on the Wiener filter in the intermediate frequency domain: (a) The Eval of T-shaped cardboard,(b) the Eval of puppet model with arms down,(c) the Eval of house-shaped cardboard,(d) the Eval of C-shaped cardboard,(e) the Eval of puppet model with arms up.

5 結 論

共焦光路中的非視域成像可以看作是一個反卷積問題,在維納濾波反卷積的過程中參數K值的選取對成像的速度以及重建的質量有很大的影響,非視域成像中每種隱藏場景的最優K值都不一樣,每次實驗都需人為調節.針對該問題本文引入了基于中頻域的非視域成像算法,可以使用瞬態圖像的中頻域信息來估計K值.相比其他取K值的算法,該算法沒有迭代和矩陣運算,計算復雜度低,能夠快速的確定K值.為了驗證本文算法的有效性,本文進行共焦非視域成像實驗,將一系列手動設置的K值的重建結果與基于中頻域的非視域成像算法得到的重建結果進行了對比.實驗結果表明,基于中頻域的維納濾波算法取得的K值能落在成像效果接近最佳的量級上.綜上,該算法具有快速、準確和參數少的特點,有效提升了共焦非視域成像的重建質量和實時性.