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