利用数据低秩性和稀疏性的位场分离

朱 丹 刘天佑* 李宏伟

(①中国地质大学(武汉)地球物理与空间信息学院，湖北武汉 430074；②中国地质大学(武汉)数学与物理学院，湖北武汉 430074)

0 引言

在应用科学和工程领域，观测的数据集往往是高维的，提取高维数据的低秩结构尤为重要[1]。主成分分析(Principal Component Analysis,PCA)作为一种经典的低秩矩阵逼近方法被广泛应用。对于被独立同分布的高斯噪声污染的观测数据，经典PCA往往能够给出低秩矩阵逼近问题的最优解，然而对于被高度污染的数据(如有色噪声)却难以准确估计[2]。

为了解决数据被高度污染的低秩矩阵逼近问题，Candes等[3-5]提出鲁棒主成分分析(Robust Principal Component Analysis,RPCA)。该方法在经典PCA的基础上，将原来求解噪声矩阵Frobenius范数极小化的数学模型，拓展到求解低秩成分矩阵的秩和噪声矩阵的L0范数极小化的数学模型。基于合理的数学模型，RPCA迅速成为图像处理和数据挖掘领域的研究热点。针对RPCA数学模型的优化问题，许多学者进行了研究。史加荣等[6]采用迭代阈值算法(Iterative thresholding,IT)求解其提出的数学模型。该算法具有迭代形式简单且收敛的优点，但是由于收敛速度慢及难以选取步长，所以应用范围有限。Lin等[2]采用加速近端梯度算法(Accelerated proximal gradient,APG)求解RPCA的数学模型。与IT算法相比，APG算法迭代次数大大降低。Lin等[2]提出精确拉格朗日乘子(Exact augmented Lagrange multipliers,EALM)法和不精确拉格朗日乘子(Inexact augmented Lagrange multipliers,IALM)法。实践证明，EALM和IALM法具有更高的计算效率和精度[6]，因而逐渐成为求解RPCA数学模型的主流算法。本文采用EALM求解RPCA数学模型。

重磁场是由具有不同密度和磁性的地质体的共同响应，即不同尺度、不同幅值异常的叠加。如何从混叠重磁场中分离出目标地质体引起的异常，是重磁数据处理的重点之一。

位场的识别与分离方法包括解析延拓、滑动平均、趋势分析、匹配滤波、插值切割以及近年发展起来的小波分析、最小曲率等[7-8]。Spector等[9]推导了航磁异常的能谱公式，提出利用能谱分析粗略估计块状体埋深的方法，并在1975年提出匹配滤波法分离重磁场；程方道等[10]提出插值切割法分离磁异常；Mickus等[11]首次将最小曲率法应用于美国某地区的布格重力异常分离；王万银等[12-13]研究了位场数据最小曲率扩边和补空方法；随着小波分析法的提出，Fedu等[14]、侯遵泽等[15]和杨文采等[16]提出一种基于离散小波变换的重磁位场分离方法。

现有的位场分离方法往往会出现参数难选择、分离效果差、虚假异常等情况，这不仅会造成主观上地质认识的错误，还给后续反演解释带来偏差。因此，精度高、自适应性强、应用简单的位场分离方法是最佳选择。朱丹等[17]提出基于奇异谱分析的位场分离方法，该方法先将位场数据表示成Hankel矩阵(原始数据是一维的情况)或块Hankel矩阵(原始数据是二维的情况)，然后利用PCA或EOF(经验正交函数分解，empirical orthogonal functions)实现数据降维，从而达到区域场与局部场分离的目的。

本文RPCA是信号处理领域的一种新方法。利用RPCA实现Hankel矩阵或块Hankel矩阵的降维，RPCA可有效处理低秩矩阵逼近问题，其特点是参数设置简单、结果稳健。本文主要研究RPCA在位场分离中的应用，介绍经典PCA和RPCA的数学模型，并从位场的低秩性和稀疏性的角度分析位场分离问题；利用EALM算法求解RPCA数学模型，并提出针对一维信号的RPCA处理方法；通过理论模型展示分离效果及参数的选取方式，并将本文方法应用于实际数据处理。

