MVDR(Capon)適応ビームフォーミングの理論と導出と実装

衛星地上局のアレイアンテナが、目的の通信衛星からの微弱な信号を受信しているとします。ところが同じ周波数帯に、別方向から強力な干渉波(隣接衛星や地上の妨害波)が飛び込んできました。アンテナを物理的に動かさずに、目的方向の感度はそのままに、干渉方向だけを「鈍感」にできないでしょうか。これを実現するのが適応ビームフォーミングであり、その代表格が本記事で扱う MVDR(Minimum Variance Distortionless Response、最小分散無歪み応答)ビームフォーマ、別名 Capon ビームフォーマ です。

固定ビームフォーミング(遅延和ビームフォーマ)は、あらかじめ決めた方向にビームを向けるだけで、どこに干渉があるかは一切考慮しません。一方 MVDR は、受信信号そのものから干渉と雑音の統計(相関行列)を学習し、「希望方向の利得は 1 に固定したまま、出力電力を最小化する」という方針で、干渉方向に自動的に深いヌル(感度の谷)を掘ります。この「データに適応してヌルを掘る」発想は、レーダーの妨害波除去(サイドローブキャンセラ)、5G/Beyond-5G の Massive MIMO 基地局、音響アレイによる話者分離、電波天文学の干渉除去など、応用先は枚挙にいとまがありません。

本記事では、まずアレイ受信の数理モデルとステアリングベクトルを整理し、MVDR の最適化問題

$$ \min_{\bm{w}} \ \bm{w}^H \bm{R} \bm{w} \quad \text{s.t.} \quad \bm{w}^H \bm{a}(\theta_0) = 1 $$

をラグランジュ未定乗数法で一行ずつ解いて、最適ウェイト

$$ \bm{w}_{\text{MVDR}} = \frac{\bm{R}^{-1} \bm{a}(\theta_0)}{\bm{a}(\theta_0)^H \bm{R}^{-1} \bm{a}(\theta_0)} $$

を導きます。そのうえで Capon スペクトルがなぜ固定ビームより鋭く干渉をヌルできるのかを説明し、Python で干渉波を含むアレイ受信を合成して放射パターンのヌルと出力 SINR の改善量を数値で確かめます。

本記事の内容

  • アレイ信号モデルとステアリングベクトルの直感的理解
  • 空間相関行列 $\bm{R}$ の意味と推定
  • MVDR の制約付き最小化からラグランジュ法による最適ウェイトの導出(省略なし)
  • Capon スペクトルが深いヌルを形成する理由
  • 出力 SINR の定義と、固定ビームフォーマとの比較
  • Python による放射パターン・スペクトル・SINR の可視化

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

アレイ受信のモデルとステアリングベクトル

まずは1本のアンテナから複数本へ

1本のアンテナは、空間のある一点で電界の時間変化を測るだけの「点」のセンサです。これでは波がどの方向から来たのかを区別できません。ところが、アンテナを一直線に等間隔で並べると、同じ平面波でも各素子に到達時刻のわずかなズレが生じます。海岸に斜めから押し寄せる波が、横一列に並んだ杭に少しずつ時間差で当たるのと同じです。この時間差(位相差)のパターンを読めば、波の到来方向(DOA: Direction of Arrival)を推定でき、さらに位相差を打ち消すように足し合わせれば特定方向の波だけを強められます。これがビームフォーミングの出発点です。

ここでは最も基本的な 等間隔直線アレイ(ULA: Uniform Linear Array) を考えます。素子数を $M$、素子間隔を $d$、波長を $\lambda$ とします。遠方から到来する平面波の方向を、アレイ正面(ブロードサイド)からの角度 $\theta$ で表します。

素子間の位相差

到来角 $\theta$ の平面波は、隣り合う素子に届くまでに余分な距離 $d \sin\theta$ を進みます。波長 $\lambda$ あたり $2\pi$ の位相が回ることから、隣接素子間の位相差は

$$ \phi = \frac{2\pi}{\lambda} d \sin\theta $$

です。$m$ 番目の素子($m = 0, 1, \dots, M-1$)は基準素子に対して $m\phi$ だけ位相が進みます。したがって、振幅 $1$・基準位相 $0$ の単一波源に対する各素子の受信複素振幅をまとめたベクトルがステアリングベクトルです。

$$ \bm{a}(\theta) = \begin{bmatrix} 1 \\ e^{j \frac{2\pi}{\lambda} d \sin\theta} \\ e^{j \frac{2\pi}{\lambda} d \, 2\sin\theta} \\ \vdots \\ e^{j \frac{2\pi}{\lambda} d (M-1)\sin\theta} \end{bmatrix} \in \mathbb{C}^{M} $$

