一種基于MIMO陣列的高爐雷達三維快速成像算法

閆昱君，敬漢丹，孫厚軍，邢光龍,*

(1.燕山大學 信息科學與工程學院，河北 秦皇島 066004；2.北京理工大學 毫米波與太赫茲技術北京市重點實驗室，北京 100081)

0 引言

現代大型高爐屬于耗能巨大的“黑箱”設備，操作人員在無法直接觀測內部料面的條件下進行復雜布料操作，難以實現對冶煉流程的有效優化。對布料表面的實時監測不僅可以優化冶煉流程，提高冶煉質量，同時可以實現節能減排的目標[1]。

基于激光和超聲波的探測方式極易受到高爐內粉塵、溫度、氣流復雜嚴酷環境的影響，由于毫米波的穿透特性，毫米波雷達成像技術可以克服粉塵等影響，是當前高爐料面監測的重要技術發展方向。電子科技大學的丁義元于1988年底研制成功了我國第一部監測高爐料面的毫米波雷達，通過反射鏡的運動實現對料面的機械掃描[2]；北京科技大學的魏紀東、陳先中等設計了沿高爐料面徑向的機械擺動雷達，通過天線的機械擺動實現對料面的掃描[3]。機械掃描雷達結構復雜，在高爐內部復雜環境下，其使用壽命必然大大降低。

近年來，提出了一種新體制雷達——多輸入多輸出(Multiple-Input Multiple-Output, MIMO)雷達。MIMO雷達采用多個發射天線和多個接收天線，各個發射天線同時發射相互正交的信號，各個接收天線對所有的回波信號進行接收并集中處理[4]。MIMO雷達具有實時性好、分辨能力高、抗干擾能力強等諸多優點，被廣泛應用于各個領域。但是，MIMO雷達在高爐料面監測背景下的應用還很有限，已知奧地利的Dominik Zankl等人設計的BLASTDAR雷達已應用于奧鋼聯鋼鐵有限公司5#高爐進行料面成像[5-6]；國內北京航空航天大學許鼎[7]和北京科技大學陳先中[8]設計了特種天線單元，設計了T型MIMO成像雷達用于武鋼7#高爐。已知的應用皆是對MIMO陣列的天線單元進行設計，而在算法層面上，國內外研究學者并沒有提出應用于高爐料面成像場景的針對性成像算法。通過對現有的MIMO雷達成像算法進行研究，發現這些算法應用到高爐料面成像這個大場景下其計算量大、效率低，難以實現對料面的實時成像。

為此，提出了一種基于MIMO陣列的高爐成像雷達距離-角譜域算法。該算法通過對目標的回波信號進行波程近似處理，得到x-y-k域均勻分布的回波數據，對其做收發陣列維快速傅立葉變換和距離維快速傅立葉逆變換，目標在kx-ky-r域實現聚焦，最后通過三維插值將目標圖像從kx-ky-r域轉換至空間直角坐標系。所提出的算法復雜度低并且利用了運算效率極高的快速傅立葉變換，與已有的MIMO陣列近場成像算法相比，運算效率得到了顯著提高。

本文結構安排如下：第1章對所提出的算法進行了理論推導；第2章給出了MIMO陣列的分辨率和陣元間隔的公式；第3章對所提算法進行了仿真驗證；第4章對比了所提算法與傳統的MIMO陣列近場成像算法的運算量；第5章對比分析了所提算法的成像質量；最后，第6章對全文進行了總結。

1 距離-角譜域算法三維成像原理

1.1 三維成像模型

所采用的陣列形式為L型MIMO陣列，設發射陣元為NT個，沿y軸排列，接收陣元為MR個，沿x軸排列，NT個發射陣元和MR個接收陣元組成了NT×MR的虛擬陣列，成像模型如圖1所示。

圖1 三維成像幾何模型Fig.1 The geometric model of three-dimensional imaging

設第mR個接收陣元的位置為

RRx[mR]=[mRdx,0,0],mR∈[0,MR-1]，

(1)

第nT個發射陣元的位置為

PTy[nT]=[0,nTdy,0],nT∈[0,NT-1]，

(2)

其中，dx和dy表示接收和發射陣元的陣元間隔。設點目標標距坐標原點的距離為R，該點目標與坐標原點的連線與x軸、z軸夾角分別為θ、φ，該點目標在三維空間的位置可以表示為

(3)

1.2 回波波程近似

由式(1)～(3)中各收發陣元和點目標在三維直角坐標系中的位置，可以求出點目標的回波波程為