1 PCA和RPCA的数学模型及EALM求解方法

假设有一高维的数据嵌套在低维线性子空间中，那么PCA和RPCA的目标就是高效且准确地提取这个低维子空间。

1.1 数学模型

假定已知一混叠低秩和高秩成分的观测矩阵D∈Rm×n，为了提取D的低秩成分，将D分解成矩阵A和E，这里A是低秩矩阵，E是服从独立同高斯分布的矩阵。那么通过PCA估计A和E可归结为求解最优化问题

(1)

式中‖·‖是Frobenius范数，其定义为

(2)

可利用奇异值分解D=UΣVT计算上述数学模型，确定r个主奇异值以及对应的左奇异向量。这r个主奇异值反映了信号的主要特征，并通过这r个主奇异值和奇异向量可以恢复出低秩矩阵A。

然而，若E是非独立同高斯分布矩阵，PCA方法不再适用。此时，将原数学模型修改为双目标极小的数学模型

(3)

式中‖E‖0是0范数，表示E中非零元素的个数。

从式(3)可以看出，与PCA只考虑‖E‖F极小不同，RPCA的数学模型要求E尽量稀疏、同时又保证了A的低秩性。

对于式(3)的双目标极小问题，可以引入惩罚函数进行求解

(4)

式中λ是权系数。

由于上述数学模型是一个非凸的优化问题，为便于求解，需要将模型中的目标函数进行松弛处理[6]

(5)

式中‖·‖1,1和‖·‖*分别表示(1,1)范数和核范数，分别定义为

(6)

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

UTU=Im,VTV=In]

(7)

式中trace(·)表示求矩阵的迹。

1.2 EALM算法

在介绍EALM算法之前，先引入两类优化问题及其最优解[2]。

令Q=(q1,q2,…,qn)为m×n维实矩阵，则优化问题

(8)

的最优解为P*=Sε(Q),Sε(Q)是关于参数ε、Q的函数，表示Q的第(i,j)元素经max(|qi,j|-ε,0)×sgn(qi,j)变换得到新的矩阵P*的第(i，j)元素，其中参数ε>0，且Sε(Q)=εS1(Q/ε)。

对于优化问题

(9)

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

下面使用EALM求解式(5)。构造增广拉格朗日函数

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

(10)

交替迭代A、E和μ，直到满足收敛条件或达到最大迭代次数为止，其计算流程如下。

(1)输入观测矩阵D∈Rm×n和权系数λ。

1.3 位场数据的延滞矩阵构建

秩的概念是针对矩阵而言的，对于重磁勘探一维观测数据，把它表示成空间延滞矩阵，即对一维数据Hankel矩阵化，对于二维观测数据，则构建块Hankel矩阵。对数据进行重排列的原因是延滞矩阵具有低秩结构[18]。

(11)

若有二维矩阵X=[xi,j]∈Rp×q，其中xi,j表示矩阵X的第i行第j列的元素。利用X的第l列构建Hankel矩阵

(12)

(13)

1.4 实现步骤

RPCA的位场分离方法包括下列3个步骤。

(2)利用RPCA实现矩阵降秩。选取合适的权参数构建双目标优化模型，并利用EALM算法求解低秩矩阵A和稀疏矩阵E。

(3)数据重构。步骤(2)中得到的低秩矩阵A和稀疏矩阵E是近似Hankel矩阵或近似块Hankel矩阵的结构，因此需要将其恢复成严格的Hankel矩阵或块Hankel矩阵。对于一维数据，采用逆对角线平均的方法；对于二维数据，先求式(13)的逆对角线上的Hankel矩阵的平均，再求各Hankel矩阵的逆对角线平均。最后，将Hankel矩阵或块Hankel矩阵恢复成数据向量或矩阵。

2 区域场的低秩性与局部场的稀疏性

矩阵秩定义为该矩阵中线性无关的行数或列数。考虑矩阵A∈Rm×m是一个方阵，其秩是矩阵A的非零特征值的个数。当A不是方阵时，即矩阵A∈Rm×n，其中m≠n，其奇异值分解A=UΣVT，且Σ=diag(σ1,σ2,…,σr,0,…,0)，其中σi是奇异值，那么矩阵A的秩是非零奇异值的个数。若rank(A)≪min(m,n)，那么称矩阵A是低秩的。