ステアリングベクトルは「到来角 $\theta$ という方向情報を、$M$ 個の複素位相のパターンに翻訳した指紋」だと捉えてください。方向が違えば位相パターンも違い、$\bm{a}(\theta)$ はその方向の波が素子上にどう刻まれるかを完全に記述します。

受信信号ベクトル

希望波(到来角 $\theta_0$、複素振幅 $s(t)$)に加えて、$L$ 個の干渉波(到来角 $\theta_i$、振幅 $i_\ell(t)$)と、各素子に独立に乗る熱雑音 $\bm{n}(t)$ があるとします。スナップショット(ある時刻 $t$ の $M$ 素子分の観測)をまとめると

$$ \bm{x}(t) = \bm{a}(\theta_0)\, s(t) + \sum_{\ell=1}^{L} \bm{a}(\theta_\ell)\, i_\ell(t) + \bm{n}(t) $$

と書けます。第1項が欲しい信号、第2項が消したい干渉、第3項が避けられない雑音です。ビームフォーミングとは、この $M$ 次元ベクトル $\bm{x}(t)$ に複素ウェイト $\bm{w} \in \mathbb{C}^M$ を掛けて1つのスカラー出力

$$ y(t) = \bm{w}^H \bm{x}(t) $$

を作る操作です($\bm{w}^H$ は共役転置)。問題は「どんな $\bm{w}$ を選べば、希望波を保ったまま干渉と雑音を最小化できるか」です。これに答えるために、まず信号の統計をまとめる相関行列を導入します。

空間相関行列とは

出力電力をウェイトで表す

ビームフォーマ出力 $y(t) = \bm{w}^H \bm{x}(t)$ の平均電力は、期待値 $\mathbb{E}[\cdot]$ を使って

$$ P(\bm{w}) = \mathbb{E}\big[\, |y(t)|^2 \,\big] = \mathbb{E}\big[\, \bm{w}^H \bm{x}(t)\, \bm{x}(t)^H \bm{w} \,\big] $$

と書けます。$\bm{w}$ は確率変数ではなく定数なので期待値の外に出せて

$$ P(\bm{w}) = \bm{w}^H \, \mathbb{E}\big[\, \bm{x}(t)\, \bm{x}(t)^H \,\big]\, \bm{w} = \bm{w}^H \bm{R}\, \bm{w} $$

となります。ここで現れた $M \times M$ 行列

$$ \bm{R} = \mathbb{E}\big[\, \bm{x}(t)\, \bm{x}(t)^H \,\big] $$

空間相関行列(共分散行列)です。対角成分は各素子の受信電力、非対角成分は素子間の相関を表し、「どの方向からどれだけの電力が来ているか」という空間的な電力分布の情報がすべて詰まっています。

Rの構造と性質

各波源と雑音が互いに無相関であれば、相関行列は各成分の寄与の和に分解できます。希望波電力を $\sigma_s^2$、干渉波電力を $\sigma_\ell^2$、雑音電力を $\sigma_n^2$ とすると

$$ \bm{R} = \sigma_s^2\, \bm{a}(\theta_0)\bm{a}(\theta_0)^H + \sum_{\ell=1}^{L} \sigma_\ell^2\, \bm{a}(\theta_\ell)\bm{a}(\theta_\ell)^H + \sigma_n^2 \bm{I} $$

です。重要な性質を確認しておきます。$\bm{R}$ はエルミート($\bm{R}^H = \bm{R}$)で、雑音項 $\sigma_n^2 \bm{I}$ があるため正定値であり、したがって逆行列 $\bm{R}^{-1}$ が必ず存在します。後の導出で $\bm{R}^{-1}$ を堂々と使えるのは、この雑音項のおかげです。

実運用では真の期待値は得られないので、$N$ 個のスナップショットから標本相関行列

$$ \hat{\bm{R}} = \frac{1}{N} \sum_{t=1}^{N} \bm{x}(t)\, \bm{x}(t)^H $$

で推定します。$N$ が十分大きければ $\hat{\bm{R}} \to \bm{R}$ に収束します。これで MVDR の材料がそろいました。次は「希望波を歪ませず、出力電力を最小化する」という方針を最適化問題として定式化します。

MVDRの最適化問題

なぜ「電力最小化」で干渉が消えるのか

直感が肝心です。もし無条件に出力電力 $\bm{w}^H \bm{R} \bm{w}$ を最小化したら、答えは自明に $\bm{w} = \bm{0}$(何も受信しない)になってしまい無意味です。そこで「希望方向 $\theta_0$ から来た波は、振幅も位相も一切いじらずにそのまま通す」という無歪み制約を課します。

$$ \bm{w}^H \bm{a}(\theta_0) = 1 $$

