レーダーは目標からの反射信号を検出する装置ですが、現実の環境では目標以外にも地面、海面、建築物、雨や雲などから大きな反射信号が返ってきます。これらの不要な反射信号をクラッタ（Clutter）と呼びます。クラッタは一般に目標信号よりもはるかに強く、そのままでは目標の検出が困難になります。
MTI（Moving Target Indication: 移動目標指示）は、パルス間の位相変化を利用して移動目標を検出する技術です。静止クラッタはパルス間で位相が変化しないのに対し、移動目標はドップラー効果により位相がシフトします。MTIフィルタはこの差異を利用して、クラッタを除去しつつ移動目標の信号を通過させます。
本記事では、クラッタの特性からスタートし、パルスキャンセラの伝達関数と周波数応答を丁寧に導出し、ブラインド速度の問題とその対策、MTI改善係数の定義と計算を解説し、最後にPythonで移動目標＋クラッタのシミュレーションを実装します。
本記事の内容
- クラッタの特性と問題点
- MTIの基本原理（パルス間位相変化による移動目標検出）
- 単パルスキャンセラ $H(z) = 1 – z^{-1}$ の伝達関数と周波数応答
- ブラインド速度 $v_{\mathrm{blind}} = \lambda/(2T_r)$ の問題
- 2パルス/3パルスキャンセラの周波数応答比較
- スタガードPRFによるブラインド速度改善
- MTI改善係数の定義と計算
- Pythonでのシミュレーションと可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
クラッタとは
クラッタの分類
クラッタは発生源によって以下のように分類されます。
地表面クラッタ: 地面、森林、都市構造物などからの反射。一般に最も強いクラッタ源であり、レーダーに対してほぼ静止しているため、ドップラーシフトはほぼゼロです。
海面クラッタ（シークラッタ）: 海面の波からの反射。波浪の運動によりドップラースペクトルが広がりを持ちます。風速に依存してスペクトル幅が変化します。
気象クラッタ（ウェザークラッタ）: 雨、雪、霧などの降水粒子からの反射。風に流されるため、風速に応じたドップラーシフトとスペクトル広がりを持ちます。
クラッタの統計的特性
クラッタの振幅分布は、一般に以下のモデルで記述されます。
レイリー分布: 多数のランダム散乱体からの反射の重ね合わせとして、振幅がレイリー分布に従うモデル。平地の地表面クラッタや海面クラッタに適用されます。
$$ p(a) = \frac{a}{\sigma^2} \exp\!\left(-\frac{a^2}{2\sigma^2}\right), \quad a \geq 0 $$
クラッタのドップラースペクトル: 静止クラッタのドップラースペクトルは $f_d = 0$ に集中しますが、内部運動やアンテナスキャンの影響でスペクトル幅を持ちます。ガウス型スペクトルでモデル化されることが多いです。
$$ S_c(f) = \frac{P_c}{\sqrt{2\pi}\sigma_c} \exp\!\left(-\frac{f^2}{2\sigma_c^2}\right) $$
ここで $P_c$ はクラッタ電力、$\sigma_c$ はスペクトル幅（標準偏差）です。
MTIの基本原理
パルス間位相変化
パルスレーダーでは、パルス繰り返し間隔（PRI: Pulse Repetition Interval）$T_r$ ごとにパルスを送信し、同じ距離セルからの反射信号を受信します。
距離 $R$ にある目標からの受信信号の位相は次の通りです。
$$ \phi = -\frac{4\pi R}{\lambda} $$
目標が速度 $v$ でレーダーに接近している場合、$n$ 番目のパルスにおける距離は次のようになります。
$$ R_n = R_0 – v \cdot n T_r $$
したがって、$n$ 番目のパルスでの受信位相は次の通りです。
$$ \phi_n = -\frac{4\pi}{\lambda}(R_0 – v \cdot nT_r) = -\frac{4\pi R_0}{\lambda} + \frac{4\pi v}{\lambda} \cdot nT_r $$
パルス間の位相変化は次のようになります。
$$ \Delta\phi = \phi_{n+1} – \phi_n = \frac{4\pi v T_r}{\lambda} = 2\pi f_d T_r $$
ここで $f_d = 2v/\lambda$ はドップラー周波数です。
MTIの基本アイデア
静止クラッタ: $v = 0$ なので $\Delta\phi = 0$。パルス間で位相が変化しません。
移動目標: $v \neq 0$ なので $\Delta\phi \neq 0$。パルス間で位相が変化します。
したがって、連続するパルスの差をとれば、静止クラッタは打ち消され、移動目標の信号だけが残ります。これがMTIの基本原理です。
単パルスキャンセラ（2パルスMTI）
定義
最も単純なMTIフィルタは、隣接する2パルスの差をとる単パルスキャンセラ（Single Delay Canceller）です。
$n$ 番目のパルスの受信信号を $x[n]$ とすると、出力 $y[n]$ は次の通りです。
$$ y[n] = x[n] – x[n-1] $$
Z変換による伝達関数
Z変換を適用すると、伝達関数は次のようになります。
$$ H(z) = 1 – z^{-1} $$
これは1次のFIR（有限インパルス応答）フィルタです。
周波数応答の導出
$z = e^{j\omega}$（$\omega = 2\pi f / f_{\mathrm{PRF}}$ は正規化角周波数、$f_{\mathrm{PRF}} = 1/T_r$ はPRF）を代入して周波数応答を求めます。
$$ \begin{align} H(e^{j\omega}) &= 1 – e^{-j\omega} \\ &= e^{-j\omega/2}\left(e^{j\omega/2} – e^{-j\omega/2}\right) \\ &= e^{-j\omega/2} \cdot 2j\sin\!\left(\frac{\omega}{2}\right) \end{align} $$
振幅応答は次の通りです。
$$ |H(e^{j\omega})| = 2\left|\sin\!\left(\frac{\omega}{2}\right)\right| $$
ドップラー周波数 $f_d$ で表すと、$\omega = 2\pi f_d T_r$ なので次のようになります。
$$ \boxed{|H(f_d)| = 2\left|\sin(\pi f_d T_r)\right|} $$
周波数応答の特徴
この周波数応答から、以下の特徴が読み取れます。
$f_d = 0$（静止物体）: $|H| = 0$。静止クラッタは完全に除去されます。
$f_d = f_{\mathrm{PRF}}/2$（最大応答）: $|H| = 2$。PRFの半分の周波数で最大応答を示します。
$f_d = f_{\mathrm{PRF}}$: $|H| = 0$。PRFの整数倍の周波数でも応答がゼロになります。
ブラインド速度
ブラインド速度の定義
単パルスキャンセラの周波数応答 $|H(f_d)| = 2|\sin(\pi f_d T_r)|$ は、$f_d = k \cdot f_{\mathrm{PRF}}$（$k = 0, 1, 2, \ldots$）でゼロになります。$f_d = 0$ は静止クラッタの除去に対応しますが、$f_d = k \cdot f_{\mathrm{PRF}}$（$k \geq 1$）は移動目標の見落としを意味します。
この現象に対応する速度をブラインド速度（Blind Speed）と呼びます。$f_d = 2v/\lambda$ の関係を用いると次のようになります。
$$ f_d = k \cdot f_{\mathrm{PRF}} = \frac{k}{T_r} $$
$$ \frac{2v_{\mathrm{blind}}^{(k)}}{\lambda} = \frac{k}{T_r} $$
$$ \boxed{v_{\mathrm{blind}}^{(k)} = \frac{k\lambda}{2T_r} = \frac{k\lambda f_{\mathrm{PRF}}}{2}} $$
最初のブラインド速度（$k = 1$）は次の通りです。
$$ v_{\mathrm{blind}} = \frac{\lambda}{2T_r} = \frac{\lambda f_{\mathrm{PRF}}}{2} $$
数値例
Xバンドレーダー（$f_c = 10 \, \mathrm{GHz}$, $\lambda = 0.03 \, \mathrm{m}$）でPRF $= 1{,}000 \, \mathrm{Hz}$ の場合を考えます。
$$ v_{\mathrm{blind}} = \frac{0.03 \times 1000}{2} = 15 \, \mathrm{m/s} = 54 \, \mathrm{km/h} $$
これは自動車程度の速度であり、多くの移動目標がブラインド速度付近に落ちてしまう可能性があります。
ブラインド速度の対策
ブラインド速度に対する主な対策は以下の通りです。
-
PRFを高くする: $v_{\mathrm{blind}} = \lambda f_{\mathrm{PRF}}/2$ なので、PRFを上げるとブラインド速度も上がります。ただし、最大非曖昧距離 $R_{\mathrm{ua}} = c/(2f_{\mathrm{PRF}})$ が短くなるトレードオフがあります。
-
スタガードPRF: パルス間隔を一定にせず、複数の異なるPRIを交互に使用する方法です。詳細は後述します。
2パルスキャンセラと3パルスキャンセラ
2パルスキャンセラ（ダブルキャンセラ）
単パルスキャンセラを2段直列に接続したものが2パルスキャンセラ（Double Delay Canceller）です。
$$ H_2(z) = (1 – z^{-1})^2 = 1 – 2z^{-1} + z^{-2} $$
出力は次のようになります。
$$ y[n] = x[n] – 2x[n-1] + x[n-2] $$
振幅応答を導出します。
$$ \begin{align} |H_2(e^{j\omega})| &= |1 – e^{-j\omega}|^2 \\ &= \left(2\left|\sin\!\left(\frac{\omega}{2}\right)\right|\right)^2 \\ &= 4\sin^2\!\left(\frac{\omega}{2}\right) \end{align} $$
ドップラー周波数で表すと次の通りです。
$$ |H_2(f_d)| = 4\sin^2(\pi f_d T_r) $$
3パルスキャンセラ
3パルスキャンセラは $(1 – z^{-1})^3$ に対応します。
$$ H_3(z) = 1 – 3z^{-1} + 3z^{-2} – z^{-3} $$
$$ |H_3(f_d)| = 8|\sin(\pi f_d T_r)|^3 $$
比較
| フィルタ | 伝達関数 | 振幅応答 | $f_d = 0$ 付近の減衰 | ノッチ幅 |
|---|---|---|---|---|
| 1段 | $1 – z^{-1}$ | $2|\sin(\pi f_d T_r)|$ | $\propto f_d$ | 狭い |
| 2段 | $(1-z^{-1})^2$ | $4\sin^2(\pi f_d T_r)$ | $\propto f_d^2$ | 中程度 |
| 3段 | $(1-z^{-1})^3$ | $8|\sin(\pi f_d T_r)|^3$ | $\propto f_d^3$ | 広い |
段数を増やすとゼロドップラー付近のクラッタ除去性能は向上しますが、ノッチ（応答がゼロになる帯域）が広くなり、低速目標も除去されるという代償があります。
スタガードPRFによるブラインド速度改善
原理
スタガードPRF（Staggered PRF）は、連続するパルスの間隔を一定にせず、2つ以上の異なるPRI $T_{r1}, T_{r2}$ を交互に使用する方法です。
2つのPRIを使う場合、それぞれのPRIに対するブラインド速度は次の通りです。
$$ v_{\mathrm{blind},1} = \frac{\lambda}{2T_{r1}}, \quad v_{\mathrm{blind},2} = \frac{\lambda}{2T_{r2}} $$
全体としてのブラインド速度が消えるためには、$v_{\mathrm{blind},1}$ と $v_{\mathrm{blind},2}$ が同時にゼロ応答を与える速度が十分高くなるように $T_{r1}$ と $T_{r2}$ を選びます。
スタガー比
スタガー比を $r = T_{r2}/T_{r1}$ と定義します。有効なブラインド速度の改善のために、$r$ を適切に選ぶ必要があります。一般に、$r$ は互いに素な整数の比（例: $r = 25/26$）が推奨されます。
2つのPRIのキャンセラの周波数応答を掛け合わせたものが全体の応答です。
$$ |H_{\mathrm{stag}}(f_d)| = |H_1(f_d)| \cdot |H_2(f_d)| = 4|\sin(\pi f_d T_{r1})| \cdot |\sin(\pi f_d T_{r2})| $$
$T_{r1}$ と $T_{r2}$ が異なるため、両方の $\sin$ が同時にゼロになるドップラー周波数が高くなり、実質的にブラインド速度が改善されます。
MTI改善係数
定義
MTI改善係数（Improvement Factor, $I_F$）は、MTI処理の有効性を定量化する指標です。次のように定義されます。
$$ \boxed{I_F = \frac{\mathrm{SIR}_{\mathrm{out}}}{\mathrm{SIR}_{\mathrm{in}}}} $$
ここで $\mathrm{SIR}_{\mathrm{in}}$ はMTI入力のクラッタに対する信号対干渉比（Signal-to-Interference Ratio）、$\mathrm{SIR}_{\mathrm{out}}$ はMTI出力のSIRです。
計算
入力においてクラッタ電力が $P_c$、目標信号電力が $P_s$ とすると、入力SIRは次の通りです。
$$ \mathrm{SIR}_{\mathrm{in}} = \frac{P_s}{P_c} $$
MTIフィルタの周波数応答を $H(f)$ とすると、出力のクラッタ電力は次の通りです。
$$ P_{c,\mathrm{out}} = \int_{-f_{\mathrm{PRF}}/2}^{f_{\mathrm{PRF}}/2} S_c(f) |H(f)|^2 df $$
ここで $S_c(f)$ はクラッタのドップラースペクトルです。ドップラー周波数 $f_d$ の目標に対する出力信号電力は次の通りです。
$$ P_{s,\mathrm{out}} = P_s \cdot |H(f_d)|^2 $$
したがって、MTI改善係数は次のように計算されます。
$$ \begin{align} I_F &= \frac{P_{s,\mathrm{out}} / P_{c,\mathrm{out}}}{P_s / P_c} \\ &= \frac{P_c \cdot |H(f_d)|^2}{\int_{-f_{\mathrm{PRF}}/2}^{f_{\mathrm{PRF}}/2} S_c(f) |H(f)|^2 df} \end{align} $$
理想クラッタ（単一スペクトル線）の場合
クラッタが完全に静止（$f = 0$ に集中）している理想的な場合、$S_c(f) = P_c \cdot \delta(f)$ です。単パルスキャンセラでは $H(0) = 0$ なので $P_{c,\mathrm{out}} = 0$ となり、改善係数は無限大になります。
有限スペクトル幅クラッタの場合
現実にはクラッタのスペクトルは有限幅を持ちます。ガウス型クラッタスペクトル $S_c(f) \propto \exp(-f^2/(2\sigma_c^2))$ に対する単パルスキャンセラの改善係数を計算しましょう。
$\sigma_c \ll f_{\mathrm{PRF}}$（クラッタのスペクトル幅がPRFに比べて十分狭い）という典型的な条件のもと、$\sin(\pi f T_r) \approx \pi f T_r$ の近似が使えます。
$$ |H(f)|^2 = 4\sin^2(\pi f T_r) \approx 4\pi^2 f^2 T_r^2 $$
出力クラッタ電力は次のようになります。
$$ \begin{align} P_{c,\mathrm{out}} &= \int_{-\infty}^{\infty} \frac{P_c}{\sqrt{2\pi}\sigma_c} \exp\!\left(-\frac{f^2}{2\sigma_c^2}\right) \cdot 4\pi^2 f^2 T_r^2 \, df \\ &= 4\pi^2 T_r^2 P_c \int_{-\infty}^{\infty} \frac{f^2}{\sqrt{2\pi}\sigma_c} \exp\!\left(-\frac{f^2}{2\sigma_c^2}\right) df \\ &= 4\pi^2 T_r^2 P_c \cdot \sigma_c^2 \end{align} $$
最後の等号では、ガウス分布の2次モーメント $\int f^2 \cdot \mathcal{N}(0, \sigma_c^2) df = \sigma_c^2$ を用いました。
最大応答（$f_d = f_{\mathrm{PRF}}/2$）での目標に対する改善係数は次の通りです。
$$ I_F = \frac{P_c \cdot |H(f_{\mathrm{PRF}}/2)|^2}{P_{c,\mathrm{out}}} = \frac{P_c \cdot 4}{4\pi^2 T_r^2 P_c \sigma_c^2} = \frac{1}{\pi^2 T_r^2 \sigma_c^2} $$
$$ \boxed{I_F = \frac{1}{\pi^2 T_r^2 \sigma_c^2}} \quad (\text{単パルスキャンセラ, ガウスクラッタ}) $$
この結果から、クラッタのスペクトル幅 $\sigma_c$ が狭いほど（クラッタが静止に近いほど）MTI改善係数が大きくなることがわかります。
Pythonでのシミュレーション
クラッタ＋移動目標信号の生成
import numpy as np
import matplotlib.pyplot as plt
# --- パラメータ設定 ---
fc = 10e9 # 搬送波周波数 [Hz] (Xバンド)
lam = 3e8 / fc # 波長 [m]
PRF = 1000 # パルス繰り返し周波数 [Hz]
Tr = 1.0 / PRF # パルス繰り返し間隔 [s]
N_pulses = 64 # パルス数
# ドップラー周波数軸
f_axis = np.linspace(-PRF/2, PRF/2, 1000)
# --- 単パルスキャンセラの周波数応答 ---
H1 = 2 * np.abs(np.sin(np.pi * f_axis * Tr))
# 2パルスキャンセラ
H2 = 4 * np.sin(np.pi * f_axis * Tr)**2
# 3パルスキャンセラ
H3 = 8 * np.abs(np.sin(np.pi * f_axis * Tr))**3
# --- 周波数応答の可視化 ---
fig, ax = plt.subplots(figsize=(12, 6))
ax.plot(f_axis, 20*np.log10(H1/np.max(H1) + 1e-10), 'b-', linewidth=2, label='1-stage: $1 - z^{-1}$')
ax.plot(f_axis, 20*np.log10(H2/np.max(H2) + 1e-10), 'r-', linewidth=2, label='2-stage: $(1 - z^{-1})^2$')
ax.plot(f_axis, 20*np.log10(H3/np.max(H3) + 1e-10), 'g-', linewidth=2, label='3-stage: $(1 - z^{-1})^3$')
ax.set_xlabel('Doppler Frequency [Hz]', fontsize=12)
ax.set_ylabel('Normalized Amplitude [dB]', fontsize=12)
ax.set_title('MTI Filter Frequency Response', fontsize=14)
ax.set_ylim([-60, 5])
ax.grid(True, alpha=0.3)
ax.legend(fontsize=12)
# ブラインド速度の表示
v_blind = lam * PRF / 2
ax.set_xlabel(f'Doppler Frequency [Hz]\n(Blind speed = {v_blind:.1f} m/s = {v_blind*3.6:.1f} km/h)')
plt.tight_layout()
plt.show()
print(f"波長 λ = {lam*100:.1f} cm")
print(f"ブラインド速度 = {v_blind:.1f} m/s = {v_blind*3.6:.1f} km/h")
MTIフィルタの適用シミュレーション
import numpy as np
import matplotlib.pyplot as plt
# --- パラメータ設定 ---
fc = 10e9 # 搬送波周波数 [Hz]
lam = 3e8 / fc # 波長 [m]
PRF = 1000 # パルス繰り返し周波数 [Hz]
Tr = 1.0 / PRF # パルス繰り返し間隔 [s]
N_pulses = 64 # パルス数
N_range = 200 # レンジセル数
np.random.seed(42)
# --- クラッタの生成 ---
# ガウス型ドップラースペクトルを持つクラッタ
sigma_c = 5 # クラッタのスペクトル幅 [Hz]
clutter_power_dB = 40 # クラッタ電力 [dB]
clutter_power = 10**(clutter_power_dB / 10)
# 各レンジセルでクラッタ信号を生成
# クラッタはパルス間で相関を持つ（低ドップラー）
clutter = np.zeros((N_pulses, N_range), dtype=complex)
for r in range(N_range):
# ガウス型のパルス間相関を持つクラッタ
# パルスインデックス
n = np.arange(N_pulses)
# ランダムなクラッタ位相
clutter_phase = np.random.uniform(0, 2*np.pi)
# ゆっくり変動するクラッタ（スペクトル幅 sigma_c）
clutter_doppler = sigma_c * np.random.randn() # クラッタのドップラー
amplitude = np.sqrt(clutter_power) * (0.5 + np.random.rand())
clutter[:, r] = amplitude * np.exp(
1j * (clutter_phase + 2 * np.pi * clutter_doppler * n * Tr)
)
# --- 移動目標の生成 ---
# (レンジセル, 速度 [m/s], 振幅)
target_list = [
(50, 30, 1.0), # 目標1: 30 m/s (108 km/h)
(100, 80, 0.8), # 目標2: 80 m/s (288 km/h)
(150, 5, 1.5), # 目標3: 5 m/s (18 km/h, 低速)
]
targets = np.zeros((N_pulses, N_range), dtype=complex)
for (r_cell, v_tgt, amp) in target_list:
fd = 2 * v_tgt / lam # ドップラー周波数
n = np.arange(N_pulses)
target_phase = np.random.uniform(0, 2*np.pi)
targets[:, r_cell] = amp * np.exp(1j * (target_phase + 2 * np.pi * fd * n * Tr))
# --- 雑音の生成 ---
noise_power = 0.01
noise = np.sqrt(noise_power / 2) * (
np.random.randn(N_pulses, N_range) + 1j * np.random.randn(N_pulses, N_range)
)
# --- 受信信号 ---
received = clutter + targets + noise
# --- MTIフィルタの適用 ---
# 単パルスキャンセラ: y[n] = x[n] - x[n-1]
mti_1stage = received[1:, :] - received[:-1, :]
# 2パルスキャンセラ: y[n] = x[n] - 2x[n-1] + x[n-2]
mti_2stage = received[2:, :] - 2*received[1:-1, :] + received[:-2, :]
# --- 各パルスの電力を平均してレンジプロファイルを生成 ---
power_before = np.mean(np.abs(received)**2, axis=0)
power_1stage = np.mean(np.abs(mti_1stage)**2, axis=0)
power_2stage = np.mean(np.abs(mti_2stage)**2, axis=0)
# --- 可視化 ---
fig, axes = plt.subplots(3, 1, figsize=(14, 12))
range_cells = np.arange(N_range)
# MTI処理前
axes[0].plot(range_cells, 10*np.log10(power_before + 1e-10), 'b-', linewidth=1)
axes[0].set_ylabel('Power [dB]', fontsize=11)
axes[0].set_title('Before MTI (Clutter + Targets + Noise)', fontsize=13)
axes[0].grid(True, alpha=0.3)
for (r_cell, v_tgt, _) in target_list:
axes[0].axvline(x=r_cell, color='red', linestyle='--', alpha=0.5,
label=f'Target: v={v_tgt} m/s' if r_cell == target_list[0][0] else '')
axes[0].legend()
# 単パルスキャンセラ後
axes[1].plot(range_cells, 10*np.log10(power_1stage + 1e-10), 'r-', linewidth=1)
axes[1].set_ylabel('Power [dB]', fontsize=11)
axes[1].set_title('After 1-Stage MTI: $y[n] = x[n] - x[n-1]$', fontsize=13)
axes[1].grid(True, alpha=0.3)
for (r_cell, v_tgt, _) in target_list:
axes[1].axvline(x=r_cell, color='red', linestyle='--', alpha=0.5)
# 2パルスキャンセラ後
axes[2].plot(range_cells, 10*np.log10(power_2stage + 1e-10), 'g-', linewidth=1)
axes[2].set_xlabel('Range Cell', fontsize=11)
axes[2].set_ylabel('Power [dB]', fontsize=11)
axes[2].set_title('After 2-Stage MTI: $y[n] = x[n] - 2x[n-1] + x[n-2]$', fontsize=13)
axes[2].grid(True, alpha=0.3)
for (r_cell, v_tgt, _) in target_list:
axes[2].axvline(x=r_cell, color='red', linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()
ドップラースペクトルの比較
import numpy as np
import matplotlib.pyplot as plt
# --- パラメータ設定 ---
fc = 10e9
lam = 3e8 / fc
PRF = 1000
Tr = 1.0 / PRF
N_pulses = 64
np.random.seed(42)
# --- 単一レンジセルの信号を生成 ---
# クラッタ（低ドップラー）
sigma_c = 5 # Hz
clutter_amp = 100
clutter_fd = sigma_c * np.random.randn()
n = np.arange(N_pulses)
clutter_sig = clutter_amp * np.exp(1j * 2 * np.pi * clutter_fd * n * Tr)
# 移動目標
target_v = 40 # m/s
target_fd = 2 * target_v / lam
target_amp = 1.0
target_sig = target_amp * np.exp(1j * 2 * np.pi * target_fd * n * Tr)
# 雑音
noise = 0.1 * (np.random.randn(N_pulses) + 1j * np.random.randn(N_pulses))
# 合計信号
signal = clutter_sig + target_sig + noise
# --- MTI処理 ---
signal_mti = signal[1:] - signal[:-1]
# --- ドップラースペクトル ---
N_fft = 256
doppler_before = np.fft.fftshift(np.fft.fft(signal, N_fft))
doppler_after = np.fft.fftshift(np.fft.fft(signal_mti, N_fft))
f_doppler = np.fft.fftshift(np.fft.fftfreq(N_fft, d=Tr))
# 速度軸
v_axis = f_doppler * lam / 2
# --- 可視化 ---
fig, axes = plt.subplots(1, 2, figsize=(16, 6))
# MTI前のドップラースペクトル
axes[0].plot(v_axis, 20*np.log10(np.abs(doppler_before)/np.max(np.abs(doppler_before)) + 1e-10),
'b-', linewidth=1.5)
axes[0].set_xlabel('Velocity [m/s]', fontsize=12)
axes[0].set_ylabel('Amplitude [dB]', fontsize=12)
axes[0].set_title('Doppler Spectrum (Before MTI)', fontsize=13)
axes[0].set_ylim([-60, 5])
axes[0].axvline(x=target_v, color='red', linestyle='--', label=f'Target: v={target_v} m/s')
axes[0].axvline(x=0, color='gray', linestyle=':', label='Clutter (v≈0)')
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)
# MTI後のドップラースペクトル
axes[1].plot(v_axis, 20*np.log10(np.abs(doppler_after)/np.max(np.abs(doppler_after)) + 1e-10),
'r-', linewidth=1.5)
axes[1].set_xlabel('Velocity [m/s]', fontsize=12)
axes[1].set_ylabel('Amplitude [dB]', fontsize=12)
axes[1].set_title('Doppler Spectrum (After MTI)', fontsize=13)
axes[1].set_ylim([-60, 5])
axes[1].axvline(x=target_v, color='red', linestyle='--', label=f'Target: v={target_v} m/s')
axes[1].axvline(x=0, color='gray', linestyle=':', label='Clutter (v≈0)')
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"目標のドップラー周波数: f_d = {target_fd:.1f} Hz")
print(f"目標の速度: v = {target_v:.1f} m/s = {target_v*3.6:.1f} km/h")
print(f"クラッタ対目標比 (入力): {20*np.log10(clutter_amp/target_amp):.1f} dB")
スタガードPRFの効果
import numpy as np
import matplotlib.pyplot as plt
# --- パラメータ設定 ---
fc = 10e9
lam = 3e8 / fc
PRF_nom = 1000 # 公称PRF [Hz]
Tr_nom = 1.0 / PRF_nom
# スタガー比
stagger_ratio = 25.0 / 26.0
Tr1 = Tr_nom
Tr2 = Tr_nom * stagger_ratio
f_axis = np.linspace(-PRF_nom/2, PRF_nom/2, 2000)
# --- 非スタガード（均一PRI）の応答 ---
H_uniform = 2 * np.abs(np.sin(np.pi * f_axis * Tr_nom))
# --- スタガードPRIの応答 ---
# 2つのPRIに対するキャンセラの積
H_stag1 = 2 * np.abs(np.sin(np.pi * f_axis * Tr1))
H_stag2 = 2 * np.abs(np.sin(np.pi * f_axis * Tr2))
H_staggered = H_stag1 * H_stag2
# 速度軸
v_axis = f_axis * lam / 2
# --- 可視化 ---
fig, axes = plt.subplots(2, 1, figsize=(14, 10))
# 周波数応答の比較
axes[0].plot(f_axis, 20*np.log10(H_uniform/np.max(H_uniform) + 1e-10),
'b-', linewidth=2, label='Uniform PRI')
axes[0].plot(f_axis, 20*np.log10(H_staggered/np.max(H_staggered) + 1e-10),
'r-', linewidth=2, label=f'Staggered PRI (ratio={stagger_ratio:.4f})')
axes[0].set_xlabel('Doppler Frequency [Hz]', fontsize=12)
axes[0].set_ylabel('Normalized Amplitude [dB]', fontsize=12)
axes[0].set_title('MTI Frequency Response: Uniform vs Staggered PRI', fontsize=14)
axes[0].set_ylim([-60, 5])
axes[0].legend(fontsize=12)
axes[0].grid(True, alpha=0.3)
# 速度応答の比較
axes[1].plot(v_axis, 20*np.log10(H_uniform/np.max(H_uniform) + 1e-10),
'b-', linewidth=2, label='Uniform PRI')
axes[1].plot(v_axis, 20*np.log10(H_staggered/np.max(H_staggered) + 1e-10),
'r-', linewidth=2, label=f'Staggered PRI (ratio={stagger_ratio:.4f})')
axes[1].set_xlabel('Velocity [m/s]', fontsize=12)
axes[1].set_ylabel('Normalized Amplitude [dB]', fontsize=12)
axes[1].set_title('MTI Velocity Response: Uniform vs Staggered PRI', fontsize=14)
axes[1].set_ylim([-60, 5])
axes[1].legend(fontsize=12)
axes[1].grid(True, alpha=0.3)
# ブラインド速度の表示
v_blind = lam * PRF_nom / 2
axes[1].axvline(x=v_blind, color='gray', linestyle=':', alpha=0.7,
label=f'1st Blind Speed = {v_blind:.1f} m/s')
axes[1].axvline(x=-v_blind, color='gray', linestyle=':', alpha=0.7)
axes[1].legend(fontsize=11)
plt.tight_layout()
plt.show()
print(f"均一PRI: ブラインド速度 = {v_blind:.1f} m/s = {v_blind*3.6:.1f} km/h")
print(f"スタガーPRI: T_r1 = {Tr1*1e3:.3f} ms, T_r2 = {Tr2*1e3:.3f} ms")
上記のコードを実行すると、MTIフィルタの周波数応答（1段/2段/3段の比較）、クラッタと移動目標を含む信号へのMTI処理の効果、ドップラースペクトル上でのクラッタ除去の様子、そしてスタガードPRFによるブラインド速度改善の効果が可視化されます。MTI処理によってクラッタが大幅に抑圧され、移動目標が浮かび上がる様子を確認できます。
まとめ
本記事では、MTI（移動目標指示）とクラッタ除去の理論を数式で丁寧に導出しました。
- クラッタは地面・海面・気象からの不要反射であり、目標信号よりはるかに強いことが多い
- MTIはパルス間の位相変化 $\Delta\phi = 2\pi f_d T_r$ を利用して移動目標を検出する
- 単パルスキャンセラ $H(z) = 1 – z^{-1}$ の振幅応答は $2|\sin(\pi f_d T_r)|$ であり、$f_d = 0$ でゼロ（クラッタ除去）
- ブラインド速度 $v_{\mathrm{blind}} = \lambda f_{\mathrm{PRF}}/2$ はMTIフィルタの応答がゼロになる速度であり、目標の見落としにつながる
- 多段キャンセラはゼロドップラー付近のクラッタ除去性能を向上させるが、ノッチ幅も広がる
- スタガードPRFにより、異なるPRIのキャンセラを組み合わせてブラインド速度を改善できる
- MTI改善係数 $I_F$ はクラッタのスペクトル幅 $\sigma_c$ が狭いほど大きくなる
次のステップとして、以下の記事も参考にしてください。