稀疏性描述的是矩阵中非零元素的个数，大多数元素为零的向量或者矩阵称为稀疏向量或稀疏矩阵[1]。考虑矩阵E∈Rm×n，若‖E‖0≪mn,那么称矩阵E是稀疏的。实际情况中，区域场具有低秩性，局部场具有稀疏性，即位场的分离问题符合RPCA的数学模型。

下面解释区域场的低秩特征。通常能够利用数据间的相关性说明矩阵秩的大小，因为矩阵不同列之间的相关性差异往往能够反映矩阵秩的差异。考虑一种极端的情况，当矩阵中每一列都是线性无关时，该矩阵为满秩矩阵，若其中某些列是线性相关时，经过初等变换，那么可以使其中部分列变换成零向量，这时该矩阵是低秩的。

区域场通常变化平缓，其相邻观测点之间具有较强的相关性。由此，区域场的延滞矩阵的相邻列之间的相关性更强，矩阵具有更低秩的特征。统计学上采用自相关系数描述这种自相关性

(14)

通过理论模型对比进一步说明区域场的低秩性。将由浅部规模较小场源引起的小尺度异常视为局部场(图1a)，将由深部规模较大场源引起的大尺度异常视为区域场(图1b)。对比区域场、局部场和高斯白噪声的自相关系数(图2)可以看出，区域场自相关性强，局部场次之，噪声不具有自相关性。因此区域场通过滑动窗口构成的延滞矩阵不同列之间具有较强的相关性，局部场的延滞矩阵不同列之间具有较弱的相关性，而噪声的延滞矩阵不同列之间不具有相关性。因此，从图3可以看出，区域场延滞矩阵比局部场延滞矩阵的非零奇异值少，噪声延滞矩阵的奇异值几乎全大于0。这说明区域场延滞矩阵的秩低于局部场，噪声延滞矩阵是满秩的。

与低秩性相比，稀疏性相对容易理解。实际应用中，向量或矩阵中的元素通常不会严格等于0，因此利用接近于0的元素或明显大于0的元素个数描述矩阵的稀疏性。深部规模大的场源引起的异常往往是大尺度的，一般会覆盖整条剖面或测区，因此异常值明显大于0的观测点个数较多，具有非稀疏性或较弱的稀疏性。对于浅部小规模场源，即使具有较大的磁化强度和剩余密度，随着观测点远离场源，它引起的异常会迅速衰减，因此仅很小的范围内位场数据数值的绝对值明显大于0。相比于区域场，局部场的稀疏性更强，并且两者稀疏性的差异通常是显著的。

图1 磁场模拟

图2 区域场、局部场和噪声的自相关系数

图3 区域场、局部场和噪声的奇异值

3 RPCA数学模型的地球物理意义

前面的分析说明尺度较大的异常具有低秩的结构，尺度较小的异常具有稀疏的结构，即区域场是低秩的，局部场是稀疏的。从区域场的逼近程度理解位场分离问题，在分离位场时往往会出现区域场的欠拟合和过拟合问题。欠拟合现象实质就是分离得到的区域场不够低秩而局部场过分稀疏；过拟合现象的实质是分离得到的区域场过于低秩导致局部场不够稀疏。利用传统位场分离方法(如匹配滤波或小波分析等)分离位场数据时，往往存在稳定性较差的问题，即参数选择的微小改变使分离结果变化显著。例如，估计对数功率谱低频线性段的斜率时，极小的估计误差就会带来深部场源似深度估计的巨大误差，这容易导致位场分离结果的欠拟合或过拟合。小波分析的应用难点在于选择哪一阶或哪几阶细节叠加构成局部场或区域场。然而，RPCA的数学模型要求E尽量稀疏的同时又保证了A的低秩性，这使系统具有更强的鲁棒性，即分离结果更稳定并且更容易避免欠拟合或过拟合现象。因为当分离的局部场过分稀疏(欠拟合)时，区域场低秩的结构被破坏并导致区域场的核范数增加；而当分离的区域场过分低秩(过拟合)时，局部场稀疏的结构被破坏并导致局部场的(1,1)范数增加。因此，RPCA数学模型的地球物理意义在于:在保证分离区域场尽可能低秩的前提下局部场尽可能的稀疏，并且，当区域场的低秩性结构和局部场的稀疏性结构被破坏时，其对应的核范数和(1,1)范数会快速增加。因此，采用RPCA分离位场很容易达成一种相对稳态，这种稳态体现在权系数变化时，分离结果几乎不变。处于这种稳态时的分离效果最佳。

