雙水平集算法及其在函數(shù)輪廓估計(jì)中的應(yīng)用

李曉杉，夏界寧

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

1 引言

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

水平集估計(jì)問題在函數(shù)最優(yōu)化[1]、遙感圖像處理[7]、醫(yī)學(xué)圖像處理[8]、數(shù)據(jù)聚類分析[9]等領(lǐng)域得到了廣泛應(yīng)用。

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

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

2 D＿LSE算法

當(dāng)數(shù)據(jù)點(diǎn)(觀測值)的獲取難度比較大或成本昂貴時(shí)，無法對數(shù)據(jù)進(jìn)行逐個(gè)遍歷，必須對數(shù)據(jù)的統(tǒng)計(jì)特征進(jìn)行分析。D＿LSE算法利用高斯過程對函數(shù)建模，能夠有效描述數(shù)據(jù)的統(tǒng)計(jì)特征[10]。通過對高斯過程的樣本進(jìn)行分類，分為上水平集、中水平集和次水平集，并進(jìn)行不斷更新，能夠?qū)崿F(xiàn)函數(shù)的估計(jì)(分類)。

2.1 D＿LSE算法概述

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

D＿LSE算法的原理是構(gòu)建一個(gè)三維曲面，將二維問題轉(zhuǎn)化為三維曲面的水平集估計(jì)問題。以高斯過程逼近采樣點(diǎn)與曲面的映射關(guān)系，不斷對置信邊界進(jìn)行引導(dǎo)采樣和分類。

假定估計(jì)函數(shù)f：D→；采樣點(diǎn)x∈D；兩個(gè)閾值h1∈Rd，h2∈Rd，且h1

(1)

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

圖1 雙水平集示意圖

D＿LSE算法主要有初始化、計(jì)算置信區(qū)間、分類、更新、選取查詢點(diǎn)五個(gè)步驟。

2)計(jì)算置信區(qū)間：對于未分類的樣本點(diǎn)，用當(dāng)前的高斯過程計(jì)算樣本點(diǎn)的置信區(qū)間。

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

μ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是即將觀測點(diǎn)的核矩陣

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

(5)

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

(6)

對置信區(qū)間取交集，Ct(x)表示每個(gè)x的單調(diào)遞減的置信區(qū)域

(7)

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

4)更新：依據(jù)式(2)、(4)、(7)更新未分類點(diǎn)的高斯過程與置信區(qū)間，并進(jìn)行分類。

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

D＿LSE算法為一種單調(diào)分類的算法，在分類錯(cuò)誤的情況下無法進(jìn)行修改。因此以置信區(qū)間引入模糊度：

(8)

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

圖2 模糊度示意圖

2.2 算法流程

算法1：D＿LSE算法

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

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， 計(jì)算μt(x)andσt(x)

16)t←t+1

2.3 精度與復(fù)雜度分析

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

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

(9)

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

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

(10)

高斯分布的熵定義如下

(11)

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

(12)

為了衡量分類的質(zhì)量，引入錯(cuò)誤分類損失

(13)

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

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

(14)

此時(shí)，算法有1-α的概率達(dá)到ε分類精度

(15)

3 仿真與結(jié)果分析

首先完成D＿LSE算法的MATLAB程序驗(yàn)證，接著給出其F1分?jǐn)?shù)評價(jià)，最后修改實(shí)驗(yàn)參數(shù)進(jìn)行對比，并進(jìn)行分類結(jié)果質(zhì)量評估。

3.1 目標(biāo)函數(shù)的建模

參考文獻(xiàn)[1]中的函數(shù)模型，其目標(biāo)閾值邊界具不連續(xù)、存在模糊區(qū)域、存在函數(shù)遠(yuǎn)離閾值區(qū)域、長度足夠的特點(diǎn)，因此能夠避免算法在在這些區(qū)域過采樣[1]。函數(shù)的三維模型如圖4所示，考慮的函數(shù)模型如下

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

(16)

相關(guān)參數(shù)設(shè)置如下：

3.2 D＿LSE仿真

3.3 算法分類質(zhì)量評估

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

表1 閾值h2分類結(jié)果評估方法表

表2 閾值h1分類結(jié)果評估方法表

3.4 實(shí)驗(yàn)結(jié)果

3.4.1 對比實(shí)驗(yàn)

1)實(shí)驗(yàn)參數(shù)1：h1=-0.5，h2=0.3，采樣點(diǎn)25*50=1250，精度參數(shù)ε=0.3，實(shí)驗(yàn)結(jié)果如圖3。

圖3 參數(shù)1實(shí)驗(yàn)結(jié)果

查詢點(diǎn)數(shù)t=132；運(yùn)行時(shí)間405.65秒。

2)實(shí)驗(yàn)參數(shù)2：h1=-0.5，h2=0.3，采樣點(diǎn)25*50=1250，精度參數(shù)ε=0.1，實(shí)驗(yàn)結(jié)果如圖4。

圖4 參數(shù)2實(shí)驗(yàn)結(jié)果

查詢點(diǎn)數(shù)t=532； 運(yùn)行時(shí)間9611.15秒。

3)實(shí)驗(yàn)參數(shù)3：h1=-0.5，h2=0.3，采樣點(diǎn)10*20=200，精度參數(shù)ε=0.3，實(shí)驗(yàn)結(jié)果如圖5。

圖5 參數(shù)3實(shí)驗(yàn)結(jié)果

查詢點(diǎn)數(shù)t=91； 運(yùn)行時(shí)間39.51秒。

3.4.2 函數(shù)最優(yōu)化問題實(shí)驗(yàn)

不斷地增大和減小閾值，能夠不斷接近實(shí)際函數(shù)的最大值和最小值。相比于傳統(tǒng)求函數(shù)最值的方法，此方法并沒有對函數(shù)凹凸性的限制，即可得到全局最優(yōu)結(jié)果。

設(shè)定參數(shù)為：采樣點(diǎn)25*50=1250，精度參數(shù)ε=0.1，實(shí)驗(yàn)結(jié)果如圖6。

圖6 函數(shù)最優(yōu)化問題實(shí)驗(yàn)結(jié)果

3.5 實(shí)驗(yàn)結(jié)果分析

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

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

4 結(jié)束語

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