定偽警報率 (CFAR) 検出

ライブスクリプトを開く

この例では、定偽警報率 (CFAR) 検出について紹介し、Phased Array System Toolbox™ の CFARDetector と CFARDetector2D を使用してセル平均化 CFAR 検出を実行する方法を示します。

はじめに

レーダーシステムが実行する重要なタスクの 1 つは、ターゲットの検出です。検出自体は非常に簡単です。信号をしきい値と比較します。したがって、検出における本質的な課題は、適切なしきい値を見つけることです。一般に、しきい値は検出と偽警報の確率の両方の関数となります。

多くのフェーズドアレイシステムでは、誤検出に関連するコストが発生するため、検出確率を最大化するだけでなく、偽警報の確率を事前設定されたレベル未満に保つ検出しきい値をもつことが望まれます。

検出しきい値を決定する方法については、多くの文献が存在します。よく知られたいくつかの結果としては、Signal Detection in White Gaussian NoiseおよびSignal Detection Using Multiple Samplesの例が参考になるでしょう。ただし、これらの古典的な結果はすべて理論的な確率に基づいており、既知の分散 (パワー) をもつホワイトガウスノイズに限定されています。実際の用途においては、カラードノイズである場合が多く、そのパワーも未知です。

CFAR テクノロジーは、これらの問題に対処します。CFAR において、テスト対象セル (CUT) と呼ばれる特定のセルに対して検出が必要な場合、隣接するセルからノイズパワーが推定されます。検出しきい値 $T$ は次で与えられます。

$T = α P_{n}$

ここで、 $P_{n}$ はノイズパワー推定値であり、 $α$ はしきい値の係数と呼ばれるスケーリング係数です。

この方程式から、しきい値がデータに適応することが明らかです。適切なしきい値係数 $α$ を使用すれば、結果として得られる偽警報の確率を一定に保つことができることが示されており、これが CFAR (定偽警報率) という名前の由来です。

セル平均化 CFAR 検出

セル平均化 CFAR 検出器は、おそらく最も広く使用されている CFAR 検出器です。他の CFAR 手法のベースライン比較としても使用されます。セル平均化 CFAR 検出器では、CUT 前後の先行するセルと後続するセル (トレーニングセルと呼ばれる) の両方からノイズサンプルを抽出します。ノイズ推定値は次のように計算できます。[1]

$P_{n} = \frac{1}{N} \sum_{m = 1}^{N} x_{m}$

ここで、 $N$ はトレーニングセルの数、 $x_{m}$ は各トレーニングセル内のサンプルです。 $x_{m}$ が二乗検波器の出力である場合、 $P_{n}$ は推定されたノイズパワーを表します。一般に、先行するトレーニングセルと後続するトレーニングセルの数は同じにします。ガードセルは、CUT の前方と後方の両方に隣接して配置されます。これらのガードセルの目的は、信号成分がトレーニングセルに漏れてノイズ推定に悪影響を与えるのを防ぐことです。

次の図は、1 次元の場合におけるこれらのセルの関係を示しています。

上記のセル平均化 CFAR 検出器において、検出器に渡されるデータが単一のパルスからのものである (つまり、パルス積分が含まれない) と仮定すると、しきい値係数は次のように記述できます。[1]

$α = N (P_{f a}^{- 1 / N} - 1)$

ここで、 $P_{f a}$ は目的の偽警報率です。

自動しきい値係数を使用した CFAR 検出

この例の残りの部分では、Phased Array System Toolbox を使用してセル平均化 CFAR 検出を実行する方法を示します。簡潔さと汎用性を損なわないために、ここでもノイズはホワイトガウスであると仮定します。これにより、CFAR と古典的な検出理論を比較できるようになります。

次のコマンドを使用して CFAR 検出器をインスタンス化できます。

cfar = phased.CFARDetector('NumTrainingCells',20,'NumGuardCells',2);

この検出器では、合計 20 個のトレーニングセルと 2 個のガードセルを使用します。つまり、CUT の両側それぞれに 10 個のトレーニングセルと 1 個のガードセルがあることになります。前述のように、信号がパルス積分のない二乗検波器からのものであると仮定すれば、しきい値はトレーニングセルの数と目的の偽警報の確率に基づいて計算できます。目的の偽警報率が 0.001 であると仮定する場合、この計算を実行できるように CFAR 検出器を次のように構成できます。

