基于快速鲁棒主成分分析的航磁异常识别与提取:以沉积变质型铁矿为例

Identification and extraction of aeromagnetic anomalies based on fast and robust principal component analysis: A case study of sedimentary metamorphic iron deposits

摘要

鲁棒主成分分析(RPCA)作为无监督特征学习多元统计分析的数据降维方法,广泛地应用于机器学习与模式识别领域。随着找矿难度的增加,为提高航空物探快速识别与提取弱异常信息的能力,本文采用基于非精确增广拉格朗日乘子(IALM)算法的鲁棒主成分分析方法,对航磁异常数据进行区域背景场与局部场分离。在保证计算结果鲁棒的同时,该方法更适用于大数据量航磁处理。理论模型试验结果表明,该空间域信息提取与识别方法分离出的区域磁场和局部磁场,与模型正演场形态一致,并且计算效率更高。对实测数据的单次运算时间远小于精确增广拉格朗日乘子(EALM)算法,计算速度可提升约100倍。最后,将该算法用于鞍山—本溪沉积变质岩型铁矿成矿集中区的航磁异常处理与解释,分离出的局部高磁异常与磁铁矿空间位置分布具有较强的相关性。由此圈定出有潜力的磁铁矿勘查区域,鞍山、本溪、弓长岭、灯塔附近的局部高磁异常区,仍存在有待勘查的隐伏磁铁矿。

Abstract

Robust Principal Component Analysis (RPCA) is a dimensionality reduction method widely used in machine learning and pattern recognition for unsupervised feature learning in multivariate statistical analysis. In response to increasing challenges in mineral exploration, this study employed an RPCA mathematical model based on the Inexact Augmented Lagrange Multiplier (IALM) algorithm to enhance the capability of airborne geophysical exploration for the rapid identification and extraction of weak anomaly information. This approach is not only robust but also well-suited for processing large volumes of data. Specifically, it was applied to separate the regional background field from local anomalies within aeromagnetic anomaly data. Theoretical model experiments demonstrated that this spatial-domain information extraction method effectively separates regional and local magnetic fields. The results are consistent with the model’s forward modeling outcomes, and the method achieves high computational efficiency. For measured data, the computation time per dataset was significantly shorter than that of the Exact Augmented Lagrange Multiplier (EALM) algorithm, representing a speed improvement of approximately two orders of magnitude (i.e., about 100 times faster). Finally, the algorithm was applied to process and interpret aeromagnetic anomalies in the Anshan-Benxi area, a sedimentary-metamorphic iron ore metallogenic district. The extracted local high magnetic anomalies exhibit strong correlations with the spatial distribution of known magnetite deposits, successfully delineating prospective exploration targets. Significant exploration potential remains within the local high magnetic anomaly areas near Anshan, Benxi, Gongchangling, and Dengta, where hidden magnetite deposits are likely to be discovered.

Graphical abstract

关键词

鲁棒主成分分析 / 航磁异常 / IALM算法 / 沉积变质岩型铁矿

Key words

RPCA / aeromagnetic anomaly / IALM algorithm / sedimentary-metamorphic iron ore

0 引言

主成分分析(Principal Component Analysis,PCA)作为一种高效且准确的高维数据处理分析方法,在提高可解释性的同时,最大程度地减少信息损失,因此,在工程科学领域被广泛应用^[1]。由于PCA的原理前提假设噪声数据为高斯分布,导致对异常噪声十分敏感,处理严重噪声数据无法得出最优解^[2],因此Candès等^[3]在PCA基础上提出鲁棒主成分分析法(Robust Principal Component Analysis,RPCA),也称作鲁棒张量分解算法。该算法广泛应用于机器学习图像处理中,其原理是通过求解凸矩阵分解目标函数最优化问题,将数学模型扩展到求解低秩成分矩阵的秩和噪声矩阵的 范数极小化的过程,从而解决严重噪声数据的低秩矩阵逼近问题。针对RPCA数学模型计算过程的优化问题,许多学者提出不同的算法。Beck和Teboulle^[4]提出迭代阈值算法(Iterative theresholoding,IT)求解数学模型,该算法在经典梯度算法基础上加以改进,迭代形式简单且保证收敛,但迭代速度较慢;Lin等^[5]提出加速近端梯度算法(Accelerated proximal gradient,APG),以减少计算中的迭代次数;Lin等^[2]后续提出精确拉格朗日乘子(Exact augmented Lagrange multiplier,EALM)法和不精确拉格朗日乘子(Inexact augmented Lagrange multiplier,IALM)法,该方法被证明计算结果稳定、计算效率更高^[6]。本文将采用EALM和IALM两种算法求解RPCA数学模型,由于在实际工作过程中,航磁原始数据具有分辨率高、数据量大的特点,因此需要提升航磁弱信息提取的计算效率与效果。