4 权系数的选择

图4 权系数与分离误差的关系

5 理论模型计算

匹配滤波和小波分析是位场分离中的常用方法。匹配滤波是一种地质意义明确的位场分离方法，其基本原理是埋深不同的两种场源在对数功率谱曲线上表现出不同梯度的线性下降特征，利用这种特征的差异构制滤波器，能够分离上下叠加的地质体的场[7,19]。小波分析能够将重磁场精细地分解到不同尺度的细节上，这些不同尺度的细节反映不同尺度和深度的场源，常被用于位场的分解和分析[15-16,20-21]。本文通过理论模型计算对RPCA与这两种方法进行对比。

5.1 RPCA与匹配滤波法的对比

通过理论模型1(图5d)对比RPCA与匹配滤波分离法。模型正演设定256个观测点，分别对总场(图5a)、理论区域场(图5b)和局部场(图5c)进行离散傅里叶变换并求取其对数功率谱。图6是非负频率对数功率谱，每条曲线共129个离散对数功率点。根据Spector等[9]提出的方法，分别在对数功率谱的低频段和中高频段的近似线性下降部分取切线，并通过两条切线的斜率及在纵轴上的截距构建区域场和局部场滤波器，最终得到的分离结果如图7所示。与理论场相比，匹配滤波分离的局部场多余了一部分低频成分而区域场缺失了一部分低频成分(图7a和图7b中蓝色虚线部分)，说明分离结果与理论值不吻合。

分离效果较差的原因主要有以下几点： ①区域场频谱的能量过于集中于低频段，导致只能利用很少的频率点判断切线位置。②由于局部场在低频段同样具有较强能量，并且位场具有频率叠加性，因此局部场的低频数据叠加到总场的低频段，导致总场低频段能量高于实际区域场的能量，使得估计的切线斜率大于实际值。③匹配滤波法稳定性较差，当低频段切线的首、尾端点位置选取变化较小时，其计算的斜率值变化大，导致分离结果变化大。

与匹配滤波相比，RPCA的参数估计更为简单。利用RPCA分离位场时，只需要估计权系数，又由于权系数的取值区间宽松，且当其位于最佳区间时分离结果十分稳定，因此容易选取合适的权系数λ。利用RPCA法分离图5中理论场，当λ取值为0.03～0.10时，分离结果稳定。因此选取中间值即0.06时的分离结果(图7红色虚线)与匹配滤波法分离结果进行对比，发现RPCA分离效果更优。

图5 理论模型1及磁异常

图6 理论模型1的磁场总场、区域场和局部场的对数功率谱

图7 理论模型1的磁场匹配滤波和RPCA法分离结果

5.2 RPCA与小波分析对比

通过小波分析法将模型1引起的总场异常分成5阶细节(图8a)。可以看出，随着分解的阶数增加，细节尺度逐渐加大。将小波细节1～5阶叠加构成局部场，把4阶小波逼近作为区域场，得到分离结果如图8b和图8c。小波分析法的优势在于尺度不变性[15-16]，因此应用相对简单，缺点主要是分解的细节没有明确的地球物理意义，难以判断哪几阶细节由局部场引起，哪阶逼近由区域场引起。从理论模型1的分离效果来看，RPCA要优于小波分析。