,

將根號內的二次項展開、化簡可以得到

,

將上式按泰勒級數展開，由于目標到天線陣列的距離較遠，省略展開后分母中含有R的高次項，可以得到

RRx+RTy≈2R-sinθcosφmRdx-
sinθsinφnTdy

。

(4)

1.3 三維成像處理

采用的信號形式為步進頻率信號，設信號起始頻率為f0，脈寬為τ，脈沖數(即步數)為N，脈沖間隔(即脈沖重復周期PRT)為T，相鄰脈沖的頻率階躍(即步長)為Δf，第n個脈沖載頻的初始相位為θn，則該發射信號的時域表達式為

exp{j[2π(f0+nΔf)t+θn]}，

設目標的散射系數為σ，則回波信號為

其中，r為目標散射點到雷達天線的距離。回波信號s0(t)和本振信號混頻后，得到視頻回波信號

。

式中，假設目標的散射系數σ為1，相位θn-θ0n表示發射信號和本振信號的初始相位差，對于該雷達成像系統，該項應為常數，可假設為零，因此可簡寫為

(5)

根據式(4)和式(5)可知，在該成像模型下，第(mR,nT)組陣元接收到的回波信號，經過與本振信號混頻處理以及回波波程近似后，其回波信號表達式為

(6)

定義空間頻率(波數)矢量、空間頻率矢量分量以及收發陣元坐標值:

其中，kx和ky可以看作空間頻率矢量在x和y方向上的分量，x和y分別代表接收陣元和發射陣元在x軸和y軸上的坐標值，將上述定義的變量代入式(6)可得

s(x,y,k)=exp(-j2kR)·

exp(jksinθcosφx)exp(jksinθsinφy)。

由于發射陣元和接收陣元在x和y平面上是均勻采樣，步進頻率信號的各個頻點在矢量k方向也是均勻的，因此通過對上式進行分析，可以發現上式的回波數據在x-y-k域是均勻分布的，對該回波數據沿發射陣列維和接收陣列維做快速傅立葉變換(FFT)，將回波數據在kx維和ky維壓縮至目標點位置：

FFT2D[s(x,y,k)]=

δ(kx-ksinθcosφ,ky-ksinθsinφ)exp(-j2kR)≈

δ(kx-k0sinθcosφ,ky-k0sinθsinφ)exp(-j2kR)，

(7)

式中，k0代表中心頻率所對應的波數矢量。經過對回波數據的二維快速傅立葉變換后，在kx-ky二維截面中，數據被壓縮至(k0sinθcosφ,k0sinθsinφ)處呈一沖擊響應函數δ。以位于(0，0，R)處的目標為例，沿收發陣列維FFT后，其(x，y，k)域和(kx，ky，k)域分布如圖2所示。

圖2 沿收發陣列維FFTFig.2 FFT along the transceiver array dimension

經過收發陣列維快速傅立葉變換后，得到沿著k方向的一串沖擊響應函數，在實際的三維空間中k與kx、ky并不是正交關系，而是在目標與原點之間距離的R方向上，因此通過對快速傅立葉變換之后的回波數據沿距離向進行快速傅立葉逆變換就可以將回波數據壓縮至目標點，從而得到目標在kx-ky-r域中的圖像：

IFFT1D{FFT2D[s(x,y,k)]}=P(kx,ky,r)≈

δ(kx-k0sinθcosφ,ky-k0sinθsinφ,r-R),

式中，P代表目標點的坐標函數，由此可知通過對回波數據的二維快速傅立葉變換和一維快速傅立葉變換即可得到目標在kx-ky-r域中的坐標，目標的位置由θ、φ、R決定。但是，目標在kx-ky-r域中的圖像并不能直觀地看出目標點在三維空間直角坐標系中的形狀分布和真實位置，因此對處理后的回波數據進行三維插值處理，將其插值到空間直角坐標系：

。

綜上所述，對波程近似后的理想回波基帶信號做收發陣列維快速傅立葉變換、距離維快速傅立葉逆變換和三維插值即可得到目標點在直角坐標系中的圖像。

上述對算法的推導是在回波波程近似后的理想回波模型下進行的，該模型下的回波波程做了相關近似處理，省略掉了分母中含有R的高次項。實際中回波模型的回波波程未做泰勒展開以及相應近似處理，含有分母中與R有關的高次項。在高爐料面這個大成像場景下，天線陣列與料面的距離較遠，因此1.2節中省略掉的分母中含有R的高次項的值是非常小的，對目標成像結果的影響可以忽略不計，因此在實際回波模型下，所提出的距離-角譜域成像算法依然可以對目標進行準確的快速三維成像。

