羅重陽，張玉潔,2*

1 中國地質(zhì)大學(xué)（武漢）數(shù)學(xué)與物理學(xué)院，武漢 430074 2 地質(zhì)探測與評估教育部重點(diǎn)實(shí)驗(yàn)室，武漢 430074

0 引言

在地球物理勘探中，磁法勘探是以地殼中巖礦石等介質(zhì)磁性差異為基礎(chǔ)，通過觀測磁場的變化規(guī)律來探查地質(zhì)構(gòu)造、尋找礦產(chǎn)等的一種物探方法（劉天佑,2007）.它主要應(yīng)用于研究區(qū)域地質(zhì)構(gòu)造、尋找磁鐵礦、工程環(huán)境問題勘察等方面.在目前的磁法勘探中，磁異常物性反演是磁異常數(shù)據(jù)定量解釋的有效手段和研究熱點(diǎn)，它不僅能夠描述目標(biāo)地質(zhì)體的幾何和物理屬性，還能對地質(zhì)體的物性參數(shù)進(jìn)行定量計(jì)算.但是，磁異常物性反演也存在著難以攻克的關(guān)鍵性問題:由于地下不同體積的場源會產(chǎn)生相同的磁異常，導(dǎo)致反演結(jié)果存在多解性問題.

為了解決反演結(jié)果的多解性這一關(guān)鍵問題，地球物理學(xué)家將正則化技術(shù)應(yīng)用于反演中.簡單來說，正則化技術(shù)是一種在模型解空間尋找特定解的方法，常常用于求解欠定性問題.對于磁異常反演或者其他地球物理場反演而言，一種非常經(jīng)典的正則化技術(shù)是吉洪諾夫正則化（Tikhonov and Arsenin,1977），利用二次函數(shù)作為懲罰項(xiàng)加入到反演的目標(biāo)函數(shù)中，從而求解出欠定問題的特定解.Li和Oldenburg（1996,1998）將吉洪諾夫正則化應(yīng)用于重力異常與磁異常反演中，極小化模型的L2范數(shù)，以得到連續(xù)，光滑的反演結(jié)果.

為了獲得聚焦和稀疏的反演模型，學(xué)者們將稀疏約束引入目標(biāo)函數(shù)中.Last和Kubic（1983）將最小體積約束引入到重力數(shù)據(jù)反演中，提出了聚焦反演來恢復(fù)具有清晰密度邊界的地下結(jié)構(gòu)，其構(gòu)建的范數(shù)被稱為最小支撐泛函.這種方法被Portniaguine和Zhdanow（2002）進(jìn)一步用于模型梯度的情況.Pilkington（2009）利用柯西范數(shù)證明了模型參數(shù)的稀疏性，對合成數(shù)據(jù)的測試表明，基于柯西范數(shù)的稀疏反演與最小二乘反演解相比，產(chǎn)生一個(gè)更加集中的效果.Meng等（2018）等建立了零范數(shù)反演目標(biāo)函數(shù)，由于極小化零范數(shù)屬于NP難問題，現(xiàn)有技術(shù)難以直接求解，他利用光滑函數(shù)序列對L0范數(shù)進(jìn)行求解逼近求解（Mohimani et al.,2009）采用這種反演方法獲得了具有清晰邊界信息的結(jié)果.

近些年，隨著壓縮感知的飛速發(fā)展，壓縮感知在各個(gè)領(lǐng)域發(fā)揮著越來越重要的作用，包括：圖像重建（高雪,2015）、通信、地球物理（唐剛,2010;馬堅(jiān)偉,2018）等.壓縮感知作為一種尋找欠定系統(tǒng)稀疏解的技術(shù)，近些年來被一些學(xué)者用于位場數(shù)據(jù)反演中.Vatankhah等（2017）討論了基于L1范數(shù)正則化的重力數(shù)據(jù)稀疏反演，該方法在實(shí)際數(shù)據(jù)重建中，得到與已知信息基本一致的反演結(jié)果.Li等（2018）、李澤林等（2019）在Li和Oldenburg（2003）的基礎(chǔ)上，將反演目標(biāo)函數(shù)中模型的L2范數(shù)擴(kuò)展為模型的Lp（0≤p≤1）范數(shù)，提出了物性上下界約束Lp（0≤p≤1）范數(shù)反演算法，獲得了具有尖銳邊界的反演結(jié)果，并進(jìn)一步將這種方法應(yīng)用于重力數(shù)據(jù)反演中.Utsugi（2019）在L1范數(shù)正則化的基礎(chǔ)上，使用了L1-L2（即彈性網(wǎng)模型）混合正則化用于磁異常反演，這種方法克服了L1范數(shù)正則化得到過于集中的反演結(jié)果的缺點(diǎn).