この制約のもとで出力電力を最小化すると、何が起きるでしょうか。出力電力は「希望波+干渉+雑音」のすべての寄与の和です。希望波の寄与は制約で固定されている(必ず通る)ので、最小化の余地があるのは干渉と雑音の寄与だけです。つまりビームフォーマは、希望波を保ったまま、自由に使える残りの自由度を総動員して干渉方向の利得を下げる(ヌルを掘る)ことで電力を下げにいきます。これが MVDR の核心的なアイデアです。

最適化問題の定式化

以上をまとめると、MVDR ウェイトは次の制約付き最小化問題の解として定義されます。

$$ \bm{w}_{\text{MVDR}} = \arg\min_{\bm{w}} \ \bm{w}^H \bm{R}\, \bm{w} \quad \text{subject to} \quad \bm{w}^H \bm{a}(\theta_0) = 1 $$

目的関数 $\bm{w}^H \bm{R} \bm{w}$ は出力電力(分散)、制約 $\bm{w}^H \bm{a}(\theta_0) = 1$ は無歪み応答です。名前 “Minimum Variance Distortionless Response” がそのまま問題の中身を表しているわけです。

この問題は、$\bm{R}$ が正定値なので目的関数は $\bm{w}$ に関して凸であり、制約は線形ですから、解は一意に定まります。次節でラグランジュ未定乗数法を使い、この $\bm{w}_{\text{MVDR}}$ の閉形式を一行ずつ導出します。

ラグランジュ未定乗数法による導出

ラグランジアンの構成

制約付き最小化を解く定石は、制約を未定乗数で目的関数に組み込んだラグランジアンを作り、その停留点を求めることです。複素ベクトルに対する制約 $\bm{w}^H \bm{a}(\theta_0) = 1$ を扱うため、複素乗数 $\lambda$ を導入し、制約とその複素共役の両方を加えた実数値ラグランジアンを構成します。

$$ \mathcal{L}(\bm{w}, \lambda) = \bm{w}^H \bm{R}\, \bm{w} – \lambda^{*}\big(\bm{w}^H \bm{a}(\theta_0) – 1\big) – \lambda\big(\bm{a}(\theta_0)^H \bm{w} – 1\big) $$

第1項が最小化したい出力電力、第2項と第3項が制約項です。$\bm{w}^H \bm{R} \bm{w}$ は実数($\bm{R}$ がエルミートなため)であり、制約項も互いに共役な対を加えているので $\mathcal{L}$ 全体が実数になります。これにより、後で出てくる複素勾配を矛盾なく扱えます。

複素勾配(Wirtinger微分)でゼロとおく

複素ベクトル $\bm{w}$ に関する最小化では、$\bm{w}$ と $\bm{w}^H$ を独立変数とみなして $\bm{w}^H$ について偏微分し、それをゼロとおく(Wirtinger微分)のが標準的な手法です。各項を $\bm{w}^H$ で微分します。

まず第1項 $\bm{w}^H \bm{R} \bm{w}$ を $\bm{w}^H$ で微分すると $\bm{R}\bm{w}$ になります。第2項 $-\lambda^*(\bm{w}^H \bm{a}(\theta_0) – 1)$ を $\bm{w}^H$ で微分すると $-\lambda^* \bm{a}(\theta_0)$ です。第3項 $-\lambda(\bm{a}(\theta_0)^H \bm{w} – 1)$ には $\bm{w}^H$ が含まれないので微分は $\bm{0}$ になります。以上をまとめて勾配をゼロとおくと

$$ \frac{\partial \mathcal{L}}{\partial \bm{w}^H} = \bm{R}\, \bm{w} – \lambda^{*} \bm{a}(\theta_0) = \bm{0} $$

を得ます。ここから停留条件として

$$ \bm{R}\, \bm{w} = \lambda^{*} \bm{a}(\theta_0) $$

が導かれます。$\bm{R}$ は正定値で逆行列を持つので、両辺に左から $\bm{R}^{-1}$ を掛けて $\bm{w}$ について解くと

$$ \bm{w} = \lambda^{*} \bm{R}^{-1} \bm{a}(\theta_0) $$

となります。これでウェイトの「形」は $\bm{R}^{-1}\bm{a}(\theta_0)$ に比例することが判明しました。あとは比例係数 $\lambda^*$ を制約から決めるだけです。

制約から未定乗数を決定する

得られた $\bm{w} = \lambda^* \bm{R}^{-1} \bm{a}(\theta_0)$ を、無歪み制約 $\bm{w}^H \bm{a}(\theta_0) = 1$ に代入します。$\bm{w}^H = (\lambda^* \bm{R}^{-1} \bm{a}(\theta_0))^H = \lambda\, \bm{a}(\theta_0)^H \bm{R}^{-H}$ ですが、$\bm{R}$ はエルミートなので $\bm{R}^{-H} = \bm{R}^{-1}$ が使えて

$$ \bm{w}^H \bm{a}(\theta_0) = \lambda\, \bm{a}(\theta_0)^H \bm{R}^{-1} \bm{a}(\theta_0) = 1 $$