航空物探勘查金属矿产时,提取局部异常弱信息对矿产资源快速勘查评价至关重要,将目标地质体或矿床产生的关键异常从观测的叠加重磁场中分离出来,是重磁数据处理的关键问题。为解决重磁位场分离关键技术,前人提出了许多方法,总体可分为两类,即空间域位场分离与频率域位场分离^[7]。常见的方法包括解析延拓法^[8],趋势分析法^[9],插值切割法^[10],结构逻辑法^[11],匹配滤波法^[12],正则化滤波法^[13],小波多尺度分析法^[14],最小曲率法^[15],以及互相关滤波法^[16]等。随着计算科学、信号处理、人工智能领域技术的进步,地球物理位场数据处理方法也随之变化与发展,其中有许多经典方法。Mickus等^[17]于1991年提出最小曲率位场分离方法,其特点是采用差分迭代格式,马杰等^[18]研究表明,最小曲率法在实际资料处理中,对于识别盆地构造和断裂特征分布具有一定优势,主要应用于地质构造研究;熊盛青^[11]基于统计学随机结构函数滑动窗口计算原理,于1997年提出航空物探局部弱信息提取的结构逻辑法,实现金属矿产的航空物探快速勘查评价,其主要应用于剖面数据处理;刘天佑^[19]和宋双等^[20]采用小波多尺度分解法将信号分解成不同尺度,于2008年识别出渤海湾凹陷区火山岩弱重磁异常,但该方法基于不同小波基计算结果的差异性较大,选取最优小波基较复杂;吴成平等^[21]采用自相关滤波法压制区域磁场,于2022年在山东齐河—禹城夕卡岩铁矿勘查中提取航磁异常,该方法对干扰噪声的压制较弱。由于每种重磁位场分离方法的数学原理、假设条件不同,上述方法均有各自的优缺点及应用前提^[22-23]。在实际找矿过程中,地质工作者需要将多种方法分离的结果与已知先验地质体或矿体进行对比,选择对应较好的位场分离方法。

随着人工智能机器学习数据、挖掘算法近年的快速发展,PCA作为无监督特征学习方法广泛应用于特征选择和稀疏学习^[24-27],这也为地球物理重磁位场分离提供新的解决思路。朱丹等^[28-29]为实现重磁区域场与局部场分离,先后提出基于奇异谱分析以及基于EALM算法求解RPCA数学模型的数据降维处理方法,并证明其效果优于匹配滤波和小波分析。然而,在处理实测大比例尺航磁数据时,仍存在计算耗时严重的缺点。本文在前人研究成果的基础上,利用IALM算法解算RPCA数学模型,分离局部航磁异常,提取航磁信息。首先阐述EALM和IALM算法的具体实现过程,针对理论磁异常模型进行数据试验,总结对比两种算法计算效果与效率的优缺点。由于沉积变质岩型铁矿石具有强磁性,磁异常主要由不同规模及埋深的铁矿反映,计算提取的局部航磁异常是此类矿床找矿的直接标志。且以往开展过大量研究工作,认为在已知矿深部或边部仍有未探明的铁矿资源。因此,选取鞍本沉积变质岩型铁矿成矿集中区开展理论方法验证,通过对比已知矿床的空间分布特征,既可以验证该方法识别提取效果的可靠性,也可以对全区范围内局部磁异常进行再次筛选,有利于隐伏磁铁矿的进一步勘查。

1 方法原理

假设给定的高维数据位于低维线性子空间附近,采用PCA计算的目标就是高效且准确地估计该低维子空间。假设观测数据包含高斯分布的矩阵E,观测的低秩与高秩混叠数据可以排列为大型矩阵D∈R^m×n,为将矩阵D分解为低秩矩阵A和高秩矩阵E,且使矩阵A和D间的差异最小化,其PCA数学模型归纳为求解最优化问题:

(1)$ \min _{\boldsymbol{A}, \boldsymbol{E}}\|\boldsymbol{E}\|_{F} \quad \operatorname{rank}(\boldsymbol{A}) \leqslant r, \boldsymbol{D}=\boldsymbol{A}+\boldsymbol{E}$