本文采用基于Lp（0

圖1 磁異常反演網(wǎng)格剖分模型Fig.1 Mesh division model of magnetic anomaly inversion

1 正演問題

如圖1所示，將地下二維異常區(qū)域進(jìn)行離散化處理，劃分為多行多列的模型單元，每行每列模型單元都緊密相鄰，每個(gè)模型單元大小形狀一致，磁化率均勻；地表有一條橫跨異常區(qū)域的測線，測線上均勻地分布著觀測點(diǎn).假設(shè)模型單元的總數(shù)為N，觀測點(diǎn)的總數(shù)為M，則記第j個(gè)模型單元在第i觀測點(diǎn)產(chǎn)生的磁異常為Δdij，公式為：

Δdij=Gijmj，

（1）

其中，mj表示第j的模型單元的磁化率，Gij表示第j個(gè)模型單元的單位磁化率在第i觀測點(diǎn)產(chǎn)生的磁異常，利用磁異常正演公式（劉天佑,2007）對Gij進(jìn)行定量計(jì)算，根據(jù)位場疊加性，第i個(gè)觀測點(diǎn)的總磁異常等于場源區(qū)域所有模型單元在該觀測點(diǎn)產(chǎn)生磁異常之和，記第i個(gè)觀測點(diǎn)的總磁異常為di：

（2）

將正演公式寫為矩陣形式：

d=Gm，

（3）

其中，d=（d1,…,dM）為觀測向量，G=（Gij）i=1,…,M ,j=1,…,N為核矩陣，m=（m1,…,mN）為模型磁化率向量.

2 深度加權(quán)

在反演過程中，由于核函數(shù)會隨著位場數(shù)據(jù)的深度而衰減，使得反演結(jié)果出現(xiàn)“趨膚效應(yīng)”，即反演結(jié)果集中于地表，為了避免這種現(xiàn)象的產(chǎn)生，Li和Oldenburg（1998）提出了深度加權(quán)函數(shù)，對反演的物性參數(shù)進(jìn)行加權(quán)，用來平衡不同深度的場源引起的場的衰減，其形式為：

（4）

其中，W是一個(gè)對角矩陣，第i個(gè)對角元素為第i個(gè)模型單元的深度wi，β為參數(shù)，不同的β參數(shù)導(dǎo)致的反演結(jié)果不一樣，Oldenburg和Li（2005）認(rèn)為反演問題可以選擇不同的β值進(jìn)行嘗試，最后選出最符合地質(zhì)意義上的解，這種加權(quán)方式在各種反演算法中得到了廣泛地應(yīng)用.

3 反演

地球反演的本質(zhì)是求解一個(gè)欠定性問題，通過向量d和核矩陣G，求解向量m，由于N?M，直接求解存在嚴(yán)重多解性問題.常用的方法是正則化技術(shù)，通過施加正則項(xiàng)來減少解的非唯一性.

針對式（3）的求解，利用正則化方法，構(gòu)造反演目標(biāo)函數(shù)：

（5）

（6）

