低秩稀疏重建分析的邊緣檢測方法

劉 軍 宋維琪 陳俊安 譚 明 胡建林 董 林

(①中國石化西北油田分公司，新疆烏魯木齊 830011；②中國石油大學(xué)(華東)，山東青島 266580)

0 引言

碳酸鹽巖地區(qū)由于溶蝕作用，斷層往往伴有縫洞。研究斷層及其伴生的縫洞，對儲層預(yù)測具有重要意義。邊緣檢測方法眾多，如方差體方法、螞蟻追蹤方法[1]、曲率屬性方法及其他方法[2-4]、多代相干體算法[5-7]、Canny算子邊緣檢測方法、sobel算子邊緣檢測算法等，最新的方法有相似屬性方法、Likelihood屬性方法、混亂屬性相干方法等。其中相似屬性方法、Likelihood屬性方法及多代相干體方法實質(zhì)上也是相干方法。上述方法是非常有效的斷層檢測方法，取得了很好的應(yīng)用效果[8-14]，但也有其自身的不足和邊緣檢測能力的限制，特別是對噪聲干擾、多邊緣干涉及弱小目標(biāo)邊緣的檢測效果不理想。曲率屬性方法、Canny算子和sobel算子邊緣檢測方法涉及梯度計算，嚴(yán)重放大了噪聲，檢測結(jié)果的可靠性差。

壓縮感知理論是近年來新的信號處理理論，在數(shù)據(jù)壓縮、傳輸領(lǐng)域得到廣泛應(yīng)用。壓縮感知理論的稀疏重建理論[15-17]為信號重建提供了新的理論方法，在圖像重建領(lǐng)域普遍獲得較好效果[18-21]；在地震勘探領(lǐng)域，人們利用壓縮感知理論去噪[22]、數(shù)據(jù)重建[23]、拓頻處理，取得了較好效果。

針對邊緣檢測存在的問題，本文首先分析斷層邊緣和縫洞邊緣的空間分布特征，根據(jù)斷層邊緣和縫洞邊緣的地震響應(yīng)特征，把低秩稀疏分析理論引入邊緣檢測，研究邊緣信息、背景信息及噪聲信息的低秩稀疏分解與重建；為了提高邊緣檢測能力和分辨率，在壓縮感知稀疏表示基礎(chǔ)上，對地震資料進(jìn)行深度稀疏化表示，結(jié)合向量稀疏表示和矩陣稀疏表示，通過低秩稀疏分析理論[24-25]，形成一種全新的邊緣檢測方法。

1 基本框架理論

1.1 信號分析低秩表示模型

低秩表示思想是構(gòu)建一個能夠表示所有樣本能力的字典，并且每個樣本都可以用字典中的原子線性表示，同類樣本的表示系數(shù)具有相關(guān)性，故整個表示系數(shù)矩陣是一個低秩矩陣。該模型可以更好地挖掘數(shù)據(jù)的結(jié)構(gòu)信息，較好地解決多線性子空間聚類和分離問題。由于矩陣滿足低秩條件，所以可以保證同類數(shù)據(jù)樣本聚集到相同的子空間。

在復(fù)雜的地震資料邊緣檢測中，相似背景和稀疏目標(biāo)組成了地震響應(yīng)的所有成分，背景又由多種子背景構(gòu)成。背景和局部目標(biāo)也是相對的，若研究目標(biāo)是斷層，則地層相對于斷層為背景；若研究目標(biāo)是斷層中的縫洞，則斷層相對于縫洞為背景。

1.2 低秩稀疏表示

根據(jù)低秩稀疏分解理論，假設(shè)地震信號S由區(qū)域相似信息和局部信息構(gòu)成，則可以分解為低秩矩陣B和稀疏矩陣D兩部分，即：S=B+D，其中B為背景部分(低秩部分或相似部分)，D為高秩部分(稀疏部分)。可利用信息主成分分析求解最優(yōu)低秩矩陣B，即求解下列優(yōu)化問題

minB,D‖S‖F(xiàn)滿足 rank(B)≤r,S=B+D

(1)

式中：r為矩陣的秩；‖S‖F(xiàn)為Frobenius范數(shù)。

通過此約束問題，可以找到矩陣S在一個R維線性子空間的投影。當(dāng)矩陣D的元素服從獨立高斯分布時，對S進(jìn)行奇異值分解，便可得到式(1)的最優(yōu)解。如果D不滿足以上條件，式(1)無解或誤差較大。于是利用穩(wěn)健的主成分分析方法解決這個問題。引入正則化參數(shù)μ，得到新的優(yōu)化問題