exp_pfa = 1e-3;
cfar.ThresholdFactor = 'Auto';
cfar.ProbabilityFalseAlarm = exp_pfa;

構成された CFAR 検出器を以下に示します。

cfar

cfar = 
  phased.CFARDetector with properties:

                   Method: 'CA'
            NumGuardCells: 2
         NumTrainingCells: 20
          ThresholdFactor: 'Auto'
    ProbabilityFalseAlarm: 1.0000e-03
             OutputFormat: 'CUT result'
      ThresholdOutputPort: false
     NoisePowerOutputPort: false

ここで、入力データをシミュレーションします。今回の焦点は、CFAR 検出器が偽警報率を特定の値未満に維持できることを示すことにあるため、それらのセル内のノイズサンプルのみをシミュレーションします。設定は次のとおりです。

データシーケンスは 23 サンプルの長さで、CUT はセル 12。これにより、CUT の両側それぞれに 10 個のトレーニングセルと 1 個のガードセルが配置される。
10 万回のモンテカルロ試行を使用して、偽警報率を計算する。

rs = RandStream('mt19937ar','Seed',2010);
npower = db2pow(-10);  % Assume 10dB SNR ratio

Ntrials = 1e5;
Ncells = 23;
CUTIdx = 12;

% Noise samples after a square law detector
rsamp = randn(rs,Ncells,Ntrials)+1i*randn(rs,Ncells,Ntrials);   
x = abs(sqrt(npower/2)*rsamp).^2;

検出を実行するには、データを検出器に渡します。この例では、CUT は 1 つだけなので、出力はすべての試行の検出結果を含む logical ベクトルになります。結果が true の場合、対応する試行にターゲットが存在することを意味します。この例では、ノイズのみを渡しているため、すべての検出結果が偽警報になります。結果として得られる偽警報率は、偽警報の数と試行回数に基づいて計算できます。

x_detected = cfar(x,CUTIdx);
act_pfa = sum(x_detected)/Ntrials

act_pfa = 
9.4000e-04

この結果は、指定したとおり、結果として得られる偽警報の確率が 0.001 未満であることを示しています。

カスタムしきい値係数を使用した CFAR 検出

この例の前半で説明したように、CFAR 検出器が適切なしきい値係数を自動的に計算できるケースはごくわずかです。たとえば、前のシナリオを使用しても、データが検出器に入る前に 10 パルスの非コヒーレント積分を行った場合、自動しきい値では目的の偽警報率を提供できなくなります。

npower = db2pow(-10);  % Assume 10dB SNR ratio
xn = 0;
for m = 1:10
    rsamp = randn(rs,Ncells,Ntrials)+1i*randn(rs,Ncells,Ntrials);
    xn = xn + abs(sqrt(npower/2)*rsamp).^2;   % noncoherent integration
end
x_detected = cfar(xn,CUTIdx);
act_pfa = sum(x_detected)/Ntrials

act_pfa = 
0

結果として得られる偽警報率が 0 であるのに、それが偽警報率が 0.001 であることよりも悪いと考える理由について疑問に思うかもしれません。結局のところ、偽警報率 0 は優れた結果ではないのでしょうか。この疑問に対する答えは、偽警報の確率が下がると検出の確率も下がるという事実にあります。この場合、真の偽警報率が許容値を大幅に下回っているということは、検出しきい値の設定が高すぎることになります。目的の偽警報の確率を使用することで、より低いコスト (たとえば、より低い送信機電力) で同じ検出確率を達成できます。

ほとんどの場合、しきい値係数は、特定の環境とシステム構成に基づいて推定する必要があります。以下に示すように、カスタムしきい値係数を使用するように CFAR 検出器を構成できます。

release(cfar);
cfar.ThresholdFactor = 'Custom';

パルス積分の例を続行し、経験的データを使用したところ、目的の偽警報率を達成するためにカスタムしきい値係数として 2.35 を使用できることがわかりました。このしきい値を使用することで、結果として得られる偽警報率が期待値と一致することが確認できます。