式中r<<min(m,n)代表子空间的目标维度,‖·‖_F是Frobenius范数,该问题可以通过首先计算D的奇异值分解,而r个主奇异值反应信号的主要特征。若E是非独立同高斯分布矩阵时,PCA将修改为RPCA双目标极小数学模型^[2]:

(2)$ \min _{\boldsymbol{A}, \boldsymbol{E}}\left[\operatorname{rank}(\boldsymbol{A}),\|\boldsymbol{E}\|_{0}\right] \quad \boldsymbol{D}=\boldsymbol{A}+\boldsymbol{E}$

式中‖E‖₀表示0范数,RPCA数学模型在保证E尽量稀疏的同时,确保了A的低秩性。引入惩罚函数进行求解,λ表示权重系数:

(3)$ \min _{\boldsymbol{A}, \boldsymbol{E}}\left[\operatorname{rank}(\boldsymbol{A})+\lambda\|\boldsymbol{E}\|_{0}\right] \quad \boldsymbol{D}=\boldsymbol{A}+\boldsymbol{E}$

为解决上述非凸的优化问题,将目标函数进行松弛处理^[30]:

(4)$ \min _{\boldsymbol{A}, \boldsymbol{E}}\left[\|\boldsymbol{A}\|_{*}+\lambda\|\boldsymbol{E}\|_{1,1}\right] \quad \boldsymbol{D}=\boldsymbol{A}+\boldsymbol{E}$

式中trace(·)表示计算矩阵的迹,‖·‖_1,1和‖·‖_*分别定义为:

(5)$ \|\boldsymbol{E}\|_{1,1}=\sum_{i=1}^{m} \sum_{j=1}^{n}\left|a_{i, j}\right|$

(6)$ \|\boldsymbol{A}\|_{*}=\max \left[\operatorname{trace}\left(\boldsymbol{U}^{T} \boldsymbol{A} \boldsymbol{V}\right): \boldsymbol{U}^{T} \boldsymbol{U}=\boldsymbol{I}_{m}, \boldsymbol{V}^{T} \boldsymbol{V}=\boldsymbol{I}_{n}\right]$

1.1 EALM与IALM算法原理

对于(4)式的凸松弛RPCA问题,需要构造增光拉格朗日函数^[31]:

(7)$ \begin{array}{c}\boldsymbol{X}=(\boldsymbol{A}, \boldsymbol{E}), f(\boldsymbol{X})=\|\boldsymbol{A}\|_{*}+\lambda\|\boldsymbol{E}\|_{1,1}, \\h(\boldsymbol{X})=\boldsymbol{D}-\boldsymbol{A}-\boldsymbol{E}\end{array}$

则拉格朗日函数为:

(8)$ \begin{array}{l}L(\boldsymbol{A}, \boldsymbol{E}, \boldsymbol{Y}, \mu)=\|\boldsymbol{A}\|_{*}+\lambda\|\boldsymbol{E}\|_{1,1}+ \\\langle\boldsymbol{Y}, \boldsymbol{D}-\boldsymbol{A}-\boldsymbol{E}\rangle+\frac{\mu}{2}\|\boldsymbol{D}-\boldsymbol{A}-E\|_{F}^{2}\end{array}$

交替迭代A、E和μ,直到满足收敛条件或达到最大迭代次数,其计算流程又可以分为EALM和IALM两种。

EALM算法流程步骤如下:

输入观测数据D∈R^m×n,权系数λ。

(1)计算:${Y}_{0}^{*}$=sgn(D)/J[sgn(D)],且令μ₀>0,ρ>1,k=0。其中sgn(D)表示对矩阵D每个元素求符号后的矩阵,J[·]=max(‖·‖₂,λ^-1‖·‖_∞)。

(2)解优化问题:(${A}_{k+1}^{*}$,${E}_{k+1}^{*}$)=arg$\underset{A,E}{min}$L(A,E,${Y}_{k}^{*}$,μ_k),数据初始化:${A}_{k+1}^{0}$=${A}_{k}^{*}$, ${E}_{k+1}^{0}$=${E}_{k}^{*}$, j=0。

(3)解优化问题:${A}_{k+1}^{j+1}$=arg$\underset{A}{min}$L(A,${E}_{k+1}^{j}$,${Y}_{k}^{*}$,μ_k),其迭代公式:

$ \boldsymbol{A}_{k+1}^{j+1}=\boldsymbol{D}_{\mu_{k}}^{-1}\left(\boldsymbol{D}-\boldsymbol{E}_{k+1}^{j}+\mu_{k}^{-1} \boldsymbol{Y}_{k}^{*}\right)$

