基于迭代自適應方法的柱面共形陣2D DOA估計

張狀和 韓東 劉德亮

摘? 要： 針對柱面共形陣的二維波達方向（2D DOA）估計問題，提出一種基于迭代自適應方法的柱面共形陣2D DOA估計算法。首先，通過歐拉旋轉變換公式將柱面共形陣全局直角坐標系變換為單元局部直角坐標系，建立柱面共形天線陣列流型;然后，將感興趣區域內信號源所有可能存在的位置通過劃分二維網格的形式來表示，利用迭代自適應算法（IAA）估計出每個潛在位置對應的信號源能量;最后，繪制能量譜圖，譜峰對應的方位角和俯仰角即為二維波達方向的估計值。理論分析與仿真實驗證明，所提算法不需要額外的參數匹配，簡化了估算步驟，在短快拍的情況下可以實現柱面共形陣的2D DOA估計。

關鍵詞： 柱面共形陣; 2D DOA估計; 坐標系變換; 信號源能量估計; 迭代自適應算法; 能量譜繪制; 理論分析

中圖分類號： TN911.7?34? ? ? ? ? ? ? ? ? ? ? ?文獻標識碼： A? ? ? ? ? ? ? ? ? ? ? ? ? 文章編號： 1004?373X（2020）11?0006?04

Cylindrical conformal array′s 2D DOA estimation based on iterative adaptive approach

ZHANG Zhuanghe， HAN Dong， LIU Deliang

（Department of Missile Engineering， Army Engineering University of PLA （Shijiazhuang Campus）， Shijiazhuang 050003， China）

Abstract： An iterative adaptive approach （IAA） based 2D DOA （two dimensional direction of arrival） estimation algorithm for cylindrical conformal array is proposed to deal with the problem of 2D DOA estimation of cylindrical conformal array. The global rectangular coordinate system of the cylindrical conformal array is transformed into the unit local rectangular coordinate system by Euler rotation transformation formula for the establishment of array manifold of the cylindrical conformal antenna array. All possible locations of signal source in the ROI （region of interest） are represented by dividing 2D grids， and the IAA is used to estimate the signal source energy corresponding to each potential location. An energy spectrum map is drawn， and the azimuth angle and pitch angle corresponding to the spectrum peak are the estimated values of the 2D DOA. Theoretical analysis and simulation experiments prove that it is unnecessary to match additional parameters in the application of the proposed algorithm， so the steps of estimation is simplified and 2D DOA estimation of cylindrical conformal array can be realized under the condition of short snapshot.

Keywords： cylindrical conformal array; 2D DOA estimation; coordinate system transformation; signal source energy estimation; iterative adaptive approach （IAA）; energy spectrum drawing; theoretical analysis

0? 引? 言

共形陣天線是一種和載體外形保持一致的天線[1]，其能夠和物體結構相融合，可以滿足飛機、導彈、衛星等飛行器的隱身和空氣動力學需求，有效減少了荷載的大小和重量。相對于傳統的線陣和平面陣，共形陣天線還具有覆蓋范圍大、易偽裝、抗干擾能力較強等特點。近些年來，學者們圍繞著共形陣列展開了多方面的研究[2?7]，其中，由于波達方向（Direction of Arrival，DOA）估計技術在通信、雷達、聲吶等領域的重要作用，因此如何實現共形陣列DOA估計受到廣泛關注。

