基于高階邊界散射變換的二面角快速成像算法

楊俊慧 諸葛曉棟

（北京航空航天大學電子信息工程學院 微波感知與安防應用北京市重點實驗室, 北京 100091）

0 引 言

毫米波與光波相比具有更長的波長，可以輕易穿透許多光學上不透明的物體，如普通的衣服材料、茂密的植被、煙霧、霧氣、墻壁等.因此，毫米波成像通常被用于穿墻成像、安全檢查、地面穿透雷達等.

為了對目標進行精確成像，我們通常需要高質量的圖像重構算法.一系列的毫米波主動成像算法已經被提出來用于目標檢測、識別和成像.后向投影(back propagation, BP)算法[1]是一種經典的時域算法，也是最常用的成像技術之一，通過建立數據域到成像域的投影并進行相干求和，可以實現最準確的成像，但存在計算量大和效率低的問題.波數域(ω-k)算法[2]通過傅里葉變換加速成像，但在頻域進行的Stolt 插值會影響成像質量.距離-多普勒(range-Doppler, RD)算法[3]也需要插值，以消除由距離徙動引起的縱向和橫向間的相互耦合.Chirp Scaling(CS)算法[4]無須進行插值操作，但在該算法中回波信號的距離譜會發生變化.基于以上傳統成像算法，一些人體成像應用如SafeView 和QPS[5-6]在檢測人體中隱藏的違禁品方面也取得了很好的效果.但它們通常面臨陣列過于龐大、計算量大和成像速度慢等問題.

由于人體安檢成像的最終目標是檢測、提取目標輪廓并最終識別隱藏的違禁品，因此，如何快速準確地獲取目標輪廓是我們重點解決的問題.邊界散射變換(boundary scattering transform, BST)算法通過建立接收波前和目標形狀之間的一一對應關系，可以直接從接收波前恢復目標的輪廓，而不需要一個完整的圖像重構過程，極大地提升了目標檢測的速度[7].SEABED (Shape Estimation Algorithm based on BST and Extraction of Directly scattered waves)[8]是一種基于邊界散射變換（boundary scattering transform,BST）的實時三維成像算法，該算法通過求導操作建立了接收波前與目標輪廓間的關系，但也正因利用了接收波前的導數，該算法對誤差相當敏感.在嘈雜的環境中，SEABED 獲得的圖像會變差.Envelope[9]使用圓的包絡來建立對應關系，這種方法不利用接收波前的導數，可以對任意簡單形狀的目標實現穩健成像，但它需要對觀測距離進行精確連接以保證成像質量.對于復雜形狀的目標或多個目標，因為每個天線從目標表面的多個散射點接收多個回波，對于一個復雜的表面，這種連接往往很困難，因為每個天線觀察到多個回波，有太多的候選連接.距離點徙動（range points migration，RPM）算法[10]通過對接收點到達方向(direction of arrival, DOA)的精確估計，實現了從接收點到目標點的直接映射，消除了距離連接過程，即使在復雜的邊界提取中也能顯著提高精度.然而，若目標表面變化尺度小于脈沖寬度，當面對復雜目標（如凹凸結構）時，該算法仍然面臨由多次散射信號引起的圖像失真問題.

需要注意的是，傳統的基于BST 的外形重構算法往往把所有的反射當作單次反射，而對于復雜的目標，特別是凹腔目標，其內部往往存在多次反射，目標的某些部分只有多次反射信號能照射到，只考慮單次反射將導致重構目標形狀的部分缺失和額外的成像偽影.

本文中，我們提出了一種基于高階邊界散射變換(high-order boundary scattering transform, HBST)的外形重構算法.與傳統BST 算法中采用單站線掃描不同，我們選擇了一個具有二面角結構的近場單站圓周掃描場景進行研究.通過建立正向散射模型，我們發現出射邊駐定相位點所在路徑的散射對回波起主要貢獻，也即在二面角內的各種散射中，鏡面反射起主導作用.基于此，我們在后續的推導中只考慮鏡面反射對散射場的貢獻，并推導了出射邊的駐定相位點與回波波前之間的對應關系，并建立了可快速而精確處理二面角內多次反射回波的外形重構算法.該算法有效地消除了傳統BST 算法中出現的偽影，提高了重構精度.仿真和實驗結果都證明了該算法對二面角結構的有效性.