通常根據(jù)反演的需求設(shè)置不同的正則項(xiàng)，根據(jù)正則項(xiàng)的構(gòu)造，式（6）的優(yōu)化問題可分為凸優(yōu)化和非凸優(yōu)化.在凸優(yōu)化反演中，具有代表性的正則項(xiàng)是L2和L1正則項(xiàng).基于L2正則化的目標(biāo)函數(shù)處處可導(dǎo)，求解相對容易.其反演結(jié)果是光滑的，分辨率不高，常用的算法有梯度下降（Li and Oldenburg,2003;Pilkington,1997）等；基于L1正則化的目標(biāo)函數(shù)在零點(diǎn)處不可導(dǎo)，求解難度相較L2正則化大幅提高.其反演結(jié)果是稀疏的，分辨率較高，常用的算法有坐標(biāo)下降算法（Friedman et al.,2007），線性規(guī)劃（Van Zon and Roy-Chowdhury,2006），稀疏梯度投影算法（Figueiredo et al.,2007），伯格曼算法（Dong et al.,2010）等.在非凸優(yōu)化反演中，具有代表性的正則項(xiàng)是Lp（0

盡管前人在凸優(yōu)化反演和非凸優(yōu)化反演中展現(xiàn)了豐富的算法，但這些算法具有特異性，即不同的目標(biāo)函數(shù)需要選擇不同的算法進(jìn)行求解，所以尋找一種適應(yīng)性強(qiáng)的算法在反演中顯得尤為必要.針對這個(gè)問題，本文將介紹ADMM算法在反演中的應(yīng)用.

3.1 ADMM算法

ADMM算法最早分別由Glowinski和Marrocco（1975）及Gabay和Mercier（1976）年提出，近些年來，ADMM算法被重新整理和研究，并且證明其適用于大規(guī)模的分布式優(yōu)化問題（Boyd et al.,2010）.ADMM算法的核心思想是在原目標(biāo)函數(shù)基礎(chǔ)上，引入新的變量，將規(guī)模較大的目標(biāo)函數(shù)分解為一系列規(guī)模較小的子問題求解，降低了求解難度，通過協(xié)調(diào)子問題的求解，從而獲得原目標(biāo)函數(shù)的全局解.考慮如下優(yōu)化問題：

minf（x）+g（x），

（7）

其中，f、g為凸函數(shù)，變量x∈Rn.引入變量z∈Rn，將式（7）改寫為：

（8）

其中，變量z∈Rm，矩陣A∈Rn×n，矩陣B∈Rm×m以及向量c∈Rp.構(gòu)造增廣拉格朗日函數(shù)：

L（x,z,y）=f（x）+g（z）+yT（Ax+Bz-c）

（9）

其中，y∈Rn是拉格朗日乘子，ρ>0是懲罰因子.

算法的迭代過程為：

（10）

令u=（1/ρ）y，上述迭代過程寫成如下形式：

（11）

上述過程表明，ADMM算法具有適應(yīng)性很強(qiáng)的特點(diǎn)：凡具備式（7）結(jié)構(gòu)，并且函數(shù)f、g為凸函數(shù)的目標(biāo)函數(shù)優(yōu)化問題，均可通過ADMM算法進(jìn)行求解，并且可以保證其收斂性.這就為后續(xù)研究不同的反演目標(biāo)函數(shù)提供了通用算法.當(dāng)函數(shù)f、g中至少一個(gè)為非凸函數(shù)時(shí)，ADMM算法可以視為一種局部優(yōu)化算法，其最終結(jié)果與起始點(diǎn)和參數(shù)有關(guān).本文將從凸優(yōu)化反演中的L1正則化和非凸優(yōu)化中的Lp（0

3.2 基于ADMM算法的凸優(yōu)化反演

當(dāng)目標(biāo)函數(shù)的正則項(xiàng)為L1范數(shù)時(shí)，利用ADMM算法框架將式（6）改寫為：

（12）

進(jìn)一步，將式（11）轉(zhuǎn)化為極小化一系列子問題進(jìn)行求解,迭代步驟：

（13）

原問題為：

（14）

屬于凸優(yōu)化問題，對變量mw求導(dǎo)，并使其等于零，得到如下等式：

（15）

（16）