文獻[2]研究了多重信號分類（Multiple Signal Classification，MUSIC）算法在共形天線中的應用。文獻[3]通過虛擬陣列變換技術，應用旋轉不變性技術估計信號參數（Estimation of Signal Parameters Via Rotational Invariance Techniques，ESPRIT）和MUSIC算法，克服了共形載體的遮蔽效應，實現了基于共形陣的2D DOA（Two Dimensional Direction of Arrival，2D DOA）估計。但是以上算法忽略了載體曲率對方向圖綜合的影響，簡化了陣列流型的數據模型。文獻[4?5]通過運用歐拉旋轉變換公式，提出了共形天線陣列流型的通用建模方法。文獻[6]通過合理設置陣元，并基于ESPRIT，解決了柱面共形陣的盲極化2D DOA估計。文獻[7]則引入了虛擬期望信號，利用最大信干噪比準則下的最優權矢量對接收信號進行加權處理，實現了共形陣DOA估計。上述文獻均將經典的子空間類算法應用到了共形陣列DOA估計問題上，但是這些算法都需要大量的快拍數據支撐，以保證對樣本協方差矩陣進行特征分解時，信號子空間與噪聲子空間不發生混疊，同時獲得高精度、高分辨率的估計結果。但是在實際的工程應用中，信號經常處于不能長時間穩定或快速時變的環境中，如高速目標追蹤、水下信號處理、城市無線通信等，采集大量的快拍數據會導致處理時間過長，增大與真實樣本的誤差或降低運算速度等問題。而近些年來，隨著稀疏信號處理理論的興起，尤其以Petre Stoica等提出的振幅和相位估計的迭代自適應方法（Iterative Adaptive Approach for Amplitude and Phase Estimation， IAA?APES）[8?9]、基于協方差稀疏迭代的譜估計方法（Sparse Iterative Covariance?based Estimation Method， SPICE）[10]為代表，使得短快拍DOA估計的實現成為可能。

為了能在短快拍條件下實現共形陣列天線的高精度2D DOA估計，本文提出了一種基于迭代自適應方法的柱面共形陣2D DOA估計算法。首先，利用歐拉旋轉變換公式，實現柱面共形陣天線全局直角坐標系到單元局部直角坐標系的變換，建立柱面共形天線陣列流型;然后，基于加權最小二乘法估計出入射信號源的能量，通過循環迭代對估計結果進行更新直至收斂;最后，繪制出三維能量譜圖，譜峰對應方位角和俯仰角即為入射信源的2D DOA。仿真實驗表明，在只有2個快拍的情況下，所提算法經過適量迭代就能實現對入射信號源2D DOA的高精度、高分辨率估計。

下面對文中的符號和算子進行說明：[[?]*]表示復數共軛;[[?]T]表示轉置;[[?]H]表示共軛轉置;[?]表示克羅內克積;[x]表示[x]的估計值;[diag（? ）]表示對角化;[E[?]]表示期望值;[?]表示Euclidean范數。

1? 柱面共形陣數據模型

柱面共形陣列天線結構如圖1所示。[M×N] 個陣元均勻地分布在載體表面，定義[n]為由下至上的陣元序號，[m]為每一圓環上逆時針方向的陣元序號，[R]表示圓柱半徑，[H]表示圓柱的高，[d]為上下兩陣元間的距離，[β] 為相鄰兩陣元間的夾角，[?]為圓環末陣元與[x]軸的夾角;并定義坐標系[[x，y，z]]表示陣列全局直角坐標系，[[x，y，z]] 表示陣元局部直角坐標系，[[r，θ，φ]]表示陣列全局極坐標，[[r，θ，φ]]表示陣元局部極坐標。

給出陣元坐標以及歐拉旋轉變換角公式[4]：

[xnm=Rcos[（m-1）β-?]] （1）

[ynm=Rsin[（m-1）β-?]] （2）

[znm=nd-（N+1）d/2] （3）

[Dnm=π-?+（M-1）βEnm=-π2Fnm=0] （4）

假設有[K]個遠場窄帶信號源入射到圖1所示柱面共形陣上，[θi]，[?i]分別表示第[i]個入射信號源的方位角與俯仰角，結合文獻[4?5]，柱面共形陣的導向矢量模型可以表示為：

[a（θ，?）=[r1e-j2πp1?uλ，r2e-j2πp2?uλ，…，rNMe-j2πpnm?uλ]T] （5）

[ri=（g2iθ+g2i?）12（k2iθ+k2i?）12cos θigk=giplcos θigk=gi?pl=giθkθ+gi?k?] （6）

式中：[pnm=[xnm，ynm，znm]]（[n=1，2，…，N]，[m=1，2，…，M]）表示第[nm] 個陣元的坐標矢量;如圖2所示，[u=[cos θsin ?，sin θsin ?，cos ?]T]定義為傳播矢量;[ri]表示第[i]個陣元對單位信號的響應;[gi]為陣元方向圖表達式;[giθ，gi?]分別由信源構成的矢量基[uθ，u?]所構成;[θigk]為[gi]和[pl]向量之間的夾角。

