基于交替方向乘子算法的二維磁異常稀疏反演

羅重陽，張玉潔,2*

1 中國地質大學（武漢）數學與物理學院，武漢 430074 2 地質探測與評估教育部重點實驗室，武漢 430074

0 引言

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

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

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

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

本文采用基于Lp（0

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

1 正演問題

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

Δdij=Gijmj，

（1）

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

（2）

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

d=Gm，

（3）

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

2 深度加權

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

（4）

其中，W是一個對角矩陣，第i個對角元素為第i個模型單元的深度wi，β為參數，不同的β參數導致的反演結果不一樣，Oldenburg和Li（2005）認為反演問題可以選擇不同的β值進行嘗試，最后選出最符合地質意義上的解，這種加權方式在各種反演算法中得到了廣泛地應用.

3 反演

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

針對式（3）的求解，利用正則化方法，構造反演目標函數：

（5）

（6）

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

盡管前人在凸優化反演和非凸優化反演中展現了豐富的算法，但這些算法具有特異性，即不同的目標函數需要選擇不同的算法進行求解，所以尋找一種適應性強的算法在反演中顯得尤為必要.針對這個問題，本文將介紹ADMM算法在反演中的應用.

3.1 ADMM算法

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

minf（x）+g（x），

（7）

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

（8）

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

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

（9）

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

算法的迭代過程為：

（10）

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

（11）

上述過程表明，ADMM算法具有適應性很強的特點：凡具備式（7）結構，并且函數f、g為凸函數的目標函數優化問題，均可通過ADMM算法進行求解，并且可以保證其收斂性.這就為后續研究不同的反演目標函數提供了通用算法.當函數f、g中至少一個為非凸函數時，ADMM算法可以視為一種局部優化算法，其最終結果與起始點和參數有關.本文將從凸優化反演中的L1正則化和非凸優化中的Lp（0

3.2 基于ADMM算法的凸優化反演

當目標函數的正則項為L1范數時，利用ADMM算法框架將式（6）改寫為：

（12）

進一步，將式（11）轉化為極小化一系列子問題進行求解,迭代步驟：

（13）

原問題為：

（14）

屬于凸優化問題，對變量mw求導，并使其等于零，得到如下等式：

（15）

（16）

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

（17）

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

（2）zk+1的求解

原問題為：

（18）

根據式（18）可以看出，向量z的每一個分量都可以獨立求解，只需對向量z的分量zi（i=1…N）進行優化求解，優化問題寫為：

（19）

（20）

繼而對解析解進行閾值處理：

（21）

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

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

3.3 基于ADMM算法的非凸優化反演

當目標函數為Lp（0

當0

（22）

同理，對向量z的分量zi（i=1,…,N） 進行優化求解，化簡：

（23）

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

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

算法2:廣義軟閾值（GST）算法輸入:si,λ,p,J.輸出:TGSTp（si;λ）.1.計算τ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;循環k=0,1,…,J;zk+1i=si-λp（zk）p-1;k=k+1;TGSTp（si;λ）=sgn（si）zk;結束.

利用此方法，從而得到優化式（23）的解：

（24）

同式（21），在每次迭代后的結果進行閾值處理：

（25）

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

算法3:基于Lp（0

4 實驗

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

4.1 二維模擬實驗

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

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

從圖2—圖4中的模擬實驗可以看出，本文提出的基于ADMM算法的L1范數和Lp（0

針對基于Lp（0

（26）

為了進一步評價本文算法的精度，利用L1范數定義了反演結果的相對誤差，用E來表示：

（27）

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

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

從表1中可以發現，在直立板狀體模型和向斜模型中，基于Lp（p=0.4,p=0.1）范數反演效果優于基于L1范數反演，而基于Lp（p=0.7）范數反演與基于L1范數的反演效果相當；在傾斜板狀體模型中，基于L1范數反演與基于Lp（p=0.7,p=0.4,p=0.1）范數反演效果相當；在向斜模型中，基于Lp（p=0.7,p=0.4,p=0.1）范數反演效果均優于基于L1范數反演；此外，在這三種模型中，本文的兩種算法都遠遠優于基于L2范數反演算法效果.

