朱 丹 劉天佑* 李宏偉

(①中國(guó)地質(zhì)大學(xué)(武漢)地球物理與空間信息學(xué)院，湖北武漢 430074；②中國(guó)地質(zhì)大學(xué)(武漢)數(shù)學(xué)與物理學(xué)院，湖北武漢 430074)

0 引言

在應(yīng)用科學(xué)和工程領(lǐng)域，觀測(cè)的數(shù)據(jù)集往往是高維的，提取高維數(shù)據(jù)的低秩結(jié)構(gòu)尤為重要[1]。主成分分析(Principal Component Analysis,PCA)作為一種經(jīng)典的低秩矩陣逼近方法被廣泛應(yīng)用。對(duì)于被獨(dú)立同分布的高斯噪聲污染的觀測(cè)數(shù)據(jù)，經(jīng)典PCA往往能夠給出低秩矩陣逼近問(wèn)題的最優(yōu)解，然而對(duì)于被高度污染的數(shù)據(jù)(如有色噪聲)卻難以準(zhǔn)確估計(jì)[2]。

為了解決數(shù)據(jù)被高度污染的低秩矩陣逼近問(wèn)題，Candes等[3-5]提出魯棒主成分分析(Robust Principal Component Analysis,RPCA)。該方法在經(jīng)典PCA的基礎(chǔ)上，將原來(lái)求解噪聲矩陣Frobenius范數(shù)極小化的數(shù)學(xué)模型，拓展到求解低秩成分矩陣的秩和噪聲矩陣的L0范數(shù)極小化的數(shù)學(xué)模型。基于合理的數(shù)學(xué)模型，RPCA迅速成為圖像處理和數(shù)據(jù)挖掘領(lǐng)域的研究熱點(diǎn)。針對(duì)RPCA數(shù)學(xué)模型的優(yōu)化問(wèn)題，許多學(xué)者進(jìn)行了研究。史加榮等[6]采用迭代閾值算法(Iterative thresholding,IT)求解其提出的數(shù)學(xué)模型。該算法具有迭代形式簡(jiǎn)單且收斂的優(yōu)點(diǎn)，但是由于收斂速度慢及難以選取步長(zhǎng)，所以應(yīng)用范圍有限。Lin等[2]采用加速近端梯度算法(Accelerated proximal gradient,APG)求解RPCA的數(shù)學(xué)模型。與IT算法相比，APG算法迭代次數(shù)大大降低。Lin等[2]提出精確拉格朗日乘子(Exact augmented Lagrange multipliers,EALM)法和不精確拉格朗日乘子(Inexact augmented Lagrange multipliers,IALM)法。實(shí)踐證明，EALM和IALM法具有更高的計(jì)算效率和精度[6]，因而逐漸成為求解RPCA數(shù)學(xué)模型的主流算法。本文采用EALM求解RPCA數(shù)學(xué)模型。

重磁場(chǎng)是由具有不同密度和磁性的地質(zhì)體的共同響應(yīng)，即不同尺度、不同幅值異常的疊加。如何從混疊重磁場(chǎng)中分離出目標(biāo)地質(zhì)體引起的異常，是重磁數(shù)據(jù)處理的重點(diǎn)之一。

位場(chǎng)的識(shí)別與分離方法包括解析延拓、滑動(dòng)平均、趨勢(shì)分析、匹配濾波、插值切割以及近年發(fā)展起來(lái)的小波分析、最小曲率等[7-8]。Spector等[9]推導(dǎo)了航磁異常的能譜公式，提出利用能譜分析粗略估計(jì)塊狀體埋深的方法，并在1975年提出匹配濾波法分離重磁場(chǎng)；程方道等[10]提出插值切割法分離磁異常；Mickus等[11]首次將最小曲率法應(yīng)用于美國(guó)某地區(qū)的布格重力異常分離；王萬(wàn)銀等[12-13]研究了位場(chǎng)數(shù)據(jù)最小曲率擴(kuò)邊和補(bǔ)空方法；隨著小波分析法的提出，F(xiàn)edu等[14]、侯遵澤等[15]和楊文采等[16]提出一種基于離散小波變換的重磁位場(chǎng)分離方法。

現(xiàn)有的位場(chǎng)分離方法往往會(huì)出現(xiàn)參數(shù)難選擇、分離效果差、虛假異常等情況，這不僅會(huì)造成主觀上地質(zhì)認(rèn)識(shí)的錯(cuò)誤，還給后續(xù)反演解釋帶來(lái)偏差。因此，精度高、自適應(yīng)性強(qiáng)、應(yīng)用簡(jiǎn)單的位場(chǎng)分離方法是最佳選擇。朱丹等[17]提出基于奇異譜分析的位場(chǎng)分離方法，該方法先將位場(chǎng)數(shù)據(jù)表示成Hankel矩陣(原始數(shù)據(jù)是一維的情況)或塊Hankel矩陣(原始數(shù)據(jù)是二維的情況)，然后利用PCA或EOF(經(jīng)驗(yàn)正交函數(shù)分解，empirical orthogonal functions)實(shí)現(xiàn)數(shù)據(jù)降維，從而達(dá)到區(qū)域場(chǎng)與局部場(chǎng)分離的目的。