1 基本原理

1.1 二面角回波分析

為研究二面角的散射特性，我們首先對其回波進行分析，對二面角采用單站圓周掃描，將整個模型放在極坐標系下研究.如圖1 所示，二面角左右對稱，張角為 2φ(即左右邊所在角分別為?φ和+φ)，角平分線位于0°觀測角，頂點位于掃描中心圓點(0,0)，天線收發點為P0(R,θ0).設置二面角的邊長L為0.3 m，圓周掃描半徑R為0.7 m.由于我們只關注二面角內部的散射回波，因此將收發天線的掃描范圍 θ0限制為?60° ～ 60°，掃描間隔為0.5°.

圖1 單站圓周掃描幾何示意圖Fig.1 Schematic geometry of monostatic cylindrical scanning

基于電磁仿真軟件Ansys HFSS 中的高頻近似算法彈跳射線(shooting and bouncing rays, SBR)法，我們對不同張角(90°，70°，55°，40°，36°)的二面角進行仿真(頻率范圍為2～40 GHz)，相應的散射回波如圖2 所示，回波圖案顯示了在不同觀測角接收到的反射信號在二面角內傳播的路徑長度.由于二面角頂點位于掃描中心且左右對稱，因此回波圖案也在0°觀測角附近對稱.隨著二面角的張角變小，散射回波中也有更高次反射出現，且同一反射次數(reflection times, RT)被接收到的觀測角位置也在不斷變化.二面角內出現的最大反射次數(maximum number of reflection times, MRT)和其張角的關系為[11-12]：

圖2 不同張角(90°, 70°, 55°, 40°和36°)二面角的時域回波Fig.2 Time-domain echo signal of dihedrals with different opening angles (90°,70°, 55°,40° and 36°)

式中，ceil 表示向上取整.

表1 展示了常見張角內出現的MRT，由于偶次反射不能正確重構二面角邊長[11]，基于圓極化波每反射一次極化方向則會改變這一規律[13-14]，使用圓極化收發器過濾偶次反射，只使用奇次反射進行目標重構.在后文中我們也僅研究二面角的奇次反射回波，不再對其偶次反射回波進行進一步分析.因此，圖2 中張角為90°、70°、55°、40°、36°的二面角對應的最大奇次反射次數為1RT、3RT、3RT、5RT、5RT.其中1RT 信號所在觀測角距0°最遠，傳播距離最短，最高次反射回波則在0°觀測角附近被接收，其傳播距離接近掃描半徑R的兩倍.

表1 常見的張角和相應的MRTTab.1 Common opening angles and corresponding MRT

1.2 駐定相位點對回波的影響

在1.1 節中我們展示了不同張角二面角的回波分布，但對二面角內的多次反射機理卻仍不明了，接下來對這一問題進行研究.假設二面角為電大目標，滿足高頻近似條件ka≥10，其中k為波數，a為目標孔徑[15-16].在此近似下，可以將電磁波的傳播近似為射線傳播，基于SBR 法假設除最后一次反射為全向散射外，電磁波在二面角內的反射均為鏡面反射，圖3為二面角內3RT 信號的傳播路徑示意圖.基于鏡像反射原理，我們可以求出2n?1次反射的鏡像反射點P′(R,θ2n?1)，其中鏡像反射點P′的角度為

圖3 二面角內三次反射回波的傳播路徑示意圖Fig.3 Schematic of the propagation path of triple reflections within the dihedral

需要注意的是， θ為首次入射邊的角度.由于在P0點發射的電磁波在二面角內反射后到達最終的出射點P所走的路徑長度與2n?1次鏡像反射點P′到P的路徑長度相同，因此可以求出該段路徑長度為

同時，可以求出出射點P到接收點P0的距離：

因此，根據總傳播距離Rall=d1+d2可以建立接收場的表達式：