となります。ここで分母に現れる $\bm{a}(\theta_0)^H \bm{R}^{-1} \bm{a}(\theta_0)$ はスカラーで、$\bm{R}^{-1}$ が正定値なので必ず正の実数です。したがって $\lambda$ について解くと

$$ \lambda = \frac{1}{\bm{a}(\theta_0)^H \bm{R}^{-1} \bm{a}(\theta_0)} $$

を得ます。$\lambda$ が実数なので $\lambda^* = \lambda$ です。これを $\bm{w} = \lambda^* \bm{R}^{-1}\bm{a}(\theta_0)$ に戻します。

最適ウェイト

以上から、MVDR ビームフォーマの最適ウェイトが閉じた形で求まりました。

$$ \boxed{\ \bm{w}_{\text{MVDR}} = \frac{\bm{R}^{-1} \bm{a}(\theta_0)}{\bm{a}(\theta_0)^H \bm{R}^{-1} \bm{a}(\theta_0)}\ } $$

分子 $\bm{R}^{-1}\bm{a}(\theta_0)$ がウェイトの方向を決め、分母のスカラーが無歪み制約 $\bm{w}^H \bm{a}(\theta_0) = 1$ を満たすように全体を正規化しています。

この式の意味を噛みしめましょう。もし干渉も指向性のある雑音もなく $\bm{R} = \sigma_n^2 \bm{I}$ なら、$\bm{R}^{-1} = \sigma_n^{-2}\bm{I}$ となって $\bm{w}_{\text{MVDR}} \propto \bm{a}(\theta_0)$、すなわち固定ビームフォーマ(遅延和ビームフォーマ)に一致します。逆に $\bm{R}$ に強い干渉成分が含まれると、$\bm{R}^{-1}$ がその方向の成分を抑え込むように働き、ウェイトが干渉方向の利得を下げる方向に「歪み」ます。MVDR が固定ビームと決定的に違うのは、この $\bm{R}^{-1}$ を通じてデータ(環境)に適応する点です。

最小出力電力

ついでに、得られた最小出力電力も計算しておきます。$\bm{w}_{\text{MVDR}}$ を $P = \bm{w}^H \bm{R} \bm{w}$ に代入し、$\bm{R}^{-1}\bm{R}\bm{R}^{-1} = \bm{R}^{-1}$ を使うと

$$ P_{\min} = \frac{\bm{a}(\theta_0)^H \bm{R}^{-1} \bm{R}\, \bm{R}^{-1} \bm{a}(\theta_0)}{\big(\bm{a}(\theta_0)^H \bm{R}^{-1} \bm{a}(\theta_0)\big)^2} = \frac{\bm{a}(\theta_0)^H \bm{R}^{-1} \bm{a}(\theta_0)}{\big(\bm{a}(\theta_0)^H \bm{R}^{-1} \bm{a}(\theta_0)\big)^2} = \frac{1}{\bm{a}(\theta_0)^H \bm{R}^{-1} \bm{a}(\theta_0)} $$

となります。この $P_{\min}$ を、$\theta_0$ を走査角 $\theta$ に置き換えて方向の関数とみなしたものが、次に説明する Capon スペクトルです。導出が一段落したので、このスペクトルがなぜ鋭いヌルを作るのかを見ていきます。

Caponスペクトルと深いヌルの理由

Caponスペクトルの定義

MVDR の最小出力電力を、注目方向 $\theta$ をパラメータとして掃引したものが Capon スペクトル(最小分散スペクトル)です。

$$ P_{\text{Capon}}(\theta) = \frac{1}{\bm{a}(\theta)^H \bm{R}^{-1} \bm{a}(\theta)} $$

これは「方向 $\theta$ に無歪みビームを向けたとき、出力に残る最小電力」を表します。$\theta$ に本当に強い波源があれば、その方向に向けたビームには大きな電力が残るのでスペクトルはピークになります。波源のない方向ではビームが自由に他の信号を抑え込めるので残留電力は小さくなります。こうして $P_{\text{Capon}}(\theta)$ は到来方向のスペクトル推定にもなります。

なぜ固定ビームより深いヌルが掘れるのか

固定ビームフォーマ(遅延和、Bartlett ビームフォーマ)のスペクトルは

$$ P_{\text{Bartlett}}(\theta) = \frac{\bm{a}(\theta)^H \bm{R}\, \bm{a}(\theta)}{\bm{a}(\theta)^H \bm{a}(\theta)} $$

で、ウェイトは常に $\bm{w} \propto \bm{a}(\theta)$ に固定されます。このウェイトはアレイ開口で決まる固定の指向性しか持たず、ヌルの位置も深さも開口長だけで決まってしまいます。$M$ 素子のアレイの放射パターンは $M-1$ 個のヌルを持ちますが、それらは等間隔に並んだ「決め打ち」のヌルで、干渉が来ている方向に都合よく一致してくれる保証はありません。