距離-角譜域算法的實現流程如圖3所示。

圖3 距離-角譜域算法流程Fig.3 Flowchart of distance-angle spectral domain algorithm

2 分辨率與陣元間隔

2.1 分辨率

在MIMO陣列成像體制中，距離向的分辨能力取決于發射信號的帶寬，方位向和高度向的分辨能力是由MIMO陣列在方位向和高度向進行采樣從而實現的。令δz表示成像系統在距離向的分辨率，δx與δy分別表示成像系統在方位向和高度向的分辨率，各方向上分辨率的計算公式為

(8)

其中，c代表光速，B代表步進頻率信號的等效帶寬，N代表步進頻率信號的步數，Δf代表步進頻率信號的步長，λc代表步進頻率信號中心頻率所對應的波長，R代表目標點到陣列平面的垂直距離，LT_x和LT_y是發射陣列分別沿x方向和y方向的陣列孔徑寬度，LR_x和LR_y是接收陣列分別沿x方向和y方向的陣列孔徑寬度。通過觀察上式，得知距離向的分辨率由信號的帶寬決定，帶寬越大，則在距離向的分辨性能越好；在系統中心頻率已定的情況下，方位向的分辨率由陣列在該向的孔徑寬度決定，陣列在該向的孔徑寬度越大，則其分辨性能越好。

2.2 陣元間隔

MIMO雷達成像性能不僅僅取決于各維度的分辨率，其在很大程度上取決于MIMO陣列的陣元間隔。如果陣元間隔過大，將會導致成像混疊，無法準確地確定目標位置；如果陣元間隔過小，將會增加系統的復雜度，從而使成像所需要的時間以及實際中的系統成本增加。因此，MIMO陣列應該滿足一定的陣元間隔條件，使得成像結果在不混疊的前提下最大限度地減少陣元數。下面給出陣元間隔所需要滿足的條件。

接收陣元間距dx滿足

，

(9)

發射陣元間距dy要滿足

，

(10)

其中，目標區域寬度為D，發射陣列孔徑寬度為LT，接收陣列孔徑寬度為LR，R為目標到MIMO陣列的距離，λmin為信號頻帶內所對應的最小波長。式(9)和式(10)即為陣元間隔應滿足的條件。由此可知，收發陣元間隔與信號波長、目標區域寬度、陣列孔徑大小和目標到MIMO陣列的距離有關。

3 仿真結果與分析

采用MATLAB對所提算法在理想回波模型和實際回波模型下分別進行仿真驗證。建立點目標，在理想回波模型下對其進行成像仿真，驗證理想回波模型下上述算法具有良好的成像性能；對同樣的點目標在實際回波模型下對其進行成像仿真并與理想回波模型下的圖像做對比，驗證實際應用中所提算法的成像性能和運算效率；最后建立料面模型，在實際回波模型下對其進行成像仿真，驗證算法的圖形重建性能和運算效率。所獲得的仿真圖像均是通過使用MATLAB 2018a進行仿真得到的。考慮到毫米波對高爐內部粉塵具有良好的穿透效果，選取系統中心頻率為77 GHz；考慮到爐內焦炭顆粒尺寸，選取分辨率為0.15 m；考慮到現代高爐料面尺寸，選取成像區域大小為10 m×10 m×4 m。根據上述指標，在圖像不混疊的前提下，根據式(17)～(19)可計算得到系統仿真的各項參數如表1所示。

表1 系統仿真參數Tab.1 Parameters of system simulation

3.1 理想回波模型下算法成像仿真

建立3個散射點目標，目標在空間直角坐標系中的位置分別為(-2，2，5)、(0，0，6)、(2，-2，7)。在近似后的理想回波模型下，對3個散射點目標進行成像仿真，驗證算法在理想回波模型下的成像性能，目標的成像結果如圖4所示。

由圖4(a)可以清晰地看出3個散射點目標在空間直角坐標系中的位置圖像。由圖4(b)、(c)、(d)點目標在各維度的最大值投影圖，可以直觀地看出在理想回波模型下，所提算法在3個維度均具有良好的目標圖像分辨性能和較低的旁瓣水平，驗證了第1章對所提算法在理想回波模型下實現的理論推導是正確的。

圖4 點目標成像仿真(理想的)Fig.4 The imaging simulation of point target(ideal)

3.2 實際回波模型下算法成像仿真