由于矩陣Gw維度比較大，式（15）直接求逆效率比較低，為了提升（15）的求解效率，采用預(yù)處理共軛梯度法（Preconditioned Conjugate Gradient,PCG），PCG算法的有優(yōu)點(diǎn)在于收斂快，計(jì)算過程中不用存儲矩陣，在重磁反演中得到廣泛應(yīng)用（Pilkington,2009；陳召曦等，2012）.在每次迭代求解后對參數(shù)進(jìn)行閾值處理，將參數(shù)限制在一定范圍.參數(shù)的約束對反演而言具有重要意義，一方面能使得參數(shù)分布在期望范圍內(nèi)，更重要的是可以提高反演的分辨率（Sun and Li,2016）.參數(shù)的約束采用上下界約束（Meng et al.,2018;Portniaguine and Zhdanov,2002；趙洋洋等，2022），形式為：

（17）

其中，0和bi分別為反演參數(shù)的上下界.

（2）zk+1的求解

原問題為：

（18）

根據(jù)式（18）可以看出，向量z的每一個(gè)分量都可以獨(dú)立求解，只需對向量z的分量zi（i=1…N）進(jìn)行優(yōu)化求解，優(yōu)化問題寫為：

（19）

（20）

繼而對解析解進(jìn)行閾值處理：

（21）

一方面，從反演目標(biāo)函數(shù)（12）的約束條件來看，變量mw和z是相等的，沒有必要對其都進(jìn)行閾值處理；另一方面，從目標(biāo)函數(shù)的正則項(xiàng)來看，求解mw時(shí)候是極小化L2范數(shù)，L2范數(shù)特點(diǎn)是對mw的每個(gè)分量給與本身的權(quán)重，即較大的分量具有較大權(quán)重，較小的分量具有較小的權(quán)重，這就使得最終求解的值往中間范圍靠攏，不容易出現(xiàn)越界值；反觀z的求解是極小化L1范數(shù)，L1范數(shù)的特點(diǎn)是對z的每個(gè)分量給與相同的權(quán)重，這就使得求解允許出現(xiàn)較大值和較小值，所以需要對最終的求解結(jié)果進(jìn)行閾值處理.因此，在算法設(shè)計(jì)中，僅對變量z進(jìn)行閾值處理.具體算法如下：

算法1:基于L1正則化的ADMM反演算法輸入:核矩陣Gw,磁異常向量d,正則化參數(shù)λ,懲罰因子ρ,深度加權(quán)矩陣W.輸出:模型參數(shù)向量m?.初始化:m0w=0,z0=0,u0=0.循環(huán):k=k+1;1.利用PCG算法求解（15）,更新mk+1w;2.利用（20）和（21）更新zk+1;3.更新uk+1;直到收斂.5.更新模型參數(shù)m?=W-1zk+1.

3.3 基于ADMM算法的非凸優(yōu)化反演

當(dāng)目標(biāo)函數(shù)為Lp（0

當(dāng)0

（22）

同理，對向量z的分量zi（i=1,…,N） 進(jìn)行優(yōu)化求解，化簡：

（23）

針對式（23）非凸優(yōu)化子問題的求解，前人做了許多工作，例如：利用凸松弛技術(shù)求解，代表算法有迭代重加權(quán)最小二乘（Chartrand and Yin,2008；Daubechies et al.,2010），迭代重加權(quán)L1（Candès et al.,2008）等，此類算法會使得式（23）的求解收斂到局部極小值，子問題的求解不精確將會進(jìn)一步影響ADMM算法在非凸問題的求解，用在此處顯然不合適，所以該子問題的求解應(yīng)該采用能夠收斂到全局極小值的算法.

當(dāng)ρ/2=1時(shí)，對于式（23）的全局極小值的求解，Zuo等（2013）在軟閾值函數(shù)基礎(chǔ)上提出了廣義軟閾值算法，具體算法如下：

算法2:廣義軟閾值（GST）算法輸入:si,λ,p,J.輸出:TGSTp（si;λ）.1.計(jì)算τGSTp（λ）=（2λ（1-p））1/（2-p）+λp（2λ（1-p））（p-1）/（2-p）;2.如果si<τGSTp（λ）,則TGSTp（si;λ）=0;3.否則,令k=0,zk=si;循環(huán)k=0,1,…,J;zk+1i=si-λp（zk）p-1;k=k+1;TGSTp（si;λ）=sgn（si）zk;結(jié)束.

