Bayes理论与多传感器数据融合

Bayes理论与多传感器数据融合

1. 背景介绍

多传感器数据融合是一种处理多源异构信源信息的方法，而Bayes理论是一种概率推理方法。

为了更好地讨论多传感器数据融合方法在具体问题中的应用，我们这里引入“单信源二元信号统计检测问题”作为问题场景，目的是更好地阐述多传感器数据融合技术是如何解决传统算法在面对多个信源时的挑战的。

0x1：单信源二元信号统计检测的基本模型

所谓单信源，是指只存在一个信源，即特征维度为1，只有1个特征。

所谓二元信号，是指信源信号只存在两种状态（二元状态），即该特征维度只有2个可选状态，是一种最简单的情况。

而统计信号处理的理论基础之一是信号的统计检测。信号的统计检测理论主要研究在受噪声干扰的随机信号中信号的有/无问题、或信号属于哪个状态的最佳判决的概念、方法和性能等问题，其数学基础就是统计检测理论，又称假设检验理论。

统计检测理论的基本模型是二元信号检测问题：

二元信号统计检测理论模型

• 信源：检测的信号发生源
• 概率转移机构：它是在信源输出其中一个为真的基础上，把噪声干扰背景中的假设H_j(j=0,1)为真的信号以一定概率关系映射到观测空间。

• 如果没有噪声干扰，信源输出的某一种确知信号将映射到观测空间中的某一确定点
• 但在噪声干扰的情况下，它将以一定的概率映射到整个观测空间

• 观测空间：观测空间 R 是在信源输出不同信号状态下，在噪声干扰背景中，由概率转移机构所生成的全部可能观测量的集合
• 判决规则：观测量落入观测空间后，就可以用来推断哪一个假设成立是合理的，即判决信号属于哪种状态。为此需要建立一个判决规则，以便使观测空间中的每一个点对应着一个相应的假设H_i(i＝0,1)。判决结果就是选择假设 H₀ 成立，还是选择假设 H₁ 成立。例如最大后验估计就是一种判决规则。所谓判决规则，本质上是对观测空间R的划分问题。

1、信源

2、概率转移机构

3、观测空间

在二元信号检测中，是把整个观测空间 R 划分为 R₀ 和 R₁ 两个子空间（判决域），并满足：

如果观测空间 R 中的某个观测量落入 R₀ 域，就判决 H₀ 成立，否则就判决假设 H₁ 成立。

4、判决规则 

用于判决的准则可以是多种多样的，但信号检测的基本工具却大多利用了数理统计的假设检测方法。

在目标检测中的决策任务实际上属于“二元假设检验问题”。所谓的假设检验问题是指给定一组假设：

通过对已有的数据集 y 进行处理，来确定当前哪一个假设 H_j 是成立的。

如何根据已有的观测数据 z 得到最优决策结果，属于最优决策策略问题。

在目标检测中的二元假设检验，例如在监视区域里检测目标的存在与否（传感器网、雷达等），如果假设 H₀ 表示无目标，H₁ 表示有目标。则 P₀ 和 P₁ 分别为 H₀ 和 H₁ 出现的先验概率。

单个传感器的决策值 u 为二元值，一般定义为： 

单个传感器的观测结果为 z：

二元假设检验问题有两种可能的错误： 

• H₀ 为真、判决 u=1 为第一类错误：虚警
• H1 为真、判决 u=0 为第二类错误：漏检

对于二元假设检验问题，数据集 z 空间可划分为两个区域：

• R₀：如果数据点 z 位于 R₀ 区，则认为假设 H₀ 成立而做出决策 u=0；
• R₁：如果数据点 z 位于 R₁ 区，则认为假设 H₁ 成立而做出决策 u=1；

笔者插入：

这里使用假设检验的词语是按照概率统计推断的理论框架来说的，从概率论的角度来看，这等同于贝叶斯估计和最大后验概率估计。它们的本质是一样的。

0x2：单信源二元信号检测的数学描述

假设 p₀(z) 和 p₁(z) 分别为 z 关于 H₀ 和 H₁ 的条件概率密度函数，下图中斜线部分红色和天蓝色表示了第一类错误概率 P_f 和第二类错误概率 P_m。

2. 多源检测融合的系统架构

