雙水平集算法及其在函數輪廓估計中的應用

李曉杉，夏界寧

(1. 中國地震局地震研究所地震預警湖北省重點實驗室，湖北 武漢 430071；2. 湖北省重大工程地震安全監測與預警處置工程技術研究中心，湖北 武漢 430071)

1 引言

水平集估計(Level Set Estimation， LSE)問題，即在某個空間上尋找信號大于某個門限值的點集。在信號采集過程中，只有非常有限的信號樣本點可用，且這些樣本含有噪聲，這就使得無法逐個“遍歷”每個樣本點，而必須使用有限且含有噪聲的樣本進行分析[1-5]。在進行數據采樣及分析時，如果觀測樣本有限或者觀測難度較大時，需要關注樣本點的統計特性，高斯過程具有十分良好的特性，因此，利用高斯過程進行目標函數的建模[6]。LSE算法本質上與二元分類算法相似，是一種單調分類的算法[1]。利用LSE算法，能夠在觀測值有限的情況下有效地進行數據分類與預測。

水平集估計問題在函數最優化[1]、遙感圖像處理[7]、醫學圖像處理[8]、數據聚類分析[9]等領域得到了廣泛應用。

文獻[1]中B. Bryan等人對小樣本情況下的函數閾值邊界問題進行了分析，利用高斯過程對目標函數進行建模，對比了不同方法的效果，并進行了數值實驗，將參數估計問題所需的樣本數和計算量減少了一個數量級[1]。文獻[2]中A. Gotovos等人建立了LSE算法框架，并進行了數據實驗，利用LSE算法對實際問題進行了處理，實現了小樣本空間下的數據分析，同時對LSE算法進行了擴展，針對隱含式閾值邊界及批量選擇樣本情況進行了延伸，并且進行了收斂性理論分析[2]。文獻[3]中Y. Ma等人使用高斯過程對函數建模，并使用貝葉斯積分來推斷其在感興趣區域上的平均值，針對有效區域搜索(Active Area Search， AAS)問題進行了分析，比以前提出的相關方法更快地識別出正區域[3]。文獻[4]中I. Bogunovic等人利用高斯過程處理貝葉斯優化和水平集估計問題，提出了截斷方差減少(Truncated Variance Reduction， TRUVAR)算法，可以納入近視算法。與執行超前的其它貝葉斯優化(Bayesian Optimization，BO)算法相比，TRUVAR算法避免了計算上復雜的對后驗和測量進行平均的步驟，并且具有嚴格的理論保證[4]。文獻[5]中Andrea Zanette等人提出了水平集估計的魯棒最大改進(Robust Maximum Improvement for Level-set Estimation， RMILE)算法，對LSE算法進行了改進，在一步先行程序中選擇最大化高于閾值點域的預期體積的下一個查詢點，并導出分析公式以封閉形式計算該數量，能夠更有效地使用從少量樣本中獲得的信息，使其適用于非常有限的樣本點，且保證了漸近收斂性，使其在特定模型的情況下得到了有效結果[5]。

與其它分類算法相比，D＿LSE算法的查詢點少、受噪聲影響小，在解決有限且含噪聲數據的分類問題時具有優勢。與傳統的LSE算法相比，D＿LSE算法考慮了兩個閾值，更適用于小樣本空間下的雙水平集估計問題。這一研究具有一定工程應用價值，例如在重力測量方面，絕對重力的測量成本高、難度大，重力數據有限且含有噪聲。相比于獲得每個測量點的具體重力值，某個區域的重力分部更加重要，希望將重力值劃分為多個等級，此時單一閾值無法滿足研究需求，需考慮多閾值分類的算法。相比于之前的研究，考慮了分類精度參數ε對算法的影響。本文首先給出了D＿LSE算法的理論依據，其次進行數值實驗，最后利用D＿LSE算法解決函數最優化問題，展示D＿LSE算法的實用性。

2 D＿LSE算法

當數據點(觀測值)的獲取難度比較大或成本昂貴時，無法對數據進行逐個遍歷，必須對數據的統計特征進行分析。D＿LSE算法利用高斯過程對函數建模，能夠有效描述數據的統計特征[10]。通過對高斯過程的樣本進行分類，分為上水平集、中水平集和次水平集，并進行不斷更新，能夠實現函數的估計(分類)。

2.1 D＿LSE算法概述

