基于非凸包的灌注區(qū)包絡(luò)自動(dòng)擬合方法

謝家靖 ，滕奇志 ，何小海 ，龔 劍

(1.四川大學(xué) 電子信息學(xué)院 圖像信息研究所，四川 成都 610065；2.成都西圖科技有限公司，四川 成都 610065)

0 引言

微觀驅(qū)替實(shí)驗(yàn)[1-2]是指利用真實(shí)巖心或者微觀仿真模型進(jìn)行驅(qū)替實(shí)驗(yàn)，通過顯微鏡觀察微觀狀態(tài)下油水在孔隙中的滲流情況，從而分析水驅(qū)油、聚合物驅(qū)油等驅(qū)替方式的驅(qū)替效果以及剩余油的分布情況。由于良好的可視性和可重復(fù)利用性，石油地質(zhì)工作者往往采用玻璃刻蝕模型來(lái)代替真實(shí)巖心進(jìn)行仿真實(shí)驗(yàn)，其實(shí)驗(yàn)過程如圖1 所示[3-4]。通過玻璃刻蝕模型可以在驅(qū)替過程中觀察到油、水、驅(qū)替劑在孔隙中的運(yùn)移過程，再通過連接攝像裝置，采集模型的動(dòng)態(tài)圖像，利用圖像處理技術(shù)對(duì)采集到的圖像進(jìn)行定量分析，從而研究不同驅(qū)替劑的驅(qū)替效率和機(jī)理[5]。

圖1 微觀驅(qū)替仿真實(shí)驗(yàn)過程

為分析驅(qū)替劑在模型中的波及程度，需要擬合灌注區(qū)包絡(luò)。灌注區(qū)顧名思義是驅(qū)替劑灌注的區(qū)域。通過包絡(luò)曲線擬合灌注入的驅(qū)替劑所波及到的范圍，由此計(jì)算出驅(qū)替劑在該時(shí)刻下的波及系數(shù)，從而為評(píng)價(jià)驅(qū)替劑的驅(qū)替效率提供重要依據(jù)。由于灌注區(qū)形態(tài)復(fù)雜且不規(guī)則，難以得到準(zhǔn)確的波及系數(shù)，目前大多系統(tǒng)中提取灌注區(qū)包絡(luò)是通過手動(dòng)擬合的方式。但手動(dòng)擬合方法步驟繁瑣，特別是對(duì)動(dòng)態(tài)序列圖進(jìn)行處理時(shí)，要進(jìn)行多次的重復(fù)操作，這給使用者帶來(lái)較大的負(fù)擔(dān)，不適合微觀驅(qū)替序列圖的批量處理。為實(shí)現(xiàn)提取灌注區(qū)包絡(luò)的自動(dòng)化，同時(shí)保證灌注區(qū)擬合的精確度，本文提出了一種基于非凸包的灌注區(qū)包絡(luò)自動(dòng)擬合方法。

1 灌注區(qū)圖像介紹

以水驅(qū)油為例，圖 2(a)、(b)分別為玻璃模型中充滿油時(shí)的飽和油圖像和驅(qū)替實(shí)驗(yàn)過程中水驅(qū)14 min39 s 的圖像，通過將驅(qū)替過程中的圖像與飽和油圖像進(jìn)行比較，可以很明顯地觀察到驅(qū)替劑灌注進(jìn)的區(qū)域。

圖2 微觀驅(qū)替實(shí)驗(yàn)圖像

為了提取出灌注區(qū)的包絡(luò)，首先需要對(duì)圖片進(jìn)行預(yù)處理，目的是為了將玻璃模型中充滿原油的孔隙部分作為目標(biāo)與仿真模型背景區(qū)分出來(lái)。目標(biāo)區(qū)域的顏色單一，且顏色與背景顏色相區(qū)分度明顯。利用這個(gè)特點(diǎn)，本文采用最大方差法，對(duì)圖像進(jìn)行分割。

以飽和油時(shí)刻的分割結(jié)果作為基準(zhǔn)圖，并將水驅(qū)過程中序列圖的分割結(jié)果與基準(zhǔn)圖進(jìn)行比較，就可以得到驅(qū)替過程中各個(gè)時(shí)刻的灌注區(qū)。以水驅(qū)14 min 39 s 時(shí)刻的圖像為例，圖3(c)紅色的區(qū)域即為注入的水所填充的區(qū)域，也就是灌注區(qū)。本文的目的是自動(dòng)計(jì)算上述預(yù)處理得到的灌注區(qū)的包絡(luò)，也就是用比較平滑的曲線將灌注的區(qū)域包住形成一個(gè)封閉的形狀，并保證包絡(luò)的精確性，從而能夠計(jì)算得到較為準(zhǔn)確的波及系數(shù)。