分布式数据融合是每个传感器都先对原始观测数据进行初步分析处理，做出本地判决结论，只把这种本地判决结论及其有关信息，或经初步分析认定可能存在某种结论但又不完全可靠的结论及其有关信息，向数据融合中心呈报。

对于单个传感器来说，它本质上还是一个单信源信号检测模型，可以继续应用上一章节所述的贝叶斯方法进行假设检验。

3. 多源信号融合的准则

0x1：基于似然函数比的融合准则

统计分析中的似然函数（likelihood function)，是指在某种非随机参数条件下观测到数据 z 的概率。

似然函数就是条件概率密度函数p(z∣θ)。

单个传感器检测中，在某种假设下，单传感器的判定结果的条件概率函数为似然函数，即p(z∣H_i)。

在融合中心，多个传感器的判定结果会得到汇总，在多元融合检测中，在某种假设下，各个传感器判定结果输出的概念为似然函数，即 p(u|H_i)，其中 u=(u₁，u₂，...，u_N)。

0x2：最大后验概率融合检测准则

最大后验概率准则是根据已有的观测数据 z，分析给定全局决策之 u 下各种假设的概率，选择最可能产生该全局决策值的假设。

在二元假设目标检测中就是取两个似然概率中较大者所对应的假设，即：

• 若 P(H₁|u) > P(H₀|u)，则选择 H₁
• 反之选择 H₀

应用贝叶斯法则得：

可得似然函数比：

以 N 个传感器为例，下面讨论 N 个传感器的最大后验概率融合准则。

假设 u={u₁，u₂，...，u_N} 中判定为“有目标”的传感器个数为 k 个，即：

似然函数比求对数：

计算似然函数比，并求对数：

其中:

一般记为：

• u₀=1表示有目标
• 否则视为无目标

最大后验概率融合检测准则的缺点：

在最大后验概率融合检测准则中，实现最佳融合检测准则必须知道假设 H₁ 或 H₀ 的先验概率 P₁ 或 P₀=1-P₁，以及各个局部检测器的虚警概率 P_fⁱ 和漏检概率 P_mⁱ。 

0x3：贝叶斯融合检测基本准则

1、准则形式化描述

在最大后验概率融合检测准则中，虚警和漏检两类错误都没有特殊加权，这样相当于假定两类错误是同等危险的。

在工程应用中，各类错误的后果并非同等严重，即不同类型的错误所造成的损失或者说要付出的代价是不同的。为了反映这些差别，应对各类错误概率分别规定不同的代价，即所谓代价函数来反映这种损失的不同。

贝叶斯融合检测准则对每一个决策结果分配相应的代价值，基于假设概率得到平均总代价，检测策略是使平均总代价最小。

以单传感器二元假设检验问题为例定义代价： 

C_ij 为假设 Hi 成立时做出 u_j 决策的代价，平均总代价为：

结合前面二元信号检测的概率密度曲线图，我们还有：

其中，p₀(z)为没有目标的量测分布概率密度函数，p₁(z)为有目标的量测分布概率密度函数。

所以平均总代价为：

我们的目标是选择区域 R₁ 最小，即使平均代价最小，即保证上式积分项小于零。  

由此得到二元假设检测的贝叶斯检测准则：

对比最大后验概率准则： 

显然，按贝叶斯准则与按最大后验概率准则得到的检测系统只是门限不同，其公式形式都是类似的。

另外，注意表达式中概率密度函数与概率的区别。而当代价的选取满足下式时：

最大后验概率准则就是贝叶斯判决准则的特例。 

2、贝叶斯融合检测准则讨论

• C₁₀ - C₀₀ 大：即虚警引起的损失大，门限应取大一些使虚警出现的可能性小一点；反之亦然。
• C₀₁ - C₁₁ 大：即漏检引起的损失大，门限应取小一些使漏检出现的可能性小一点；反之亦然。

在概率论理论体系中，贝叶斯融合检测准则是多传感器系统优化决策的主流技术，是迄今为止理论上最完整的信息融合方法。

在各种先验概率及各种错误决策代价已知的情况下，贝叶斯方法是最优的方法，但是如何获得所需先验概率知识及各种错误决策的代价是应用该方法的一个关键问题。

0x4：最小误差概率准则

在有一些应用场合，对两类错误没有什么特殊的区别，那么令所有误差的代价函数最小也是一个合理的准则，即令下式成立：

那么代价函数式为： 