本文RPCA是信號(hào)處理領(lǐng)域的一種新方法。利用RPCA實(shí)現(xiàn)Hankel矩陣或塊Hankel矩陣的降維，RPCA可有效處理低秩矩陣逼近問(wèn)題，其特點(diǎn)是參數(shù)設(shè)置簡(jiǎn)單、結(jié)果穩(wěn)健。本文主要研究RPCA在位場(chǎng)分離中的應(yīng)用，介紹經(jīng)典PCA和RPCA的數(shù)學(xué)模型，并從位場(chǎng)的低秩性和稀疏性的角度分析位場(chǎng)分離問(wèn)題；利用EALM算法求解RPCA數(shù)學(xué)模型，并提出針對(duì)一維信號(hào)的RPCA處理方法；通過(guò)理論模型展示分離效果及參數(shù)的選取方式，并將本文方法應(yīng)用于實(shí)際數(shù)據(jù)處理。

1 PCA和RPCA的數(shù)學(xué)模型及EALM求解方法

假設(shè)有一高維的數(shù)據(jù)嵌套在低維線(xiàn)性子空間中，那么PCA和RPCA的目標(biāo)就是高效且準(zhǔn)確地提取這個(gè)低維子空間。

1.1 數(shù)學(xué)模型

假定已知一混疊低秩和高秩成分的觀測(cè)矩陣D∈Rm×n，為了提取D的低秩成分，將D分解成矩陣A和E，這里A是低秩矩陣，E是服從獨(dú)立同高斯分布的矩陣。那么通過(guò)PCA估計(jì)A和E可歸結(jié)為求解最優(yōu)化問(wèn)題

(1)

式中‖·‖是Frobenius范數(shù)，其定義為

(2)

可利用奇異值分解D=UΣVT計(jì)算上述數(shù)學(xué)模型，確定r個(gè)主奇異值以及對(duì)應(yīng)的左奇異向量。這r個(gè)主奇異值反映了信號(hào)的主要特征，并通過(guò)這r個(gè)主奇異值和奇異向量可以恢復(fù)出低秩矩陣A。

然而，若E是非獨(dú)立同高斯分布矩陣，PCA方法不再適用。此時(shí)，將原數(shù)學(xué)模型修改為雙目標(biāo)極小的數(shù)學(xué)模型

(3)

式中‖E‖0是0范數(shù)，表示E中非零元素的個(gè)數(shù)。

從式(3)可以看出，與PCA只考慮‖E‖F(xiàn)極小不同，RPCA的數(shù)學(xué)模型要求E盡量稀疏、同時(shí)又保證了A的低秩性。

對(duì)于式(3)的雙目標(biāo)極小問(wèn)題，可以引入懲罰函數(shù)進(jìn)行求解

(4)

式中λ是權(quán)系數(shù)。

由于上述數(shù)學(xué)模型是一個(gè)非凸的優(yōu)化問(wèn)題，為便于求解，需要將模型中的目標(biāo)函數(shù)進(jìn)行松弛處理[6]

(5)

式中‖·‖1,1和‖·‖*分別表示(1,1)范數(shù)和核范數(shù)，分別定義為

(6)

‖A‖*=max[trace(UTAV)∶

UTU=Im,VTV=In]

(7)

式中trace(·)表示求矩陣的跡。

1.2 EALM算法

在介紹EALM算法之前，先引入兩類(lèi)優(yōu)化問(wèn)題及其最優(yōu)解[2]。

令Q=(q1,q2,…,qn)為m×n維實(shí)矩陣，則優(yōu)化問(wèn)題

(8)

的最優(yōu)解為P*=Sε(Q),Sε(Q)是關(guān)于參數(shù)ε、Q的函數(shù)，表示Q的第(i,j)元素經(jīng)max(|qi,j|-ε,0)×sgn(qi,j)變換得到新的矩陣P*的第(i，j)元素，其中參數(shù)ε>0，且Sε(Q)=εS1(Q/ε)。

對(duì)于優(yōu)化問(wèn)題

(9)

有閉形式解P*=Dε(Q)，其中Dε(Q)=USε(Σ)VT，UΣVT是矩陣Q的奇異值分解，并且Dε(Q)=εD1(Q/ε)。

下面使用EALM求解式(5)。構(gòu)造增廣拉格朗日函數(shù)

L(A,E,Y,μ)=‖A‖*+λ‖E‖1,1+