min [rank(B)+μ‖D‖F(xiàn)]滿足S=B+D

(2)

上述優(yōu)化問題是一個NP-hard問題，因此對目標(biāo)函數(shù)進(jìn)行松弛，即將非凸函數(shù)問題轉(zhuǎn)化為凸函數(shù)問題，利用核范數(shù)代替矩陣的秩，如L1范數(shù)代替L0范數(shù)。此時的凸優(yōu)化問題為

min[‖B‖1+μ‖D‖1]滿足S=B+D

(3)

利用壓縮感知理論，將式(3)改為

min‖S-A(B+D)‖F(xiàn)滿足‖B‖2≤r,‖D‖0≤k

(4)

式中：A為字典矩陣；k為稀疏目標(biāo)的稀疏度，取決于稀疏目標(biāo)的大小和數(shù)目，而低秩化程度估計的準(zhǔn)確度直接影響背景的抑制效果。由式(4)可知，模型在矩陣的秩和向量的稀疏度約束下求目標(biāo)函數(shù)的最優(yōu)解。通過求解以上最優(yōu)問題，使背景相似信息和局部稀疏信息最佳分離，有效分離地質(zhì)體邊緣信息、地層信息及噪聲，從而進(jìn)行邊緣檢測。

1.3 向量信號稀疏與矩陣低秩稀疏特征分析

信號的稀疏分解可以是一維向量，如一維小波分解，也可以是二維矩陣，如曲波分解。矩陣的秩刻畫了矩陣不相關(guān)行或列的數(shù)目，一個矩陣中的不相關(guān)的行或列越少，則該矩陣的秩越小，反之亦然。矩陣奇異值分解能夠提取高維數(shù)據(jù)的局部特征及約減維數(shù)，其本質(zhì)是將分類能力弱的數(shù)據(jù)集合(矩陣)的類間方差增大，提高分類能力。矩陣的奇異值個數(shù)越少，奇異值分解越稀疏，因此稀疏和低秩是一致的，反映數(shù)據(jù)的壓縮能力。

2 地震資料低秩稀疏邊緣檢測方法

基于壓縮感知原理和假設(shè)，可以將邊緣檢測問題轉(zhuǎn)化為壓縮與矩陣的低秩稀疏分解問題。通過低秩稀疏分解方法壓制地層信息，從而分離邊緣信息。通過數(shù)據(jù)的矩陣結(jié)構(gòu)變換，把地層信息變換為低秩部分，把邊緣信息變換為稀疏部分，利用奇異值分解壓制地層信息進(jìn)而分離邊緣信息。信號重建的可靠性、精度及分辨率都和信號的稀疏化分解程度有關(guān)。由低秩稀疏分解理論[24-25]可知，要得到唯一的分解結(jié)果，要求稀疏部分獨立和高度稀疏。實際地震記錄的邊緣信息、地層信息及噪聲等遠(yuǎn)不能滿足高度稀疏與獨立的條件。

2.1 地震信號向量稀疏表示及去噪優(yōu)化

地震資料稀疏表示是提高邊緣檢測能力的關(guān)鍵環(huán)節(jié)之一。利用壓縮感知理論對地震信號向量稀疏表示，采用平穩(wěn)小波分解方法對地震信號進(jìn)行多尺度小波分解，建立小波稀疏表示字典。有效信號壓縮重建就是稀疏系數(shù)優(yōu)化和稀疏系數(shù)加權(quán)組合的過程。含噪信號通過某種基分解，使信號和噪聲的特征不同而分離。通常情況下，噪聲系數(shù)呈均勻分布，信號系數(shù)呈冪律分布。對于能量較大信號，稀疏分解使有效信號的能量強(qiáng)于噪聲，可根據(jù)壓縮感知理論的優(yōu)化方法壓制噪聲；對于能量較小信號，稀疏分解使有效信號與噪聲能量水平相當(dāng)，此時壓縮感知理論優(yōu)化方法對噪聲的壓制效果不理想。假設(shè)在橫向上信號稀疏分解后有效信號系數(shù)相關(guān)，而噪聲系數(shù)不相關(guān)，則可采用橫向相關(guān)法對稀疏系數(shù)去噪處理。

2.2 矩陣低秩稀疏表示及優(yōu)化方法

以上分析了信號在某種正交基下的稀疏表示方法。為了充分挖掘數(shù)據(jù)結(jié)構(gòu)方面的信息，以下從矩陣秩角度討論稀疏化問題。