(4)解优化问题:${E}_{k+1}^{j+1}$=arg$\underset{E}{min}$L(${A}_{k+1}^{j+1}$,E,${Y}_{k}^{*}$,μ_k),其迭代公式:

$ \boldsymbol{E}_{k+1}^{j+1}=\boldsymbol{S}_{\lambda \mu_{k}}^{-1}\left(\boldsymbol{D}-\boldsymbol{A}_{k+1}^{j+1}+\mu_{k}^{-1} \boldsymbol{Y}_{k}^{*}\right)$

(5)若${A}_{k+1}^{j+1}$或${E}_{k+1}^{j+1}$不收敛,则j=j+1,并返回步骤(3);若满足收敛条件,则

$ \begin{array}{l}\boldsymbol{A}_{k}^{*}=\boldsymbol{A}_{k+1}^{j+1} \\\boldsymbol{E}_{k}^{*}=\boldsymbol{E}_{k+1}^{j+1}\end{array}$

(6)计算:${Y}_{k+1}^{*}$=${Y}_{k}^{*}$+μ_k(D-${A}_{k+1}^{*}$-${E}_{k+1}^{*}$),且μ_k+1=ρμ_k。

(7)若满足收敛条件,则${A}_{k+1}^{*}$和${E}_{k+1}^{*}$为拉格朗日函数公式(8)的最优解;若${A}_{k+1}^{*}$、${E}_{k+1}^{*}$和${Y}_{k+1}^{*}$不收敛,则k=k+1,并返回步骤(2)。

实际上,奇异值分解(SVD)消耗大量时间,在计算过程中可以不必解决子问题:

$ \left(\boldsymbol{A}_{k+1}^{*}, \boldsymbol{E}_{k+1}^{*}\right)=\underset{A, E}{\operatorname{argmin}} L\left(\boldsymbol{A}, \boldsymbol{E}, \boldsymbol{Y}_{k}^{*}, \mu_{k}\right)$

在求解该子问题时更新A_k和E_k一次就足以使其收敛到RPCA问题的最优解。

IALM算法流程具体实现步骤如下:

输入观测数据D∈R^m×n,权系数λ。

(1)计算Y₀=D/J(D),且令E₀=0,μ₀>0,ρ>1,k=0。

(2)解优化问题:A_k+1=arg$\underset{A}{min}$L(A,E_k,Y_k,μ_k),其迭代公式:

$ \boldsymbol{A}_{k+1}=\boldsymbol{D}_{\mu_{k}}^{-1}\left(\boldsymbol{D}-\boldsymbol{E}_{k}+\mu_{k}^{-1} \boldsymbol{Y}_{k}\right)$

(3)解优化问题:E_k+1=arg$\underset{E}{min}$L(A_k+1,E,Y_k,μ_k),其迭代公式:

$ \boldsymbol{E}_{k+1}=\boldsymbol{S}_{\lambda \mu_{k}}^{-1}\left(\boldsymbol{D}-\boldsymbol{A}_{k+1}+\mu_{k}^{-1} \boldsymbol{Y}_{k}\right)$

(4)计算:Y_k+1=Y_k+μ_k(D-A_k+1-E_k+1),且μ_k+1=ρμ_k。

(5)若满足收敛条件,则A_k+1和E_k+1为拉格朗日函数公式(8)的最优解;若A_k+1、E_k+1和Y_k+1不收敛,则k=k+1,并返回步骤(2)。

参数ρ和μ₀不影响计算结果,且与收敛速度相关,本文中涉及的EALM和IALM算法均采用参数ρ=1.5,μ=1.25/‖D‖₂,大量的数值实验表明^[2],对于几何级增长的μ_k,IALM算法仍然是以Q-linear的方式收敛。

1.2 航磁异常识别提取步骤

基于IALM或EALM实现航磁异常识别提取,主要包括以下两个步骤:

(1)基于研究区域化极磁异常或原始磁异常数据,选取一系列权参数λ,并构建双目标优化模型,利用IALM或EALM算法求解低秩矩阵A和稀疏矩阵E。其中最为重要的是权系数λ的选取^[31],前人的研究结果建议,λ的最佳取值范围为0~1/$\sqrt{max(m,n)}$时,按照较小间隔取值,得到的系列结果通常较为稳态^[29]。