(10)

交替迭代A、E和μ，直到滿(mǎn)足收斂條件或達(dá)到最大迭代次數(shù)為止，其計(jì)算流程如下。

(1)輸入觀測(cè)矩陣D∈Rm×n和權(quán)系數(shù)λ。

1.3 位場(chǎng)數(shù)據(jù)的延滯矩陣構(gòu)建

秩的概念是針對(duì)矩陣而言的，對(duì)于重磁勘探一維觀測(cè)數(shù)據(jù)，把它表示成空間延滯矩陣，即對(duì)一維數(shù)據(jù)Hankel矩陣化，對(duì)于二維觀測(cè)數(shù)據(jù)，則構(gòu)建塊Hankel矩陣。對(duì)數(shù)據(jù)進(jìn)行重排列的原因是延滯矩陣具有低秩結(jié)構(gòu)[18]。

(11)

若有二維矩陣X=[xi,j]∈Rp×q，其中xi,j表示矩陣X的第i行第j列的元素。利用X的第l列構(gòu)建Hankel矩陣

(12)

(13)

1.4 實(shí)現(xiàn)步驟

RPCA的位場(chǎng)分離方法包括下列3個(gè)步驟。

(2)利用RPCA實(shí)現(xiàn)矩陣降秩。選取合適的權(quán)參數(shù)構(gòu)建雙目標(biāo)優(yōu)化模型，并利用EALM算法求解低秩矩陣A和稀疏矩陣E。

(3)數(shù)據(jù)重構(gòu)。步驟(2)中得到的低秩矩陣A和稀疏矩陣E是近似Hankel矩陣或近似塊Hankel矩陣的結(jié)構(gòu)，因此需要將其恢復(fù)成嚴(yán)格的Hankel矩陣或塊Hankel矩陣。對(duì)于一維數(shù)據(jù)，采用逆對(duì)角線(xiàn)平均的方法；對(duì)于二維數(shù)據(jù)，先求式(13)的逆對(duì)角線(xiàn)上的Hankel矩陣的平均，再求各Hankel矩陣的逆對(duì)角線(xiàn)平均。最后，將Hankel矩陣或塊Hankel矩陣恢復(fù)成數(shù)據(jù)向量或矩陣。

2 區(qū)域場(chǎng)的低秩性與局部場(chǎng)的稀疏性

矩陣秩定義為該矩陣中線(xiàn)性無(wú)關(guān)的行數(shù)或列數(shù)?？紤]矩陣A∈Rm×m是一個(gè)方陣，其秩是矩陣A的非零特征值的個(gè)數(shù)。當(dāng)A不是方陣時(shí)，即矩陣A∈Rm×n，其中m≠n，其奇異值分解A=UΣVT，且Σ=diag(σ1,σ2,…,σr,0,…,0)，其中σi是奇異值，那么矩陣A的秩是非零奇異值的個(gè)數(shù)。若rank(A)?min(m,n)，那么稱(chēng)矩陣A是低秩的。

稀疏性描述的是矩陣中非零元素的個(gè)數(shù)，大多數(shù)元素為零的向量或者矩陣稱(chēng)為稀疏向量或稀疏矩陣[1]?？紤]矩陣E∈Rm×n，若‖E‖0?mn,那么稱(chēng)矩陣E是稀疏的。實(shí)際情況中，區(qū)域場(chǎng)具有低秩性，局部場(chǎng)具有稀疏性，即位場(chǎng)的分離問(wèn)題符合RPCA的數(shù)學(xué)模型。

下面解釋區(qū)域場(chǎng)的低秩特征。通常能夠利用數(shù)據(jù)間的相關(guān)性說(shuō)明矩陣秩的大小，因?yàn)榫仃嚥煌兄g的相關(guān)性差異往往能夠反映矩陣秩的差異。考慮一種極端的情況，當(dāng)矩陣中每一列都是線(xiàn)性無(wú)關(guān)時(shí)，該矩陣為滿(mǎn)秩矩陣，若其中某些列是線(xiàn)性相關(guān)時(shí)，經(jīng)過(guò)初等變換，那么可以使其中部分列變換成零向量，這時(shí)該矩陣是低秩的。

區(qū)域場(chǎng)通常變化平緩，其相鄰觀測(cè)點(diǎn)之間具有較強(qiáng)的相關(guān)性。由此，區(qū)域場(chǎng)的延滯矩陣的相鄰列之間的相關(guān)性更強(qiáng)，矩陣具有更低秩的特征。統(tǒng)計(jì)學(xué)上采用自相關(guān)系數(shù)描述這種自相關(guān)性

(14)