2.2.1 張量奇異值分解

矩陣稀疏分解有很多方式，為了更好地分析低秩問題，本文采用奇異值分解稀疏表示方法。數(shù)據(jù)矩陣的秩既與數(shù)據(jù)本身尺度有關(guān)，還與數(shù)據(jù)的排列方式(數(shù)據(jù)結(jié)構(gòu))有關(guān)。地質(zhì)體是三維的，地質(zhì)體地震響應(yīng)信息也是三維的，把三維數(shù)據(jù)分解成二維矩陣進(jìn)行奇異值稀疏分解，涉及張量分析及高階奇異值分解。一般奇異值分解針對2階張量矩陣，對于高階張量，奇異值分解是高階奇異值分解。高階奇異值分解難以分析其低秩稀疏特征，為了研究簡便，這里采用張量分解的方法，即把一個子數(shù)據(jù)體塊分解成一系列二維矩陣，對于不同的分解模型，通過分析每個矩陣的低秩稀疏特征，構(gòu)建具有最優(yōu)低秩稀疏結(jié)構(gòu)信息的張量分解模型。

2.2.2 矩陣低秩優(yōu)化

(1)矩陣張量構(gòu)建與優(yōu)化

對于一個數(shù)據(jù)集合，不同矩陣組建方式有不同的矩陣分解特征，從而有不同的信息表示能力。根據(jù)多線性代數(shù)理論，張量是一個N維數(shù)組，則張量的階數(shù)也為N。向量是一階張量，矩陣是2階張量，立體矩陣是3階張量。數(shù)據(jù)信息包括能量信息、結(jié)構(gòu)信息等。立體矩陣包含了全面的結(jié)構(gòu)信息，為了充分挖掘結(jié)構(gòu)信息，對三維地震數(shù)據(jù)進(jìn)行立體矩陣分析建模。將子塊三維數(shù)據(jù)分解成二維矩陣的模式很多，如果把所有分解模式都列出計算，計算量太大而不適合于實際分析。在斷層和縫洞檢測中，斷層在二維情況下具有線狀特征，洞邊緣在二維情況下具有等軸(近似圓)特性。根據(jù)檢測的兩類邊緣信息，分別構(gòu)建行、列數(shù)不同的矩陣和行、列數(shù)相同的矩陣。對于斷層檢測問題，在縱向構(gòu)建非方陣，即線、道及45°線方向共四個二維矩陣，在橫向構(gòu)建方陣，即過計算點的一個方陣；對于縫洞檢測問題，在縱、橫向均構(gòu)建方陣。

(2)奇異值分解去噪優(yōu)化

矩陣低秩優(yōu)化包括矩陣建模和奇異值分解去噪，以下討論矩陣去噪優(yōu)化。矩陣低秩優(yōu)化就是通過某種方式盡可能獲得具有最小秩的表示矩陣，也就是尋找一個最佳的低秩表示，達(dá)到矩陣低秩稀疏優(yōu)化的目的，從而提高算法的準(zhǔn)確性和邊緣檢測的分辨率。無論采用何種矩陣建模，都不能直接消除數(shù)據(jù)集合中噪聲的影響。矩陣分析理論表明，矩陣中個別元素的突變會使矩陣的數(shù)目和大小發(fā)生明顯變化，這種變化與該元素和矩陣的方差呈線性關(guān)系。矩陣奇異值分解去噪低秩優(yōu)化問題就是使矩陣秩最小化問題。

S=AZ+T+P

(5)

式中：B=AZ是分解字典矩陣，Z是系數(shù)矩陣；T為稀疏目標(biāo)；P為噪聲；λ、β為兩個待定參數(shù)。

將式(5)寫為

S=AZ+T+P

(6)

通過增廣拉格朗日乘子方法，將式(6)的約束優(yōu)化問題轉(zhuǎn)化為一個無約束優(yōu)化問題，優(yōu)化目標(biāo)函數(shù)為

(7)

式中：等式右邊第1項代表奇異值分解；第2項代表壓縮感知稀疏分解重建；第3項代表噪聲，主要為奇異值分解的噪聲，壓縮感知的噪聲已去除；第4項和第5項為信號自身的約束，γ、ρ為懲罰參數(shù)。后三項都是2范數(shù)，以保證計算過程既要壓制噪聲，又要保證新的重建信號是期望信號(去除噪聲的原始信號)。式(7)綜合了小波向量稀疏優(yōu)化、矩陣稀疏優(yōu)化及各種約束問題，利用現(xiàn)有的固定變量迭代方法求解，根據(jù)研究問題的特點和算法實現(xiàn)要求，采用閾值迭代方法進(jìn)行優(yōu)化。

