付青青, 楊立華, 付建棟, 吳愛平

(1.長江大學(xué)電信學(xué)院, 湖北 荊州 434023; 2.中國石油集團(tuán)測井有限公司, 陜西 西安 710201)

0 引 言

超聲成像測井不僅可以在裸眼井中反映井眼幾何形狀,識別裂縫、孔洞、層理等地層非均質(zhì)性,還能在套管井中檢查射孔質(zhì)量、分析套管損壞等。測井過程中,由于超聲換能器的形狀和工作頻率、井孔的大小和形狀、井壁的傾角和表面結(jié)構(gòu)、井液泥漿密度等因素會影響超聲波的衰減[1-2],因此,一幅超聲成像測井圖像的灰度分布通常會集中在某一個很窄的動態(tài)范圍內(nèi),使得圖像對比度低,導(dǎo)致圖像中的細(xì)節(jié)常?？床磺宄?給目標(biāo)地質(zhì)體特征的識別和探測帶來困難。在成像測井圖像處理中,常采用直方圖修正的方法增強(qiáng)圖像質(zhì)量[3-6]。胡剛等[3]采用全局直方圖均衡對測井圖像進(jìn)行處理[3],雖增強(qiáng)了圖像的整體信息,但對局部細(xì)節(jié)信息增強(qiáng)能力不明顯。閆建平等[4]在利用圖像直方圖局部均衡化實現(xiàn)測井圖像動態(tài)增強(qiáng)研究的基礎(chǔ)上,采用Morphing技術(shù)對顏色進(jìn)行線性插值,消除動態(tài)增強(qiáng)產(chǎn)生的臺階問題,效果明顯。涂繼輝等[5]利用某井段一個局部區(qū)域數(shù)據(jù)統(tǒng)計其直方圖,從而實現(xiàn)局部圖像的最佳均衡,該方法對局部細(xì)節(jié)的增強(qiáng)效果優(yōu)于全局直方圖均衡算法,但這種動態(tài)直方圖增強(qiáng)方法會引入平坦區(qū)域少量噪聲放大和過度增強(qiáng)等問題。

本文針對超聲測井圖像對比度低、細(xì)節(jié)模糊的問題,借鑒一種限定對比度直方圖均衡方法對超聲測井圖像進(jìn)行增強(qiáng),達(dá)到了突出局部細(xì)節(jié),限制平坦區(qū)域過度增強(qiáng)的目的,對進(jìn)一步提高成像測井圖像解釋精度具有重要意義。

1 直方圖均衡化

直方圖均衡又稱為全局直方圖均衡,它是最基本的直方圖均衡方法[7],其基本思想是通過對原圖像進(jìn)行某種灰度變換,使變換后圖像的直方圖能均勻分布,這樣就能使原圖像中具有相近灰度且占有大量像素點的區(qū)域灰度范圍展寬,使大區(qū)域中的微小變化顯現(xiàn)出來,使圖像更清晰。對于連續(xù)圖像,設(shè)r代表原圖像的灰度級,假定r歸一化為0≤r≤1,r=0代表黑(最暗),r=1代表白(最亮),p(r)為原始圖像灰度分布的概率分布函數(shù),直方圖均衡處理實際上就是尋找一個灰度變換函數(shù)T,使變換后的灰度值滿足s=T(r),其中,s歸一化為0≤s≤1,要求處理后圖像灰度分布的概率密度函數(shù)p(s)=1,變換函數(shù)T(r)必須滿足2個條件[8]:

(1) 在0≤r≤1區(qū)間內(nèi)是單值且單調(diào)遞增函數(shù),確保灰度變換,原始圖像的每個灰度級r都對應(yīng)產(chǎn)生一個灰度級s,且灰度變換前后不倒置(使灰度級保持從黑到白的次序)。

(2) 在0≤s≤1區(qū)間內(nèi)0≤T(r)≤1,確保輸出灰度級與輸入灰度級有同樣的范圍。

設(shè)原圖像的灰度范圍為[r,r+dr],包含的像素個數(shù)為p(r)dr,經(jīng)過單調(diào)增的一對一變換,變換后灰度范圍為[s,s+ds],包含像素個數(shù)為p(s)ds,變換前后的像素個數(shù)應(yīng)相等,即

p(r)dr=p(s)ds

(1)

(2)

式(2)右邊為直方圖p(r)的累積分布函數(shù),也是直方圖均衡化對應(yīng)的灰度變換函數(shù)。