通過(guò)理論模型對(duì)比進(jìn)一步說(shuō)明區(qū)域場(chǎng)的低秩性。將由淺部規(guī)模較小場(chǎng)源引起的小尺度異常視為局部場(chǎng)(圖1a)，將由深部規(guī)模較大場(chǎng)源引起的大尺度異常視為區(qū)域場(chǎng)(圖1b)。對(duì)比區(qū)域場(chǎng)、局部場(chǎng)和高斯白噪聲的自相關(guān)系數(shù)(圖2)可以看出，區(qū)域場(chǎng)自相關(guān)性強(qiáng)，局部場(chǎng)次之，噪聲不具有自相關(guān)性。因此區(qū)域場(chǎng)通過(guò)滑動(dòng)窗口構(gòu)成的延滯矩陣不同列之間具有較強(qiáng)的相關(guān)性，局部場(chǎng)的延滯矩陣不同列之間具有較弱的相關(guān)性，而噪聲的延滯矩陣不同列之間不具有相關(guān)性。因此，從圖3可以看出，區(qū)域場(chǎng)延滯矩陣比局部場(chǎng)延滯矩陣的非零奇異值少，噪聲延滯矩陣的奇異值幾乎全大于0。這說(shuō)明區(qū)域場(chǎng)延滯矩陣的秩低于局部場(chǎng)，噪聲延滯矩陣是滿(mǎn)秩的。

與低秩性相比，稀疏性相對(duì)容易理解。實(shí)際應(yīng)用中，向量或矩陣中的元素通常不會(huì)嚴(yán)格等于0，因此利用接近于0的元素或明顯大于0的元素個(gè)數(shù)描述矩陣的稀疏性。深部規(guī)模大的場(chǎng)源引起的異常往往是大尺度的，一般會(huì)覆蓋整條剖面或測(cè)區(qū)，因此異常值明顯大于0的觀測(cè)點(diǎn)個(gè)數(shù)較多，具有非稀疏性或較弱的稀疏性。對(duì)于淺部小規(guī)模場(chǎng)源，即使具有較大的磁化強(qiáng)度和剩余密度，隨著觀測(cè)點(diǎn)遠(yuǎn)離場(chǎng)源，它引起的異常會(huì)迅速衰減，因此僅很小的范圍內(nèi)位場(chǎng)數(shù)據(jù)數(shù)值的絕對(duì)值明顯大于0。相比于區(qū)域場(chǎng)，局部場(chǎng)的稀疏性更強(qiáng)，并且兩者稀疏性的差異通常是顯著的。

圖1 磁場(chǎng)模擬

圖2 區(qū)域場(chǎng)、局部場(chǎng)和噪聲的自相關(guān)系數(shù)

圖3 區(qū)域場(chǎng)、局部場(chǎng)和噪聲的奇異值

3 RPCA數(shù)學(xué)模型的地球物理意義

前面的分析說(shuō)明尺度較大的異常具有低秩的結(jié)構(gòu)，尺度較小的異常具有稀疏的結(jié)構(gòu)，即區(qū)域場(chǎng)是低秩的，局部場(chǎng)是稀疏的。從區(qū)域場(chǎng)的逼近程度理解位場(chǎng)分離問(wèn)題，在分離位場(chǎng)時(shí)往往會(huì)出現(xiàn)區(qū)域場(chǎng)的欠擬合和過(guò)擬合問(wèn)題。欠擬合現(xiàn)象實(shí)質(zhì)就是分離得到的區(qū)域場(chǎng)不夠低秩而局部場(chǎng)過(guò)分稀疏；過(guò)擬合現(xiàn)象的實(shí)質(zhì)是分離得到的區(qū)域場(chǎng)過(guò)于低秩導(dǎo)致局部場(chǎng)不夠稀疏。利用傳統(tǒng)位場(chǎng)分離方法(如匹配濾波或小波分析等)分離位場(chǎng)數(shù)據(jù)時(shí)，往往存在穩(wěn)定性較差的問(wèn)題，即參數(shù)選擇的微小改變使分離結(jié)果變化顯著。例如，估計(jì)對(duì)數(shù)功率譜低頻線(xiàn)性段的斜率時(shí)，極小的估計(jì)誤差就會(huì)帶來(lái)深部場(chǎng)源似深度估計(jì)的巨大誤差，這容易導(dǎo)致位場(chǎng)分離結(jié)果的欠擬合或過(guò)擬合。小波分析的應(yīng)用難點(diǎn)在于選擇哪一階或哪幾階細(xì)節(jié)疊加構(gòu)成局部場(chǎng)或區(qū)域場(chǎng)。然而，RPCA的數(shù)學(xué)模型要求E盡量稀疏的同時(shí)又保證了A的低秩性，這使系統(tǒng)具有更強(qiáng)的魯棒性，即分離結(jié)果更穩(wěn)定并且更容易避免欠擬合或過(guò)擬合現(xiàn)象。因?yàn)楫?dāng)分離的局部場(chǎng)過(guò)分稀疏(欠擬合)時(shí)，區(qū)域場(chǎng)低秩的結(jié)構(gòu)被破壞并導(dǎo)致區(qū)域場(chǎng)的核范數(shù)增加；而當(dāng)分離的區(qū)域場(chǎng)過(guò)分低秩(過(guò)擬合)時(shí)，局部場(chǎng)稀疏的結(jié)構(gòu)被破壞并導(dǎo)致局部場(chǎng)的(1,1)范數(shù)增加。因此，RPCA數(shù)學(xué)模型的地球物理意義在于:在保證分離區(qū)域場(chǎng)盡可能低秩的前提下局部場(chǎng)盡可能的稀疏，并且，當(dāng)區(qū)域場(chǎng)的低秩性結(jié)構(gòu)和局部場(chǎng)的稀疏性結(jié)構(gòu)被破壞時(shí)，其對(duì)應(yīng)的核范數(shù)和(1,1)范數(shù)會(huì)快速增加。因此，采用RPCA分離位場(chǎng)很容易達(dá)成一種相對(duì)穩(wěn)態(tài)，這種穩(wěn)態(tài)體現(xiàn)在權(quán)系數(shù)變化時(shí)，分離結(jié)果幾乎不變。處于這種穩(wěn)態(tài)時(shí)的分離效果最佳。