cfar.CustomThresholdFactor = 2.35;
x_detected = cfar(xn,CUTIdx);
act_pfa = sum(x_detected)/Ntrials

act_pfa = 
9.6000e-04

CFAR 検出しきい値

CFAR 検出は、セル内の入力信号レベルがしきい値レベルを超えたときに発生します。各セルのしきい値レベルは、しきい値係数と、トレーニングセルから推定されたノイズパワーによって決まります。定偽警報率を維持するために、検出しきい値はトレーニングセル内のノイズパワーに比例して増減します。ThresholdOutputPort プロパティを使用して、各検出に使用されたしきい値を出力するように CFAR 検出器を構成します。自動しきい値係数と 200 個のトレーニングセルを使用します。

release(cfar);
cfar.ThresholdOutputPort = true;
cfar.ThresholdFactor = 'Auto';
cfar.NumTrainingCells = 200;

次に、ノイズパワーが増加する二乗則入力信号を作成します。

rs = RandStream('mt19937ar','Seed',2010);
Npoints = 1e4;
rsamp = randn(rs,Npoints,1)+1i*randn(rs,Npoints,1);
ramp = linspace(1,10,Npoints)';
xRamp = abs(sqrt(npower*ramp./2).*rsamp).^2;

信号内のすべてのセルについて検出結果としきい値を計算します。

[x_detected,th] = cfar(xRamp,1:length(xRamp));

次に、CFAR しきい値を入力信号と比較します。

plot(1:length(xRamp),xRamp,1:length(xRamp),th,...
  find(x_detected),xRamp(x_detected),'o')
legend('Signal','Threshold','Detections','Location','Northwest')
xlabel('Time Index')
ylabel('Level')

Figure contains an axes object. The axes object with xlabel Time Index, ylabel Level contains 3 objects of type line. One or more of the lines displays its values using only markers These objects represent Signal, Threshold, Detections.

ここで、しきい値は信号のノイズパワーに応じて増加し、定偽警報率が維持されています。検出は、信号レベルがしきい値を超えたところで行われます。

CFAR と古典的なネイマン・ピアソン検出器の比較

このセクションでは、CFAR 検出器の性能を、ネイマン・ピアソンの原理を使用した古典的な検出理論と比較します。最初の例に戻り、真のノイズパワーが既知であると仮定すると、理論的なしきい値は次のように計算できます。

T_ideal = npower*db2pow(npwgnthresh(exp_pfa));

この古典的なネイマン・ピアソン検出器の偽警報率は、この理論的なしきい値を使用して計算できます。

act_Pfa_np = sum(x(CUTIdx,:)>T_ideal)/Ntrials

act_Pfa_np = 
9.5000e-04

ノイズパワーが既知であるため、古典的な検出理論によっても目的の偽警報率が達成されます。CFAR 検出器によって達成される偽警報率も同様です。

release(cfar);
cfar.ThresholdOutputPort = false;
cfar.NumTrainingCells = 20;
x_detected = cfar(x,CUTIdx);
act_pfa = sum(x_detected)/Ntrials

act_pfa = 
9.4000e-04

次に、両方の検出器が現場に展開されたと仮定します。ここで、ノイズパワーが予想よりも 1 dB 大きいとします。この場合、理論上のしきい値を使用すると、結果として得られる偽警報の確率は目的とする値の 4 倍に跳ね上がります。

npower = db2pow(-9);  % Assume 9dB SNR ratio
rsamp = randn(rs,Ncells,Ntrials)+1i*randn(rs,Ncells,Ntrials);
x = abs(sqrt(npower/2)*rsamp).^2;   
act_Pfa_np = sum(x(CUTIdx,:)>T_ideal)/Ntrials

act_Pfa_np = 
0.0041

これに対し、CFAR 検出器のパフォーマンスには影響がありません。

x_detected = cfar(x,CUTIdx);
act_pfa = sum(x_detected)/Ntrials

act_pfa = 
0.0011

したがって、CFAR 検出器はノイズパワーの不確実性に対してロバストであり、現場での用途に適しています。

最後に、カラードノイズが発生している状態で CFAR 検出を使用します。まず、古典的な検出しきい値をデータに適用します。