利用此方法，從而得到優(yōu)化式（23）的解：

（24）

同式（21），在每次迭代后的結(jié)果進(jìn)行閾值處理：

（25）

其中，0和bi分別為反演參數(shù)的上下界.具體算法如下：

算法3:基于Lp（0

4 實(shí)驗(yàn)

為了驗(yàn)證本文提出的算法的有效性，我們采用了磁異常反演中常用的三種二維模型進(jìn)行模擬實(shí)驗(yàn)，同時(shí)，與傳統(tǒng)的基于L2范數(shù)的反演方法（Li and Oldenburg,2003）進(jìn)行了對比.

4.1 二維模擬實(shí)驗(yàn)

三種模型分別為直立板狀體模型、傾斜板狀體模型以及向斜模型.首先將地下區(qū)域進(jìn)行離散化處理，剖分為20行40列，一共800個(gè)形狀大小相同的網(wǎng)格模型單元，每個(gè)模型單元是邊長為25 m的正方形，并且每個(gè)模型體單元認(rèn)為是均勻磁化的，磁化強(qiáng)度為100 A·m-1，磁化傾角為I=45°，測線磁方位角為A=0°，即有效磁化角為45°，觀測點(diǎn)個(gè)數(shù)設(shè)置為51.

圖2、圖3、圖4是三種不同模型分別用本文提出基于的Lp范數(shù)（其中p分別取1,0.7,0.4,0.1）以及L2范數(shù)的反演結(jié)果，在這些圖中，加粗框線的位置代表場源的真實(shí)位置，顏色深淺代表每個(gè)模型單元的磁化率的大小.

從圖2—圖4中的模擬實(shí)驗(yàn)可以看出，本文提出的基于ADMM算法的L1范數(shù)和Lp（0

針對基于Lp（0

（26）

為了進(jìn)一步評價(jià)本文算法的精度，利用L1范數(shù)定義了反演結(jié)果的相對誤差，用E來表示：

（27）

其中，m代表原始結(jié)果，m*代表反演結(jié)果.每個(gè)模型的反演相對誤差見表1.

表1 反演相對誤差Table 1 Relative error of inversion

從表1中可以發(fā)現(xiàn)，在直立板狀體模型和向斜模型中，基于Lp（p=0.4,p=0.1）范數(shù)反演效果優(yōu)于基于L1范數(shù)反演，而基于Lp（p=0.7）范數(shù)反演與基于L1范數(shù)的反演效果相當(dāng)；在傾斜板狀體模型中，基于L1范數(shù)反演與基于Lp（p=0.7,p=0.4,p=0.1）范數(shù)反演效果相當(dāng)；在向斜模型中，基于Lp（p=0.7,p=0.4,p=0.1）范數(shù)反演效果均優(yōu)于基于L1范數(shù)反演；此外，在這三種模型中，本文的兩種算法都遠(yuǎn)遠(yuǎn)優(yōu)于基于L2范數(shù)反演算法效果.

圖2 直立板狀體模型反演結(jié)果（a） L1范數(shù)反演結(jié)果；（b） L1范數(shù)反演數(shù)據(jù)擬合；（c） Lp（p=0.7）范數(shù)反演結(jié)果；（d） Lp（p=0.7）范數(shù)反演數(shù)據(jù)擬合；（e） Lp（p=0.4）范數(shù)反演結(jié)果；（f） Lp（p=0.4）范數(shù)反演數(shù)據(jù)擬合；（g） Lp（p=0.1）范數(shù)反演結(jié)果；（h） Lp（p=0.1）范數(shù)反演數(shù)據(jù)擬合；（i） L2范數(shù)反演結(jié)果；（j） L2范數(shù)反演數(shù)據(jù)擬合.Fig.2 Inversion results of vertical plate model（a） L1 norm inversion results;（b） L1 norm inversion data fitting;（c） Lp（p=0.7） norm inversion results;（d） Lp（p=0.7） norm inversion data fitting;（e） Lp（p=0.4） norm inversion results;（f） Lp（p=0.4） norm inversion data fitting;（g） Lp（p=0.1） norm inversion results;（h） Lp（p=0.7） norm inversion results;（i） L2 norm inversion results;（j） L2 norm inversion data fitting.