4 權(quán)系數(shù)的選擇

圖4 權(quán)系數(shù)與分離誤差的關(guān)系

5 理論模型計(jì)算

匹配濾波和小波分析是位場(chǎng)分離中的常用方法。匹配濾波是一種地質(zhì)意義明確的位場(chǎng)分離方法，其基本原理是埋深不同的兩種場(chǎng)源在對(duì)數(shù)功率譜曲線(xiàn)上表現(xiàn)出不同梯度的線(xiàn)性下降特征，利用這種特征的差異構(gòu)制濾波器，能夠分離上下疊加的地質(zhì)體的場(chǎng)[7,19]。小波分析能夠?qū)⒅卮艌?chǎng)精細(xì)地分解到不同尺度的細(xì)節(jié)上，這些不同尺度的細(xì)節(jié)反映不同尺度和深度的場(chǎng)源，常被用于位場(chǎng)的分解和分析[15-16,20-21]。本文通過(guò)理論模型計(jì)算對(duì)RPCA與這兩種方法進(jìn)行對(duì)比。

5.1 RPCA與匹配濾波法的對(duì)比

通過(guò)理論模型1(圖5d)對(duì)比RPCA與匹配濾波分離法。模型正演設(shè)定256個(gè)觀測(cè)點(diǎn)，分別對(duì)總場(chǎng)(圖5a)、理論區(qū)域場(chǎng)(圖5b)和局部場(chǎng)(圖5c)進(jìn)行離散傅里葉變換并求取其對(duì)數(shù)功率譜。圖6是非負(fù)頻率對(duì)數(shù)功率譜，每條曲線(xiàn)共129個(gè)離散對(duì)數(shù)功率點(diǎn)。根據(jù)Spector等[9]提出的方法，分別在對(duì)數(shù)功率譜的低頻段和中高頻段的近似線(xiàn)性下降部分取切線(xiàn)，并通過(guò)兩條切線(xiàn)的斜率及在縱軸上的截距構(gòu)建區(qū)域場(chǎng)和局部場(chǎng)濾波器，最終得到的分離結(jié)果如圖7所示。與理論場(chǎng)相比，匹配濾波分離的局部場(chǎng)多余了一部分低頻成分而區(qū)域場(chǎng)缺失了一部分低頻成分(圖7a和圖7b中藍(lán)色虛線(xiàn)部分)，說(shuō)明分離結(jié)果與理論值不吻合。

分離效果較差的原因主要有以下幾點(diǎn)： ①區(qū)域場(chǎng)頻譜的能量過(guò)于集中于低頻段，導(dǎo)致只能利用很少的頻率點(diǎn)判斷切線(xiàn)位置。②由于局部場(chǎng)在低頻段同樣具有較強(qiáng)能量，并且位場(chǎng)具有頻率疊加性，因此局部場(chǎng)的低頻數(shù)據(jù)疊加到總場(chǎng)的低頻段，導(dǎo)致總場(chǎng)低頻段能量高于實(shí)際區(qū)域場(chǎng)的能量，使得估計(jì)的切線(xiàn)斜率大于實(shí)際值。③匹配濾波法穩(wěn)定性較差，當(dāng)?shù)皖l段切線(xiàn)的首、尾端點(diǎn)位置選取變化較小時(shí)，其計(jì)算的斜率值變化大，導(dǎo)致分離結(jié)果變化大。