这里Pe为误差概率，使 C 最小的策略就是使误差概率 Pe 最小。为使 Pe 最小，则应使上式积分项尽可能小，于是有：

等价于：

因此最小误差概率准则与最大后验概率准则完全相同。准则为：

4. 基于贝叶斯融合理论的多源信号融合算法

0x1：多源信号贝叶斯融合公式

前面说到，在概率论理论体系中，贝叶斯融合检测准则是多传感器系统优化决策的主流技术。

这个章节我们来讨论一下如何基于贝叶斯融合理论，实现多元信号融合算法。

以分布式并行结构为例，对于N个传感器融合检测贝叶斯风险：

其中：

上面贝叶斯风险是为整个系统的贝叶斯风险，首先由各个传感器本身作出判决后将判决结果送入融合中心，融合中心的判决结果由 u₀ 表示，每个传感器的判决结果统一由 u=(u₁，u₂，....，u_N) 表示。

融合系统的发现概率和虚警概率的计算公式如下：

由上述公式合并可推得融合系统贝叶斯风险为： 

由融合系统贝叶斯风险公式可知，为使贝叶斯风险达到最小的最优分布检测必要条件为： 

由以上必要条件可进一步推得： 

或者：

通过上式可以看出，要想运用该理论到工程实践中，就必须知道每个传感器的似然函数。而一般在传感器技术文档中，是不会给出似然函数概率的，通常以概率密度的方式给出。

下面就推导将上式转化为概率密度表达的方式。 

设每个传感器有条件概率密度：p_i(z_i | H₁)，p_i(z_i | H₀)，以及定义：

在解决了似然函数用条件概率密度替换的问题后，以上计算还需要先验概率密度知识。

在多传感器环境下，各个传感器的先验概率密度在计算上是非常困难的。为此人们重点研究各个传感器满足检测概率大于虚警概率，且彼此统计独立的条件下的贝叶斯融合检测准则。

根据最优融合规则的单调性，设每个传感器满足：

进一步地，在各个传感器统计独立的条件下的贝叶斯融合检测准则可以简化为：

上式中各个条件概率为： 

求解最优贝叶斯融合检测准则，需要在给定 C_F 和 C_D 的条件下求出：

由于最优融合检测准则与分布的各个传感器的贝叶斯判决门限是耦合的，所以在求解最优贝叶斯融合检测准则的同时，需要求解贝叶斯融合检测准则最佳判决门限： 

接下来讨论求解最优贝叶斯融合检测准则和各个传感器判决门限的迭代算法。

0x2：求解贝叶斯融合检测准则最佳判决门限的迭代算法

1、第一步 

初始化条件。设各个传感器的观测分布(已知)：

以及假设各个传感器的贝叶斯判决门限：

2、第二步

设置循环变量 n，和终止控制量 ξ，准备循环。估算贝叶斯融合检测准则：

计算对应的贝叶斯融合风险R_B⁽¹⁾：

3、第三步

以上一次计算的各个传感器的贝叶斯判决门限为基础计算新的判决门限：

 根据各个传感器新的门限计算对应传感器的

，以计算的各个传感器的贝叶斯判决门限为基础：

。

估算贝叶斯融合检测准则：

4、第四步

计算 k+1 此次迭代的贝叶斯风险 R_B(K+1)：

5. 仿真实验

0x1：实验背景

clc;
j = 100;
V1 = 2.2;
V2 = 2.5;
V3 = 2.8;
delta1 = 1;
delta2 = 1.5;
delta3 = 1.2;

%初始化三个传感器的门限值
gama = zeros(j,3);

%代价初始化
C00 = 0;
C11 = 0;
C01 = 1;
C10 = 1;


%传感器1概率初始化
pf1 = zeros(j,1);
pm1 = zeros(j,1);
pd1 = zeros(j,1);
p001 = zeros(j,1);

%传感器2概率初始化
pf2 = zeros(j,1);
pm2 = zeros(j,1);
pd2 = zeros(j,1);
p002 = zeros(j,1);

%传感器3概率初始化
pf3 = zeros(j,1);
pm3 = zeros(j,1);
pd3 = zeros(j,1);
p003 = zeros(j,1);

Rb = zeros(j,1);

%先验概率初始化
P0=1/j:(1/j):1;
P1 = 1 - P0;