圖3 傾斜板狀體模型反演結(jié)果（a） L1范數(shù)反演結(jié)果；（b） L1范數(shù)反演數(shù)據(jù)擬合；（c） Lp（p=0.7）范數(shù)反演結(jié)果；（d） Lp（p=0.7）范數(shù)反演數(shù)據(jù)擬合；（e） Lp（p=0.4）范數(shù)反演結(jié)果；（f） Lp（p=0.4）范數(shù)反演數(shù)據(jù)擬合；（g） Lp（p=0.1）范數(shù)反演結(jié)果；（h） Lp（p=0.1）范數(shù)反演數(shù)據(jù)擬合；（i） L2范數(shù)反演結(jié)果；（j） L2范數(shù)反演數(shù)據(jù)擬合.Fig.3 Inversion results of inclined plate model（a） L1 norm inversion results;（b） L1 norm inversion data fitting;（c） Lp（p=0.7） norm inversion results;（d） Lp（p=0.7） norm inversion data fitting;（e） Lp（p=0.4） norm inversion results;（f） Lp（p=0.4） norm inversion data fitting;（g） Lp（p=0.1） norm inversion results;（h） Lp（p=0.7） norm inversion results;（i） L2 norm inversion results;（j） L2 norm inversion data fitting.

圖4 向斜模型反演結(jié)果（a） L1范數(shù)反演結(jié)果；（b） L1范數(shù)反演數(shù)據(jù)擬合；（c） Lp（p=0.7）范數(shù)反演結(jié)果；（d）Lp（p=0.7）范數(shù)反演數(shù)據(jù)擬合；（e） Lp（p=0.4）范數(shù)反演結(jié)果；（f） Lp（p=0.4）范數(shù)反演數(shù)據(jù)擬合；（g） Lp（p=0.1）范數(shù)反演結(jié)果；（h） Lp（p=0.1）范數(shù)反演數(shù)據(jù)擬合；（i） L2范數(shù)反演結(jié)果；（j） L2范數(shù)反演數(shù)據(jù)擬合.Fig.4 Inversion results of synclinal model（a） L1 norm inversion results;（b） L1 norm inversion data fitting;（c） Lp（p=0.7） norm inversion results;（d） Lp（p=0.7） norm inversion data fitting;（e） Lp（p=0.4） norm inversion results;（f） Lp（p=0.4） norm inversion data fitting;（g） Lp（p=0.1） norm inversion results;（h） Lp（p=0.7） norm inversion results;（i） L2 norm inversion results;（j） L2 norm inversion data fitting.

結(jié)合上述的表現(xiàn)說明ADMM算法可以應(yīng)用于凸優(yōu)化反演和非凸優(yōu)化反演中.

圖5 極小化加權(quán)L1增強(qiáng)稀疏性（a） 稀疏向量m0，可行域d=Gm，和半徑為‖m0‖的范數(shù)球；（b） 存在m≠m0，使得‖m‖<‖m0‖；（c） 加權(quán)L1球，不存在m≠m0，使得‖Wm‖≤‖Wm0‖.Fig.5 Weighted L1 minimization to enhanced sparsity（a） Sparse vector m0,feasible set d=Gm,and L1 ball of radius‖m0‖;（b） There exists an m≠m0 for which ‖m‖<‖m0‖;（c） Weighted L1 ball,There exists no m≠m0 for which ‖Wm‖≤‖Wm0‖.

為了討論ADMM算法在反演中的性能，以模擬實(shí)驗(yàn)中直立板狀體模型為例分ADMM算法反演的精度（以模型的相對誤差來衡量）與迭代次數(shù)的關(guān)系.