與匹配濾波相比，RPCA的參數(shù)估計(jì)更為簡(jiǎn)單。利用RPCA分離位場(chǎng)時(shí)，只需要估計(jì)權(quán)系數(shù)，又由于權(quán)系數(shù)的取值區(qū)間寬松，且當(dāng)其位于最佳區(qū)間時(shí)分離結(jié)果十分穩(wěn)定，因此容易選取合適的權(quán)系數(shù)λ。利用RPCA法分離圖5中理論場(chǎng)，當(dāng)λ取值為0.03～0.10時(shí)，分離結(jié)果穩(wěn)定。因此選取中間值即0.06時(shí)的分離結(jié)果(圖7紅色虛線(xiàn))與匹配濾波法分離結(jié)果進(jìn)行對(duì)比，發(fā)現(xiàn)RPCA分離效果更優(yōu)。

圖5 理論模型1及磁異常

圖6 理論模型1的磁場(chǎng)總場(chǎng)、區(qū)域場(chǎng)和局部場(chǎng)的對(duì)數(shù)功率譜

圖7 理論模型1的磁場(chǎng)匹配濾波和RPCA法分離結(jié)果

5.2 RPCA與小波分析對(duì)比

通過(guò)小波分析法將模型1引起的總場(chǎng)異常分成5階細(xì)節(jié)(圖8a)。可以看出，隨著分解的階數(shù)增加，細(xì)節(jié)尺度逐漸加大。將小波細(xì)節(jié)1～5階疊加構(gòu)成局部場(chǎng)，把4階小波逼近作為區(qū)域場(chǎng)，得到分離結(jié)果如圖8b和圖8c。小波分析法的優(yōu)勢(shì)在于尺度不變性[15-16]，因此應(yīng)用相對(duì)簡(jiǎn)單，缺點(diǎn)主要是分解的細(xì)節(jié)沒(méi)有明確的地球物理意義，難以判斷哪幾階細(xì)節(jié)由局部場(chǎng)引起，哪階逼近由區(qū)域場(chǎng)引起。從理論模型1的分離效果來(lái)看，RPCA要優(yōu)于小波分析。

為了對(duì)比RPCA與小波分析法在復(fù)雜模型情況下的分離效果，建立如圖9f所示模型。總場(chǎng)由三部分構(gòu)成，分別是深部模型的響應(yīng)、淺部模型的響應(yīng)和高幅值隨機(jī)噪聲。調(diào)節(jié)λ，發(fā)現(xiàn)λ在0.1200～0.2000時(shí)，分離結(jié)果穩(wěn)定。選取λ為0.1700分離總場(chǎng)，得到的噪點(diǎn)如圖9e所示。對(duì)剩余場(chǎng)繼續(xù)分離，當(dāng)權(quán)系數(shù)在0.0200～0.0400變化時(shí)，分離結(jié)果是穩(wěn)定的，因此選取權(quán)系數(shù)0.0300繼續(xù)分離剩余場(chǎng)，得到結(jié)果如圖9b和圖9d所示。利用小波分析分離總場(chǎng)，分別采用位場(chǎng)小波基、Daubechies小波簇、Coiflets小波簇、Symlets小波簇、Discrete Meyer小波基、Biorthogonal小波簇等15種小波基、Reverse Biorthogonal小波簇分離總場(chǎng)[22]，不同小波基的分離結(jié)果差異較大，其中Symlets小波基分離的區(qū)域場(chǎng)誤差最小(均方誤差為5.160nT)，分離的區(qū)域場(chǎng)如圖9c中藍(lán)色虛線(xiàn)所示。對(duì)于高幅值隨機(jī)噪點(diǎn)，小波分析法無(wú)法分離?？梢钥闯觯琑PCA法的分離效果(均方根誤差為3.289nT)優(yōu)于小波分析法。

圖8 理論模型1的小波分析和RPCA分離結(jié)果

圖9 理論模型及RPCA和小波分析法的分離結(jié)果

(a)理論總場(chǎng)；(b)理論區(qū)域場(chǎng)和RPCA法分離的區(qū)域場(chǎng)；(c)理論區(qū)域場(chǎng)和小波分析法分離的區(qū)域場(chǎng)；(d)理論局部場(chǎng)和RPCA分離的局部場(chǎng)；(e)噪聲和RPCA分離的噪聲；(f)理論模型

5.3 二維位場(chǎng)數(shù)據(jù)試驗(yàn)

建立三維地質(zhì)模型2(圖10)，其中淺部場(chǎng)源包括不同形狀、不同尺度的4個(gè)地質(zhì)體，深部地質(zhì)體是一“L”形柱體。此模型的正演重、磁場(chǎng)見(jiàn)圖11。