3 算法實現(xiàn)步驟

由稀疏分解、稀疏系數(shù)去噪優(yōu)化及矩陣建模等一系列方法得到多尺度的優(yōu)化小波系數(shù)，并進(jìn)行張量矩陣建模。對張量建模分解矩陣進(jìn)行奇異值分解，再對矩陣奇異值進(jìn)行低秩優(yōu)化，根據(jù)優(yōu)化后的奇異值結(jié)果重建邊緣信息。具體算法步驟如下。

(1)地震資料平穩(wěn)小波分解。利用平穩(wěn)小波變換方法對地震資料進(jìn)行多尺度小波分解，得到地震資料的小波稀疏字典。

(2)多尺度小波系數(shù)優(yōu)化。為了最大程度地壓制噪聲，利用有效信號和噪聲的分解系數(shù)分布特征，利用壓縮感知中的MP算法和系數(shù)特征擬合方法，優(yōu)選最優(yōu)小波系數(shù)去噪；利用優(yōu)選的分解系數(shù)通過橫向相鄰點進(jìn)行相關(guān)計算，再次選出相關(guān)系數(shù)大的小波系數(shù)。

(3)根據(jù)多尺度優(yōu)化小波系數(shù)建立張量矩陣并進(jìn)行建模。通過張量建模，構(gòu)建地震數(shù)據(jù)邊緣信息分類最優(yōu)的矩陣表示模式。

(4)張量矩陣奇異值分解。對三維地震數(shù)據(jù)子塊(立體矩陣)進(jìn)行二維矩陣分解，并對各個分解矩陣進(jìn)行奇異值分解。

(5)矩陣奇異值低秩優(yōu)化。根據(jù)目標(biāo)函數(shù)(式(7))，采用閾值迭代方法優(yōu)化計算，以提高矩陣的低秩稀疏優(yōu)化程度，從而壓制奇異值分解噪聲。

(6)雙稀疏和雙優(yōu)化結(jié)果融合與重建。將小波分解稀疏、矩陣奇異值分解稀疏和小波系數(shù)優(yōu)化、矩陣低秩稀疏優(yōu)化結(jié)果進(jìn)行雙稀疏、雙優(yōu)化與重建。

4 模型分析

為了測試、分析低秩稀疏邊緣檢測方法(下文簡稱本文方法)對邊緣信息的識別能力及方法適應(yīng)性，分別對小斷層理論模型和實際模型進(jìn)行測試、分析。圖1為小斷層反射系數(shù)模型及其邊緣檢測結(jié)果。由圖可見，本文方法(圖1e、圖1f)清晰地識別了三條小斷層，對小斷層的識別效果明顯好于相干方法(圖1c、圖1d)，其中對主頻為20Hz雷克子波合成的含噪記錄(圖中沒有展示)的相干檢測結(jié)果(圖1d)中沒有斷層顯示。上述結(jié)果說明本文方法的抗噪性較強(qiáng)。

圖1 小斷層反射系數(shù)模型及其邊緣檢測結(jié)果(a)小斷層反射系數(shù)模型；(b)主頻為25Hz雷克子波與圖a褶積形成的含噪記錄；(c)相干(25Hz)；(d)相干(20Hz)；(e)本文方法(25Hz)；(f)本文方法(20Hz)。模型中存在三條小斷層，斷距分別為3個采樣點(30m)、2個采樣點(20m)和1和個采樣點(10m)

圖2為實際斷層模型及其邊緣檢測結(jié)果。由圖可見，本文方法很好地指示了常規(guī)斷層(圖2c)，同時在斷層模型剖面(圖2a)上未能指示應(yīng)出現(xiàn)斷層的部位，相干結(jié)果較模糊、散亂(圖2b)，但本文方法較清晰地揭示了斷層并且斷層連續(xù)性好(圖2c)。上述結(jié)果說明本文方法的適用性較強(qiáng)。

圖2 實際斷層模型及其邊緣檢測結(jié)果(a)斷層模型剖面；(b)相干方法；(c)本文方法

5 實際資料應(yīng)用效果分析