從圖6中可以看出，對于一個(gè)固定的正則化參數(shù)，ADMM算法前期迭代過程中收斂較快，后期迭代過程中收斂較慢，反演精度要求很高的情況下，需要大量的迭代，反演精度要求不高的情況下，該算法可以以較少的迭代次數(shù)達(dá)到要求.同時(shí)，由于Lp（p=0.1）范數(shù)反演是非凸問題，較L1范數(shù)反演更復(fù)雜，并且，在求解過程中聯(lián)合了廣義軟閾值算法，使得過程更復(fù)雜，所以收斂速度更慢.

4.2 二維真實(shí)數(shù)據(jù)實(shí)驗(yàn)

實(shí)際數(shù)據(jù)采用從青海省尕林格鐵礦保護(hù)區(qū)獲取的剖面磁異常數(shù)據(jù)，用本文方法對礦區(qū)兩條沿測線的剖面數(shù)據(jù)進(jìn)行反演,測線圖7是青海省尕林格礦區(qū)磁異常平面等值線圖，該地區(qū)的主磁鐵礦被大面積的沙礫層圍巖所覆蓋，磁異常的最大值為1600 nT，磁鐵礦的平均磁化強(qiáng)度為40 A·m-1.對于該地區(qū)的磁異常反演，前人已經(jīng)做了許多工作（劉雙等,2012;Liu et al.,2014,2018）.212和196為兩條平行的穿過主要的礦體暴露區(qū)域的測線，每條測線1200 m，分別有61個(gè)觀測點(diǎn)，礦體區(qū)域劃為為20行48列的模型單元，每個(gè)模型單元為邊長25 m的正方形.

圖6 ADMM算法迭代誤差圖（a） L1范數(shù)反演迭代過程；（b） Lp（p=0.1）范數(shù)反演迭代過程.Fig.6 Iteration error diagram of ADMM algorithm（a） L1 norm inversion iterative process；（b） Lp（p=0.1） norm inversion iterative process.

圖7 青海省尕林格礦區(qū)磁異常平面等值線圖（Liu et al.,2014）Fig.7 Plane contour map of magnetic anomaly in Galinge mining area,Qinghai province（Liu et al.,2014）

在圖8和圖9中，加粗框線的位置代表場源的真實(shí)位置，顏色深淺代表每個(gè)模型單元的磁化率的大小.針對真實(shí)數(shù)據(jù)，本文的兩種算法均能大致反映出場源的位置和走向，同時(shí)，物性參數(shù)均分布在合理范圍內(nèi)，數(shù)值大小符合實(shí)際.不同之處在于基于L1范數(shù)的反演異常曲線的擬合程度較高；基于Lp（0

圖8 196測線反演結(jié)果（a） L1范數(shù)反演結(jié)果；（b） L1范數(shù)反演數(shù)據(jù)擬合；（c） Lp（p=0.7）范數(shù)反演結(jié)果；（d） Lp（p=0.7）范數(shù)反演數(shù)據(jù)擬合；（e） Lp（p=0.4）范數(shù)反演結(jié)果；（f）Lp（p=0.4）范數(shù)反演數(shù)據(jù)擬合；（g） Lp（p=0.1）范數(shù)反演結(jié)果；（h） Lp（p=0.1）范數(shù)反演數(shù)據(jù)擬合.Fig.8 Inversion results of 196 survey line（a） L1 norm inversion results;（b） L1 norm inversion data fitting;（c） Lp（p=0.7） norm inversion results;（d） Lp（p=0.7） norm inversion data fitting;（e） Lp（p=0.4） norm inversion results;（f） Lp（p=0.4） norm inversion data fitting;（g） Lp（p=0.1） norm inversion data fitting;（h） Lp（p=0.1） norm inversion data fitting.