深部地質(zhì)體正演區(qū)域重、磁異常分別見(jiàn)圖12a、圖12c，淺部地質(zhì)體正演的局部重、磁異常如圖12b、圖12d所示。采用RPCA算法分離重磁場(chǎng)，當(dāng)λ取值在0.0020附近變化時(shí)，所得的重力分離結(jié)果穩(wěn)定，當(dāng)λ取值在0.0030附近變化時(shí)，所得的重力分離結(jié)果穩(wěn)定。小波分析和匹配濾波的參數(shù)選擇原則與實(shí)驗(yàn)1一樣，圖12i～圖12p為兩種方法選擇不同參數(shù)得到的最優(yōu)結(jié)果。

重力場(chǎng)的分離中， RPCA、小波分析、匹配濾波的均方根誤差分別是0.016、0.020、0.028mGal。 其中RPCA的分離誤差最小，小波分析殘留了較多深部地質(zhì)體引起的異常(幅值約為0.05mGal)，匹配濾波分離的局部場(chǎng)的異常幅值遠(yuǎn)小于實(shí)際局部場(chǎng)異常的幅值。磁場(chǎng)的分離中，RPCA、小波分析、匹配濾波的均方根誤差分別是45.24、60.25、46.94nT，其中RPCA的分離誤差略小于匹配濾波，但小波分析和匹配濾波分離的局部場(chǎng)殘留了部分深部地質(zhì)體引起的異常(幅值分別是120nT和110nT)?？梢?jiàn)復(fù)雜模型的重磁場(chǎng)分離中，RPCA的誤差相對(duì)較小。

圖10 模型2示意圖

圖11 模型2重(a)、磁(b)異常場(chǎng)正演結(jié)果

圖12 模型2重、磁場(chǎng)異常圖

(a)～(d)分別為理論計(jì)算的區(qū)域重力異常、局部重力異常、區(qū)域磁異常和局部磁異常；(e)～(h)分別為RPCA法分離的區(qū)域重力異常、局部重力異常、區(qū)域磁異常和局部磁異常；(i)～(l)分別為小波分析法分離的區(qū)域重力異常、局部重力異常、區(qū)域磁異常和局部磁異常；(m)～(p)分別為匹配濾波法分離的區(qū)域重力異常、局部重力異常、區(qū)域磁異常和局部磁異常

6 RPCA應(yīng)用實(shí)例

衛(wèi)寧北山—香山礦集區(qū)位于寧夏西部，大地構(gòu)造位于北祁連走廊過(guò)渡帶，東段位于阿拉善地塊南緣、鄂爾多斯地塊西緣的轉(zhuǎn)換交界處[23]，是寧夏金屬礦產(chǎn)勘查重要地區(qū)之一。衛(wèi)寧北山位于礦集區(qū)的北部，經(jīng)多年勘查工作，發(fā)現(xiàn)多處以金、鉛、銀等為主的金屬礦床，如金場(chǎng)子金礦、照璧山鐵礦、大通溝銅礦等[24]。為此在該地區(qū)開(kāi)展重磁勘探，目的是研究隱伏中酸性巖體的空間侵位特征。

研究區(qū)的基底層由早古生代次深海斜坡相陸源碎屑—泥質(zhì)沉積為主的濁積巖系構(gòu)成。劉志堅(jiān)[25]推斷在衛(wèi)寧北山西部的單梁山西段(二人山—金場(chǎng)子北)及其東部一帶可能存在隱伏的中酸性巖體，與出露地表的閃長(zhǎng)玢巖脈很可能為同源、同期形成。

由于區(qū)內(nèi)出露主要以沉積巖為主，因此位場(chǎng)地球物理方法是尋找深部隱伏巖體的有效方法。宋新華等[23]指出，區(qū)內(nèi)基底具有一定磁性，沉積地層(石炭—奧陶系)巖石標(biāo)本多為無(wú)磁性或弱磁性，密度約為2.68g/cm3；巖漿巖體具有磁性，密度為2.60g/cm3。結(jié)合物性特征可知，區(qū)內(nèi)巖漿巖體相對(duì)圍巖具有低密、高磁的物性特征，由此推斷局部低重、高磁異常主要由隱伏巖漿巖體引起。

研究區(qū)布格重力異常和化極磁異常[26]分別如圖13和圖14所示?？梢?jiàn)區(qū)內(nèi)布格重力異常為-20～35mGal，呈南低北高趨勢(shì)；區(qū)內(nèi)化極磁異常為-60～120nT，北部分布大面積正磁異常。緯度37.55°以北的磁異常均大于20nT，其中以二人山地區(qū)磁異常最強(qiáng)，磁異常極大值為120nT；二人山北西部磁異常明顯，極大值約為90nT；二人山東部、北部磁異常較弱，極大值約為50nT。結(jié)合物性資料推斷，布格重力異常的區(qū)域性變化與致密基底起伏有關(guān)，基底的隆起會(huì)引起布格重力的升高，隱伏巖漿巖體會(huì)引起局部低重異常?；瘶O磁異常主要由磁性基底和巖漿巖體引起，其中磁性基底的磁性較弱且埋深較大，主要引起幅值較低且變化平緩的磁異常，巖漿巖體磁性較強(qiáng)且埋深較淺，主要引起幅值較高且變化快的磁異常。因此，本次研究重點(diǎn)是從布格重力異常和化極磁異常中提取出局部低重、高磁異常，具有局部低重、高磁異常組合特征的區(qū)域即是隱伏巖漿巖體的遠(yuǎn)景區(qū)。