2 非凸包擬合算法

凸包算法通常用來(lái)計(jì)算有限點(diǎn)集的閉合包絡(luò)，已被廣泛應(yīng)用于圖像處理、數(shù)據(jù)分析等領(lǐng)域。所謂凸包(Convex Hull)，如圖 4(a)、圖 4(b)所示，即對(duì)于平面上離散的點(diǎn)集Q，尋找包含該點(diǎn)集的面積最小的凸多邊形P，凸包也就是構(gòu)成 P 的所有頂點(diǎn)，本文用 CH(Q)來(lái)表示 Q 的凸包。常用的凸包算法有Graham 掃 描 算 法 、Jarvis 步 進(jìn) 法 等[6-7]。

對(duì)于較為規(guī)則的點(diǎn)集，如圖4(a)所示，使用凸包算法能夠較好地?cái)M合點(diǎn)集的形狀；但對(duì)于類似圖4(b)的不規(guī)則點(diǎn)集，凸包算法并不能準(zhǔn)確描述點(diǎn)集的包絡(luò)。在微觀驅(qū)替實(shí)驗(yàn)中，灌注區(qū)通常是不規(guī)則的，采用傳統(tǒng)的凸包算法去擬合灌注區(qū)，將引入過多的無(wú)效區(qū)域，為后續(xù)的波及系數(shù)計(jì)算帶來(lái)明顯的誤差。對(duì)于這種形態(tài)的點(diǎn)集，用非凸包(Non-convex Hull)或者凹包(Concave Hull)來(lái)描述更加準(zhǔn)確[8-9]，如圖 4(c)、圖4(d)所示，非凸包能夠更好擬合不規(guī)則點(diǎn)集的邊界。本文采用 NCH(Q)來(lái)表示 Q 的非凸包。

圖3 目標(biāo)分割圖像

圖4 凸包和非凸包

本文采用漸進(jìn)迭代的思想，在凸包算法的基礎(chǔ)上挖掘非凸包點(diǎn)，最終形成非凸包。目前對(duì)于非凸包數(shù)學(xué)上并沒有明確的定義。對(duì)于同一個(gè)點(diǎn)集，采用不同的挖掘深度D，獲取的非凸包形態(tài)也不相同[10]。挖掘深度顧名思義， 即挖掘非凸包點(diǎn)的程度，本文定義為在凸包點(diǎn)之間尋找非凸包點(diǎn)的次數(shù)。如圖 4 所示，圖(c)和圖(d)采用同樣的非凸包算法，但采用不同的挖掘深度， 得到的包絡(luò)也有所不同。若 D 值過小，就無(wú)法精確地描述點(diǎn)集 Q 的形狀；D 值過大，則非凸包所構(gòu)成的多邊形邊界就不夠平滑。因此，對(duì)于不同的點(diǎn)集，需要采用合適的挖掘深度D。綜上所述，本文對(duì)傳統(tǒng)的凸包算法進(jìn)行改進(jìn)，算法具體如下：

(1)使用傳統(tǒng)凸包算法提取點(diǎn)集Q 的凸包CH(Q)，構(gòu)成凸多邊形 S={s0，s1，s2…sN-1}，其中 si是凸包點(diǎn)，也就是凸多邊形的頂點(diǎn)；

(2)設(shè)置挖掘深度 D，邊長(zhǎng)閾值 L；

(3)在每?jī)蓚€(gè)相鄰?fù)拱c(diǎn)si和si+1之間進(jìn)行挖掘，以 si和 si+1構(gòu)成的邊長(zhǎng)為搜索直徑，計(jì)算待挖掘點(diǎn)與 si和 si+1構(gòu)成的夾角 θ，并選取 θ 最大的點(diǎn)為非凸包點(diǎn)，存入NCH(Q)中；

(4)如果挖掘次數(shù)大于 D，或者 si和 si+1的距離小于L，則停止挖掘，尋找下兩個(gè)相鄰?fù)拱c(diǎn)之間的非凸包點(diǎn)，重復(fù)步驟(3)；

(5)最終生成的包絡(luò)點(diǎn)集可表示為：

以 s0和 s1為例，以 s0和 s1構(gòu)成的邊長(zhǎng)為搜索直徑，計(jì)算待挖掘點(diǎn)與 s0和 s1所構(gòu)成的夾角，并選取夾角最大的點(diǎn)為挖掘點(diǎn)，如圖 5(b)所示，選取 p 作為非凸包點(diǎn)，并存入 NCH(Q)。通過不斷迭代，逐步縮小范圍，直到滿足停止挖掘的條件，再繼續(xù)搜索下兩個(gè)凸包點(diǎn)之間的疑似非凸包點(diǎn)。通過這樣的方法，可以明顯去除凸包中絕大多數(shù)的無(wú)效區(qū)域。