也即在某一觀測點P0(R,θ0)接收到的散射回波為多條傳播路徑在出射邊上各點的散射回波之和.如圖4 所示，出射邊上每個點P(r,θ)均向接收點P0(R,θ0)散射回波，回波傳播路徑d=|PP′|+|PP0|.根據式(5)中的接收場表達式可以得到圖4 所示的回波.需要注意的是，當收發點位于大觀測角度時(|θ0|＞φ)，由于遮擋效應，可能存在接收點無法接收到部分邊散射回波的情況.如圖5 所示，由于左側邊遮擋，此時只能接收到(rmin,L)區域的散射回波，其中 θ表示出射點所在邊的角度(由于本文僅討論奇次反射，因此和入射點所在邊的角度相同)；當|θ0|＞φ時，rmin為觀測點P0與遮擋邊端點P1的連線和出射點所在邊OP3的交點P2到原點O的長度，可以表示為

圖4 右側邊出射的散射波在二面角內的等效傳播路徑Fig.4 Schematic of the equivalent propagation path of the scattered wave emitted from the right side of the dihedral

圖5 二面角的遮擋效應示意圖Fig.5 Schematic of the shading effect of the dihedral

但在實際散射中，圖4 中各條路徑的散射回波對總接收場的貢獻并不相同.根據駐定相位原理[17]，我們對式(5)中的指數項求導：

可求得對式(5)中接收場積分貢獻最大的點(駐定相位點)，也即在觀測角 θ0接收到的散射回波中貢獻最大的散射點在出射邊的位置r：

不難求證，該散射點所對應傳播路徑在二面角內均按照鏡面反射來傳播[18].

圖6（a）和圖7(a)分別展示了在70°二面角右側邊出射的1RT 和3RT 回波信號圖，在兩條白虛線之間的觀測角內存在強回波信號；圖6（b）和圖7 (b)展示了在不同觀測角下各出射點r∈(0,L)所在路徑相應的距離分布圖，其中強度表示距離信息.圖中紅色點為根據公式(8)求出的各觀測角下駐定相位點的位置，由于在大部分觀測角下駐定相位點都不在出射邊上(也即r?(0,L)，此時L=0.3 m)，因此只有白色虛線附近才出現駐定相位點.這意味著當駐定相位點位于出射邊上時，在該觀測角下能接收到強散射回波，這也側面說明在所有出射點中，駐定相位點對回波貢獻最大.需要注意的是，雖然在圖6 (b)中左邊虛線附近的觀測角度下也存在駐定相位點，但此時由于左側邊的遮擋，在這些觀測角下并不能接收到散射回波.

圖6 70°二面角右側邊出射的1RT 信號的回波和不同觀測角下各出射點所在路徑的距離分布Fig.6 Echo and distance distribution of the path of all exit points under different observation angles (of 1RT signal emitted from the right side of a 70° dihedral)

圖7 70°二面角右側邊出射的3RT 信號的回波和不同觀測角下各出射點所在路徑的距離分布Fig.7 Echo and distance distribution of the path of all exit points under different observation angles (of 3RT signal emitted from the right side of a 70° dihedral)

我們將公式(8)代入公式(3)和公式(4)，可以求出駐定相位點所在路徑的總傳播距離為

同時，根據式(8)可以推導出觀測角 θ0和駐定相位點位置r的關系：

根據r∈(rmin,L)這一條件可以求出 θ0的范圍，再結合式(9)和式(10)即可求出各觀測角 θ0下駐定相位點所在路徑的總傳播距離Rall.如圖8 所示，我們給出了90°、70°和40°二面角駐定相位點所在路徑的總傳播距離Rall和觀測角 θ0之間的關系，其中黃色虛線由公式（9）得到，由于只在某些觀測角度能接收到相應回波，因此有效區域如黑色實線所示.可以看出，該曲線與圖2中相應二面角回波波前位置的一致性，因此可以通過駐定相位點所在路徑的傳播距離表示回波波前.

圖8 在各觀測角下不同張角二面角出射點駐定相位點所在路徑的總傳播距離Fig.8 Total propagation distance of the path where the stationary phase point stands of different opening angles of the dihedral at each observation angle

1.3 HBST 算法

從1.2 節可知，駐定相位點所在的全鏡面反射路徑對回波貢獻最大，因此只考慮鏡面反射對回波的貢獻，并基于此建立駐定相位點(r,θ)與回波波前點(Rall,θ0)間的一一對應關系.鏡面反射路徑的傳播規律為最終垂直入射二面角邊后沿原路返回[18]，圖9(a)～(c)分別對應二面角內1RT、3RT、5RT 的傳播路徑，同時也展示了各次反射中的極限情況，紅色實線表示各次反射能照射到的區域，紫色實線表示能接收到相應反射信號的角度范圍.

