做陣列測向這幾年我見過太多方法在仿真里漂亮一到實測就崩盤。MUSIC的分辨率確實讓人著迷但現(xiàn)實中的相干源、單快拍、低信噪比隨便來一個都能讓子空間類方法啞火。后來接觸到了IAA迭代自適應(yīng)方法才真正理解了什么叫“非參數(shù)、不需要先驗信息也能穩(wěn)著出譜”。不過IAA論文里那套推導(dǎo)符號極其勸退我當(dāng)年也是抄完代碼后好一陣子才想明白它本質(zhì)上是一條從加權(quán)最小二乘WLS出發(fā)通過迭代更新權(quán)重來逼近DOA估計的路線。這篇文章就把這條路線完整走一遍適合正在啃IAA論文、或者想在DOA估計中引入IAA的同學(xué)作為參考。1. 從問題說起DOA估計到底在解什么1.1 陣列信號模型先把符號統(tǒng)一DOA估計說白了就是回答一個問題來波方向在哪。我們假定接收端是一個由M個陣元組成的陣列以均勻線陣為例陣元間距通常是半個波長。當(dāng)某個遠(yuǎn)場窄帶信號以角度 $\theta$ 入射時第m個陣元相對參考陣元會有一個相位差把所有陣元的相位差寫成向量就是這個角度對應(yīng)的導(dǎo)向矢量$$a(\theta) \begin{bmatrix} 1 \ e^{j2\pi\fracgv5xshg9wyt{\lambda}\sin\theta} \ \vdots \ e^{j2\pi(M-1)\fracgv5xshg9wyt{\lambda}\sin\theta} \end{bmatrix}$$如果有K個信號同時入射每個快拍的接收數(shù)據(jù)可以寫成$$y(n) A(\theta)s(n) e(n)$$其中 $A(\theta)$ 是 $M \times K$ 的陣列流形矩陣它把K個導(dǎo)向矢量按列拼接起來$s(n)$ 是K個信號的復(fù)包絡(luò)$e(n)$ 是噪聲通常建模為高斯白噪聲。整個DOA估計問題就是從一段觀測 $y(1), \dots, y(N_t)$ 中把 $\theta_1, \dots, \theta_K$ 找出來。注意這里我故意沒有說“先估計信源個數(shù)再估計角度”因為實際工程里信源個數(shù)這個參數(shù)本身就不容易給定。IAA的高明之處就在于它跳過了“先定階、再測向”這個傳統(tǒng)套路。1.2 傳統(tǒng)方法的翻車現(xiàn)場最樸素的測向方法是延遲求和也叫常規(guī)波束成形。它就是把陣列數(shù)據(jù)往每個候選方向上投影哪個方向投影能量大就認(rèn)為哪個方向有信號$$P_{CBF}(\theta) a^H(\theta)\hat{R}_y a(\theta)$$這里 $\hat{R}_y \frac{1}{N_t}\sum_n y(n)y^H(n)$ 是采樣協(xié)方差矩陣。這個方法的優(yōu)點是穩(wěn)健、計算量小缺點是分辨率被“瑞利限”卡死。陣列孔徑就那么長波束主瓣寬度擺在那里兩個角度距離小于波束寬度時這個方法是分辨不出來的。接著大家就想到了Capon波束形成也就是最小方差無失真響應(yīng)。它在保持期望方向增益為1的前提下最小化輸出功率相當(dāng)于把干擾方向“壓”出零點。這個思路在理論上是漂亮的但問題也明顯它需要求協(xié)方差矩陣的逆快拍不足或者協(xié)方差矩陣病態(tài)時結(jié)果會劇烈抖動對角加載那個參數(shù)又得反復(fù)試。然后是MUSIC。MUSIC把接收數(shù)據(jù)分解成信號子空間和噪聲子空間利用兩個子空間的正交性構(gòu)造譜峰。這個方法分辨率高但有幾個致命問題必須先知道信源個數(shù)相干信號會導(dǎo)致信號子空間“降秩”譜峰直接消失低信噪比下噪聲子空間估計不準(zhǔn)性能掉得很快。我在實測里最頭疼的就是相干源問題兩個目標(biāo)一旦相干MUSIC幾乎必然只剩一個峰。1.3 IAA想解決的三件事IAA的出現(xiàn)其實是沖著傳統(tǒng)方法的三個痛點去的第一不需要預(yù)先知道信源個數(shù)。IAA是在整個候選角度網(wǎng)格上估計功率譜有能量的角度自然形成峰值信源個數(shù)是“事后數(shù)峰”得到的。第二能處理相干信號。因為IAA不是子空間分解類方法它不依賴信號子空間的秩所以相干源對它來說只是兩個功率點。第三對快拍數(shù)不敏感。推導(dǎo)到后面你就會發(fā)現(xiàn)IAA的核心更新公式在單快拍下也能成立。這在實際場景里非常有用比如機載雷達(dá)、跳頻通信里往往沒有足夠的時間積累多個快拍。當(dāng)然代價也有那就是計算量。IAA需要迭代每次迭代都要對協(xié)方差矩陣求逆網(wǎng)格點數(shù)一多計算負(fù)擔(dān)確實上去了。這塊我在后面實操部分會展開講。2. WLS框架IAA的地基2.1 最小二乘估計的局限要理解IAA得先回到最小二乘這層地基。假設(shè)我們有線性觀測模型$$y As e$$其中 $y$ 是 $M \times 1$ 的觀測向量$A$ 是 $M \times K$ 的已知字典矩陣$s$ 是待估計的信號向量$e$ 是噪聲。經(jīng)典LS估計就是最小化誤差的2范數(shù)平方$$\hat{s}_{LS} \arg\min_s |y - As|_2^2$$解是 $\hat{s}_{LS} (A^H A)^{-1} A^H y$。這個解在噪聲是獨立同分布的白噪聲時是很好的因為每個觀測分量的可信度相同最小化平方和很公平。問題在于實際陣列信號里的噪聲和干擾并不是白化的。換個說法$e$ 的協(xié)方差矩陣 $Q E[ee^H]$ 不一定是單位陣的倍數(shù)。比如存在一個強干擾源時干擾會在某些方向上注入很大的能量此時如果還對所有觀測分量一視同仁地做最小二乘估計結(jié)果會被強干擾帶偏方差非常大。這就是LS在復(fù)雜電磁環(huán)境下的先天不足。2.2 加權(quán)最小二乘的求解過程WLS的思路很直觀既然不同觀測分量的可信度不一樣那就給可信度高的分量更大的權(quán)重給可信度低的分量更小的權(quán)重。數(shù)學(xué)上就是最小化加權(quán)范數(shù)$$\hat{s}_{WLS} \arg\min_s (y - As)^H Q^{-1} (y - As)$$這里的 $Q^{-1}$ 就是權(quán)重矩陣。為什么要用協(xié)方差矩陣的逆而不是直接用協(xié)方差矩陣因為你希望當(dāng)某個方向噪聲功率大時對應(yīng)分量的權(quán)重變小。協(xié)方差矩陣 $Q$ 度量了噪聲的“大小”它的逆天然就把大噪聲壓下去了。對 $s$ 求導(dǎo)并令導(dǎo)數(shù)為零可以得到$$\hat{s}_{WLS} (A^H Q^{-1} A)^{-1} A^H Q^{-1} y$$這個公式是WLS的標(biāo)準(zhǔn)形式。你可以把它理解成先用 $Q^{-1/2}$ 對觀測向量和字典矩陣做一次“白化”然后在白化后的空間里做普通LS。白化讓各個分量的噪聲變成等方差LS在這個空間里又是最優(yōu)的。從直覺上講$Q^{-1}$ 的作用有兩層。一層是噪聲白化讓不同陣元的噪聲功率一樣另一層是干擾抑制如果某個方向有強干擾而這份干擾又反映在 $Q$ 里那么 $Q^{-1}$ 會對準(zhǔn)那個方向形成一個凹陷。理論上當(dāng)干擾的協(xié)方差建模得足夠準(zhǔn)確加權(quán)最小二乘就能把干擾的影響基本消除。2.3 把WLS放進(jìn)DOA估計的語境現(xiàn)在把WLS拿來做DOA估計。我們并不想一次性把所有方向的信號幅度都解出來而是想逐個角度地判斷“這個方向有沒有信號”。所以對每個候選角度 $\theta_k$構(gòu)造一個“單參數(shù)模型”$$y a_k s_k \text{其他信號} \text{噪聲}$$如果“其他信號噪聲”的協(xié)方差矩陣能夠被估計出來比如記為 $R$那么對 $s_k$ 的WLS估計就是$$\hat{s}_k \frac{a_k^H R^{-1} y}{a_k^H R^{-1} a_k}$$這個式子其實是單參數(shù)WLS的特例因為 $A$ 退化成一個向量 $a_k$。同時它也跟Capon波束成形器的輸出完全一致。換句話說一旦你用某個協(xié)方差矩陣 $R$ 定義了權(quán)重那么“自適應(yīng)測向”的核心運算就是上面這個式子。這就是IAA和WLS之間最關(guān)鍵的橋IAA無非是把這個 $R$ 從一個固定常量變成隨迭代不斷更新的量。每一次迭代都在用當(dāng)前對信號功率的估計結(jié)果重新構(gòu)造協(xié)方差矩陣然后再對每個角度做一次WLS估計。權(quán)重矩陣不再是拍腦袋定的而是從數(shù)據(jù)里迭代學(xué)出來的。3. IAA核心推導(dǎo)迭代自適應(yīng)公式手把手推一遍3.1 初始化匹配濾波給的起點IAA的第一步是用匹配濾波初始化每個候選角度的功率。匹配濾波的思路很樸素用每個方向的導(dǎo)向矢量去跟接收數(shù)據(jù)做相關(guān)相關(guān)能量大的方向就可能是信號方向。初始信號幅度估計為$$\hat{s}_k^{(0)} \frac{a_k^H y}{a_k^H a_k}$$初始功率就取模平方$$p_k^{(0)} |\hat{s}_k^{(0)}|^2$$對均勻線陣來說$a_k^H a_k M$所以匹配濾波輸出其實就是常規(guī)波束成形的輸出。這一步雖然分辨率不高但能給出一個大致的功率分布足以支撐第一輪協(xié)方差矩陣的構(gòu)建。如果你處理的是多快拍數(shù)據(jù)初始功率也可以寫成所有快拍平均的結(jié)果$$p_k^{(0)} \frac{1}{N_t}\sum_{n1}^{N_t}|\hat{s}_k^{(0)}(n)|^2$$初始化的質(zhì)量會影響收斂速度但I(xiàn)AA對初始化不算太敏感。我試過用全零以外的多種初始化方式包括直接把 $p_k^{(0)}$ 設(shè)成同樣的常數(shù)迭代十幾輪之后基本都能收斂到相近的結(jié)果。當(dāng)然用匹配濾波初始化是最穩(wěn)、最省事的做法。3.2 協(xié)方差矩陣建模從功率到干擾抑制有了每個角度的初始功率就可以構(gòu)造全面的協(xié)方差矩陣。假設(shè)我們把整個角度域離散成 $K_g$ 個候選格點那么信號協(xié)方差矩陣可以寫成$$R \sum_{k1}^{K_g} p_k a_k a_k^H \sigma I$$其中 $\sigma I$ 是噪聲項實際實現(xiàn)里通常用對角加載來代替保證矩陣可逆。這個 $R$ 的物理含義非常清楚它表示在當(dāng)前功率估計下陣列接收數(shù)據(jù)的協(xié)方差結(jié)構(gòu)。如果某個格點 $k$ 的功率 $p_k$ 很大說明那里大概率有一個信號那么 $R$ 里就包含了來自這個“信號”的貢獻(xiàn)。為什么IAA要用包含所有候選角度的 $R$而不是把當(dāng)前估計的角度 $k$ 本身也去掉嚴(yán)格來說估計第 $k$ 個角度的信號幅度時更“干凈”的做法是用刪除了第 $k$ 個方向貢獻(xiàn)的協(xié)方差矩陣。但那樣每個角度都要單獨構(gòu)造一個不同的逆矩陣計算量不可接受。IAA選擇用統(tǒng)一的 $R$ 來近似代價是當(dāng)前角度自身的功率會混在干擾協(xié)方差里但在高分辨率網(wǎng)格和迭代收斂后這種影響會變得很小。這是IAA在計算量和精確性之間做的巧妙折中。3.3 單角度WLS求解核心公式誕生接下來是整篇文章最核心的推導(dǎo)。對第 $k$ 個候選角度我們希望估計信號幅度 $s_k$。把其他所有格點的貢獻(xiàn)都當(dāng)作“干擾”用當(dāng)前協(xié)方差矩陣 $R$ 來白化它。于是構(gòu)造如下加權(quán)最小二乘問題$$\min_{s_k} \left(y - a_k s_k\right)^H R^{-1} \left(y - a_k s_k\right)$$展開目標(biāo)函數(shù)$$J(s_k) y^H R^{-1} y - y^H R^{-1} a_k s_k - s_k^H a_k^H R^{-1} y s_k^H a_k^H R^{-1} a_k s_k$$令 $\alpha -y^H R^{-1} a_k$$\beta a_k^H R^{-1} a_k$注意 $\beta$ 是一個正實數(shù)因為 $R^{-1}$ 是Hermitian正定矩陣。于是$$J y^H R^{-1} y \alpha s_k \alpha^* s_k^* \beta |s_k|^2$$把 $s_k u jv$ 拆成實部和虛部分別對 $u$ 和 $v$ 求導(dǎo)并令其為零。經(jīng)過整理可以得到$$\hat{s}_k \frac{a_k^H R^{-1} y}{a_k^H R^{-1} a_k}$$這個式子太重要了值得停下來多看兩眼。它的結(jié)構(gòu)是一個歸一化的匹配濾波但匹配空間不是原始的觀測空間而是經(jīng)過 $R^{-1}$ “白化”之后的空間。$R^{-1}$ 在這里同時起到了兩個作用一是壓制其他方向強信號帶來的干擾二是把非白噪聲白化使得最終估計結(jié)果近似最優(yōu)。說白了這就是“自適應(yīng)”二字的來源。第一次迭代時 $R$ 由匹配濾波初始化得到里面的干擾信息還不準(zhǔn)但每迭代一輪$R$ 變得更準(zhǔn)確$R^{-1}$ 對干擾的抑制能力也更強于是 $p_k$ 的估計就更準(zhǔn)反過來又讓下一輪 $R$ 更準(zhǔn)。這就是一個典型的期望最大化式循環(huán)。3.4 功率更新與迭代循環(huán)得到信號幅度估計后功率更新很簡單$$p_k |\hat{s}_k|^2$$多快拍情況下則是把所有快拍的估計結(jié)果取平均$$p_k \frac{1}{N_t}\sum_{n1}^{N_t}|\hat{s}_k(n)|^2$$至此一次完整的迭代就結(jié)束了。整個IAA算法可以濃縮成下面這個循環(huán)初始化用匹配濾波得到每個候選角度的初始功率 $p_k^{(0)}$。構(gòu)建協(xié)方差矩陣$R \sum_k p_k a_k a_k^H \sigma I$。對每個候選角度計算信號幅度$\hat{s}_k \frac{a_k^H R^{-1} y}{a_k^H R^{-1} a_k}$。更新功率$p_k |\hat{s}_k|^2$。檢查收斂如果 $\frac{|p^{(i)} - p^{(i-1)}|_2}{|p^{(i-1)}|_2} \epsilon$停止否則回到第2步。收斂之后把 $p_k$ 按角度畫出來就是IAA的功率譜。譜峰對應(yīng)的角度就是DOA估計結(jié)果。這里有一個值得琢磨的細(xì)節(jié)第3步和第4步之間其實是有內(nèi)在一致性的。如果信噪比很高、干擾抑制得很干凈那么 $\hat{s}_k$ 會非常接近真實信號幅度功率更新自然準(zhǔn)確反過來如果 $R$ 里錯誤地把某個沒有信號的角度的功率設(shè)得很大那么 $R^{-1}$ 就會在那個方向形成一個“坑”下一輪這個方向的功率就會被壓下去。這種負(fù)反饋機制保證了算法不太容易發(fā)散也是IAA穩(wěn)健性的核心保障。4. 從公式到代碼實現(xiàn)中的關(guān)鍵細(xì)節(jié)4.1 初始化與迭代停止條件先說初始化。理論上可以用任何非負(fù)的功率向量啟動但匹配濾波初始化有三個好處計算量小只需要做一次矩陣向量乘物理意義明確等價于常規(guī)波束成形的輸出迭代收斂快因為初始值已經(jīng)離真實功率分布不遠(yuǎn)了。迭代停止條件有兩種常見做法。一種是嚴(yán)格檢查收斂比如設(shè)置 $\epsilon 10^{-3}$ 或 $10^{-4}$每次迭代后計算功率向量變化的相對范數(shù)。另一種更工程化的做法是固定迭代次數(shù)比如統(tǒng)一迭代10到15次。論文里一般認(rèn)為10次左右已經(jīng)能得到很穩(wěn)定的結(jié)果15次以上基本沒有肉眼可見的變化。我自己在MATLAB里跑仿真時通常設(shè)固定15次省去每次判斷收斂的開銷。需要提醒的是收斂閾值不要設(shè)得太苛刻。IAA在迭代后期功率譜的變化幅度非常小但嚴(yán)格收斂可能需要更多輪次計算收益卻不明顯。工程上講與其多跑5輪去追求千分之一的譜變化不如把省下來的算力放到提高網(wǎng)格密度上去。4.2 噪聲項處理與對角加載$R$ 的構(gòu)造里如果完全沒有噪聲項當(dāng)候選格點數(shù) $K_g$ 小于陣元數(shù) $M$ 的時候$A P A^H$ 很可能不滿秩求逆直接出問題。就算 $K_g \geq M$數(shù)值上也可能接近奇異。所以實踐中幾乎都會做對角加載也就是在 $R$ 上加一個 $\sigma I$。$\sigma$ 怎么選這直接決定算法穩(wěn)定性。我的經(jīng)驗是用當(dāng)前 $R$ 的跡取一個比例$$\sigma \delta \cdot \frac{\text{trace}(A P A^H)}{M}$$其中 $\delta$ 在 $10^{-3}$ 到 $10^{-2}$ 之間比較合適。這個取值思路是加載量跟信號總功率保持一個固定的相對水平這樣在信噪比變化時能自適應(yīng)地調(diào)整加載強度而不是用一個絕對常數(shù)。另一個思路是利用采樣協(xié)方差矩陣的底噪水平比如把 $\hat{R}_y$ 的最小特征值放大若干倍作為 $\sigma$。這種方法在信噪比較低時更精準(zhǔn)但需要額外做一次特征分解計算量稍大。對普通仿真和大多數(shù)實測場景前者已經(jīng)夠用。4.3 計算復(fù)雜度與工程加速IAA的復(fù)雜度大頭在后三行第2步要對 $M \times M$ 矩陣求逆第3步要對每個候選角度計算兩個二次型 $a_k^H R^{-1} y$ 和 $a_k^H R^{-1} a_k$。如果網(wǎng)格點數(shù) $K_g 181$陣元數(shù) $M 8$迭代15次整體計算量大概是 $15 \times (181 \times 8^2 8^3)$這在現(xiàn)代CPU上算毫秒級完全不是問題。但如果陣元數(shù)漲到64、網(wǎng)格點到361或者需要對幾百個快拍逐個處理時就要考慮加速手段。這里分享三個實測有效的做法。第一個做法是用Cholesky分解替代顯式求逆。對正定矩陣 $R$先做Cholesky分解 $R L L^H$然后把 $a_k^H R^{-1} y$ 拆成 $a_k^H (L^H)^{-1} L^{-1} y$用兩次前代/回代求解線性方程組避免顯式計算逆矩陣。數(shù)值穩(wěn)定性更好速度也更快。第二個做法是預(yù)計算那些與迭代無關(guān)的部分。所有候選角度的導(dǎo)向矢量可以事先存成 $M \times K_g$ 的矩陣每次迭代中反復(fù)用的 $a_k^H$ 和 $a_k$ 也都是現(xiàn)成的。真正需要每次計算的只有 $R^{-1}$以及它跟 $y$、$a_k$ 的乘積。第三個做法是并行化。第3步里每個角度 $k$ 的計算是相互獨立的天然適合用MATLAB的parfor或者Python的多進(jìn)程并行。當(dāng)網(wǎng)格點很多時并行效率可以逼近線性加速。5. 常見問題與實測經(jīng)驗5.1 問題排查速查表我把實際調(diào)試中遇到過的問題整理成了一張表遇到類似現(xiàn)象可以直接對照排查?，F(xiàn)象可能原因解決方案譜峰不明顯整個譜都很平對角加載量過大把 $\delta$ 調(diào)小或改用特征值底噪估計譜峰位置在幾輪迭代中漂移網(wǎng)格太粗加密角度網(wǎng)格或?qū)ψV峰附近做二次插值相干源只出一個峰網(wǎng)格失配或功率初始化偏差適當(dāng)提高迭代次數(shù)并檢查兩個源是否落在相鄰格點矩陣求逆報錯或出現(xiàn)NaN$R$ 接近奇異檢查是否忘了加對角加載項多快拍時功率譜有毛刺快拍間信號有起伏先對快拍做歸一化再進(jìn)入IAA這些情況里最隱蔽的是第二個。IAA本身對網(wǎng)格失配比MUSIC要不敏感一些但如果你把網(wǎng)格設(shè)得太粗比如3度一個格點兩個真實角度落在兩個格點之間譜峰就可能在相鄰格點之間跳來跳去。解決辦法一是加密網(wǎng)格二是對最終譜峰用拋物線擬合把小數(shù)級的角度差補回來。5.2 幾個容易忽略的經(jīng)驗提醒第一WLS的權(quán)重矩陣 $R^{-1}$ 不是“越白越好”。有些同學(xué)看到 $R^{-1}$ 就以為是對所有干擾做白化想讓噪聲完全變成白噪聲。但I(xiàn)AA的 $R$ 里包含信號本身的功率這個信號功率在對角線上會抬高 $R$ 的跡從而削弱 $R^{-1}$ 對信號方向的放大作用。換句話說IAA的權(quán)重矩陣實際上是一種“溫和”的白化它壓制干擾但并不徹底抵消信號。這個特性讓IAA在低信噪比下比Capon更穩(wěn)定。第二不要追求迭代到完美收斂。我見過有人把收斂閾值設(shè)成 $10^{-8}$結(jié)果跑了50輪還沒停下來。IAA的本質(zhì)是迭代加權(quán)最小二乘它沒有全局最優(yōu)解的那種“保證”但它的譜峰位置在早期迭代里就已經(jīng)基本穩(wěn)定。后面那些迭代更多是在微調(diào)譜峰的銳度和旁瓣電平。固定10到15次迭代既省時間又足夠準(zhǔn)。第三如果你把IAA的結(jié)果當(dāng)成下一步處理的基礎(chǔ)比如交給跟蹤濾波器時建議把譜峰旁邊的兩個格點功率也保留下來。IAA的譜峰不是純粹的脈沖它會帶有一定的展寬直接丟掉這些信息可能會造成角度估計偏差。我做過一次對比用譜峰加左右兩點做加權(quán)平均得到的角度比單純用譜峰位置要穩(wěn)定得多。第四關(guān)于深度學(xué)習(xí)與IAA結(jié)合的方向?，F(xiàn)在DOA估計領(lǐng)域里像SubspaceNet這類數(shù)據(jù)驅(qū)動方法越來越多但I(xiàn)AA這種“可解釋的迭代優(yōu)化”并不會過時。很多混合方案是讓網(wǎng)絡(luò)先給出一個粗略的目標(biāo)數(shù)目和角度先驗然后用IAA做精細(xì)估計。如果你有精力可以往這個方向試試我個人認(rèn)為這是工程落地價值很高的一個分支。最后說一點我自己的體會IAA這套推導(dǎo)表面上看是一堆矩陣公式本質(zhì)上其實就是“迭代加權(quán)最小二乘”這六個字。權(quán)重矩陣不是固定不變的而是隨著當(dāng)前對信號功率估計的更新不斷自我修正。想清楚這一點你就不會再被論文里那些符號繞暈。我在實際項目里用得最多的是把IAA當(dāng)做一個“穩(wěn)定器”——在MUSIC因為相干源失效、Capon因為協(xié)方差病態(tài)發(fā)抖的時候用IAA兜底出角度初值。雖然它比常規(guī)方法多跑好幾輪矩陣求逆但換來的是在復(fù)雜電磁環(huán)境下不用提心吊膽地調(diào)參數(shù)。測向這個領(lǐng)域穩(wěn)定壓倒一切。