为了对比RPCA与小波分析法在复杂模型情况下的分离效果，建立如图9f所示模型。总场由三部分构成，分别是深部模型的响应、浅部模型的响应和高幅值随机噪声。调节λ，发现λ在0.1200～0.2000时，分离结果稳定。选取λ为0.1700分离总场，得到的噪点如图9e所示。对剩余场继续分离，当权系数在0.0200～0.0400变化时，分离结果是稳定的，因此选取权系数0.0300继续分离剩余场，得到结果如图9b和图9d所示。利用小波分析分离总场，分别采用位场小波基、Daubechies小波簇、Coiflets小波簇、Symlets小波簇、Discrete Meyer小波基、Biorthogonal小波簇等15种小波基、Reverse Biorthogonal小波簇分离总场[22]，不同小波基的分离结果差异较大，其中Symlets小波基分离的区域场误差最小(均方误差为5.160nT)，分离的区域场如图9c中蓝色虚线所示。对于高幅值随机噪点，小波分析法无法分离。可以看出，RPCA法的分离效果(均方根误差为3.289nT)优于小波分析法。

图8 理论模型1的小波分析和RPCA分离结果

图9 理论模型及RPCA和小波分析法的分离结果

(a)理论总场；(b)理论区域场和RPCA法分离的区域场；(c)理论区域场和小波分析法分离的区域场；(d)理论局部场和RPCA分离的局部场；(e)噪声和RPCA分离的噪声；(f)理论模型

5.3 二维位场数据试验

建立三维地质模型2(图10)，其中浅部场源包括不同形状、不同尺度的4个地质体，深部地质体是一“L”形柱体。此模型的正演重、磁场见图11。

深部地质体正演区域重、磁异常分别见图12a、图12c，浅部地质体正演的局部重、磁异常如图12b、图12d所示。采用RPCA算法分离重磁场，当λ取值在0.0020附近变化时，所得的重力分离结果稳定，当λ取值在0.0030附近变化时，所得的重力分离结果稳定。小波分析和匹配滤波的参数选择原则与实验1一样，图12i～图12p为两种方法选择不同参数得到的最优结果。

重力场的分离中， RPCA、小波分析、匹配滤波的均方根误差分别是0.016、0.020、0.028mGal。 其中RPCA的分离误差最小，小波分析残留了较多深部地质体引起的异常(幅值约为0.05mGal)，匹配滤波分离的局部场的异常幅值远小于实际局部场异常的幅值。磁场的分离中，RPCA、小波分析、匹配滤波的均方根误差分别是45.24、60.25、46.94nT，其中RPCA的分离误差略小于匹配滤波，但小波分析和匹配滤波分离的局部场残留了部分深部地质体引起的异常(幅值分别是120nT和110nT)。可见复杂模型的重磁场分离中，RPCA的误差相对较小。

图10 模型2示意图

图11 模型2重(a)、磁(b)异常场正演结果

图12 模型2重、磁场异常图

(a)～(d)分别为理论计算的区域重力异常、局部重力异常、区域磁异常和局部磁异常；(e)～(h)分别为RPCA法分离的区域重力异常、局部重力异常、区域磁异常和局部磁异常；(i)～(l)分别为小波分析法分离的区域重力异常、局部重力异常、区域磁异常和局部磁异常；(m)～(p)分别为匹配滤波法分离的区域重力异常、局部重力异常、区域磁异常和局部磁异常

6 RPCA应用实例

卫宁北山—香山矿集区位于宁夏西部，大地构造位于北祁连走廊过渡带，东段位于阿拉善地块南缘、鄂尔多斯地块西缘的转换交界处[23]，是宁夏金属矿产勘查重要地区之一。卫宁北山位于矿集区的北部，经多年勘查工作，发现多处以金、铅、银等为主的金属矿床，如金场子金矿、照璧山铁矿、大通沟铜矿等[24]。为此在该地区开展重磁勘探，目的是研究隐伏中酸性岩体的空间侵位特征。

研究区的基底层由早古生代次深海斜坡相陆源碎屑—泥质沉积为主的浊积岩系构成。刘志坚[25]推断在卫宁北山西部的单梁山西段(二人山—金场子北)及其东部一带可能存在隐伏的中酸性岩体，与出露地表的闪长玢岩脉很可能为同源、同期形成。