為了驗證距離-角譜域算法在實際回波模型下成像準確性并與理想回波模型下作對比，建立3個與3.1中位置相同的散射點目標。在實際回波模型下，采用所提算法對其進行成像仿真，使用MATLAB中的計時函數，給出距離-角譜域算法在實際回波模型下的算法處理時間，驗證該成像算法處理具有較高的運算效率。目標的成像結果如圖5所示。

圖5 點目標成像仿真(實際的)Fig.5 The imaging simulation of point target(actual)

將圖5的仿真結果與圖4的仿真結果做比較，發現所提出的距離-角譜域算法在實際回波模型下的成像結果同理想回波模型下的成像結果幾乎完全一致，說明在實際回波模型下，所提出的距離-角譜域算法在該大場景下可以準確地對目標點進行成像。通過使用MATLAB中的計時函數，得到該算法對10 m×10 m×4 m的大場景成像所花費的時間約為0.397 s。

為了更加直觀地展現該成像算法在實際回波模型下對高爐料面的成像性能和運算效率，在這建立一個高爐料面模型，采用距離-角譜域算法對其仿真成像。通過查閱文獻得知高爐煉鐵時料線的形狀為V型，其料面在三維空間中可近似看成一個錐面。在10 m×10 m×4 m的三維空間中，建立一個料面模型，其在空間直角坐標系中的位置如圖6所示。

圖6 料面目標三維空間坐標Fig.6 3D coordinate map of the surface target

使用多個散射點目標構建如上圖所示的料面目標，對其仿真成像，得到成像結果如圖7所示。

由圖7可以清晰直觀地看到料面目標在三維空間中的形狀與位置，與圖6所建立的目標模型對比，說明該算法在大場景下對料面具有良好的圖像重建性能。通過使用MATLAB中的計時函數，得到使用該算法在10 m×10 m×4 m的空間中對所建立的料面目標模型進行成像處理所需要的時間僅為0.438 s，說明該算法在對料面目標模型進行三維成像時具有較高的運算效率。

圖7 料面目標成像仿真Fig.7 The imaging simulation of surface target

4 運算量與分析

為了分析所提出的距離-角譜域算法的運算量并與MIMO雷達反向投影算法(BP算法)[9]和距離徙動算法(wk算法)[10]做比較，首先對BP算法、wk算法和距離-角譜域算法進行逐步分解，分析每一步所需要的浮點運算量，并分別給出3種算法的總運算量，然后計算了具體參數下3種算法的總運算量，直觀地看出3種算法的運算量大小，最后定性分析所提算法與其余兩種算法相比的優勢。

4.1 BP算法、wk算法和所提算法運算量

3種算法大部分的運算量都來自于以下運算：

? FFT和IFFT：長度為N的FFT或IFFT的浮點運算量為5Nlog2(N)[11]。

? 相位相乘：一次復數相位相乘需要6次浮點運算[11]。

? 插值：一次一維線性插值需要4次浮點運算，一次三維線性插值需要28次浮點運算。

根據上述運算開銷，分析并給出了3種算法每一步處理所需要的浮點運算量以及總的運算量。設收發陣元數分別為M、N，步進頻率信號步數為k。為了便于3種算法運算量的對比，需要對數據量和三維成像空間網格進行統一。對于wk算法，為了使該算法的成像區域擴充到大場景，需要在3個維度上按陣元間隔將回波數據補零至所需成像場景大小。為了便于3種算法運算量的對比，在這里將BP算法的像素點和所提算法的三維插值點均按陣元間隔大小劃分，即設BP算法三維像素點數、wk算法補零至所需成像場景后的回波數據點數和距離-角譜域算法的三維插值點數均為x、y、z。

BP算法運算量如表2所示。

表2 BP算法運算量Tab.2 Computation quantity of BP algorithm

wk算法運算量如表3所示。

表3 wk算法運算量Tab.3 Computation quantity of wk algorithm

距離-角譜域算法近似實現運算量如表4所示。

表4 距離-角譜域算法運算量Tab.4 Computation quantity of distance-angle spectral domain algorithm

4.2 具體參數下三種算法總運算量

為了更加直觀地看出3種算法的運算量大小，令收發陣元數M和N均為80，步數k為100，x、y、z均為2 000，通過計算可以得知3種算法的總運算量如表5所示。

表5 3種算法總運算量Tab.5 Total computation quantity of three algorithms

由表5可以直觀地看出，所提出的距離-角譜域算法的總運算量少于傳統的BP算法和wk算法的總運算量。