就A. Gotovos的原LSE算法做出擴展，原LSE算法考慮單一閾值問題，應用單一水平集進行區域的分類，把區域分為超水平集(Super-Level Set)與次水平集兩部分。將此問題進行擴展，考慮多閾值問題，應用多個水平集，將觀測值分為多個區域。

D＿LSE算法的原理是構建一個三維曲面，將二維問題轉化為三維曲面的水平集估計問題。以高斯過程逼近采樣點與曲面的映射關系，不斷對置信邊界進行引導采樣和分類。

假定估計函數f：D→；采樣點x∈D；兩個閾值h1∈Rd，h2∈Rd，且h1

(1)

雙水平集的示意圖如圖1。將這種方法擴展到多水平集的情況下，可用來處理多閾值問題。

圖1 雙水平集示意圖

D＿LSE算法主要有初始化、計算置信區間、分類、更新、選取查詢點五個步驟。

2)計算置信區間：對于未分類的樣本點，用當前的高斯過程計算樣本點的置信區間。

假定每次只分類單個數據點，則對于任意t≥ 2，數據點xt∈D，噪聲測量值yt=f(xt)+nt，其中nt為零均值高斯白噪聲，方差為σ2。若GP的先驗g(0，k(x，x′)) 超過f；并給出At={x1，…，xt}中t點的噪聲測量值yt=[y1， …，yt]T，其中yi=f(xi)+ni，且ni～N(0，σ2)，i=1， …，t。根據高斯過程的性質，后驗仍是高斯過程，因此，能夠通過如下公式對所擬合的高斯分布進行更新

μt(x)=kt(x)T·(Kt+σ2I)-1·yt

(2)

kt(x，x′)=k(x，x′)-kt(x)T·(Kt+σ2I)-1·kt(x)

(3)

(4)

其中kt(x)=[k(x1，x)， …，k(xt，x)]T，且Kt是即將觀測點的核矩陣

Kt=[k(x，x′)]；x，x′∈At

(5)

由式(2)、(4)可推斷出置信區間如下

(6)

對置信區間取交集，Ct(x)表示每個x的單調遞減的置信區域

(7)

3)分類：通過置信區間與上閾值、下閾值的關系對該點進行分類。定義集合Ut表示未分類的點集；Ht表示超水平集，包含高于上閾值h2水平的點；Mt表示中水平集，包含高于下閾值h1且低于上閾值h2水平的點；Lt表示次水平集，包含低于下閾值h1水平的點。

4)更新：依據式(2)、(4)、(7)更新未分類點的高斯過程與置信區間，并進行分類。

5)選取查詢點：若存在未分類的點，則需要依據最大模糊度準則選擇新的查詢點來獲取信息，直至所有點均完成分類。

D＿LSE算法為一種單調分類的算法，在分類錯誤的情況下無法進行修改。因此以置信區間引入模糊度：

(8)

模糊度的示意圖如圖2。模糊度越大的點越接近閾值，包含更多信息，因此以模糊度最大準則選取查詢點。

圖2 模糊度示意圖

2.2 算法流程

算法1：D＿LSE算法

輸入： 采樣集D，高斯先驗參數(μ0=0，k，δ0)，上閾值h2，下閾值h1，精度參數ε。

1)H0←?，M0←?，L0←?，U0←?，

C0(X)←，forallx∈D

2)t←1

3)whileUt-1≠?do

4)Ht←Ht-1，Mt←Mt-1，Lt←Lt-1，

Ut←Ut-1

5)forallx∈Ut-1do

6)Ct(x)←Ct-1(x)∩Qt(x)

7)ifmin(Ct(x))+ε>h2then

8)Ut←Ut{x}，Ht←Ht∪{x}

9)elseifmax(C0(x))-ε≤h2and

min(Ct(x))+ε>h1then

10)Ut←Ut{x}，Mt←Mt∪{x}

11)elseifmax(C0(x))-ε≤h1then

12)Ut←Ut{x}，Lt←Lt∪{x}

13)xt←arg maxx∈Ut(at(x))

14)yt←f(x)+nt

15) 對所有x∈Ut， 計算μt(x)andσt(x)

16)t←t+1

2.3 精度與復雜度分析

互信息I(X：Y)可用來表示X中包含Y中信息的多少[10]。在D＿LSE算法中，考慮t個觀測值yt=(yi)1≤i≤t和估計函數f的互信息如下

I(yt：f)=H(yt)-H(yt|f)

(9)

其中H(yt)為觀測值的熵，H(yt|f)為條件熵。