對于數(shù)字圖像離散情況,直方圖均衡化的計算過程:

(1) 根據(jù)原始圖像的統(tǒng)計信息,算出原始直方圖,p(i)=ni/n。其中,i=0,1,…,L-1;L為灰度級的個數(shù);ni表示各灰度級的像素數(shù),n表示像素的總數(shù)。

(3) 借助累積直方圖,映射各像素點的灰度值,運(yùn)用變換得到公式(3),其中的INT()是取整符號。

j=INT[(L-1)pj+0.5]

(3)

(4) 確定新的映射i→j,據(jù)此,將原圖像的灰度值f(m,n)=i修正為g(m,n)=j。

(5) 統(tǒng)計變換后各灰度級的像素個數(shù)nj,并由此給出處理后圖像的直方圖p(j)=nj/n。

超聲測井圖像進(jìn)行直方圖均衡化后,能夠擴(kuò)展圖像的動態(tài)范圍,增加圖像的整體對比度,質(zhì)量得到一定程度的提升,使感興趣的特征更容易被檢測或識別。如果超聲測井圖像中表示某些特征的像素點較少,當(dāng)小到一定程度后,在式(3)中運(yùn)用累積直方圖映射各像素點的灰度值時,不同的pj有可能會得到相同的j,這樣會造成灰度值合并現(xiàn)象,導(dǎo)致測井圖像中細(xì)節(jié)丟失。另外,對于長達(dá)幾千m測井?dāng)?shù)據(jù)形成的測井圖像,不同的深度有不同的特征,使用全局直方圖均衡算法,不可避免地對少數(shù)噪聲點也進(jìn)行了放大,如果存在局部區(qū)域過亮的部分,會有過度增強(qiáng)的問題[9-10],反而會使這部分圖像包含的細(xì)節(jié)更加不清晰。因此,超聲測井圖像的增強(qiáng),需要采用突出局部細(xì)節(jié),限制平坦區(qū)過度增強(qiáng)的方法。

2 限定對比度直方圖均衡化方法

限定對比度的自適應(yīng)直方圖均衡化算法由Karel Zuiderveld提出[11],基本思想是通過限制局部直方圖的高度,從而可以限制噪聲放大和削弱局部對比度的過增強(qiáng),已在許多領(lǐng)域得到了廣泛應(yīng)用[10,12-13]。針對彩色的測井圖像,本文將RGB圖像變換到HSI(Hue Saturation Intensity)空間,僅對亮度分量進(jìn)行限定對比度直方圖均衡化處理,再將結(jié)果作RGB輸出,這樣保留了圖像的顏色信息,處理前后圖像的整體色調(diào)保持不變。

限定對比度直方圖均衡化的基本流程如下[9]:

(1) 將輸入原圖像劃分為若干大小相等的子塊,每個子塊相互連續(xù)并且互不重疊。

(2) 統(tǒng)計各子塊內(nèi)的局部直方圖。

(3) 計算子塊內(nèi)所有像素數(shù)平均分配到的各灰度級值,并用參數(shù)ANum表示,若用GNum表示子塊可能的灰度級,則以上過程可表示成

(4)

式中,XP、YP分別表示子塊X、Y方向的像素數(shù)。

(4) 設(shè)定一個剪切系數(shù)CV,并假設(shè)其具有[0,1]的取值范圍。針對不同的圖像,通過仿真結(jié)果調(diào)整到最佳值,則實際剪切極限值NV可用公式(5)表示,其中round是四舍五入函數(shù)。

NV=ANum+round[CV(XP×YP-ANum)]

(5)

(5) 利用剪切極限,對局部直方圖各個灰度級上的像素進(jìn)行剪切,多余下來的像素數(shù)目重新分配到各直方圖的各灰度級中,設(shè)已被剪切的像素總數(shù)為NClip,于是得到每個灰度級應(yīng)該分到的剪切像素數(shù)NAcp。

NClip=∑{max[H(i)-NV,0]}

(6)

式中,H(i),i=i0,i1,…,in表示子塊的局部直方圖。

(7)

如果用CH表示重新分配后的直方圖,則有

(8)

(6) 設(shè)經(jīng)過以上分配后,像素點剩余數(shù)目為NLeft,則分配像素的步長為

(9)

從最小灰度級開始按步長搜索,遇到像素數(shù)小于剪切閾值的位置,則分配1個像素,完成從最小灰度級到最大灰度級的循環(huán),直到NLeft為0,直方圖分配才算完成,最終得到新的直方圖。