圖9 二面角內一、三、五次反射的傳播路徑和極限情況Fig.9 Propagation paths and the limit cases of 1RT, 3RT and 5RT within the dihedral

不難發現，隨著反射次數增加，電磁波在二面角內能照射到的區域也越來越廣，最高次反射能照射到二面角內所有區域.理論上可以通過最高次反射回波重構二面角外形，但由于二面角存在缺口時僅用最高次反射無法完整重構目標，因此在我們的算法中所有反射回波均會被用來進行重構.

由于電磁波在二面角內最終垂直入射后沿原路返回，故圖3 中的觀測點P0、出射點P和鏡像反射點P′在一條直線上，如圖10 所示.因此，電磁波由觀測點P0出射并在二面角內反射后再回到觀測點P0的總傳播距離Rall等效于觀測點P0到鏡像反射點P′的距離|P0P′|，其中等效反射點Pe位于|P0P′|中點，也即|P0Pe|=Rall/2，觀測點P0發射的電磁波到等效反射點Pe的往返距離等于真實傳播距離Rall.若將該2n?1次反射當作一次反射處理，則重構點位于等效反射點Pe(re,θe)處.當n=1時，如圖10(a)所示，等效反射點Pe與駐定相位點P重合；當n＞1時，如圖10(b)所示，等效反射點Pe位于駐定相位點P外側，其中駐定相位點P所在角度θ 與等效反射點Pe所在角度θe的關系為：θ=θe/(2n?1)，此時若將2n?1次反射當作一次反射處理，則重構點位于角度 θe上，在二面角外側出現成像偽影.

圖10 一次反射和高次反射路徑中等效反射點和駐相點的關系Fig.10 Relationship between equivalent reflection points and standing phase points in 1RT and high-order RT paths

根據|θe?θ0|=∠PeOP0，我們可以求出駐定相位點P(r,θ)所在角度：

根據式(8)和式(11)，我們可以建立由回波波前點(RAB,θ0)到駐定相位點(r,θ)的重構關系，即針對多次反射回波信號的HBST 算法.當處理2n?1次反射回波信號時，將式(11)中的分母設置為相應的2n?1，即可得到正確的目標角度.此外，根據電磁波在二面角內的反射路徑可知，駐定相位點(r,θ)可以視為2n?1次反射中的首次反射點r1和第2n?1次反射點r2n?1，通過以下公式可以不斷迭代求出其他反射點：

1.4 算法應用流程

基于式(11)可知，為了正確進行外形重構，有必要區分不同反射次數和不同入射方向上的波前點(Rall,θ0).然而，若二面角內的m1RT信號由正確入射方向的m2RT算法(即式(11)中分母取為m2)處理，那么處理后的單邊角θ′=m1/m2·θ.當m1≠m2時, θ′仍為常數但不等于正確角度 θ.然而，如果將首先入射到二面角右側邊或左側邊的信號視為首先入射到左側邊或右側邊的信號，則會出現角度偏差|θ ?θ′|=2ψ/(2n?1).由于 ψ不是一個常數，所以 θ′的值也是變化的.在此基礎上，我們提出一種方法，以避免手動劃分屬于不同RTs 和入射方向的波前，從而實現成像的自動化.

對于MRT=2nmax?1的波前，我們可以通過各次反射算法(即式(11)分母中n=1,...,nmax)來分別處理所有回波波前點，同時不區分入射方向，得到nmax組結果.通過設置一個適當的角度閾值，保留在所有nmax組結果中均存在的角度范圍對應的點，舍棄其他的點，可以根據所有保留點的位置(r,θ)來重構目標形狀.