考慮到入射信號源均相互獨立，噪聲為加性高斯白噪聲，因此快拍數據模型可以表示為：

[X（n）=AS（n）+N（n）] （7）

[S（n）=[s1，s2，…，sn]T] （8）

[N（n）=[n1，n2，…，nn]T] （9）

[A（θ，?）=[a（θ1，?1），a（θ2，?2），…，a（θNM，?NM）]] （10）

式中：[A（θ，?）]表示導向矩陣，由導向矢量構成;[S]表示信號矢量;[N]表示加性高斯白噪聲矢量。通常情況下，實際感興趣的目標信號源數大于掃描區域內信號源的個數，導致[S（n）]中只有一小部分的數值是非零的，因此稀疏類算法適用于DOA估計。

2? 柱面共形陣2D DOA估計算法

首先，構造[K×K]的能量矩陣[P]，其對角線元素包含了掃描網格上每個角度的信號能量。能量矩陣的結構如下：

[P=diag{p1，p2，…，pi}] （11）

[pi=1Nn=1Nsi（n）2，? ? i=1，2，…，K]? （12）

定義噪聲和干擾的協方差矩陣如下：

[Q（θi，?i）=R-pia（θi，?i）aH（θi，?i）]? ?（13）

式中[R=A（θ，?）PAH（θ，?）]。

加權最小二乘法的代價函數由下式給出：

[n=1Nx（n）-si（n）a（θi，?i）2W] （14）

式中：[y2W=yHWy];[si（n）]是第[i]個遠場信號源在快拍數[n]處的復幅值;[Q-1（θi，?i）]是噪聲的協方差矩陣，由文獻[7]可知，當[W=Q-1（θi，?i）]時[si（n）]的值取最小，因此，對式（14）求關于[si（n）]的最小化，表達形式為：

[si（n）=aH（θi，?i）Q-1（θi，?i）x（n）aH（θi，?i）Q-1（θi，?i）a（θi，?i）， n=1，2，…，N] （15）

式中：[Q-1（θi，?i）]可由[R-1]代替，式（15）可以精簡為：

[si（n）=aH（θi，?i）R-1x（n）aH（θi，?i）R-1a（θi，?i），? ?n=1，2，…，N] （16）

式（16）得到了[（θi，?i）]處對應信源復幅值的估計值。將[si（n）]代入式（12）即可計算第[i]個掃描點的能量。

算法流程總結如下：

初始化：

[s（0）i（n）=aH（θi，?i）x（n）a（θi，?i）2，? ? i=1，2，…，K] （17）

[p（0）i=1Nn=1Ns（0）i（n）2，? ? i=1，2，…，K] （18）

[P（0）=diag{p（0）1，p（0）2，…，p（0）K}] （19）

迭代過程如下：

1） 構造能量矩陣[p（l）i=1Nn=1Ns（l-1）i（n）2]，[i=1，][2，…，K];[p（l）i]是第[i]個網格點對應的信號源在第[l]次迭代的能量估計值。

2） 計算信號協方差矩陣[R=A（θ，?）P（l）AH（θ，?）]。

3） 更新[s（l）i（n）=aH（θi，?i）R-1x（n）aH（θi，?i）R-1a（θi，?i）]。

4） 當迭代收斂，則停止迭代，否則返回式（1），根據經驗，一般迭代10～15次就會得到較好的估計值。

計算能量值：迭代結束后，以最后一次迭代得到的[p（end）i]作為第[i]個掃描點對應的譜峰函數繪制出三維譜峰圖，譜峰對應的方位角、俯仰角即為真實信源的二維波達角度。

3? 仿真實驗與分析

為驗證本文所提方法的有效性，對比2D?MUSIC算法展開實驗。為了不失一般性，設置如圖1所示的接收陣列，[M=5]，[N=10]，相鄰兩圓環的間距為入射信號波長的一半，即[d=λ/2];噪聲為加性高斯白噪聲，與入射信號源相互獨立;信號源個數[K=2]，且相互獨立;[θ1=55°]，[?1=30°]，[θ2=60°]，[?2=35°];假設在DOA估計之前已對感興趣的區域有一定認知，所以提前設定方位角度掃描范圍為[[50°，65°]]，俯仰角度掃描范圍為[[25°，40°]]，掃描間隔為[0.5°];迭代次數預設[l=10]。單元方向圖設置為[4]：