一方 MVDR は $\bm{R}^{-1}$ を通じて干渉の位置を「知って」います。$\bm{R}$ の中で干渉方向の固有値は雑音より遥かに大きく、$\bm{R}^{-1}$ ではその方向の成分が $1/\sigma_\ell^2$ のように極端に小さくなります。最適ウェイト $\bm{w}_{\text{MVDR}} \propto \bm{R}^{-1}\bm{a}(\theta_0)$ はこの抑圧された成分を引き継ぐため、干渉方向に正確にヌルを配置します。干渉が強い($\sigma_\ell^2$ が大きい)ほど $\bm{R}^{-1}$ での抑圧が強まり、ヌルは深くなります。これが「データ適応型ビームフォーマは干渉方向に深いヌルを掘れる」ことの正体です。

直感的にまとめると、固定ビームは「地図を見ずに決まった形のビームを向ける」のに対し、MVDR は「受信データから干渉の地図($\bm{R}$)を作り、その地図を見て干渉だけを狙い撃ちで消す」のです。応用先としては、レーダーのサイドローブ妨害除去や、衛星通信における同一チャネル干渉の抑圧が典型的です。では実際に数値で確かめましょう。

Pythonでの実装

ステアリングベクトルとシナリオ設定

まず、ULA のステアリングベクトルを生成する関数と、希望波・干渉波・雑音を含む受信スナップショットを合成する部分を実装します。素子間隔は半波長($d = \lambda/2$)とします。半波長間隔はグレーティングローブ(不要な大ピーク)を避けつつ最大の角度分解能を得る標準的な選択です。

import numpy as np
import matplotlib.pyplot as plt

# 乱数シードを固定して再現性を確保
rng = np.random.default_rng(0)

# ===== アレイとシナリオの設定 =====
M = 16            # 素子数
d_over_lambda = 0.5   # 素子間隔 / 波長(半波長間隔)
N = 2000          # スナップショット数(相関行列推定用)

theta0 = 0.0      # 希望波の到来角 [deg](ブロードサイド)
theta_intf = [-30.0, 25.0]   # 干渉波の到来角 [deg]
SNR_dB = 0.0      # 希望波のSNR [dB](素子あたり)
INR_dB = [30.0, 25.0]        # 各干渉波のINR(干渉対雑音比)[dB]

def steering_vector(theta_deg, M, d_over_lambda):
    """到来角theta_deg [deg]に対するULAのステアリングベクトルを返す"""
    theta = np.deg2rad(theta_deg)
    m = np.arange(M)
    # 各素子の位相 exp(j 2pi d/lambda m sin(theta))
    return np.exp(1j * 2 * np.pi * d_over_lambda * m * np.sin(theta))

ここでは 16 素子の ULA に対し、ブロードサイド($0^\circ$)に希望波、$-30^\circ$ と $+25^\circ$ に強力な干渉波(INR 30 dB と 25 dB)を置きました。希望波は SNR 0 dB と、干渉に比べてはるかに弱い設定です。この「弱い希望波が強い干渉に埋もれている」状況こそ、適応ビームフォーミングの真価が問われる場面です。

受信信号の合成と相関行列の推定

次に、設定した電力に従って各波源を生成し、$N$ スナップショット分の受信ベクトルを合成して標本相関行列 $\hat{\bm{R}}$ を作ります。

# 各信号の振幅(電力の平方根)。雑音電力は1に正規化
sigma_n = 1.0
amp_s = sigma_n * 10**(SNR_dB / 20)
amp_i = [sigma_n * 10**(inr / 20) for inr in INR_dB]

a0 = steering_vector(theta0, M, d_over_lambda)

# 受信スナップショット行列 X (M x N) を合成
X = np.zeros((M, N), dtype=complex)
# 希望波(複素ガウス信号)
s = amp_s * (rng.standard_normal(N) + 1j * rng.standard_normal(N)) / np.sqrt(2)
X += np.outer(a0, s)
# 干渉波
for th, amp in zip(theta_intf, amp_i):
    a = steering_vector(th, M, d_over_lambda)
    ig = amp * (rng.standard_normal(N) + 1j * rng.standard_normal(N)) / np.sqrt(2)
    X += np.outer(a, ig)
# 雑音
noise = sigma_n * (rng.standard_normal((M, N)) + 1j * rng.standard_normal((M, N))) / np.sqrt(2)
X += noise

# 標本相関行列の推定
R = (X @ X.conj().T) / N
print("相関行列のサイズ:", R.shape)
print("最大固有値 / 最小固有値:", 
      np.max(np.linalg.eigvalsh(R)).real, "/", np.min(np.linalg.eigvalsh(R)).real)