圖13 研究區(qū)布格重力異常疊加剩余重力異常圖

利用RPCA法提取局部低重、高磁異常。首先利用RPCA法分離布格重力異常，得到區(qū)域重力異常和局部重力異常。當(dāng)λ取約0.0020時(shí)，分離得到的區(qū)域重力異常和局部重力異常穩(wěn)定，其中局部重力異常負(fù)值部分如圖13所示。然后利用RPCA法分離化極磁異常，得到區(qū)域磁異常和局部磁異常。當(dāng)λ約為0.0040時(shí)，分離的區(qū)域磁異常和局部磁異常穩(wěn)定，其中局部磁異常正值部分如圖14所示。最后，將分離的局部低重異常和局部高磁異常疊合(圖15)，得到7個(gè)局部低重高磁異常疊合區(qū)(A1～A7)，推斷為隱伏巖體可能賦存的區(qū)域。從疊合圖可以看出，二人山地區(qū)已知出露的閃長(zhǎng)玢巖巖脈與分離的局部低重高磁異常對(duì)應(yīng)關(guān)系很好，其中局部低重、高磁異常向西仍有延伸(A7區(qū))，預(yù)測(cè)該區(qū)為隱伏巖漿巖體分布區(qū)。此外，A1～A6區(qū)也具有局部低重、高磁異常的組合特征，地表未見(jiàn)巖漿巖體或巖脈出露，將其解釋為隱伏巖漿巖體區(qū)。

圖14 研究區(qū)化極磁異常疊加剩余磁異常圖

圖15 局部低重、高磁異常疊合圖及巖漿巖預(yù)測(cè)區(qū)域

7 結(jié)論和討論

本文主要研究魯棒主成分分析在位場(chǎng)分離中的應(yīng)用。以模型數(shù)據(jù)為基礎(chǔ)，討論了位場(chǎng)信號(hào)的低秩性和稀疏性特征，從低秩性和稀疏性的角度論述了RPCA在位場(chǎng)分離問(wèn)題中的地球物理意義，并結(jié)合理論模型提出適用于位場(chǎng)分離的參數(shù)選取方法。

區(qū)域場(chǎng)具有低秩性和局部場(chǎng)具有稀疏性是利用RPCA法實(shí)現(xiàn)位場(chǎng)區(qū)域場(chǎng)與局部場(chǎng)分離的基礎(chǔ)。從RPCA數(shù)學(xué)模型出發(fā)，該方法同時(shí)要求分離出來(lái)的區(qū)域場(chǎng)低秩且局部場(chǎng)稀疏，因此對(duì)分離出來(lái)的場(chǎng)具有更強(qiáng)的約束力，進(jìn)而分離的結(jié)果更準(zhǔn)確。由于核范數(shù)和(1,1)范數(shù)的共同約束，采用RPCA分離位場(chǎng)很容易達(dá)成一種相對(duì)穩(wěn)態(tài)，而這一系列的相似分離結(jié)果對(duì)應(yīng)位場(chǎng)分離的最佳結(jié)果。因此，與傳統(tǒng)位場(chǎng)分離方法相比，RPCA法的精度更高，并且在一定程度上更容易避免欠擬合或過(guò)擬合現(xiàn)象。同時(shí)對(duì)于稀疏且大幅值的噪聲，RPCA法能夠有效分離。

最后，將RPCA法用于寧夏衛(wèi)寧北山地區(qū)的隱伏巖漿巖體預(yù)測(cè)，將分離得到的局部低重、高磁組合異常作為判斷巖漿巖發(fā)育區(qū)的標(biāo)志，預(yù)測(cè)出7個(gè)火成巖發(fā)育區(qū)。

本方法的不足之處是需要對(duì)原數(shù)據(jù)進(jìn)行變換得到延滯矩陣，即對(duì)一維數(shù)據(jù)Hankel矩陣化及求兩類(lèi)優(yōu)化問(wèn)題最優(yōu)解需要占用一定的內(nèi)存，計(jì)算速度也沒(méi)有頻率域快。但是隨著工作站、并行機(jī)等硬件的普及算法的改進(jìn)，這些問(wèn)題最終可以得到解決。