(7) 對剪切后的每個子區(qū)域的灰度直方圖分別進(jìn)行直方圖均衡化。

(8) 把每個子塊的中心點作為參考點,獲取其灰度值,對圖像中的每一個像素進(jìn)行灰度線性插值,采用雙線性插值的方法,每個像素點的映射由其相鄰的4個參考點對應(yīng)區(qū)域的映射確定,如圖1中黑色小矩形(x,y)代表任意目標(biāo)點,f(x,y)表示該點待計算的灰度值,與其相鄰的4個區(qū)域的中心點用黑色圓點表示,分別記為A(x1,y1)、B(x2,y2)、C(x3,y3)、D(x4,y4),其中x1=x3,x2=x4,y1=y2,y3=y4,則目標(biāo)點f(x,y)的灰度值可由A、B、C、D灰度值的線性組合表示成公式(10),在邊界區(qū)域的像素,其灰度計算用鄰近2個樣本點進(jìn)行線性插值;在圖像4個角的點用鄰近1個樣本點進(jìn)行計算。

f(x,y)=a[bf(x1,y1)+(1-b)f(x2,y2)]+

(1-a)[bf(x3,y3)+(1-b)f(x4,y4)]

(10)

圖1 雙線性插值示意圖

3 實驗結(jié)果及分析

3.1 主觀評價

圖2 長慶宜×井圖像增強(qiáng)效果圖

圖3 長慶宜×井圖像的直方圖

采用限定對比度的自適應(yīng)直方圖均衡化方法,使用版本為R2008b的Matlab軟件對超聲測井圖像的增強(qiáng)效果進(jìn)行驗證。圖2(a)是一幅長慶宜×井超聲測井原始彩色幅度圖的一段截圖;圖2(b)是使用全局直方圖均衡算法處理的結(jié)果;圖2(c)是使用限定對比度直方圖均衡處理的結(jié)果。根據(jù)仿真結(jié)果,調(diào)整參數(shù)至最佳,其中分割的子塊數(shù)目為8×8個,剪切系數(shù)CV為0.01。從圖2中可以看出,2種方法都一定程度上增強(qiáng)了圖像的對比度,但是圖2(b)中對應(yīng)原始圖像高亮區(qū)域明顯過度發(fā)亮,淹沒了細(xì)節(jié)信息,造成圖像局部區(qū)域過度增強(qiáng),有一定程度失真。圖2(c)中對原始圖像中高亮局部區(qū)域的細(xì)節(jié)增強(qiáng)效果明顯好于圖2(b)中的效果,沒有出現(xiàn)過度增強(qiáng)的問題。圖3分別給出了圖2(a)、(b)、(c)對應(yīng)的灰度直方圖信息,橫坐標(biāo)為像素灰度級,縱坐標(biāo)為像素在某灰度級出現(xiàn)的個數(shù)。從圖3(a)可以看出原始圖像像素集中分布在灰度級50～240的范圍,對比度偏低。圖3(b)是經(jīng)過全局直方圖均衡算法處理后的直方圖,其像素近似分布在0～255之間,同時在低灰度級和高灰度級像素比較集中,而中間部分比較稀疏,整體不均勻。從直方圖還可以看出100左右兩邊附近的有些灰度級的像素個數(shù)為0,這是由于較少像素的灰度被合并了,淹沒了圖像的細(xì)節(jié)信息,同時在255灰度級附近像素個數(shù)為500左右,直方圖有尖狀化現(xiàn)象,這是導(dǎo)致亮區(qū)過度增強(qiáng)的原因。圖3(c)經(jīng)過限定對比度直方圖均衡后,其直方圖的動態(tài)范圍被均勻擴(kuò)展到10～255灰度范圍內(nèi),沒有出現(xiàn)細(xì)節(jié)退化和過度增強(qiáng)的問題。

3.2 客觀評價

(11)

RPSN為最大信號量和噪聲強(qiáng)度的比值,它的表達(dá)式為

(12)

對于常見的256級灰度圖像,fmax=255。

(13)

把信息論中熵值的概率應(yīng)用到圖像信息源,假如圖像的灰度級為[1,L],可以通過直方圖得到各灰度級的概率ps(sk),k=1,2,…,L,這時,圖像的熵為

(14)

圖像處理后均方誤差(EMS)越小、峰值信噪比(RPSN)越大說明處理效果好。信息熵越大,信息的無序度越高,包含的信息量越大。