npower = db2pow(-10);
fcoeff = maxflat(10,'sym',0.2);
x = abs(sqrt(npower/2)*filter(fcoeff,1,rsamp)).^2;   % colored noise
act_Pfa_np = sum(x(CUTIdx,:)>T_ideal)/Ntrials

act_Pfa_np = 
0

結果として得られる偽警報率が要件を満たせていないことに注意してください。しかし、CFAR 検出器をカスタムしきい値係数と共に使用すれば、目的の偽警報率を得ることができます。

release(cfar);
cfar.ThresholdFactor = 'Custom';
cfar.CustomThresholdFactor = 12.85;
x_detected = cfar(x,CUTIdx);
act_pfa = sum(x_detected)/Ntrials

act_pfa = 
0.0010

レンジ-ドップラーイメージに対する CFAR 検出

前のセクションでは、CUT の先行するトレーニングセルと後続するトレーニングセルを単一次元で計算してノイズ推定を行いました。イメージに対して CFAR 検出を実行することもできます。セルはイメージ内のピクセルに対応し、ガードセルとトレーニングセルは CUT の周囲に帯状に配置されます。検出しきい値は、CUT の周囲にある四角形のトレーニングバンド内のセルから計算されます。

上図において、ガードバンドのサイズは [2 2]、トレーニングバンドのサイズは [4 3] です。サイズインデックスは、それぞれ行と列の次元における CUT の各側にあるセルの数を示しています。ガードバンドのサイズは、行と列の次元に沿って同じであるため、2 として定義することもできます。

次に、2 次元 CFAR 検出器を作成します。偽警報の確率として 1e-5 を使用し、ガードバンドサイズを 5 セル、トレーニングバンドサイズを 10 セルに指定します。

cfar2D = phased.CFARDetector2D('GuardBandSize',5,'TrainingBandSize',10,...
  'ProbabilityFalseAlarm',1e-5);

次に、レンジ-ドップラーイメージを読み込んでプロットします。イメージには、2 つの静止したターゲットからのリターンと、レーダーから遠ざかる 1 つのターゲットからのリターンが含まれています。

[resp,rngGrid,dopGrid] = helperRangeDoppler;

Figure contains an axes object. The axes object with title Range Doppler Map, xlabel Doppler (Hz), ylabel Range (m) contains an object of type image.

CFAR を使用して、レンジ-ドップラー空間でオブジェクトを検索し、検出マップをプロットします。-10 ～ 10 kHz、1000 ～ 4000 m の範囲を検索します。まず、この領域のテスト対象セルを定義します。

[~,rangeIndx] = min(abs(rngGrid-[1000 4000]));
[~,dopplerIndx] = min(abs(dopGrid-[-1e4 1e4]));
[columnInds,rowInds] = meshgrid(dopplerIndx(1):dopplerIndx(2),...
  rangeIndx(1):rangeIndx(2));
CUTIdx = [rowInds(:) columnInds(:)]';

各テスト対象セルについて検出結果を計算します。この例では、検索領域内の各ピクセルがセルになります。レンジ-ドップラーイメージの検出結果のマップをプロットします。

detections = cfar2D(resp,CUTIdx);
helperDetectionsMap(resp,rngGrid,dopGrid,rangeIndx,dopplerIndx,detections)

Figure contains an axes object. The axes object with title Range Doppler CFAR Detections, xlabel Doppler (Hz), ylabel Range (m) contains an object of type image.

3 つのオブジェクトが検出されています。時間経過に伴うレンジ-ドップラーイメージのデータキューブも同様に、cfar2D への入力信号として提供することができ、1 ステップで検出結果が計算されます。

まとめ

この例では、CFAR 検出器の基本的な概念を紹介しました。特に、Phased Array System Toolbox を使用し、信号とレンジ-ドップラーイメージに対してセル平均化 CFAR 検出を実行する方法について説明しました。セル平均化 CFAR 検出器によって提供されるパフォーマンスと、理論的に計算されたしきい値を備えた検出器によって提供されるパフォーマンスの比較は、CFAR 検出器が実際の現場での用途により適していることを明確に示しています。

参考文献

[1] Mark Richards, Fundamentals of Radar Signal Processing, McGraw Hill, 2005