圖2 直立板狀體模型反演結果（a） L1范數反演結果；（b） L1范數反演數據擬合；（c） Lp（p=0.7）范數反演結果；（d） Lp（p=0.7）范數反演數據擬合；（e） Lp（p=0.4）范數反演結果；（f） Lp（p=0.4）范數反演數據擬合；（g） Lp（p=0.1）范數反演結果；（h） Lp（p=0.1）范數反演數據擬合；（i） L2范數反演結果；（j） L2范數反演數據擬合.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 傾斜板狀體模型反演結果（a） L1范數反演結果；（b） L1范數反演數據擬合；（c） Lp（p=0.7）范數反演結果；（d） Lp（p=0.7）范數反演數據擬合；（e） Lp（p=0.4）范數反演結果；（f） Lp（p=0.4）范數反演數據擬合；（g） Lp（p=0.1）范數反演結果；（h） Lp（p=0.1）范數反演數據擬合；（i） L2范數反演結果；（j） L2范數反演數據擬合.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 向斜模型反演結果（a） L1范數反演結果；（b） L1范數反演數據擬合；（c） Lp（p=0.7）范數反演結果；（d）Lp（p=0.7）范數反演數據擬合；（e） Lp（p=0.4）范數反演結果；（f） Lp（p=0.4）范數反演數據擬合；（g） Lp（p=0.1）范數反演結果；（h） Lp（p=0.1）范數反演數據擬合；（i） L2范數反演結果；（j） L2范數反演數據擬合.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.

結合上述的表現說明ADMM算法可以應用于凸優化反演和非凸優化反演中.

圖5 極小化加權L1增強稀疏性（a） 稀疏向量m0，可行域d=Gm，和半徑為‖m0‖的范數球；（b） 存在m≠m0，使得‖m‖<‖m0‖；（c） 加權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算法在反演中的性能，以模擬實驗中直立板狀體模型為例分ADMM算法反演的精度（以模型的相對誤差來衡量）與迭代次數的關系.

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

4.2 二維真實數據實驗

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

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

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

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

圖8 196測線反演結果（a） L1范數反演結果；（b） L1范數反演數據擬合；（c） Lp（p=0.7）范數反演結果；（d） Lp（p=0.7）范數反演數據擬合；（e） Lp（p=0.4）范數反演結果；（f）Lp（p=0.4）范數反演數據擬合；（g） Lp（p=0.1）范數反演結果；（h） Lp（p=0.1）范數反演數據擬合.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測線反演結果（a） L1范數反演結果；（b） L1范數反演數據擬合；（c） Lp（p=0.7）范數反演結果；（d） Lp（p=0.7）范數反演數據擬合；（e） Lp（p=0.4）范數反演結果；（f） Lp（p=0.4）范數反演數據擬合；（g） Lp（p=0.1）范數反演結果；（h） Lp（p=0.1）范數反演數據擬合.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） 在反演過程中，當正則化參數較大時，正則項的權重相對較大，一方面會導致目標函數的誤差項過大，另一方面，會使得反演結果的模型單元呈現二值性（即0和40），但是模型單元的實際情況并非全部是均勻的，可能會導致擬合曲線局部的變化；當正則化參數較小的時候，目標函數的誤差平方項的權重相對較大，正則項的影響較小，目標函數誤差項較小，同時，反演結果也具備稀疏性，不會呈現出二值性的特點.

5 展望

ADMM算法在近幾年受到了較高的關注度，其主要原因除了上文提到的可以分離凸問題，降低目標函數的求解難度之外，另外一個原因是是可用于解決大規模優化問題.

考慮式（6）的一般化的反演問題，當核矩陣Gw過大時，改寫為：

（28）

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

（29）

式（28）被稱為一致性優化問題，利用ADMM算法將式（29）改寫為一系列子問題進行求解（Boyd et al.,2010），每一個局部變量（mw）i可以單獨求解，結合并行計算的技術,使局部變量（mw）i的求解同時進行，從而達到處理大規模問題的能力.

6 結論

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

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