出力された最大固有値と最小固有値の比に注目してください。最小固有値は雑音電力(約 1)程度であるのに対し、最大固有値は数百〜千のオーダーになります。これは強い干渉が相関行列に大きな固有値として現れていることを示します。MVDR の $\bm{R}^{-1}$ は、まさにこの大きな固有値方向を強く抑え込むことでヌルを作ります。

MVDR・固定ビームフォーマのウェイト計算

導出した式に従って MVDR ウェイトを計算し、比較対象として固定ビーム(遅延和)ウェイトも用意します。

# Rの逆行列(数値安定のため微小な対角ローディングを加える)
R_inv = np.linalg.inv(R + 1e-6 * np.eye(M))

# MVDRウェイト: w = R^{-1} a0 / (a0^H R^{-1} a0)
num = R_inv @ a0
w_mvdr = num / (a0.conj().T @ num)

# 固定ビーム(遅延和)ウェイト: w = a0 / M (無歪み制約 w^H a0 = 1 を満たす正規化)
w_cbf = a0 / M

# 無歪み制約の確認(どちらも 1 になるはず)
print("MVDR  w^H a0 =", np.round(w_mvdr.conj().T @ a0, 4))
print("固定  w^H a0 =", np.round(w_cbf.conj().T @ a0, 4))

両方のウェイトについて $\bm{w}^H \bm{a}(\theta_0)$ がほぼ $1$ になっていることが確認できます。これは両者ともに希望方向 $\theta_0$ では同じ無歪み応答を持つことを意味します。違いは「希望方向以外」での振る舞いであり、それを次に放射パターンで可視化します。

放射パターンの可視化

ウェイトを固定したまま、走査角 $\theta$ を $-90^\circ$ から $90^\circ$ まで変えて応答 $|\bm{w}^H \bm{a}(\theta)|^2$ を計算したものが放射パターン(ビームパターン)です。

# 走査角の範囲
angles = np.linspace(-90, 90, 1801)
A = np.array([steering_vector(th, M, d_over_lambda) for th in angles]).T  # M x len(angles)

# 各角度での応答(dB)
resp_mvdr = 20 * np.log10(np.abs(w_mvdr.conj().T @ A) + 1e-12)
resp_cbf  = 20 * np.log10(np.abs(w_cbf.conj().T @ A) + 1e-12)

plt.figure(figsize=(11, 6))
plt.plot(angles, resp_cbf, 'b-', lw=1.5, label='Conventional (delay-and-sum)')
plt.plot(angles, resp_mvdr, 'r-', lw=1.8, label='MVDR (Capon)')
plt.axvline(theta0, color='green', ls='--', alpha=0.7, label='Desired ($\\theta_0$)')
for th in theta_intf:
    plt.axvline(th, color='orange', ls=':', alpha=0.8)
plt.xlabel('Angle [deg]')
plt.ylabel('Array response [dB]')
plt.title('Beam pattern: MVDR vs Conventional beamformer')
plt.ylim([-80, 5])
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.savefig('mvdr_beampattern.png', dpi=150, bbox_inches='tight')
plt.show()

このグラフから、MVDR の最大の特長がはっきり読み取れます。希望方向(緑の破線、$0^\circ$)では両者とも応答が $0$ dB(利得 1)で一致しています。ところが干渉方向(オレンジの点線、$-30^\circ$ と $+25^\circ$)を見ると、固定ビーム(青)は通常のサイドローブ程度($-15$ dB 前後)しか抑圧できていないのに対し、MVDR(赤)は $-50$ dB を下回る鋭く深いヌルを正確にその方向に掘っています。これは前節で説明した「$\bm{R}^{-1}$ が干渉方向の強い成分を抑え込む」効果が現れた結果です。固定ビームのヌルは開口で決まる固定位置にあるため干渉方向には一致しませんが、MVDR は干渉の位置に適応してヌルを置いています。

Caponスペクトルによる到来方向推定

次に Capon スペクトル $P_{\text{Capon}}(\theta) = 1/(\bm{a}(\theta)^H \bm{R}^{-1}\bm{a}(\theta))$ を計算し、固定ビーム(Bartlett)スペクトルと比較して、空間スペクトルとしての分解能の違いを見ます。

# Caponスペクトルと Bartlett スペクトル
P_capon = np.zeros(len(angles))
P_bart  = np.zeros(len(angles))
for k in range(len(angles)):
    a = A[:, k]
    P_capon[k] = 1.0 / np.real(a.conj().T @ R_inv @ a)
    P_bart[k]  = np.real(a.conj().T @ R @ a) / M

# dB正規化
P_capon_dB = 10 * np.log10(P_capon / np.max(P_capon))
P_bart_dB  = 10 * np.log10(P_bart  / np.max(P_bart))