圖11(a)和(b)顯示了不區分入射方向，使用1RT 和3RT 算法對接收到的60°二面角回波波前點(Rall,θ0)進行重構的重構點云圖.由于在60°二面角中存在1RT 和3RT 回波，因此當采用正確入射方向的重構算法時會出現恒定的重構角(圖11(a)中的60°和180°角以及圖11(b)中的20°和60°角)，表2 解釋了這些角度的對應關系，而其他非恒定角度的重構點則是由錯誤入射方向算法處理得到的.保留用1RT 和3RT 算法處理得到的兩組結果中都包含的角度范圍的點(也即圖11(a)～(b)中黑色虛線覆蓋的點)，舍棄其他的點，則可以得到最終的重構目標形狀，如圖11(c)所示.需要注意的是，由于BST 算法是基于點之間的關系，我們通常使用點云圖來表示接收到的波前和重構目標的形狀.傳統的BST 算法中所有的重構點都用黑色點表示，因此無法判斷散射強度，這里我們用不同顏色的點來表示相應區域的功率強度.

表2 用m2RT算法處理60°二面角m1RT波前時重構目標的開口角度（正確入射方向）Tab.2 Reconstructed opening angle of 60° dihedral when processing the m1RT wavefront with m2RT algorithm in correct incidence direction

圖11 用HBST 算法處理60°二面角的波前點得到的點云重構結果Fig.11 The point cloud reconstruction results obtained by processing extracted wavefront points of 60° dihedral with HBST algorithm

圖12 為HBST 算法的流程圖，圖13 為對該算法的圖例說明.HBST 算法具體流程下：

圖12 基于HBST 的重構算法流程圖Fig.12 Flowchart of the proposed HBST based reconstruction algorithm

圖13 基于HBST 的重構算法示意圖Fig.13 Illustration of the proposed HBST based reconstruction algorithm

1)選擇合適的功率閾值，并獲得功率大于該閾值的強散射回波的波前點信息(功率值、觀測角 θ0和傳播距離Rall).

圖13(a)為60°二面角的反射回波，圖13(b)則是從圖13(a)中提取的波前點的點云圖，其中點的顏色越深，對應的功率值越大.圖13(a)是一幅圖像，圖13(b)是所有提取的波前點的集合.

2)人為地確定nmax.首先用傳統的BST 算法重構目標，重構結果中最小的張角就是目標的真實角度2φ.根據公式(1)求出MRT，再根據nmax=ceil(MRT/2)確定nmax，然后用各(2n?1)RT (n∈[1,nmax])算法分別處理波前點.如果重構點 (r,θ)由同一RT 算法進行處理，不區分入射方向，將其視為一組數據，則能得到nmax組結果.

3)每組結果中都有許多角度，篩選出每組結果中點數最多的前nmax個角度范圍.由于用正確的入射方向處理所得到的角度是恒定的，所以這些角度的點的數量最多.例如，圖11(a)中角度為60°和180°、圖11(b)中角度為20°和60°的點最多.需要注意的是，在處理過程中，實際角度和理論角度之間存在一定的偏差，因此我們需要設置角度閾值(例如1°)來限制角度范圍.此外，在某一組結果中，對應正確角度的重構點的數量不一定是最多的，因此不能假設包含最多的點的角度是目標角度，它可能是nmax個恒定角度中的任何一個.因此須找到每組結果中點數最多的前nmax個角度范圍.通過這種方式，可以在nmax組結果中得到nmax·nmax個角度范圍，然后計算每個角度范圍的平均值，以便在下一步中更好地進行比較.

4)根據獲得的nmax·nmax個角度平均值，找到所有nmax組結果中均存在的目標角度.可以通過在點云圖中保留目標角度下的點而放棄其他點來獲得目標形狀，如圖13(c)所示.

5）根據所選點的信息，可以通過式(13)將點云圖轉換成圖像，計算出圖像中每個像素的平均功率值Pij即為圖像中每個像素點的強度：

式中：(i,j)為該像素點在圖像中的位置；n為像素點(i,j)中包含的重構點數量；Pt為每個重構點的功率.

2 仿真和實驗結果

我們使用與圖2 相同的掃描結構和頻率，并采用BP 算法、傳統的BST 算法和本文提出的基于HBST 的重構算法分別對張角為70°、 55°和45°二面角進行重構，結果如圖14 (a)～(c)所示.此時重構二面角的張角分別為69.56°、54.98°、45.04°，可見重構角度具有較小的誤差.

圖14 不同張角二面角的仿真重構結果Fig.14 Simulation reconstruction results of dihedrals with different opening angles