為了從觀測值yt中獲得更多估計函數f的信息，最大化互信息I(yt：f)

(10)

高斯分布的熵定義如下

(11)

假定信號服從高斯分布，且噪聲為高斯白噪聲，則由式(13)、(15)可得

(12)

為了衡量分類的質量，引入錯誤分類損失

(13)

lh(x)的最大值反映了分類的質量，max(lh(x))越小則分類質量越高。對于準確度為ε的分類，若max(lh(x))≤ε，則可以認為分類在ε精度的條件下是完全正確的，即分類是ε精確的。

考慮D＿LSE算法的迭代次數T來分析算法的復雜度。對于閾值h∈，精度參數ε>0，α∈(0，1)，D＿LSE算法的最大迭代次數T滿足如下不等式

(14)

此時，算法有1-α的概率達到ε分類精度

(15)

3 仿真與結果分析

首先完成D＿LSE算法的MATLAB程序驗證，接著給出其F1分數評價，最后修改實驗參數進行對比，并進行分類結果質量評估。

3.1 目標函數的建模

參考文獻[1]中的函數模型，其目標閾值邊界具不連續、存在模糊區域、存在函數遠離閾值區域、長度足夠的特點，因此能夠避免算法在在這些區域過采樣[1]。函數的三維模型如圖4所示，考慮的函數模型如下

f(x，y)=sin(10x)+cos(4y)-cos(3xy)

(16)

相關參數設置如下：

3.2 D＿LSE仿真

3.3 算法分類質量評估

對于雙水平集問題，可以看作是兩層二元分類問題。若僅僅使用分類準確率這一指標評估分類質量，則在出現數據不平衡問題時效果很差，例如有100個樣本點，其中99個為正，1個為負，而直接預測所有點都為正，這樣雖然準確率達到99%，但預測是不具有高質量的。對于二分類問題，F1＿score方法能夠很好地評估其分類質量[11-12]。對于閾值h2，假定超水平集H為正，中水平集M為負；對于閾值h1，假定中水平集M為正，次水平集L為負。精確率與召回率的定義見表1和表2。

表1 閾值h2分類結果評估方法表

表2 閾值h1分類結果評估方法表

3.4 實驗結果

3.4.1 對比實驗

1)實驗參數1：h1=-0.5，h2=0.3，采樣點25*50=1250，精度參數ε=0.3，實驗結果如圖3。

圖3 參數1實驗結果

查詢點數t=132；運行時間405.65秒。

2)實驗參數2：h1=-0.5，h2=0.3，采樣點25*50=1250，精度參數ε=0.1，實驗結果如圖4。

圖4 參數2實驗結果

查詢點數t=532； 運行時間9611.15秒。

3)實驗參數3：h1=-0.5，h2=0.3，采樣點10*20=200，精度參數ε=0.3，實驗結果如圖5。

圖5 參數3實驗結果

查詢點數t=91； 運行時間39.51秒。

3.4.2 函數最優化問題實驗

不斷地增大和減小閾值，能夠不斷接近實際函數的最大值和最小值。相比于傳統求函數最值的方法，此方法并沒有對函數凹凸性的限制，即可得到全局最優結果。

設定參數為：采樣點25*50=1250，精度參數ε=0.1，實驗結果如圖6。

圖6 函數最優化問題實驗結果

3.5 實驗結果分析

D＿LSE水平集估計結果表明其能夠有效估計函數的真實輪廓。最終分類結果的F1得分接近于1，表明分類質量較高。

通過改變相關參數：采樣點數、精度參數進行對比。精度參數ε越小，最終F1得分越高，表明分類質量越高，但查詢點數越多，程序運行時間越長；采樣點數越多，分類質量越高，但查詢點數越多，程序運行時間越長。不斷增大/減小閾值，能夠有效逼近函數的全局最大/最小值，解決函數全局最優化問題。

4 結束語

本文提出了小樣本空間下的雙水平集算法，主要結論有：1)考慮了兩個不同的閾值，以高斯過程描述信號的統計特性，采用構造三維隱曲面的方法完成二維曲線的估計，提出了一種雙水平集估計算法；2)將D＿LSE算法應用于函數輪廓估計問題中，且以數值實驗驗證了算法的有效性；3)討論了采樣點數和分類精度參數ε的對算法的影響，并進行了對比實驗，采樣點數越多、精度參數ε越小則算法質量和復雜度越高；4)將D＿LSE算法應用于函數最優化問題，證明了算法在最優化問題上的可行性。