朱 丹 劉天佑* 李宏偉

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

0 引言

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

為了解決數(shù)據(jù)被高度污染的低秩矩陣逼近問題，Candes等[3-5]提出魯棒主成分分析(Robust Principal Component Analysis,RPCA)。該方法在經(jīng)典PCA的基礎(chǔ)上，將原來求解噪聲矩陣Frobenius范數(shù)極小化的數(shù)學(xué)模型，拓展到求解低秩成分矩陣的秩和噪聲矩陣的L0范數(shù)極小化的數(shù)學(xué)模型?；诤侠淼臄?shù)學(xué)模型，RPCA迅速成為圖像處理和數(shù)據(jù)挖掘領(lǐng)域的研究熱點(diǎn)。針對RPCA數(shù)學(xué)模型的優(yōu)化問題，許多學(xué)者進(jìn)行了研究。史加榮等[6]采用迭代閾值算法(Iterative thresholding,IT)求解其提出的數(shù)學(xué)模型。該算法具有迭代形式簡單且收斂的優(yōu)點(diǎn)，但是由于收斂速度慢及難以選取步長，所以應(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é)模型。

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

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

現(xiàn)有的位場分離方法往往會(huì)出現(xiàn)參數(shù)難選擇、分離效果差、虛假異常等情況，這不僅會(huì)造成主觀上地質(zhì)認(rèn)識(shí)的錯(cuò)誤，還給后續(xù)反演解釋帶來偏差。因此，精度高、自適應(yīng)性強(qiáng)、應(yīng)用簡單的位場分離方法是最佳選擇。朱丹等[17]提出基于奇異譜分析的位場分離方法，該方法先將位場數(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ū)域場與局部場分離的目的。

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

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

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

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

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

(1)

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

(2)

可利用奇異值分解D=UΣVT計(jì)算上述數(shù)學(xué)模型，確定r個(gè)主奇異值以及對應(yīng)的左奇異向量。這r個(gè)主奇異值反映了信號(hào)的主要特征，并通過這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的低秩性。

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

(4)

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

由于上述數(shù)學(xué)模型是一個(gè)非凸的優(yōu)化問題，為便于求解，需要將模型中的目標(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算法之前，先引入兩類優(yōu)化問題及其最優(yōu)解[2]。

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

(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/ε)。

對于優(yōu)化問題

(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和μ，直到滿足收斂條件或達(dá)到最大迭代次數(shù)為止，其計(jì)算流程如下。

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

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

秩的概念是針對矩陣而言的，對于重磁勘探一維觀測數(shù)據(jù)，把它表示成空間延滯矩陣，即對一維數(shù)據(jù)Hankel矩陣化，對于二維觀測數(shù)據(jù)，則構(gòu)建塊Hankel矩陣。對數(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的位場分離方法包括下列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矩陣。對于一維數(shù)據(jù)，采用逆對角線平均的方法；對于二維數(shù)據(jù)，先求式(13)的逆對角線上的Hankel矩陣的平均，再求各Hankel矩陣的逆對角線平均。最后，將Hankel矩陣或塊Hankel矩陣恢復(fù)成數(shù)據(jù)向量或矩陣。

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

矩陣秩定義為該矩陣中線性無關(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)，那么稱矩陣A是低秩的。

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

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

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

(14)

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

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

圖1 磁場模擬

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

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

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

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

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

5.1 RPCA與匹配濾波法的對比

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

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

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

圖5 理論模型1及磁異常

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

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

5.2 RPCA與小波分析對比

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

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

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

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

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

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

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

深部地質(zhì)體正演區(qū)域重、磁異常分別見圖12a、圖12c，淺部地質(zhì)體正演的局部重、磁異常如圖12b、圖12d所示。采用RPCA算法分離重磁場，當(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é)果。

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

圖10 模型2示意圖

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

圖12 模型2重、磁場異常圖

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

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

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

研究區(qū)布格重力異常和化極磁異常[26]分別如圖13和圖14所示?？梢妳^(qū)內(nèi)布格重力異常為-20～35mGal，呈南低北高趨勢；區(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ū)已知出露的閃長玢巖巖脈與分離的局部低重高磁異常對應(yīng)關(guān)系很好，其中局部低重、高磁異常向西仍有延伸(A7區(qū))，預(yù)測該區(qū)為隱伏巖漿巖體分布區(qū)。此外，A1～A6區(qū)也具有局部低重、高磁異常的組合特征，地表未見巖漿巖體或巖脈出露，將其解釋為隱伏巖漿巖體區(qū)。

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

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

7 結(jié)論和討論

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

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

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

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