由于区内出露主要以沉积岩为主，因此位场地球物理方法是寻找深部隐伏岩体的有效方法。宋新华等[23]指出，区内基底具有一定磁性，沉积地层(石炭—奥陶系)岩石标本多为无磁性或弱磁性，密度约为2.68g/cm3；岩浆岩体具有磁性，密度为2.60g/cm3。结合物性特征可知，区内岩浆岩体相对围岩具有低密、高磁的物性特征，由此推断局部低重、高磁异常主要由隐伏岩浆岩体引起。

研究区布格重力异常和化极磁异常[26]分别如图13和图14所示。可见区内布格重力异常为-20～35mGal，呈南低北高趋势；区内化极磁异常为-60～120nT，北部分布大面积正磁异常。纬度37.55°以北的磁异常均大于20nT，其中以二人山地区磁异常最强，磁异常极大值为120nT；二人山北西部磁异常明显，极大值约为90nT；二人山东部、北部磁异常较弱，极大值约为50nT。结合物性资料推断，布格重力异常的区域性变化与致密基底起伏有关，基底的隆起会引起布格重力的升高，隐伏岩浆岩体会引起局部低重异常。化极磁异常主要由磁性基底和岩浆岩体引起，其中磁性基底的磁性较弱且埋深较大，主要引起幅值较低且变化平缓的磁异常，岩浆岩体磁性较强且埋深较浅，主要引起幅值较高且变化快的磁异常。因此，本次研究重点是从布格重力异常和化极磁异常中提取出局部低重、高磁异常，具有局部低重、高磁异常组合特征的区域即是隐伏岩浆岩体的远景区。

图13 研究区布格重力异常叠加剩余重力异常图

利用RPCA法提取局部低重、高磁异常。首先利用RPCA法分离布格重力异常，得到区域重力异常和局部重力异常。当λ取约0.0020时，分离得到的区域重力异常和局部重力异常稳定，其中局部重力异常负值部分如图13所示。然后利用RPCA法分离化极磁异常，得到区域磁异常和局部磁异常。当λ约为0.0040时，分离的区域磁异常和局部磁异常稳定，其中局部磁异常正值部分如图14所示。最后，将分离的局部低重异常和局部高磁异常叠合(图15)，得到7个局部低重高磁异常叠合区(A1～A7)，推断为隐伏岩体可能赋存的区域。从叠合图可以看出，二人山地区已知出露的闪长玢岩岩脉与分离的局部低重高磁异常对应关系很好，其中局部低重、高磁异常向西仍有延伸(A7区)，预测该区为隐伏岩浆岩体分布区。此外，A1～A6区也具有局部低重、高磁异常的组合特征，地表未见岩浆岩体或岩脉出露，将其解释为隐伏岩浆岩体区。

图14 研究区化极磁异常叠加剩余磁异常图

图15 局部低重、高磁异常叠合图及岩浆岩预测区域

7 结论和讨论

本文主要研究鲁棒主成分分析在位场分离中的应用。以模型数据为基础，讨论了位场信号的低秩性和稀疏性特征，从低秩性和稀疏性的角度论述了RPCA在位场分离问题中的地球物理意义，并结合理论模型提出适用于位场分离的参数选取方法。

区域场具有低秩性和局部场具有稀疏性是利用RPCA法实现位场区域场与局部场分离的基础。从RPCA数学模型出发，该方法同时要求分离出来的区域场低秩且局部场稀疏，因此对分离出来的场具有更强的约束力，进而分离的结果更准确。由于核范数和(1,1)范数的共同约束，采用RPCA分离位场很容易达成一种相对稳态，而这一系列的相似分离结果对应位场分离的最佳结果。因此，与传统位场分离方法相比，RPCA法的精度更高，并且在一定程度上更容易避免欠拟合或过拟合现象。同时对于稀疏且大幅值的噪声，RPCA法能够有效分离。

最后，将RPCA法用于宁夏卫宁北山地区的隐伏岩浆岩体预测，将分离得到的局部低重、高磁组合异常作为判断岩浆岩发育区的标志，预测出7个火成岩发育区。

本方法的不足之处是需要对原数据进行变换得到延滞矩阵，即对一维数据Hankel矩阵化及求两类优化问题最优解需要占用一定的内存，计算速度也没有频率域快。但是随着工作站、并行机等硬件的普及算法的改进，这些问题最终可以得到解决。