(2)选取不同权重系数λ计算时,分析分离后的局部磁异常与区域场的稳定性,避免计算结果发生欠拟合与过拟合现象。在具体运算过程中,可以对比已知强磁性矿体、地质体的空间分布规律,分析选取的权重系数是否合理。

2 模型数据试验

为分析基于IALM和EALM算法的航磁异常识别与提取效果的区别,本文设计如下三维磁性体理论模型试验(图1a),网格大小为100×100×100,网格间距为10。模型中包含埋深不同的场源体,浅部场源形状为“一”“L”“十”,以及球形体,四种几何形状。深部场源为“L”形柱体,其磁化率存在一定差异,浅部近地表中心球体的磁化率最弱,设定磁化率值为10 000×10^-5SI,深部地质体场源磁化率最强,设定磁化率值为30 000×10^-5SI。

2.1 理论模型正演

组合磁性体模型的正演计算磁场结果如图1b所示,深部孤立地质体正演的区域磁场见图1c,浅部多个组合异常体正演的局部磁异常见图1d。在正演计算磁场结果中,埋深不同的地质体场源产生上下叠加的磁场,将浅部弱信息进行有效提取与分离,是重磁位场数据处理的关键难点。

2.2 磁异常分离计算

本文采用EALM和IALM两种算法求解RPCA数学模型,分别对化极磁异常进行分离。基于前文所述异常识别提取的步骤原则,λ值范围应在0~0.1间逐步变化,当λ值取0.04时,分离结果见图2。对比前文中深、浅源磁性体理论模型正演计算结果,基于两种计算方法均可分离出深部磁性体的形态与幅值,见图2a与图2c。在实际找矿过程中,地质工作者更关注局部异常,两种算法均可对浅部异常弱信息有效提取和分离,分离结果见图2b与图2c。由此可以看出,基于EALM和IALM算法的RPCA方法,对航磁异常的分离效果一致,且均符合理论模型的正演磁场响应规律。

组合磁性体模型的正演计算磁场结果如图1b所示,深部孤立地质体正演的区域磁场见图1c,浅部多个组合异常体正演的局部磁异常见图1d。在正演计算磁场结果中,埋深不同的地质体场源产生上下叠加的磁场,将浅部弱信息进行有效提取与分离,是重磁位场数据处理的关键难点。

2.3 计算效率对比

当EALM和IALM算法求解RPCA数学模型计算效果相同时,两者的计算效率将决定所选算法的实用性。本文分别对理论模型正演磁场结果与鞍山—本溪铁矿成矿集中区实测磁异常数据分别进行计算,正演理论模型计算得到的磁异常数据大小为100×100网格,实测数据大小为800×1 060网格。两种试验算法均采用相同工作站处理器运算,计算效率对比结果见图3。当权重系数λ取值增大时,计算耗时整体增加,当λ过大时,两种算法均出现欠拟合现象,曲线最终回落。由于采集的实测数据中混叠的高秩异常或关键特征信息远大于理论模型,因此计算更为耗时。通过对比试验可以得出,对实测航磁数据的RPCA数学模型求解过程中,采用IALM算法的计算效率远高于EALM算法,该算法可以极大程度减少重磁位场数据处理的时间。

3 应用实例

我国铁矿床类型有沉积变质型、岩浆型、接触交代-热液型(夕卡岩型)、火山岩型、沉积型和风化淋滤型,其中以沉积变质型规模最大,是铁矿石的最重要来源^[32-34]。沉积变质岩型、岩浆型、夕卡岩型、火山岩型铁矿均有磁异常显示,风化淋滤型、沉积型异常反映较弱,该方法主要适用于沉积变质岩型、岩浆型、夕卡岩型、火山岩型矿床的位场分离计算。本文仅以鞍山—本溪典型沉积变质岩型铁矿集中区进行方法验证,目前已知区内分布有弓长岭、鞍山、南芬、歪头山和齐大山等大型鞍山式典型铁矿床^[35],鞍山深部铁矿资源远大于地表及浅部已探明的规模^[36-37]。区内磁铁矿的磁性最强,磁化率值最高约为 1 SI,是唯一可以引起强磁异常的矿体,其他岩矿石的磁性均较弱^[38-39]。鞍山—本溪地区以鞍山式沉积变质型铁矿为代表,曾开展过系列航空磁测工作,发现多处隐伏低缓航磁异常,圈闭磁铁矿引起的磁异常是直接找矿的重要标志。

3.1 地质成矿条件

