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₁ 出現的先驗概率。

，且 P0 + P1 = 1

單個傳感器的決策值 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貝葉斯風險曲線','融合系統貝葉斯風險曲線');
 

Relevant Link:      

https://blog.csdn.net/haxiongha/article/details/82682358