C = C01*P1 + C00*P0;
CF = P0*(C10-C00);
CD = P1*(C01-C11);
t = (CF./CD);

%初始判决门限，因为初始情况下，门限公式后面的无法计算

gama(1,1) = t(1);
gama(1,2) = t(1);
gama(1,3) = t(1);


for i = 1:j    
%根据门限计算各个传感器的四类概率
p001(i) = normcdf(gama(i,1),0,delta1);
pm1(i) = normcdf(gama(i,1),V1,delta1);
pd1(i) =1 -  pm1(i);
pf1(i) = 1 - p001(i);


p002(i) = normcdf(gama(i,2),0,delta2);
pm2(i) = normcdf(gama(i,2),V2,delta2);
pd2(i) =1 -  pm2(i);
pf2(i) = 1 - p002(i);

p003(i) = normcdf(gama(i,3),0,delta3);
pm3(i) = normcdf(gama(i,3),V3,delta3);
pd3(i) =1 -  pm3(i);
pf3(i) = 1 - p003(i);

%根据传感器的检测概率和虚假概率，计算融合中心融合准则
Pu000_H1 = pm1(i)*pm2(i)*pm3(i);
Pu000_H0 = p001(i)*p002(i)*p003(i);

Pu001_H1 = pm1(i)*pm2(i)*pd3(i);
Pu001_H0 = p001(i)*p002(i)*pf3(i);

Pu010_H1 = pm1(i)*pd2(i)*pm3(i);
Pu010_H0 = p001(i)*pf2(i)*p003(i);

Pu011_H1 = pm1(i)*pd2(i)*pd3(i);
Pu011_H0 = p001(i)*pf2(i)*pf3(i);

Pu100_H1 = pd1(i)*pm2(i)*pm3(i);
Pu100_H0 = pf1(i)*p002(i)*p003(i);

Pu101_H1 = pd1(i)*pm2(i)*pd3(i);
Pu101_H0 = pf1(i)*p002(i)*pf3(i);

Pu110_H1 = pd1(i)*pd2(i)*pm3(i);
Pu110_H0 = pf1(i)*pf2(i)*p003(i);

Pu111_H1 = pd1(i)*pd2(i)*pd3(i);
Pu111_H0 = pf1(i)*pf2(i)*pf3(i);

Pu0_1_u000 = (Pu000_H1/Pu000_H0) > t(i);
Pu0_1_u001 = (Pu001_H1/Pu001_H0) > t(i);
Pu0_1_u010 = (Pu010_H1/Pu010_H0) > t(i);
Pu0_1_u011 = (Pu011_H1/Pu011_H0) > t(i);
Pu0_1_u100 = (Pu100_H1/Pu100_H0) > t(i);
Pu0_1_u101 = (Pu101_H1/Pu101_H0) > t(i);
Pu0_1_u110 = (Pu110_H1/Pu110_H0) > t(i);
Pu0_1_u111 = (Pu111_H1/Pu111_H0) > t(i);

%计算融合系统的贝叶斯风险

Rb(i) = C(i) + Pu0_1_u000*(CF(i)*Pu000_H0 - CD(i)*Pu000_H1) + ...
    Pu0_1_u001*(CF(i)*Pu001_H0 - CD(i)*Pu001_H1)+...
    Pu0_1_u010*(CF(i)*Pu010_H0 - CD(i)*Pu010_H1)+...
    Pu0_1_u011*(CF(i)*Pu011_H0 - CD(i)*Pu011_H1)+...
    Pu0_1_u100*(CF(i)*Pu100_H0 - CD(i)*Pu100_H1)+...
    Pu0_1_u101*(CF(i)*Pu101_H0 - CD(i)*Pu101_H1)+...
    Pu0_1_u110*(CF(i)*Pu110_H0 - CD(i)*Pu110_H1)+...
    Pu0_1_u111*(CF(i)*Pu111_H0 - CD(i)*Pu111_H1);

%计算A(u)，确定某一个传感器的值后组合两外两个传感器的值
%确定u1
Au1_00 = Pu0_1_u100 - Pu0_1_u000;
Au1_01 = Pu0_1_u101 - Pu0_1_u001;
Au1_10 = Pu0_1_u110 - Pu0_1_u010;
Au1_11 = Pu0_1_u111 - Pu0_1_u011;