研究区位于中国东北辽宁地区,该地区沉积变质型铁矿形成于新太古代,一般具有黑白相间的条带状构造,白色条带由石英燧石组成,黑色条带主要成分是铁矿物,矿石主要是磁铁矿和假象赤铁矿等,脉石主要是石英,铁矿石具有强磁性^[40-41]。区内构造发育,断裂构造主要为近EW向、NE向和NW向,地层类别出露较全,见图4。地层主要包括新生代第四系,中生界三叠系、侏罗系和白垩系,主要岩性有中酸性火山岩、火山碎屑岩等;古生界有寒武系、奥陶系、石炭系、二叠系,主要岩性有灰岩、页岩、砂岩等;新太古界主要岩性为变粒岩、浅粒岩、斜长角闪岩、片岩和磁铁石英岩,其原岩是一套以火山岩为主的火山-沉积岩建造,茨沟岩组和樱桃园岩组被认为是鞍山式铁矿的主要赋矿层^[42-44](表1)。

3.2 航磁异常特征

岩石感应磁化强度受到地球磁场的影响,磁化方向会随位置的改变而变化,磁异常的峰值与磁源体的位置存在差异,从而增加物探解释的难度,为方便对异常特征进行解释,一般要对磁异常进行化极处理^[45]。本文得到的化极磁异常表明,研究区内地质体磁性差异较大,浅源磁异常场与深源磁性地质体产生的磁异常场相叠加,见图5所示。区内北部反映为大面积的正磁异常特征,局部呈现明显强磁特征,磁场正异常幅值最大超过23 000 nT,负异常最小值约为-2 000 nT,主要分布于鞍山南部、辽阳南部区域,以及弓长岭至本溪间零星分布。

一般认为,由于沉积变质岩型铁矿自身具有较强磁性,航磁异常与铁矿床空间分布关系密切。控矿的条带状铁建造可以形成区域性的条带状磁异常,沉积变质岩型铁矿常位于条带状磁异常边部,或负背景场中的局部高磁异常^[46],从磁异常中有效分离矿致信息是寻找沉积变质岩型铁矿的关键。

3.3 矿致异常分析

采用基于IALM快速降维求解RPCA数学模型的方法分离与提取鞍本地区1∶5万航磁化极异常,λ在0~0.03间变化取值,经多次试验计算,确定当权系数λ取值为0.01时,得到区域磁异常和局部磁异常最为合理,见图6。计算的区域磁异常和局部磁异常结果稳健,局部异常极值点、梯度带与已知矿点位置相吻合。

区域背景磁异常反映埋深较大的变质地层或岩浆岩体,磁性较强的地质体在深部具有一定延伸所引起的磁异常幅值通常较高,局部磁异常(图6)主要由浅部矿床或磁性较强变质岩所引起。将局部磁异常与研究区已知铁矿点位置相叠合,可以发现,分离后的局部磁异常与铁矿床位置具有强相关性。

依据分离后的局部磁异常叠加磁铁矿矿床分布位置推断,鞍山、弓长岭、本溪、灯塔局部高磁异常区附近具有寻找隐伏磁铁矿的潜力。鞍山附近铁矿床处于局部高磁异常极大值处,以及磁梯度带边缘,其主体局部异常形态完整,呈北西向条带状分布。矿体南侧北东向分布的串珠状磁异常,反映沿断裂带两侧可能发育磁铁矿,规模尚待查证,见图7a。本溪北部孤立零星分布的高磁异常点与铁矿点位置较吻合,且附近可能存在磁铁矿引起的高磁异常,具有一定成矿潜力,见图7b。弓长岭附近已知铁矿点分布于局部高磁异常带,据分离后的局部磁异常分布特征,可大致圈闭主矿体的位置与规模,其主矿体南侧紧邻区域存在4处较强高磁异常,其中1处异常已知见矿,其余3处也可能由磁铁矿所引起。

值得指出的是,无论是火山岩型(包括海相火山岩型和陆相火山岩型),还是夕卡岩型或岩浆型(包括攀枝花式和大庙式)铁矿,铁矿石均主要由磁铁矿组成^[32-33,47],其矿床反映的磁异常规模大、强度高、形态规则,因此,也可采用RPCA方法开展试验性研究。

4 结论

本文主要研究基于IALM算法鲁棒主成分分析方法在航磁异常识别与提取中的应用,具体结论与推断如下:

(1)RPCA作为一种数据降维方法,能够很好地处理稀疏噪声信号问题,通过将数据分解为低秩和稀疏部分,对噪声和异常值具有更强的鲁棒性。选择不同的权重系数λ可实现不同尺度上的分解,为尽可能地分离出浅部矿致磁异常,需要结合已知先验地质信息进行反复试验。