當[0≤θ≤π2]時：

[gθ（θ，?）=J2πdλsin θ-J0πdλsin θ×（cos ?-jsin ?）]? （20）

[g?（θ，?）=J2πdλsin θ+J0πdλsin θ×cos θ（sin ?-jcos ?）] （21）

其他時，

[gθ（θ，?）=g?（θ，?）=0]

式中[J0，J2]分別是零階和二階第一類貝塞爾函數。

3.1? 仿真實驗1

設置信噪比（Signal to Noise Ratio，SNR）為15 dB，快拍數[N1=2，N2=200]，能量譜分別如圖3，圖4所示。

圖3，圖4仿真結果表明，當快拍數足夠大時，兩種算法都能成功地把三個近場信源分辨出來;而當快拍數減小到[N1=2] 時，2D?MUSIC算法已經失效，真實信號源的譜峰混疊在偽峰之中，已無法分辨出真實入射角度，但所提算法仍然可以分辨出兩個銳利的譜峰，且旁瓣水平較低，方位角和俯仰角的估計偏差也很小。

3.2? 仿真實驗2

定義實驗成功須滿足下式：

[（θ）-θ≤3°，? ? （?）-?≤3°] （22）

圖5設置快拍數[N2=200]，信噪比從-15 dB變化到15 dB。圖6設置信噪比為10 dB，快拍數從10變化到200，其他實驗條件不變。200次蒙特卡洛實驗觀察分辨概率和均方根誤差隨快拍數的變化曲線。

從圖5的實驗結果可以看出，所提算法的估計均方根誤差要小于2D?MUSIC算法，說明所提算法的估計精度更高。從圖6的實驗結果可以看出，當快拍數大于70時，所提算法的分辨概率趨于100%，而2D?MUSIC算法的分辨概率依舊很低。

4? 結? 語

針對柱面共形陣的2D DOA估計問題，提出一種基于IAA的估計算法。本文利用歐拉旋轉變換公式，將柱面共形陣天線全局直角坐標系變換到單元局部直角坐標系，討論了稀疏類算法在DOA估計應用上的原理，使用IAA算法解決了短快拍條件下的柱面共形陣的2D DOA估計問題，計算機仿真驗證了所提算法的有效性。

參考文獻

[1] JOSEFSSO L， PERSSON P. Conformal array antenna theory and design [M]. New Jersey： John Wiley & Sons， Inc.， 2006.

[2] 孫翠珍，毛昕蓉.多重信號分類算法在共形陣列天線中的應用[J].計算機仿真，2014，31（11）：171?174.

[3] 王瑞革，王法棟.基于虛擬陣列變換的共形陣列信號DOA估計[J].雷達科學與技術，2016，14（4）：448?452.

[4] 王布宏，郭英，王永良，等.共形天線陣列流形的建模方法[J].電子學報，2009，37（3）：481?484.

[5] 李龍軍，王布宏，夏春和.稀疏共形陣列天線方向圖綜合[J].電子學報，2017，45（1）：104?111

[6] 張羚，郭英，齊子森，等.柱面共形陣列天線盲極化2D DOA估計[J].空軍工程大學學報（自然科學版），2016，17（3）：78?84.

[7] 楊群，曹祥玉，高軍，等.一種未知信源數的共形陣DOA算法[J].西安電子科技大學學報（自然科學版），2015，42（3）：122?128.

[8] YARDIBI T， LI J， STOICA P， et al. Source localization and sensing： a nonparametric iterative adaptive approach based on weighted least squares [J]. IEEE transactions on aerospace and electronic system， 2010， 46（1）： 425?443.

[9] STOICA P， LI J， HE H. Spectral analysis of nonuniformly sampled data： a new approach versus the periodogram [J]. IEEE transactions on signal processing， 2009， 57（3）： 843?858.

[10] STOICA P， BABU P， LI J. New method of sparse parameter estimation in separable models and its use for spectral analysis of irregularly sampled data [J]. IEEE transactions on signal processing， 2011， 59（1）： 35?47.