圖5 非凸包算法

3 灌注區(qū)的自動(dòng)擬合方法

不同于上述算法中離散的點(diǎn)集，由于微觀驅(qū)替圖像中的灌注區(qū)屬于塊狀點(diǎn)集，直接使用該點(diǎn)集計(jì)算包絡(luò)會(huì)帶來(lái)巨大的計(jì)算量，本文在擬合包絡(luò)前先提取點(diǎn)集的邊緣點(diǎn)，剔除點(diǎn)集的內(nèi)部點(diǎn)，通過減少點(diǎn)數(shù)量大幅提升計(jì)算速度。

另外由于孔隙的非均質(zhì)性，灌注區(qū)的分布存在不確定性，圖像中可能存在多個(gè)灌注區(qū)。以圖3(c)為例，可以明顯看到有3 個(gè)間隔較大的紅色區(qū)域。為更加準(zhǔn)確地分析灌注區(qū)的特征，進(jìn)一步排除無(wú)效區(qū)域，需要在擬合包絡(luò)前先對(duì)點(diǎn)集進(jìn)行聚類計(jì)算。通常傳統(tǒng)的聚類算法需要人工設(shè)定參數(shù)，不合理的參數(shù)會(huì)導(dǎo)致分類結(jié)果不符合預(yù)期。本文采用迭代的思想，使用一種自適應(yīng)參數(shù)的 DBSCAN 算法[11]，尋找較為合適的鄰域(Eps)和領(lǐng)域內(nèi)最少點(diǎn)集數(shù)(MinPts)。

綜上所述，本文提出了一種基于非凸包的灌注區(qū)包絡(luò)自動(dòng)擬合方法，具體步驟如下：

(1)為減少數(shù)據(jù)計(jì)算量，提取灌注區(qū)的邊界點(diǎn)，得到點(diǎn)集Q；

(2)設(shè)定 Eps 和 MinPts 的初始值，使用 DBSCAN算法對(duì)點(diǎn)集Q 進(jìn)行聚類，獲取分類數(shù)，不斷增加Eps，繼續(xù)進(jìn)行聚類，如果連續(xù)三次進(jìn)行迭代得到的分類數(shù)不變，則選取當(dāng)前的Eps 為最優(yōu)參數(shù)；

(3)使用最優(yōu)參數(shù)輸入 DBSCAN 進(jìn)行聚類，得到灌注區(qū)點(diǎn)集 Qi(0≤i

(4)采用上節(jié)中的非凸包擬合算法對(duì)各個(gè)灌注區(qū)點(diǎn)集計(jì)算，得到包絡(luò)點(diǎn)集 B(Qi)，如圖 6(a)所示；

(5)為模擬水前沿線的真實(shí)特征，使用曲線擬合來(lái)平滑B(Qi)所構(gòu)成的多邊形邊界；

(6)計(jì)算各個(gè)灌注區(qū)的面積 Ai，也即步驟(5)得到的曲線所圍成的面積，A 為模型面積，則波及系數(shù)EA可以表示為：

圖6(b)的結(jié)果顯示，該方法能夠較好地?cái)M合灌注區(qū)，排除了絕大部分無(wú)效區(qū)域，所得到的結(jié)果更精確，擬合的水前沿線在直觀上更符合實(shí)際情況。

圖6 灌注區(qū)的自動(dòng)擬合方法

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

在實(shí)際微觀驅(qū)替實(shí)驗(yàn)中，往往需要對(duì)整個(gè)實(shí)驗(yàn)過程中采集到的序列圖進(jìn)行批量處理，因此需要考慮灌注區(qū)包絡(luò)提取速度。影響灌注區(qū)包絡(luò)提取速度的兩個(gè)主要因素為點(diǎn)集數(shù)量以及挖掘深度。對(duì)于前者， 通過提取邊界點(diǎn)， 可以大幅度減少點(diǎn)集數(shù)量，從而提升計(jì)算速度；對(duì)于后者，挖掘深度的增大能夠提升灌注區(qū)包絡(luò)的擬合精度，但同時(shí)也會(huì)帶來(lái)較大計(jì)算量。為選取較為合適的挖掘深度，選取實(shí)驗(yàn)中多張不同時(shí)刻具有代表性的驅(qū)替圖像，分別計(jì)算不同挖掘深度所得到的波及系數(shù)與手動(dòng)擬合之間的誤差值。