(2)通过理论模型数据试验,证明IALM算法可以在保证准确率和稳定性的前提下,实现重磁位场分离,且计算效率远高于EALM算法,该算法应用于实测航磁数据时具有实用性。

(3)将IALM算法的鲁棒主成分分析方法应用于鞍山—本溪沉积变质岩型铁矿成矿集中区,基于分离提取的局部磁异常与已知磁铁矿成矿点叠加图预测磁铁矿潜在分布,推断在鞍山、弓长岭、本溪、灯塔局部高磁异常附近,仍存在有待重点勘查的磁铁矿隐伏区。

参考文献

[1] JOLLIFFE I T, CADIMA J. Principal component analysis: a review and recent developments[J]. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 2016, 374(2065): 20150202.

[2] LIN Z, CHEN M, MA Y. The augmented lagrange multiplier method for exact recovery of corrupted low-rank matrices[J]. arXiv preprint arXiv, 2010, 1009: 5055.

[3] CANDÈS E J, LI X, MA Y, et al. Robust principal component analysis?[J]. Journal of the ACM (JACM), 2011, 58(3): 1-37.

[4] BECK A, TEBOULLE M. A fast iterative shrinkage-thresholding algorithm for linear inverse problems[J]. SIAM Journal on Imaging Sciences, 2009, 2(1): 183-202.

[5] LIN Z C, GANESH A, WRIGHT J, et al. Fast convex optimization algorithms for exact recovery of a corrupted low-rank matrix[C]// 2009 3rd IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing. Aruba: IEEE, 2009.

[6] 史加荣, 郑秀云, 魏宗田, 等. 低秩矩阵恢复算法综述[J]. 计算机应用研究, 2013, 30(6): 1601-1605.

[7] 邢怡, 姚长利, 滕菲. 重磁异常分离方法技术的进展[J]. 地质论评, 2015, 61(增刊1): 883-884.

[8] 杨文采. 地球物理反演的理论与方法[M]. 北京: 地质出版社, 1997.

[9] 李春芳. 空间域位场分离方法研究[D]. 西安: 长安大学, 2011.

[10] 文百红, 程方道. 用于划分磁异常的新方法: 插值切割法[J]. 中南大学学报, 1990, 21(3): 229-235.

[11] 熊盛青. 航空物探勘查金属矿产时提取弱信息的解释方法[J]. 现代地质, 1997(1): 67-73.

[12] 施志群, 李保国. 匹配滤波应用于重力数据处理[J]. 物化探计算技术, 1987(3): 238-245.

[13] 何涛, 王万银, 黄金明, 等. 正则化方法在比值类位场边缘识别方法中的研究[J]. 物探与化探, 2019, 43(2): 308-319.

[14] 侯遵泽, 杨文采. 中国重力异常的小波变换与多尺度分析[J]. 地球物理学报, 1997(1): 85-95.

[15] 纪晓琳, 王万银, 邱之云. 最小曲率位场分离方法研究[J]. 地球物理学报, 2015, 58(3): 1042-1058.

[16] 周绵远, 强风. 一种提取弱异常的方法[J]. 地球科学与环境学报, 1994, 16(4): 66-69.

[17] MICKUS K L, AIKEN C L, KENNEDY W D. Regional-residual gravity anomaly separation using the minimum-curvature technique[J]. Geophysics, 1991, 56(2): 279-283.

[18] 马杰, 王万银, 纪晓琳. 利用重力场研究塞萨尔盆地及邻区构造特征[J]. 地质科技情报, 2019, 38(1): 285-294.

[19] 刘天佑. 位场勘探数据处理新方法[M]. 北京: 科学出版社, 2007.

[20] 宋双, 刘天佑, 苏俊青. 重磁小波多尺度分解及在火山岩识别中的应用[J]. 工程地球物理学报, 2008, 5(6): 675-679

[21] 吴成平, 于长春, 熊盛青, 等. 自相关滤波方法及其在山东齐河—禹城铁矿勘查应用[J]. 物探与化探, 2022, 46(4): 940-945.

[22] 谢磊磊, 蒋甫玉, 黄岩. 基于位场分离方法寻找梅山铁矿接替资源[J]. 地球物理学进展, 2017, 32(6): 2560-2570.

[23] 韩彤. 空间域位场分离方法对比研究[D]. 西安: 长安大学, 2019.