圖9 212測線反演結(jié)果（a） L1范數(shù)反演結(jié)果；（b） L1范數(shù)反演數(shù)據(jù)擬合；（c） Lp（p=0.7）范數(shù)反演結(jié)果；（d） Lp（p=0.7）范數(shù)反演數(shù)據(jù)擬合；（e） Lp（p=0.4）范數(shù)反演結(jié)果；（f） Lp（p=0.4）范數(shù)反演數(shù)據(jù)擬合；（g） Lp（p=0.1）范數(shù)反演結(jié)果；（h） Lp（p=0.1）范數(shù)反演數(shù)據(jù)擬合.Fig.9 Inversion results of 212 survey line（a） L1 norm inversion results;（b） L1 norm inversion data fitting;（c） Lp（p=0.7） norm inversion results;（d） Lp（p=0.7） norm inversion data fitting;（e） Lp（p=0.4） norm inversion results;（f） Lp（p=0.4） norm inversion data fitting;（g） Lp（p=0.1） norm inversion data fitting;（h） Lp（p=0.1） norm inversion data fitting.

在196測線中，Lp（0

針對基于Lp（0

（2） 在反演過程中，當(dāng)正則化參數(shù)較大時(shí)，正則項(xiàng)的權(quán)重相對較大，一方面會導(dǎo)致目標(biāo)函數(shù)的誤差項(xiàng)過大，另一方面，會使得反演結(jié)果的模型單元呈現(xiàn)二值性（即0和40），但是模型單元的實(shí)際情況并非全部是均勻的，可能會導(dǎo)致擬合曲線局部的變化；當(dāng)正則化參數(shù)較小的時(shí)候，目標(biāo)函數(shù)的誤差平方項(xiàng)的權(quán)重相對較大，正則項(xiàng)的影響較小，目標(biāo)函數(shù)誤差項(xiàng)較小，同時(shí)，反演結(jié)果也具備稀疏性，不會呈現(xiàn)出二值性的特點(diǎn).

5 展望

ADMM算法在近幾年受到了較高的關(guān)注度，其主要原因除了上文提到的可以分離凸問題，降低目標(biāo)函數(shù)的求解難度之外，另外一個(gè)原因是是可用于解決大規(guī)模優(yōu)化問題.

考慮式（6）的一般化的反演問題，當(dāng)核矩陣Gw過大時(shí)，改寫為：

（28）

其中，gi是核矩陣Gw中的第i行，di是磁異常向量d的第i個(gè).將數(shù)據(jù)空間進(jìn)行分割，引入局部引入局部變量（mw）i∈Rn和全局變量z，將式（28）改寫為：

（29）

式（28）被稱為一致性優(yōu)化問題，利用ADMM算法將式（29）改寫為一系列子問題進(jìn)行求解（Boyd et al.,2010），每一個(gè)局部變量（mw）i可以單獨(dú)求解，結(jié)合并行計(jì)算的技術(shù),使局部變量（mw）i的求解同時(shí)進(jìn)行，從而達(dá)到處理大規(guī)模問題的能力.

6 結(jié)論

本文利用ADMM算法框架，將比較難求解的主問題分解為一系列相對容易求解的子問題，同時(shí)，每個(gè)子問題給出了相應(yīng)的求解步驟.為了驗(yàn)證算法的有效性，分別進(jìn)行了模擬實(shí)驗(yàn)和真實(shí)數(shù)據(jù)實(shí)驗(yàn).模擬實(shí)驗(yàn)表明，本文算法均能夠較好的反映出異常場源的形狀和位置以及物性參數(shù)的大小，其中，基于Lp（0

磁異常反演除了面臨多解性問題之外，另一個(gè)需要解決大規(guī)模數(shù)據(jù)反演.傳統(tǒng)反演算法屬于串行計(jì)算，計(jì)算能力有限，面對大規(guī)模問題難以求解，而ADMM算法具有分布式特點(diǎn)，能使算法并行化求解，達(dá)到處理大規(guī)模問題的能力，并在其他的方面已經(jīng)成功應(yīng)用，例如：高維數(shù)據(jù)變量選擇（Wang et al.,2018）,這也為今后大規(guī)模重磁反演提供了新思路.下一步計(jì)劃將結(jié)合分布式計(jì)算原理，實(shí)現(xiàn)并行計(jì)算，從而在保證反演精度的情況下，解決大規(guī)模反演問題.