plt.figure(figsize=(11, 6))
plt.plot(angles, P_bart_dB, 'b-', lw=1.5, label='Bartlett (conventional)')
plt.plot(angles, P_capon_dB, 'r-', lw=1.8, label='Capon (MVDR)')
plt.axvline(theta0, color='green', ls='--', alpha=0.7, label='Desired')
for th in theta_intf:
    plt.axvline(th, color='orange', ls=':', alpha=0.8)
plt.xlabel('Angle [deg]')
plt.ylabel('Normalized spatial spectrum [dB]')
plt.title('Capon vs Bartlett spatial spectrum')
plt.ylim([-40, 3])
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.savefig('capon_spectrum.png', dpi=150, bbox_inches='tight')
plt.show()

このスペクトルから、Capon スペクトル(赤)が希望波と2つの干渉波の合計3方向に鋭いピークを立てているのが分かります。一方 Bartlett スペクトル(青)はピークが幅広く、特に近接した波源では分離できずなだらかな山になりがちです。Capon が各波源を細いピークで分離できるのは、$\bm{R}^{-1}$ により注目方向以外の電力を抑圧した「最小分散」の枠組みが、レイリー限界に縛られた固定ビームより高い角度分解能を与えるためです。ピークの高さが電力をそのまま反映していない点には注意が必要ですが、波源の「位置」の推定には Capon が圧倒的に優れています。

出力SINRの比較

最後に、適応ビームフォーミングの目的である出力 SINR(信号対干渉雑音比)の改善量を数値化します。出力 SINR は、希望波の出力電力を干渉+雑音の出力電力で割ったものです。希望波の相関行列を $\bm{R}_s = \sigma_s^2 \bm{a}(\theta_0)\bm{a}(\theta_0)^H$、干渉+雑音の相関行列を $\bm{R}_{i+n}$ とすると

$$ \text{SINR}(\bm{w}) = \frac{\bm{w}^H \bm{R}_s\, \bm{w}}{\bm{w}^H \bm{R}_{i+n}\, \bm{w}} = \frac{\sigma_s^2\, |\bm{w}^H \bm{a}(\theta_0)|^2}{\bm{w}^H \bm{R}_{i+n}\, \bm{w}} $$

で計算できます。理論モデルから $\bm{R}_{i+n}$ を作り、両ビームフォーマの出力 SINR を比較します。

# 干渉+雑音の相関行列(理論モデル)
R_in = sigma_n**2 * np.eye(M, dtype=complex)
for th, amp in zip(theta_intf, amp_i):
    a = steering_vector(th, M, d_over_lambda)
    R_in += amp**2 * np.outer(a, a.conj())

def output_sinr(w, a0, sigma_s2, R_in):
    """出力SINRを計算"""
    sig = sigma_s2 * np.abs(w.conj().T @ a0)**2
    intf_noise = np.real(w.conj().T @ R_in @ w)
    return sig / intf_noise

sigma_s2 = amp_s**2
sinr_mvdr = output_sinr(w_mvdr, a0, sigma_s2, R_in)
sinr_cbf  = output_sinr(w_cbf,  a0, sigma_s2, R_in)

# 最適SINR(理論上限): SINR_opt = sigma_s2 * a0^H R_in^{-1} a0
R_in_inv = np.linalg.inv(R_in)
sinr_opt = sigma_s2 * np.real(a0.conj().T @ R_in_inv @ a0)

print(f"固定ビーム 出力SINR : {10*np.log10(sinr_cbf):6.2f} dB")
print(f"MVDR      出力SINR : {10*np.log10(sinr_mvdr):6.2f} dB")
print(f"最適       出力SINR : {10*np.log10(sinr_opt):6.2f} dB")
print(f"MVDRによる改善量   : {10*np.log10(sinr_mvdr/sinr_cbf):6.2f} dB")

この出力から、MVDR が固定ビームに対して数十 dB のオーダーで出力 SINR を改善することが読み取れます。固定ビームでは干渉波が素子あたり 30 dB と希望波より遥かに強いため、希望方向にビームを向けてもサイドローブ経由で大量の干渉が漏れ込み、出力 SINR は大きく劣化します。MVDR は干渉方向に深いヌルを掘るため、干渉の漏れ込みが激減し、SINR が干渉のない場合に近い水準まで回復します。さらに MVDR の出力 SINR が理論上限 $\text{SINR}_{\text{opt}} = \sigma_s^2\, \bm{a}(\theta_0)^H \bm{R}_{i+n}^{-1}\bm{a}(\theta_0)$ にほぼ一致することも確認でき、MVDR が干渉雑音環境で SINR を最大化する最適ビームフォーマであることが数値的に裏付けられます。

干渉強度に対するSINRの振る舞い

最後に、干渉の強さ(INR)を変えながら両ビームフォーマの出力 SINR がどう変化するかを掃引し、MVDR の頑健性を確認します。