利用實際資料驗證本文方法的應(yīng)用效果。圖3為不同規(guī)模矩陣的直接奇異值檢測結(jié)果。由圖可見：7×5階矩陣(橫向的窗口尺度大于縱向)檢測結(jié)果對地層邊緣壓制效果差(圖3b)；5×17階矩陣檢測結(jié)果清楚顯示了斷層，地層壓制效果較好(圖3c)。因此，對于同樣的方法，不同規(guī)模矩陣的直接奇異值檢測結(jié)果不同。矩陣行、列元素(橫、縱向窗口的尺度)的關(guān)系體現(xiàn)了低秩化揭示橫、縱向方向特性結(jié)構(gòu)信息的能力，同時說明矩陣建模對于揭示結(jié)構(gòu)信息的重要性。上述分析表明，直接對地震信號進(jìn)行奇異值分解，只能揭示較大的斷層，不能揭示小斷層，并且檢測結(jié)果模糊、不清晰。

圖3 不同規(guī)模矩陣直接奇異值檢測結(jié)果(a)地震剖面；(b)7×5階；(c)5×17階

圖4為不同尺度組合雙稀疏邊緣檢測結(jié)果，其計算過程是對地震信號進(jìn)行小波稀疏分解后再進(jìn)行奇異值分解，在計算過程中對地震信號進(jìn)行5級分解，通過把多尺度信息與雙稀疏(小波分解與矩陣分解)分析融合到算法實現(xiàn)，極大地提高了斷層尤其是小斷層的檢測能力。由圖4可見：①信號向量的稀疏分解稀疏度越大，奇異值分解低秩度也越大，即實際檢測結(jié)果分辨率提高，地層信息得到壓制。②在高頻邊緣信息具有完備性，低頻邊緣信息具有冗余性；向量稀疏促進(jìn)了矩陣低秩稀疏，再次說明雙稀疏對提高邊緣檢測能力的意義。

圖4 不同尺度組合雙稀疏邊緣檢測結(jié)果(a)第2尺度；(b)第2、第3尺度組合；(c)所有尺度組合

相似背景組成的矩陣具有低秩特征，突變信號與噪聲具有稀疏特征。如果相似背景有噪聲干擾，則背景的低秩化程度大大降低，背景越“干凈”其組成的矩陣相似性越高，矩陣低秩化程度越大，稀疏化程度也越大。“干凈”背景下存在斷層或縫洞等邊緣信息時，分離效果最好。本文的研究主旨是通過稀疏分解壓制噪聲、增強(qiáng)稀疏化、最佳低秩化，從而最大程度地分離背景信息和邊緣信息，增強(qiáng)邊緣檢測能力。

圖5為優(yōu)化前、后邊緣檢測結(jié)果。由圖可見，聯(lián)合優(yōu)化結(jié)果明顯、清晰，極大地壓制了噪聲(圖5b)。圖6為本文方法與相干法邊緣檢測結(jié)果對比。由圖可見，在小斷層檢測能力、分辨率及斷層的連續(xù)性方面，本文方法的邊緣檢測效果都明顯優(yōu)于相干方法。圖7為洞邊緣層切片，由圖可見，與相干法洞邊緣層切片(圖7b)相比，本文方法洞邊緣層切片(圖7c)對洞邊緣的識別效果更好。

圖5 優(yōu)化前(a)、后(b)邊緣檢測結(jié)果

圖6 本文方法(a)與相干法(b)邊緣檢測結(jié)果

圖7 洞邊緣層切片(a)原始切片；(b)相干方法；(c)本文方法

6 結(jié)束語

斷層和縫洞邊緣檢測實質(zhì)上是利用地震資料壓制地層信息和噪聲、突出邊緣信息的過程。地層信息具有低秩特性，噪聲和邊緣信息具有稀疏特性。提高邊緣檢測與識別能力，就是提高邊緣信息的類間分類能力和噪聲壓制能力。根據(jù)噪聲和有效信號的稀疏域分布特征，利用冪律特征約束，通過擬合方法對能量相近的噪聲和有效信號系數(shù)再次優(yōu)化。為了挖掘數(shù)據(jù)結(jié)構(gòu)信息，提高邊緣檢測能力，圍繞相似背景目標(biāo)和稀疏目標(biāo)分離，利用張量矩陣建模和奇異值分解方法，通過雙稀疏和雙優(yōu)化，即小波分解稀疏、矩陣分解稀疏和小波系數(shù)優(yōu)化、矩陣低秩優(yōu)化，增強(qiáng)了尺度分析的正交性和類間分類能力，提高了重建結(jié)果分辨率和精度。與相干方法相比，本文方法對于斷層和縫洞邊緣具有更好的刻畫能力，較全面地揭示了較小的縫洞和小斷層信息。