奈曼-皮尔逊准则
如果各个假设的先验概率 P ( H i ) ( i = 0 , 1 ) P(H_i)(i=0,1) P(Hi)(i=0,1)已知,并且每种判决的代价函数已知,这时,可以用贝叶斯准则进行判决。
如果进一步假设正确判决不付出的代价,而错误判决代价相同,即 c 00 = c 11 = 0 , c 01 = c 10 = 1 c_{00}=c_{11}=0,c_{01}=c_{10}=1 c00=c11=0,c01=c10=1,这时可以用最小平均错误概率准则进行判决。
在不知道先验概率 P ( H i ) ( i = 0 , 1 ) P(H_i)(i=0,1) P(Hi)(i=0,1)的情况下,可以采用极小化极大准则进行判决。
但在很多场景中,确定先验概率及各类错误的代价都是及其困难的。例如,在雷达系统中,我们既不能预知先验概率 P ( H i ) ( i = 0 , 1 ) P(H_i)(i=0,1) P(Hi)(i=0,1),也无法对各种判决给定代价函数 c i j ( i , j = 0 , 1 ) c_{ij}(i,j=0,1) cij(i,j=0,1)。回避这个困难的最简单的方法是应用条件概率 P F P_F PF和 P D P_D PD。一般来说,我们总希望 P F P_F PF尽可能小而 P D P_D PD尽可能大。对大多数实际问题来说,这两者是互相矛盾的。我们可以限定其中一个概率而使另一个概率为最大或最小,这就是奈曼-皮尔逊准则,简称N-P准则。
一般对虚警 P F P_F PF要求较高,所以,先限定虚警 P F P_F PF,再去求最小漏报或最大检测,所以N-P准则有时也叫恒虚警检测。 恒虚警检测技术(CFAR, Constant False-Alarm Rate)是雷达系统在保持虚警概率恒定条件下对接收机输出的信号与噪声作判别以确定目标信号是否存在的技术。
定义
信号接收机输出端的信号用 x ( t ) x(t) x(t) 表示,这里存在两种情况:
-
噪声和信号同时存在: x ( t ) = s ( t ) + n ( t ) x(t)=s(t)+n(t) x(t)=s(t)+n(t)
-
只有噪声存在: x ( t ) = n ( t ) x(t)=n(t) x(t)=n(t)
用 H 0 H_{0} H0 和 H 1 H_{1} H1 分别表示接收机的无信号输入和有信号输入的假设;
用 D 0 D_{0} D0 和 D 1 D_{1} D1 分别表示检测器作出无信号和有信号的判决结果。
于是接收机的输入与检测器的判决将有四种情况: -
H 0 H_{0} H0 为真,判为 D 0 D_{0} D0 ,即接收机无信号输入,检测器判为无信号,称为正确不发现;
-
H 0 H_{0} H0 为真,判为 D 1 D_{1} D1 ,即接收机无信号输入,检测器判为有信号,称为虚警;
-
H 1 H_{1} H1 为真,判为 D 0 D_{0} D0 ,即接收机有信号输入,检测器判为无信号,称为漏警;
-
H 1 H_{1} H1 为真,判为 D 1 D_{1} D1 ,即接收机有信号输入,检测器判为有信号,称为正确检测;
其中第一种情况和第四种情况属于正确判决,其余两种属于错误判决。
用
p
(
z
∣
H
0
)
p\left(z \mid H_{0}\right)
p(z∣H0) 和
p
(
z
∣
H
1
)
p\left(z \mid H_{1}\right)
p(z∣H1) 分别表示无信号输入和有信号输入接收机时,接收机输出端的信号电平的概率密度函数;
用
Z
0
Z_{0}
Z0 和
Z
1
Z_{1}
Z1 分别表示检测器作出无信号和有信号判决的判决区域,当输入的电平在
Z
0
Z_{0}
Z0 区域判为无信号,在
Z
1
Z_{1}
Z1 区域判为有信号。
举例
数字通信系统中,备选假设为 H 1 H_{1} H1 时,信源输出 1 ;备选假设为 H 0 H_{0} H0 时,信源输出为 0 。信号在通信信道上迭加了噪声 n ( t ) n(t) n(t) ,附加噪 声 n ( t ) n(t) n(t) 为零均值,方差为1的高斯噪声。试构造一个 P F = 0.1 P_{F}=0.1 PF=0.1 的聂曼-皮尔逊接收机。
解:在 H 1 H_{1} H1 和 H 0 H_{0} H0 两种假设下,若 z z z 为接收信号, n n n 为噪声,则接收机的输出信号可以写为
- H 1 : z = 1 + n H_{1}: z=1+n H1:z=1+n
- H 0 : z = n H_{0}: z=n H0:z=n
由于 n n n 是高斯变量,均值为 0、方差为 1。在这两种假设下, z z z 的概率密度函数为
p
(
z
∣
H
0
)
=
1
2
π
exp
(
−
z
2
2
)
(1)
p\left(z \mid H_{0}\right)=\frac{1}{\sqrt{2 \pi}} \exp \left(-\frac{z^{2}}{2}\right)\tag{1}
p(z∣H0)=2π1exp(−2z2)(1)
p
(
z
∣
H
1
)
=
1
2
π
exp
[
−
(
z
−
1
)
2
2
]
(2)
p\left(z \mid H_{1}\right)=\frac{1}{\sqrt{2 \pi}} \exp \left[-\frac{(z-1)^{2}}{2}\right]\tag{2}
p(z∣H1)=2π1exp[−2(z−1)2](2)
似然比为
λ ( z ) = p ( z ∣ H 1 ) p ( z ∣ H 0 ) = exp ( z − 1 2 ) (3) \lambda(z)=\frac{p\left(z \mid H_{1}\right)}{p\left(z \mid H_{0}\right)}=\exp \left(z-\frac{1}{2}\right)\tag{3} λ(z)=p(z∣H0)p(z∣H1)=exp(z−21)(3)
那么判决规则为:若 λ ( z ) \lambda(z) λ(z) 大于 λ \lambda λ ,判决有信号输入;若 λ ( z ) \lambda(z) λ(z) 小于 λ \lambda λ 时,判决无信号输入。 用对数似然比的形式将上式写为
ln λ ( z ) = ln p ( z ∣ H 1 ) p ( z ∣ H 0 ) = z − 1 2 (4) \ln \lambda(z)=\ln \frac{p\left(z \mid H_{1}\right)}{p\left(z \mid H_{0}\right)}=z-\frac{1}{2}\tag{4} lnλ(z)=lnp(z∣H0)p(z∣H1)=z−21(4)
则判决规则变成:若信号电平 z z z 大于 ( ln λ + 1 2 ) \left(\ln \lambda+\frac{1}{2}\right) (lnλ+21) ,判决有信号输入;若 z z z 小于 ( ln λ + 1 2 ) \left(\ln \lambda+\frac{1}{2}\right) (lnλ+21) 时,判决无信号输入。 接下来求 λ \lambda λ ,设 ( ln λ + 1 2 ) = P F = γ \left(\ln \lambda+\frac{1}{2}\right) =P_F=\gamma (lnλ+21)=PF=γ ,由于 P F = 0.1 P_{F}=0.1 PF=0.1 ,即
∫ γ ∞ p ( z ∣ H 0 ) d z = ∫ γ ∞ 1 2 π exp ( − z 2 2 ) d z = P F = 0.1 (5) \int_{\gamma}^{\infty} p\left(z \mid H_{0}\right) d z=\\ \int_{\gamma}^{\infty} \frac{1}{\sqrt{2 \pi}} \exp \left(-\frac{z^{2}}{2}\right) d z=P_F=0.1\tag{5} ∫γ∞p(z∣H0)dz=∫γ∞2π1exp(−2z2)dz=PF=0.1(5)
用 n o r m i n v norminv norminv函数即可解得 γ \gamma γ的值。
norminv((1-0.1),0,sqrt(1))
解得
γ = 1.2816 \gamma = 1.2816 γ=1.2816
利用 e r f c erfc erfc函数也可以解 γ \gamma γ的值。需要将(5)式变为 e r f c erfc erfc函数的格式。 e r f c erfc erfc函数为
e r f c ( x ) = ∫ x ∞ 2 π exp ( − t 2 ) d t (6) erfc(x)=\int_{x}^{\infty}\frac{2}{\sqrt{\pi}} \exp \left(-t^2\right) dt\tag{6} erfc(x)=∫x∞π2exp(−t2)dt(6)
令
t
=
z
2
t=\displaystyle\frac{z}{\sqrt{2}}
t=2z
则(5)式变为
∫ γ ∞ p ( z ∣ H 0 ) d z = 1 2 ∫ γ 2 ∞ 2 π exp ( − t 2 ) d t = 1 2 e r f c ( γ 2 ) = 0.1 (6) \int_{\gamma}^{\infty} p\left(z \mid H_{0}\right) d z=\\\frac{1}{2}\int_{\frac{\gamma}{\sqrt{2}}}^{\infty} \frac{2}{\sqrt{ \pi}} \exp \left(-t^2\right) d t=\\\frac{1}{2}erfc(\frac{\gamma}{\sqrt{2}})=0.1\tag{6} ∫γ∞p(z∣H0)dz=21∫2γ∞π2exp(−t2)dt=21erfc(2γ)=0.1(6)
由(6)可得
γ = 2 e r f c i n v ( 0.2 ) = 1.2816 (7) \gamma=\sqrt{2} erfcinv(0.2)=1.2816\tag{7} γ=2erfcinv(0.2)=1.2816(7)
此时,检测概率为
P D = P ( D 1 ∣ H 1 ) = ∫ γ ∞ p ( z ∣ H 1 ) d z = ∫ γ ∞ 1 2 π exp [ − ( z − 1 ) 2 2 ] d z (8) P_{D}=P\left(D_{1} \mid H_{1}\right)=\int_{\gamma}^{\infty} p\left(z \mid H_{1}\right) d z=\\\int_{\gamma}^{\infty} \frac{1}{\sqrt{2 \pi}} \exp \left[-\frac{(z-1)^{2}}{2}\right] d z\tag{8} PD=P(D1∣H1)=∫γ∞p(z∣H1)dz=∫γ∞2π1exp[−2(z−1)2]dz(8)
用 n o r m c d f normcdf normcdf函数可以算出 P D P_D PD。
1-normcdf(1.2816,1,1)
算得检测概率
P D = 0.3891 P_D=0.3891 PD=0.3891
画图
1 P F P_F PF和门限之间的关系图
在噪声服从
N
(
0
,
1
)
N(0,1)
N(0,1)的情况下,将
P
F
P_F
PF作为自变量,门限
γ
\gamma
γ作为变量作图:
显然,随着 P F P_F PF增大,门限会下降。
2 P F P_F PF和 P D P_D PD之间的关系图
在噪声服从 N ( 0 , 1 ) N(0,1) N(0,1)的情况下,将 P F P_F PF作为自变量,检测概率 P D P_D PD作为变量作图:
3 仿真 P F P_F PF值
通过理论计算,得到门限,用门限对仿真数据进行检测,计算虚警的个数,进而计算虚警概率,和理论值作比较,画图。
代码如下:
%在N(0,1)噪声,PF=0.1的情况下,用蒙特卡洛仿真计算PF_sim
loopTime = 60;
simLen = 20000000;
%构造信号
x = randi([0 1],1,simLen);
sigma = 1;
s = randn(1,simLen)*sigma;
z = x + s;
pf = linspace(0,1,loopTime);
PF_sim = zeros(1,loopTime);
for i=1:loopTime
%计算门限
gamma = norminv((1-pf(i)),0,1);
%根据门限做虚警概率的计算
%检测
y = z>gamma;
%虚警的索引
pf_vector = y-x>0.5;
%虚警总数
pf_number = sum(pf_vector);
%虚警概率
pf_sim(i)=pf_number/sum(x==1);
end
figure(4)
plot(pf,'o--');
hold on
plot(pf_sim,'^--');
hold off
legend('pf','pf_sim')
可以看到,理论和仿真结果是一致的。