%确定u2
Au2_00 = Pu0_1_u010 - Pu0_1_u000;
Au2_01 = Pu0_1_u011 - Pu0_1_u001;
Au2_10 = Pu0_1_u110 - Pu0_1_u100;
Au2_11 = Pu0_1_u111 - Pu0_1_u101;

%确定u3
Au3_00 = Pu0_1_u001 - Pu0_1_u000;
Au3_01 = Pu0_1_u011 - Pu0_1_u010;
Au3_10 = Pu0_1_u101 - Pu0_1_u100;
Au3_11 = Pu0_1_u111 - Pu0_1_u110;

%计算后验概率
Pu1_00_H0 = p002(i)*p003(i);
Pu1_01_H0 = p002(i)*pf3(i);
Pu1_10_H0 = pf2(i)*p003(i);
Pu1_11_H0 = pf2(i)*pf3(i);

Pu2_00_H0 = p001(i)*p003(i);
Pu2_01_H0 = p001(i)*pf3(i);
Pu2_10_H0 = pf1(i)*p003(i);
Pu2_11_H0 = pf1(i)*pf3(i);

Pu3_00_H0 = p001(i)*p002(i);
Pu3_01_H0 = p001(i)*pf2(i);
Pu3_10_H0 = pf1(i)*p002(i);
Pu3_11_H0 = pf1(i)*pf2(i);

Pu1_00_H1 = pm2(i)*pm3(i);
Pu1_01_H1 = pm2(i)*pd3(i);
Pu1_10_H1 = pd2(i)*pm3(i);
Pu1_11_H1 = pm2(i)*pm3(i);

Pu2_00_H1 = pm1(i)*pm3(i);
Pu2_01_H1 = pm1(i)*pd3(i);
Pu2_10_H1 = pd1(i)*pm3(i);
Pu2_11_H1 = pm1(i)*pm3(i);

Pu3_00_H1 = pm1(i)*pm2(i);
Pu3_01_H1 = pm1(i)*pd2(i);
Pu3_10_H1 = pd1(i)*pm2(i);
Pu3_11_H1 = pm1(i)*pm2(i);

%计算三种传感器的融合判决门限

T1 =(CF(i)*(Au1_00*Pu1_00_H0 + Au1_01*Pu1_01_H0 + Au1_10*Pu1_10_H0 + Au1_11*Pu1_11_H0))/...
    (CD(i)*(Au1_00*Pu1_00_H1 + Au1_01*Pu1_01_H1 + Au1_10*Pu1_10_H1 + Au1_11*Pu1_11_H1));

T2 =(CF(i)*(Au2_00*Pu2_00_H0 + Au2_01*Pu2_01_H0 + Au2_10*Pu2_10_H0 + Au2_11*Pu2_11_H0))/...
    (CD(i)*(Au2_00*Pu2_00_H1 + Au2_01*Pu2_01_H1 + Au2_10*Pu2_10_H1 + Au2_11*Pu2_11_H1));

T3 =(CF(i)*(Au3_00*Pu3_00_H0 + Au3_01*Pu3_01_H0 + Au3_10*Pu3_10_H0 + Au3_11*Pu3_11_H0))/...
    (CD(i)*(Au3_00*Pu3_00_H1 + Au3_01*Pu3_01_H1 + Au3_10*Pu3_10_H1 + Au3_11*Pu3_11_H1));

gama(i+1,1) = log(T1)/V1 + V1/2;
gama(i+1,2) = log(T2)/V2 + V2/2;
gama(i+1,3) = log(T3)/V3 + V3/2;

end
%传感器1贝叶斯风险
R1 = (C00*p001+C10*pf1).*P0' + (C01*pm1+C11*pd1).*P1';
%传感器2贝叶斯风险
R2 = (C00*p002+C10*pf2).*P0' + (C01*pm2+C11*pd2).*P1';
%传感器3贝叶斯风险
R3 = (C00*p003+C10*pf3).*P0' + (C01*pm3+C11*pd3).*P1';

plot(P0,R1,'b',P0,R2,'g',P0,R3,'black',P0,Rb,'r'); 
title('贝叶斯风险曲线图');
xlabel('p0');
ylabel('贝叶斯风险值');
legend('传感器1贝叶斯风险曲线','传感器2贝叶斯风险曲线','传感器3贝叶斯风险曲线','融合系统贝叶斯风险曲线');