表1給出了圖2(a)的超聲測井圖像經(jīng)過全局直方圖均衡和限定對比度直方圖均衡算法各參數(shù)相應(yīng)的結(jié)果。從表1中可以看出,全局直方圖均衡算法處理過的圖像,其均方誤差大于限定對比度直方圖均衡算法處理過的圖像,前者的峰值信噪比和信息熵均小于后者,表明采用限定對比度直方圖均衡化算法具有較好的圖像增強(qiáng)效果。

表1 超聲測井圖像的客觀評價結(jié)果

4 結(jié) 論

(1) 超聲測井圖像灰度直方圖通常會集中在很窄的動態(tài)范圍內(nèi),導(dǎo)致圖像對比度低、細(xì)節(jié)模糊現(xiàn)象。采用限定對比度直方圖均衡方法對超聲測井圖像進(jìn)行增強(qiáng),達(dá)到了突出局部細(xì)節(jié),限制平坦區(qū)域過度增強(qiáng)的目的。

(2) 對圖像處理效果及直方圖的形態(tài)進(jìn)行分析,通過客觀評價參數(shù)均方誤差、峰值信噪比和信息熵對比2種直方圖增強(qiáng)方法,仿真結(jié)果顯示,限定對比度直方圖均衡算法對于超聲測井圖像質(zhì)量的改善比全局直方圖均衡算法效果更好,對實際的超聲測井資料評價具有一定的意義,也為測井圖像處理提供了又一種可行的方法。

參考文獻(xiàn):

[1] 李長文, 強(qiáng)毓明, 余春昊, 等. 提高超聲成像測井圖像質(zhì)量的途徑 [J]. 石油儀器, 2009, 23(5): 51-54.

[2] 付青青, 謝凱. 應(yīng)用維納濾波技術(shù)改善超聲波成像測井圖像效果 [J]. 測井技術(shù), 2012, 36(6): 581-584.

[3] 胡剛, 陳凡. 基于直方圖均衡化的成像測井彩色圖像增強(qiáng) [J]. 國外測井技術(shù), 2011, 4: 68-70.

[4] 閆建平, 首祥云, 邵在平, 等. 成像測井圖像的動態(tài)增強(qiáng)及Morphing方法 [J]. 測井技術(shù), 2006, 30(4): 364-366.

[5] 涂繼輝, 余厚全, 李長文, 等. 基于超聲測井圖像的動態(tài)直方圖均衡算法研究 [J]. 電視技術(shù), 2011, 35(19): 113-114.

[6] 王曉峰, 彭天慈, 雷剛, 等. XRMI動態(tài)增強(qiáng)及全井眼成像方法研究與應(yīng)用 [J]. 測井技術(shù), 2015, 39(4): 432-437.

[7] 曾煒赫, 楊俊煒, 邢宇翔. 動態(tài)直方圖雙向均衡化的圖像增強(qiáng)方法 [J]. 中國體視學(xué)與圖像分析, 2014, 19(2): 129-139.

[8] GONZAKEZ R C, WOODS R E. DigitalImage Processing [M]. 阮秋琦, 阮宇智, 譯. 2版. 北京: 電子工業(yè)出版社, 2007: 70-74.

[9] 趙靜. 彩色圖像增強(qiáng) [D]. 西安: 西安電子科技大學(xué), 2014.

[10] 張璞. 基于視頻交通車輛檢測系統(tǒng)的圖像增強(qiáng)的研究 [D]. 青島: 青島大學(xué), 2012.

[11] ZUIDERVELD K. Contrast Limited Adaptive Histogram Equalization [M]. Chapter VIII. 5, Graphics gems IV, San Diego: Academic Press Professional, Inc, 1994: 474-485.

[12] ALAMEEN Z, SULONG G, REHMAN A, et al. An Innovative Technique for Contrast Enhancement of Computed Tomography Images Using Normalized Gamma-corrected Contrast-limited Adaptive Histogram Equalization [J]. EURASIP Journal on Advances in Signal Processing, 2015, 2015(1): 1-12.

[13] 楊有, 李波. CLAHE和細(xì)節(jié)放大相結(jié)合的檔案圖像增強(qiáng)方法 [J]. 中國圖象圖形學(xué)報, 2011, 16(4): 522-527.

[14] 張德豐. Matlab數(shù)字圖像處理 [M]. 北京: 機(jī)械工業(yè)出版社, 2012: 213-218.