由于在張角小于90°時高階反射出現，此時BP 算法和傳統的BST 算法成像結果中出現重構偽影，只有二面角的外邊沿被正確重構.對于張角為70°、55°和45°的二面角進行重構，分別出現了210°、165°和135°的偽影，這是因為3RT 被當作1RT 處理.

利用HBST 算法、BP 算法、BST 算法和文獻[11]中算法分別對45°和70°二面角進行成像，并計算每種成像結果的相對平均誤差(relative mean error,RME)和計算時間，結果如表3 所示.所有四種算法都是用MATAB 代碼實現的，并在一臺裝有英特爾i5-8300H 中央處理器(CPU)@2.3 GHz 和24 GB 隨機存取存儲器(random-access memory, RAM)的計算機上運行.其中RME 表示成像結果與基準圖像相比的相對誤差，其定義如下：

表3 不同算法準確度和計算時間的比較Tab.3 Comparison of accuracy and calculation time of different algorithms

式中：N為成像區像素點的個數；和分別為基準圖像和目標圖像中第j個像素點的邏輯值，當像素點的功率值大于閾值功率時，其邏輯值為1，反之為0.RME 越小，說明目標圖像與基準圖像越相似，即相應的成像算法的精度越高.

從表3 可以看出，由于考慮了多次反射，文獻[11]中的算法與BP 和BST 算法相比具有更高的成像精度，但文獻[11]中的算法需要更長的計算時間，而HBST 算法由于計算量較小，所以速度更快.同時，當選擇合適的閾值時，由于HBST 算法建立了目標回波波前與目標外形間的一一映射關系，因此目標存在的區域對應高反射率，沒有目標的區域反射率為0.而文獻[11]則是將回波信號向目標空間投影，并通過相干疊加來重構目標反射率分布，目標所在位置反射率高，但沒有目標的位置反射率值雖然相對較小，但卻并不為0.因此在此情況下本文算法有更高的精度.

我們用實驗數據對HBST 算法進行進一步驗證.在暗室中進行近場圓周掃描實驗，在矢量網絡分析儀的端口上連接一個喇叭天線作為發射和接收天線.將一個電機連接到轉盤上，控制機械轉盤的旋轉.同時將一個低反射率的泡沫臺放置在轉盤上，并在轉盤上放置要測量的目標.電機和矢量網絡分析儀均由計算機控制.此外，實驗頻率范圍為4～20 GHz，頻率步長為0.02 GHz.

如圖15(a)所示，測試了由多個二面角結構組成的復雜形狀的目標.對目標的下半部分(圖16(a) 中紅線標記部分)進行180°圓周掃描，掃描半徑約為0.7 m；圖15(b)顯示了目標的下半部分的輪廓，中間60°二面角的邊長為0.09 m，而外側兩邊邊長為0.07 m.使用不同的算法對目標進行重構，結果如圖16 所示.結果顯示：BP 和BST 算法將3RT 視為1RT，存在180°偽影；相比之下，本文基于HBST 的算法可以準確地重構目標形狀.

圖15 復雜目標示意圖Fig.15 Schematic diagram of complex target

圖16 使用不同算法對復雜目標底部進行重構的實驗結果Fig.16 Experimental reconstruction results of the bottom part of the complex-shaped target using different algorithms

3 結 論

本文研究了不同張角二面角的回波分布規律，同時發現鏡面反射回波在所有散射回波中的貢獻最大.基于此，我們僅考慮鏡面反射回波對散射場的貢獻，建立了反射回波波前與目標輪廓之間的一一對應關系，提出了能夠處理二面角內多次反射的HBST 算法，實現了二面角的快速成像.但該方法的問題在于對回波波前的提取精度要求較高，通常要求距離分辨率小于1 cm.這是由于波前位置的微小誤差將導致重構點位置的較大偏差，尤其是當最高次反射波前存在偏移時.因此該方法希望成像帶寬盡可能大(大于15 GHz)，以便有更精細的距離分辨率來獲得更精確的波前位置.同時，由于目前的算法主要針對二面角或者多個二面角的組合目標，且當二面角表面較粗糙時(如存在凹陷或凸起)，該算法的重構結果將惡化.因此我們下一步的工作將提升該算法的普適性，以應用于其他存在多次反射的復雜目標.