4.3 算法運算量對比與分析

MIMO雷達BP算法是在時域上對每個像素點進行相干疊加，逐個像素點聚焦的串行處理方式導致BP算法的運算效率非常低。MIMO雷達wk算法是在波數域上對雷達回波數據進行相應成像處理，雖然使用了運算效率高的FFT和IFFT處理，但是其弊端在于在大場景下為了擴充成像區域，必須將雷達的回波數據按陣元間隔來補零至所需成像場景大小，這樣將導致最后的成像像素點數過多，存在大量的冗余。

所提出的距離-角譜域算法中的插值點是按陣元間隔大小來劃分的，這是為了與wk算法的運算量做對比，而在實際應用中，無須按照陣元間隔大小來劃分插值點，只需要將插值點間隔設為分辨率大小即可使圖像具有良好的性能，因此所提出的距離-角譜域算法在實際應用中其運算量遠遠小于BP算法和wk算法。

5 成像質量與分析

BP算法是最標準的用于成像質量對比的算法，其成像精度高，常用于陣列成像質量的標準化分析；wk算法在波數域插值時存在近似，有誤差，因此不適合用來做標準化分析。為了驗證所提出的算法具有良好的成像質量，在實際回波模型下將其與BP算法的成像質量進行對比與分析。

為了避免多個點目標之間的相互影響，采用單個點目標，對目標位于成像區域中心處和邊緣處兩種情況，使用所提出的算法和BP算法在表1的系統仿真參數下對點目標進行成像仿真，對比分析兩種情況下算法的成像精度和旁瓣水平。

5.1 目標位于成像區域中心處成像質量

在成像區域中心處建立單個點目標，其在直角坐標系中的位置為(0，0，6)，使用所提出的算法對其進行成像，各維度成像結果如圖8所示。

在成像區域中心處，建立同樣的單個散射點目標，在相同的仿真參數下，使用BP算法對其進行成像仿真，得到各個維度的成像仿真結果如圖9所示。

由圖8和圖9，可知兩種算法在成像區域中心處各維度的成像誤差和旁瓣水平，如表6所示。

表6 成像質量對比(中心)Tab.6 Comparison of imaging quality(center)

由表6可知，在成像區域中心處，所提出的算法和BP算法在各維度的成像精度都非常高，誤差為0；在旁瓣水平上，所提出的算法在x維和y維略高于BP算法，在z維略低于BP算法。

5.2 目標位于成像區域邊緣處成像質量

在成像區域邊緣處建立單個點目標，其在直角坐標系中的位置為(3，3，7)，使用所提出的算法對其進行成像，各維度成像結果如圖10所示。

使用BP算法對其進行成像仿真，得到3個維度的成像結果如圖11所示。

由圖10和圖11，可以直觀清晰地看出兩種算法在成像區域邊緣處對目標成像時，各個維度的成像誤差和旁瓣水平，對比如表7所示。

表7 成像質量對比(邊緣)Tab.7 Comparison of imaging quality(edge)

由表7可以直觀地看出，當目標在成像區域邊緣處時，所提出的距離-角譜域算法和BP算法在3個維度的成像精度都非常高，誤差為0；在旁瓣水平上，所提出的距離-角譜域算法在3個維度上略高于BP算法。

綜上所述，所提出的距離-角譜域算法在該成像場景下具有良好的成像質量。

圖8 距離-角譜域算法各維度成像仿真(中心)Fig.8 The imaging simulation of distance-angle spectraldomain algorithm in all dimensions(center)

圖9 BP算法各維度成像仿真(中心)Fig.9 The imaging simulation of BP algorithm in all dimensions(center)

圖10 距離-角譜域算法各維度成像仿真(邊緣)Fig.10 The imaging simulation of distance-angle spectraldomain algorithm in all dimensions(edge)

圖11 BP算法各維度成像仿真(邊緣)Fig.11 The imaging simulation of BP algorithm inall dimensions(edge)

6 結論

本文提出了一種基于MIMO陣列的大場景快速成像算法——距離-角譜域算法。對該算法在理想回波模型和實際回波模型下分別進行成像仿真，驗證了該算法在保證了成像性能的前提下具有較高的運算效率。文中還對比分析了距離-角譜域算法與已有的MIMO陣列近場成像算法的運算量和成像質量，體現了所提算法的高運算效率和良好的成像質量。下一步將開展MIMO陣列近場毫米波成像系統的工程化研究工作，實現所提算法的工程化應用，對算法進行進一步的驗證和改進。