[24] 周志华. 机器学习[M]. 北京: 清华大学出版社, 2016.

[25] 叶舒婉, 侯卫生, 杨玠, 等. 三维地质智能建模研究进展[J]. 地学前缘, 2025, 32(4): 182-198.

[26] 郑澳月, 费金娜, 陈永清, 等. 应用奇异值分解(SVD)-主成分分析(PCA)组合模型定量圈定与评价腾冲地块锡钨和铅锌多金属找矿靶区[J]. 地学前缘, 2025, 32(1): 283-301.

[27] 陈永清, 郑澳月, 费金娜, 等. 应用二维经验模分解(2D-EMD)-主成分分析(PCA)组合模型定量圈定与评价胶东金与关键金属找矿靶区[J]. 地学前缘, 2025, 32(1): 266-282.

[28] 朱丹, 刘天佑, 李宏伟. 基于奇异谱分析的重磁位场分离方法[J]. 地球物理学报, 2018, 61(9): 3800-3811.

[29] 朱丹, 刘天佑, 李宏伟. 利用数据低秩性和稀疏性的位场分离[J]. 石油地球物理勘探, 2019, 54(4): 925-936.

[30] WRIGHT J, GANESH A, RAO S, et al. Robust principal component analysis: exact recovery of corrupted low-rank matrices via convex optimization[C/OL]// Advances in Neural Information Processing Systems, 2009, 22: 74349. http://hdl.handle.net/2142/74349.

[31] CANDES E, RECHT B. Exact matrix completion via convex optimization[J]. Communications of the ACM, 2012, 55(6): 111-119.

[32] 张招崇, 李厚民, 李建威, 等. 我国铁矿成矿背景与富铁矿成矿机制[J]. 中国科学: 地球科学, 2021, 51(6): 827-852.

[33] ZHANG Z C, HOU T, SANTOSH M, et al. Spatio-temporal distribution and tectonic settings of the major iron deposits in China: an overview[J]. Ore Geology Reviews, 2014, 57: 247-263.

[34] 李厚民, 李延河, 李立兴, 等. 沉积变质型铁矿成矿条件及富铁矿形成机制[J]. 地质学报, 2022, 96(9): 3211-3233.

[35] 李延河, 张增杰, 侯可军, 等. 辽宁鞍本地区沉积变质型富铁矿的成因: Fe, Si, O, S 同位素证据[J]. 地质学报, 2014, 88(12): 2351-2372.

[36] 熊盛青, 范正国, 黄旭钊, 等. 全国矿产资源潜力评价中磁测资料应用研究[M]. 北京: 地质出版社, 2021.

[37] 范正国, 黄旭钊, 谭林, 等. 鞍山地区地质构造及深部铁矿[J]. 地质与勘探, 2013, 49(6): 1153-1163.

[38] 范正国, 黄旭钊, 谭林, 等. 鞍山铁矿地质地球物理特征及找矿远景[J]. 矿床地质, 2012(201): 839-840.

[39] 范正国, 周道卿, 熊盛青, 等. 航磁异常找矿案例[M]. 北京: 地质出版社, 2015.

[40] 张秋生. 中国早前寒武纪地质及成矿作用[M]. 长春: 吉林人民出版社, 1988.

[41] 阴江宁, 肖克炎, 娄德波. 中国铁矿预测模型与资源潜力分析[J]. 地学前缘, 2018, 25(3): 107-117.

[42] 付海涛. 鞍本地区EW向深大断裂地质特征[J]. 地质找矿论丛, 2014, 29(4): 471-479.

[43] 张连强, 王洪波, 陈檬, 等. 辽宁省灯塔市大达连沟铁矿地质特征[J]. 地质找矿论丛, 2018, 33(1): 54-59

[44] 刘忠元, 付海涛, 刘陆山. 鞍本地区超大型铁矿分布规律及特征[J]. 化工矿产地质, 2015, 37(2): 65-70.

[45] 熊盛青, 杨海, 丁燕云, 等. 中国航磁大地构造单元划分[J]. 中国地质, 2018, 45(4): 658-680.

[46] 杨海, 熊盛青, 杨雪, 等. 中国航磁异常特征与铁矿床空间分布关系[J]. 矿床地质, 2022, 41(5): 893-916.

[47] 赵一鸣. 中国主要富铁矿床类型及地质特征[J]. 矿床地质, 2013, 32(4): 685-704.

基金资助

国家重点研发计划项目(2022YFC2903704)