inr_sweep = np.arange(0, 41, 2)   # 干渉INRを 0〜40 dB で掃引
sinr_mvdr_list, sinr_cbf_list = [], []

for inr in inr_sweep:
    amp = sigma_n * 10**(inr / 20)
    R_in_s = sigma_n**2 * np.eye(M, dtype=complex)
    for th in theta_intf:
        a = steering_vector(th, M, d_over_lambda)
        R_in_s += amp**2 * np.outer(a, a.conj())
    R_tot = sigma_s2 * np.outer(a0, a0.conj()) + R_in_s
    Rinv = np.linalg.inv(R_tot + 1e-6*np.eye(M))
    num = Rinv @ a0
    w_m = num / (a0.conj().T @ num)
    w_c = a0 / M
    sinr_mvdr_list.append(10*np.log10(output_sinr(w_m, a0, sigma_s2, R_in_s)))
    sinr_cbf_list.append(10*np.log10(output_sinr(w_c, a0, sigma_s2, R_in_s)))

plt.figure(figsize=(10, 6))
plt.plot(inr_sweep, sinr_cbf_list, 'b-o', ms=4, label='Conventional')
plt.plot(inr_sweep, sinr_mvdr_list, 'r-s', ms=4, label='MVDR')
plt.xlabel('Interference-to-Noise Ratio (INR) [dB]')
plt.ylabel('Output SINR [dB]')
plt.title('Output SINR vs interference strength')
plt.grid(True, alpha=0.3)
plt.legend()
plt.tight_layout()
plt.savefig('mvdr_sinr_vs_inr.png', dpi=150, bbox_inches='tight')
plt.show()

このグラフから、決定的な違いが見えます。固定ビーム(青)は干渉が強くなるにつれて出力 SINR が一直線に劣化していきます。これは固定のヌルが干渉方向に一致していないため、干渉が強くなるほど漏れ込みもそのまま増えるからです。対照的に MVDR(赤)は INR が増えても出力 SINR がほぼ一定の高い水準を保ちます。干渉が強くなるほど $\bm{R}^{-1}$ がその方向を強く抑え込み、ヌルがさらに深くなるため、結果的に漏れ込みが抑えられ続けるのです。「干渉が強いほど MVDR の優位性が際立つ」というこの振る舞いは、適応ビームフォーミングが妨害波の激しい実環境で重宝される理由そのものです。

まとめ

本記事では、MVDR(Capon)適応ビームフォーミングについて、最適化問題の定式化からラグランジュ未定乗数法による導出、そして Python による検証までを一通り解説しました。

  • 問題設定: アレイ受信を $\bm{x}(t) = \bm{a}(\theta_0)s(t) + \sum_\ell \bm{a}(\theta_\ell)i_\ell(t) + \bm{n}(t)$ とモデル化し、出力電力は $\bm{w}^H \bm{R}\bm{w}$、空間相関行列は $\bm{R} = \mathbb{E}[\bm{x}\bm{x}^H]$ で表される
  • 最適化問題: 希望方向の利得を固定する無歪み制約 $\bm{w}^H\bm{a}(\theta_0)=1$ のもとで出力電力 $\bm{w}^H\bm{R}\bm{w}$ を最小化する。これにより干渉と雑音の寄与だけが最小化される
  • 最適ウェイト: ラグランジュ法により $\bm{w}_{\text{MVDR}} = \bm{R}^{-1}\bm{a}(\theta_0) / (\bm{a}(\theta_0)^H\bm{R}^{-1}\bm{a}(\theta_0))$。$\bm{R} = \sigma_n^2\bm{I}$ なら固定ビームに一致し、干渉があると $\bm{R}^{-1}$ が干渉方向の利得を下げる
  • Caponスペクトル: $P_{\text{Capon}}(\theta) = 1/(\bm{a}(\theta)^H\bm{R}^{-1}\bm{a}(\theta))$。固定ビームより高い分解能で波源を分離し、干渉方向に深いヌルを形成する
  • 性能: Python 実験で MVDR が干渉方向に $-50$ dB を超える深いヌルを掘り、固定ビームに対して出力 SINR を大幅に改善し、理論上限にほぼ一致することを確認した

MVDR は、相関行列 $\bm{R}$ にバッチで一度逆行列をかける手法でした。しかし実環境では信号環境が時々刻々と変化するため、$\bm{R}$ を逐次的に更新しながらウェイトを追従させる適応アルゴリズム(LMS、RLS、サンプル行列反転 SMI など)が必要になります。また、$\bm{R}$ の推定誤差や希望方向のステアリングベクトル誤差に対する頑健化(対角ローディングや LCMV への拡張)も実用上は重要なテーマです。

次のステップとして、以下の記事も参考にしてください。