如圖 7、圖 8 所示，D 代表挖掘深度，隨著 D 的不斷增大，波及系數(shù)逐漸逼近手動(dòng)擬合得到的波及系數(shù)值。當(dāng) D 小于 20 時(shí)，此時(shí)挖掘深度過小，并沒有完全排除無(wú)效灌注區(qū)；當(dāng)挖掘深度大于 20 時(shí)，自動(dòng)提取得到的波及系數(shù)和平均絕對(duì)誤差不再變化，趨于穩(wěn)定。因此本文選取挖掘深度為20，在盡量減小挖掘次數(shù)的同時(shí)保證提取的精度。

由于灌注區(qū)的波及范圍并沒有絕對(duì)的參考值，直觀判斷具有一定的主觀性，因此本文方法是采用多人手動(dòng)擬合灌注區(qū)所得到波及系數(shù)的平均值作為比較的參考值，在一定程度上弱化了主觀性對(duì)結(jié)果的影響，具有較高的可信度。手動(dòng)擬合方式得到的結(jié)果如圖9(b)所示。為檢驗(yàn)上述方法的準(zhǔn)確性，本文選取多張連續(xù)時(shí)間的具有代表性的驅(qū)替圖像，表 1 中 E1、E2、E3分 別 表 示 采 用 凸 包 算 法 、 本 文 方法以及手動(dòng)擬合的方式擬合灌注區(qū)包絡(luò)所計(jì)算得到的波及系數(shù)。其中圖片編號(hào)為0 的圖片為基準(zhǔn)圖片，即飽和油圖片。

圖7 不同挖掘深度的波及系數(shù)對(duì)比

圖8 不同挖掘深度的誤差分析

圖9 三種方法擬合包絡(luò)結(jié)果對(duì)比

從表1 中的數(shù)據(jù)可知，采用這種方法自動(dòng)提取灌注區(qū)包絡(luò)，計(jì)算得到的波及系數(shù)與手動(dòng)擬合的結(jié)果比較相近，可以準(zhǔn)確地反映驅(qū)替劑所波及的范圍。采用平均絕對(duì)誤差(MAE)對(duì)表1 中的數(shù)據(jù)進(jìn)行誤差分析：

MAEE1，E3=0.268 0

表1 波及系數(shù)計(jì)算

MAEE2，E3=0.016 1

根據(jù)計(jì)算結(jié)果，本文算法相比傳統(tǒng)凸包算法進(jìn)行灌注區(qū)包絡(luò)擬合有更高的精確度。將上述數(shù)據(jù)放入同一折線圖，如圖10 所示，可以看到本文算法與凸包算法相比，在各個(gè)時(shí)刻均有更高的精確度。驅(qū)替實(shí)驗(yàn)的14 min 51 s 之前，本文方法計(jì)算得到的波及系數(shù)與手動(dòng)擬合得到的波及系數(shù)幾乎一致，由于此時(shí)的灌注區(qū)面積相對(duì)較小，目標(biāo)個(gè)數(shù)較少，本文方法能夠準(zhǔn)確地?cái)M合灌注區(qū)包絡(luò)；但是隨著灌注區(qū)面積的增大，目標(biāo)個(gè)數(shù)增多，灌注區(qū)的包絡(luò)更為復(fù)雜，甚至灌注區(qū)內(nèi)部存在非灌注區(qū)，因此沒有完全排除無(wú)效灌注區(qū)域。

圖10 三種方法的波及系數(shù)變化

5 結(jié)論

由于灌注區(qū)形態(tài)復(fù)雜且不規(guī)則，采用凸包算法難以準(zhǔn)確提取灌注區(qū)的包絡(luò)，往往會(huì)引入大面積的無(wú)效區(qū)域。本文方法先對(duì)灌注區(qū)點(diǎn)集進(jìn)行邊緣提取，減少計(jì)算量，再對(duì)灌注區(qū)進(jìn)行聚類分析，最后采用非凸包算法對(duì)每個(gè)灌注區(qū)提取包絡(luò)。另外，采用不同的挖掘深度對(duì)不同時(shí)刻具有代表性的驅(qū)替圖像進(jìn)行灌注區(qū)包絡(luò)計(jì)算，通過誤差分析得到最合適的挖掘深度，在保證精確度的情況下盡可能地提升計(jì)算速度。實(shí)驗(yàn)結(jié)果證明采用本文方法得到的波及系數(shù)相比凸包算法具有更高的精確性，與參照值的相關(guān)程度更高。相比步驟繁瑣的手動(dòng)擬合方式，本文方法在保證結(jié)果精確性的同時(shí)減少了大量的人工操作，具有比較高的應(yīng)用價(jià)值。