程 三 張志勇* 周 峰② 李 曼 翟斌軍

(①東華理工大學(xué)地球物理與測控技術(shù)學(xué)院，江西南昌 330013； ②東華理工大學(xué)放射性地質(zhì)與勘探技術(shù)國防重點學(xué)科實驗室，江西南昌 330013)

0 引言

由于地球物理場的等效性、觀測數(shù)據(jù)集的有限性和數(shù)據(jù)觀測存在誤差等，地球物理反問題不可避免地存在多解性[1]。正則化技術(shù)[2]作為減小反問題多解性的有效策略之一，廣泛應(yīng)用于地球物理反演。該技術(shù)通過在數(shù)據(jù)和模型權(quán)重矩陣引入先驗信息，降低反問題的不適定性； 當(dāng)擁有足夠的先驗信息時，可有效提高反演的穩(wěn)定性。相反，當(dāng)先驗信息不足甚至錯誤時，將會增加產(chǎn)生錯誤結(jié)果的可能性； 降低反問題多解性的另一途徑是提高觀測數(shù)據(jù)集的數(shù)量和質(zhì)量[3]，但無法改變地球物理場等效性造成的反演多解性。此外，減少反演未知數(shù)數(shù)量是有效降低反問題不適定性的另一重要途徑[4]，但是以犧牲反演分辨率為代價； 另一方面，如果離散的反演網(wǎng)格無法真實構(gòu)建實際異常體，也可能會導(dǎo)致反演的多解性、影響反演分辨率[5]。高質(zhì)量的反演網(wǎng)格應(yīng)該在確保反演效果的同時減少反演單元數(shù)，進而從一定程度上降低反問題的不適定性。

提高地球物理反演離散網(wǎng)格質(zhì)量是一個亟需解決的問題，但是實際的反問題由于缺乏地下結(jié)構(gòu)邊界、尺寸和形狀等先驗信息，難以建立理想的初始反演網(wǎng)格[5]。雖然可以通過全局網(wǎng)格加密，實現(xiàn)對反演目標(biāo)體更好的逼近，但這無疑會過多地增加單元數(shù)量，對反演效率和穩(wěn)定性造成不利影響。B?hm等[6]提出了一種變網(wǎng)格反演方案，即先使用粗網(wǎng)格進行反演迭代，并將粗網(wǎng)格反演結(jié)果作為先驗信息，對粗網(wǎng)格進行加密與合并獲取更高質(zhì)量的反演網(wǎng)格，然后重新反演以提高反演分辨率。Ascher等[7]對B?hm等的方案進行了深入研究，發(fā)展了一種自適應(yīng)反演方案，并證明了粗網(wǎng)格選擇的正則化因子在細(xì)網(wǎng)格是可以延續(xù)使用的； 此外，交叉檢驗與L曲線正則化因子選取方法在自適應(yīng)反演中依然適用[8]。很多學(xué)者在地震層析成像[6，9-12]、圖像重建[13-15]、醫(yī)學(xué)成像[16]、電阻率成像[17]以及電磁法反演[18-19]、磁測數(shù)據(jù)反演[20-21]等方面開展了自適應(yīng)反演方案的研究。自適應(yīng)反演方案中，最重要的是選擇網(wǎng)格并對選定網(wǎng)格進行細(xì)化，理論上網(wǎng)格細(xì)化最有效的方案就是分析靈敏度矩陣的奇異值和海森矩陣的特征值，但這兩種方法的計算量過大，對于二維、三維反演問題難于實用。目前較為實用的網(wǎng)格細(xì)化標(biāo)準(zhǔn)有物性變化[17]、物性梯度變化[16]和物性相關(guān)的給定函數(shù)[22]等； 此外，基于小波技術(shù)的自適應(yīng)反演方案成功應(yīng)用于直流電阻率反演[23]。在自適應(yīng)反演過程中，當(dāng)網(wǎng)格細(xì)化過程中允許網(wǎng)格結(jié)點位置發(fā)生改變時，對于簡單模型的反演效果明顯[9]，但可能會產(chǎn)生局部過小的網(wǎng)格[12]，約束網(wǎng)格面積能在一定程度上避免此問題。

本文以二維MT反演為例，研究了自適應(yīng)網(wǎng)格細(xì)化反演算法。在反演初期使用粗網(wǎng)格，通過減少反演單元數(shù)的方式降低反問題的不適定性； 基于網(wǎng)格細(xì)化策略，隨反演迭代的進行，網(wǎng)格自適應(yīng)加密，以提高正則化反演效果。反演過程中，細(xì)化網(wǎng)格反演以前一重網(wǎng)格的反演結(jié)果作為參考模型與初始模型，確保反演沿正確的方向進行，同時提高反演穩(wěn)定性，逐步改善反演效果。采用四種網(wǎng)格細(xì)化方案進行自適應(yīng)反演，對比分析了四種網(wǎng)格細(xì)化方案的特點和反演效果，并通過Hessian矩陣的特征值分布分析了這四種方案對反演結(jié)果的影響。最后，開展了理論模型和實測數(shù)據(jù)試算，結(jié)果表明反演過程穩(wěn)定，自適應(yīng)反演算法具有較好的實用性。

1 漸近自適應(yīng)正則化反演

1.1 二維MT正演

對于大地電磁問題，考慮地下電阻率變化，取時諧因子為e-iωt。假設(shè)地下電性結(jié)構(gòu)為二維，取走向為z軸，電導(dǎo)率只在(x，y)平面變化，則走向方向的電場和磁場的偏微分方程可分別表示為

(1)

(2)

式中：σ為電導(dǎo)率；μ為磁導(dǎo)率；ε為介電常數(shù)；ω為角頻率； i為虛數(shù)單位；Ez為電場的z分量；Hz為磁場的z分量。采用非結(jié)構(gòu)化三角單元對研究區(qū)域進行離散，考慮空氣層的邊界條件，式(1)和式(2)對應(yīng)的變分問題可通過有限單元法[24]轉(zhuǎn)化為求解大型稀疏線性方程組

KU=P

(3)

式中：K為剛度矩陣；U為節(jié)點上的場值；P為由邊界條件形成的源項。通過求解式(3)的線性方程組可得主場分量，通過麥克斯韋方程可求得輔助場分量，最后可計算視電阻率和相位

(4)

(5)

式中：ρs為視電阻率；Z為阻抗；φ為相位； Im(·)表示取復(fù)數(shù)虛部； Re(·)表示取復(fù)數(shù)實部。

1.2 正則化反演目標(biāo)函數(shù)

建立二維MT正則化反演目標(biāo)函數(shù)，其中模型穩(wěn)定函數(shù)采用Zhdanov的一般表達(dá)式[25]，則有

Φ(m) =λΦm(m)+Φd(m)

(6)

式中：m為離散模型；Φd(m)為數(shù)據(jù)誤差函數(shù)；λ為正則化因子；Φm(m)為模型穩(wěn)定函數(shù)，本文采用最小結(jié)構(gòu)穩(wěn)定函數(shù)[26]，非結(jié)構(gòu)化網(wǎng)格模型粗糙度采用最小二乘算法計算[27]；Wm為模型權(quán)重矩陣；mapr為先驗?zāi)Ｐ?；Wd為數(shù)據(jù)權(quán)重矩陣；dobs為觀測數(shù)據(jù)；A(m)為與m相關(guān)的正演算子。自適應(yīng)反演過程中考慮到粗網(wǎng)格正演存在誤差，則Wd可表示為

(7)

式中：ξobs為觀測數(shù)據(jù)誤差；ξcal為正演計算誤差，通過網(wǎng)格細(xì)化后的正演數(shù)據(jù)評估得到；N為觀測數(shù)據(jù)個數(shù)。

1.3 目標(biāo)函數(shù)優(yōu)化

采用高斯—牛頓最優(yōu)化法求解目標(biāo)函數(shù)(式(6))，并基于迭代法求解高斯—牛頓方程。為了獲得模型參數(shù)，取mn=mn-1+Δm，其中n為迭代次數(shù)，對A(m)進行一階泰勒展開，則

(8)

對式(6)關(guān)于Δm求偏導(dǎo)并令其為0，可得高斯—牛頓迭代方程

Hn-1Δm=-?Φ(mn-1)

(9)

(10)

(11)

式中：J為靈敏度矩陣，表示正演算子對m的偏導(dǎo)數(shù)； Δm為模型改變量；H為海森矩陣。采用互換定理[28]計算J，采用穩(wěn)定雙共軛梯度法[29](BICGSTAB)求解式(9)。最后可得新的模型參數(shù)

mn=mn-1+αΔm

(12)

式中α為改進量的搜索步長[25]。

1.4 自適應(yīng)反演算法

自適應(yīng)反演算法在整個反演過程中采用可變網(wǎng)格。在反演初期使用粗網(wǎng)格，在反演迭代過程中基于網(wǎng)格優(yōu)化方案進行自適應(yīng)加密，細(xì)化網(wǎng)格的反演以前一重網(wǎng)格的反演結(jié)果作為初始模型和參考模型以確保反演的準(zhǔn)確性和穩(wěn)定性。自適應(yīng)反演流程(圖1)如下：

圖1 自適應(yīng)反演流程

(1)設(shè)置反演網(wǎng)格最大細(xì)化次數(shù)及自適應(yīng)反演中每一重網(wǎng)格最大迭代次數(shù)，根據(jù)反演需求設(shè)置每重網(wǎng)格細(xì)化的單元比例和最小單元面積約束等；

(2)采用非結(jié)構(gòu)化三角單元進行網(wǎng)格離散生成網(wǎng)格，為確保正演精度，在測量點周圍采用局部加密技術(shù)；

(3)對當(dāng)前網(wǎng)格，采用高斯—牛頓最優(yōu)化法求解目標(biāo)函數(shù)，直至達(dá)到當(dāng)前重網(wǎng)格的最大反演迭代次數(shù)；

(4)根據(jù)網(wǎng)格優(yōu)化方案計算單元優(yōu)化參考信息，并根據(jù)優(yōu)化比例選擇優(yōu)化單元，然后對網(wǎng)格自適應(yīng)加密；

(5)以細(xì)化后的網(wǎng)格作為下次反演的網(wǎng)格，設(shè)置反演的初始模型與參考模型，并返回第三步開始當(dāng)前網(wǎng)格反演，判斷是否達(dá)到反演中止條件，若達(dá)到，則終止反演。

應(yīng)用低阻體模型說明自適應(yīng)反演過程。在圍巖電阻率為100Ω·m均勻半空間，設(shè)置一個尺寸為400m×400m、電阻率為10Ω·m的低阻異常體，其頂部埋深為200m，在0.1～100Hz頻段以10為底的對數(shù)等間距取10個頻點，測點距為40m，共60個測點。在反演過程中基于模型梯度信息進行網(wǎng)格細(xì)化，每次細(xì)化比例為反演總網(wǎng)格數(shù)量的2%。圖2為前三次網(wǎng)格細(xì)化結(jié)果及基于模型梯度信息優(yōu)化方案的自適應(yīng)反演結(jié)果，每一重網(wǎng)格反演的初始正則化因子為200，迭代過程中以0.9的速率下降。為防止網(wǎng)格在某些區(qū)域過度細(xì)化，初始最小面積約束設(shè)為以測點間距為正方形的單元面積，并隨網(wǎng)格細(xì)化高于最小面積約束的單元面積以2為倍數(shù)增加。圖2中黑色正方形為實際異常體邊界，可見自適應(yīng)反演方案在反演初期使用粗網(wǎng)格，隨反演迭代進行，基于網(wǎng)格細(xì)化策略對網(wǎng)格自適應(yīng)加密，反演分辨率得到有效提高，反演目標(biāo)體的位置空間信息也更為清晰。

1.5 網(wǎng)格優(yōu)化策略

在自適應(yīng)反演過程中網(wǎng)格優(yōu)化策略對反演至關(guān)重要，不同網(wǎng)格優(yōu)化策略的網(wǎng)格加密位置有所差異，最后的反演效果會有所不同。為此，本文研究了四種網(wǎng)格優(yōu)化方案，對它們的網(wǎng)格加密特點進行了對比； 此外，為了分析四種優(yōu)化方案在網(wǎng)格加密過程中對反演的影響，對各自的Hessian矩陣進行SVD，并對歸一化特征值分布進行分析。四種網(wǎng)格優(yōu)化方案的參數(shù)分別為模型靈敏度矩陣J[29]、模型改變量Δmi、模型梯度信息η和“邊—角”檢測參數(shù)R[22]。J、Δmi、η、R具體為

(13)

(14)

(15)

R=det(M)-k[trace (M)]2

(16)

(17)

為說明不同網(wǎng)格優(yōu)化方案的加密特點，對上述低阻模型采用四種優(yōu)化方案進行自適應(yīng)反演，整個反演過程網(wǎng)格共細(xì)化4次，圖3為四種優(yōu)化方案的最后一重網(wǎng)格細(xì)化結(jié)果及反演結(jié)果。由圖3a和圖3d對比可知，基于“邊—角”檢測和梯度信息的網(wǎng)格優(yōu)化方案類似，主要在異常體邊界加密，能夠有效提高反演的效果。圖3b為基于靈敏度細(xì)化方案的反演結(jié)果，一般近地表的單元靈敏度較大，此方案會由近地表逐漸向深部加密，因此該方案具有對反演區(qū)域整體優(yōu)化的特點； 圖3c為基于模型改變量方案的反演結(jié)果，該方案加密區(qū)域集中于異常體內(nèi)部。

圖2 各重網(wǎng)格自適應(yīng)反演結(jié)果(a)第一重網(wǎng)格； (b)第二重網(wǎng)格； (c)第三重網(wǎng)格； (d)第四重網(wǎng)格

圖3 不同方案網(wǎng)格細(xì)化結(jié)果及反演結(jié)果對比(a)基于“邊—角”檢測； (b)基于靈敏度信息； (c)基于模型改變量； (d)基于梯度信息

為分析上述四種網(wǎng)格優(yōu)化方案對反演的影響，對Hessian矩陣的歸一化特征值分布進行分析，在不考慮先驗信息的情況下Hessian矩陣可近似表示為JTJ。在反問題中，Hessian矩陣的特征值分布分為主成分特征值(接近0)和離群值兩部分，其中主成分特征值與網(wǎng)格尺寸無關(guān)，網(wǎng)格數(shù)量增加只會提高主成分特征值的數(shù)量[31-33]； 此外，線性問題的病態(tài)特征與Hessian矩陣特征值衰減為0的速度有關(guān)，在相同條件下衰減速率越快說明不適定性越強[31-32，34]； Hessian矩陣特征值分布有時也會出現(xiàn)極小值，與主成分特征值相差數(shù)個數(shù)量級，這種情況被認(rèn)為是由于零空間的存在引起的[35-36]。因此本文對自適應(yīng)反演中每一重網(wǎng)格最后一次迭代的Hessian矩陣的特征值進行歸一化分析，如圖4所示，可見整個特征值分布出現(xiàn)兩部分極值，即極大和極小特征值部分，中間較平滑部分可認(rèn)為是主成分特征值。對于基于模型靈敏度的網(wǎng)格優(yōu)化方案，隨網(wǎng)格的增加，主成分特征值的衰減速率基本沒有變化，說明此方案對反演的不適定性影響較小，因此不僅可以用于進行自適應(yīng)反演，也可以用于生成高質(zhì)量的初始反演網(wǎng)格； 對于基于模型梯度信息和“邊—角”檢測的網(wǎng)格優(yōu)化方案，在第一次網(wǎng)格加密后特征值衰減速率基本不變，但隨著網(wǎng)格數(shù)持續(xù)增加特征值衰減速率增加，意味著反演的不適定性增強； 對于基于模型改變量的網(wǎng)格優(yōu)化方案，在第一次網(wǎng)格加密后特征值衰減速率基本不變，但隨著網(wǎng)格數(shù)繼續(xù)增加，特征值衰減速率明顯增加，相比其他方案反演不適定性增加速率更快。

由圖3和圖4可知，基于靈敏度信息加密方式具有全局優(yōu)化加密的特點； 而其他三種網(wǎng)格優(yōu)化方案，網(wǎng)格加密均直接與反演異常體分布區(qū)域相關(guān)，在異常體分布區(qū)域進行網(wǎng)格適量加密，對異常體逼近程度提高時，反問題的不適定性受到的影響不大，但持續(xù)增加后反問題的不適定性會快速增加，模型改變量尤為明顯?？傮w上看，反演單元數(shù)越大，反問題的不適定性越強，在自適應(yīng)反演過程中應(yīng)盡量避免冗余網(wǎng)格的產(chǎn)生。

圖4 不同網(wǎng)格細(xì)化方案的自適應(yīng)反演的Hessian矩陣特征值分布(a)基于“邊—角”檢測； (b)基于靈敏度信息； (c)基于模型改變量； (d)基于梯度信息

2 模型試算

2.1 水平地形模型

在背景電阻率為100Ω·m的均勻半空間下設(shè)置三個異常體，其中兩個為正方形低阻異常體，尺寸為400m×400m，電阻率為10Ω·m，頂面深度分別為200m和400m； 另一個為矩形高阻異常體，尺寸為600m×400m，頂面深度為200m，電阻率為1000Ω·m。異常體的實際位置在反演結(jié)果中以黑色矩形框標(biāo)識。在0.1～100Hz頻段以10為底的對數(shù)等間距取10個頻點，測點距為20m，共540個測點，反演區(qū)域為測線長度與地下2000m形成的矩形區(qū)。每一重網(wǎng)格初始正則化因子均為500，以0.95的速率下降，整個反演過程中網(wǎng)格共細(xì)化四次，每次優(yōu)化比例為反演總網(wǎng)格的2%，基于四種網(wǎng)格優(yōu)化策略進行反演； 為避免某些區(qū)域過度細(xì)化，初始最小面積約束設(shè)為以測點間距為正方形的單元面積，并隨網(wǎng)格細(xì)化，以2為倍數(shù)增加。圖5為四種網(wǎng)格優(yōu)化方案反演結(jié)果。

為驗證自適應(yīng)反演算法的反演效果，采用固定網(wǎng)格的方案進行反演，并與自適應(yīng)反演的效果進行對比。在固定網(wǎng)格的反演中，整個反演區(qū)域采用均勻網(wǎng)格離散，最大單元面積為5000m2，反演單元數(shù)為18100，反演過程中采用經(jīng)驗選取的方法選擇正則化因子，初始正則化因子為500。整個反演過程中共進行了34次迭代，圖6為固定網(wǎng)格反演結(jié)果。

圖7a為五種反演方案的數(shù)據(jù)均方根(RMS)誤差曲線，圖7b為模型參數(shù)誤差曲線； 表1為四種方案的網(wǎng)格變化。由圖7可知，四種優(yōu)化方案的自適應(yīng)反演，初始RMS均快速下降，整個自適應(yīng)反演過程中RMS誤差均呈穩(wěn)定下降趨勢并逐漸趨近于1，以及模型參數(shù)誤差的穩(wěn)定上升說明反演過程穩(wěn)定[25]； 固定網(wǎng)格方案反演RMS誤差呈穩(wěn)定下降趨勢并逐漸趨近于1，說明反演過程穩(wěn)定，但下降速率較自適應(yīng)反演略慢，說明粗網(wǎng)格較細(xì)網(wǎng)格的反演迭代收斂要快。根據(jù)表1和圖5可知，四種優(yōu)化方案最后的網(wǎng)格數(shù)量接近，網(wǎng)格數(shù)量較初始網(wǎng)格增加了約六分之一，四種網(wǎng)格細(xì)化策略都能較好地反演出異常體，其中低阻物性接近真值，證明了四種網(wǎng)格加密方案在自適應(yīng)反演中的有效性。基于“邊—角”檢測技術(shù)與物性梯度信息的網(wǎng)格加密方式都以物性兩個方向的梯度變化為基礎(chǔ)，兩者網(wǎng)格加密方式較為相似，反演結(jié)果相近； 基于模型靈敏度的加密方式，反演區(qū)域整體由上至下逐漸加密，這是由于地表觀測的地球物理方法模型靈敏度由淺到深逐漸減小； 基于模型改變量的網(wǎng)格加密能夠有效在高、低阻異常區(qū)進行網(wǎng)格加密，異常體周圍網(wǎng)格較大，能夠有效反演出異常體。

由圖5a、圖5c、圖5d和圖6可知，基于“邊—角”檢測技術(shù)、物性梯度信息和模型改變量的局部網(wǎng)格加密方案對頂面深度為200m的低阻體邊界描述較固定網(wǎng)格方案略好，頂面深度為400m的低阻體物性恢復(fù)得更好，且反演單元數(shù)約為固定網(wǎng)格的77%。說明基于局部網(wǎng)格加密的自適應(yīng)反演方案能在保證反演效果的同時減小反演單元數(shù)，從而減少內(nèi)存消耗。

圖5 四種細(xì)化方案的自適應(yīng)反演結(jié)果(a)基于“邊—角”檢測； (b)基于靈敏度信息； (c)基于模型改變量； (d)基于梯度信息

圖6 固定網(wǎng)格反演結(jié)果

圖7 水平地形模型五種方案的反演參數(shù)隨迭代次數(shù)的變化曲線(a)數(shù)據(jù)的RMS誤差； (b)模型參數(shù)誤差

表1 水平地形模型四種自適應(yīng)加密方案的單元數(shù)

2.2 起伏地形模型

為測試自適應(yīng)反演方案在帶地形反演中的效果，在起伏地形下加入三個異常體，其中兩個低阻異常體，尺寸為400m×400m，頂面深度為400m，電阻率為10Ω·m； 一個高阻異常體尺寸為400m×400m，頂面深度為400m，電阻率為1000Ω·m，取地下背景電阻率為100Ω·m。異常體的實際位置在反演結(jié)果中以黑色矩形框標(biāo)識。在0.1～100Hz頻段以10為底的對數(shù)等間距取15個頻點，測點距為20m，共240個測點。每一重網(wǎng)格初始正則化因子為300，以0.95的速率下降，整個反演過程中網(wǎng)格共細(xì)化四次，每次優(yōu)化比例為反演總網(wǎng)格的3%。限于篇幅只討論基于模型梯度信息和模型改變量兩種網(wǎng)格優(yōu)化方案的自適應(yīng)反演。

圖8為基于模型梯度信息和模型改變量網(wǎng)格優(yōu)化方案的反演結(jié)果，圖9為數(shù)據(jù)的RMS誤差和模型參數(shù)誤差曲線，表2為兩種方案的網(wǎng)格變化。由圖9可知，兩種優(yōu)化方案反演初始RMS快速下降，整個自適應(yīng)反演過程中RMS總體呈下降的趨勢并逐漸趨近于1，說明反演過程穩(wěn)定。根據(jù)表2和圖8可知：兩種優(yōu)化方案最后的網(wǎng)格數(shù)量接近，網(wǎng)格數(shù)量較初始網(wǎng)格增加了約四分之一，兩種方案都能較好地反演出目標(biāo)體，其中低阻物性接近真值，且與真實位置相一致； 高阻物體受地形影響，反演位置與實際位置略有偏移； 此外，基于模型改變量網(wǎng)格加密方案在高阻物體周圍網(wǎng)格加密相對較少，網(wǎng)格明顯偏大，這與模型改變量優(yōu)先細(xì)化異常明顯區(qū)域有關(guān)。模型梯度信息的方案可有效避免這種情況。

圖8 起伏地形模型兩種方案的自適應(yīng)反演結(jié)果(a)基于梯度信息； (b)基于模型改變量

圖9 起伏地形模型兩種方案的反演參數(shù)隨迭代變化曲線(a)數(shù)據(jù)的RMS誤差； (b)模型參數(shù)誤差

表2 起伏地形模型兩種方案的自適應(yīng)反演單元數(shù)

3 實測數(shù)據(jù)反演

選擇A花崗巖區(qū)一條MT實測剖面驗證自適應(yīng)反演方案能力。該測線長約為7km，點距約為50m，采用四種網(wǎng)格細(xì)化方案進行自適應(yīng)反演。初始正則化因子設(shè)為4000，以0.9的速率下降。每次網(wǎng)格優(yōu)化比例為反演總網(wǎng)格的5%，采用相同的網(wǎng)格面積約束方案。圖10為四種網(wǎng)格細(xì)化方案反演結(jié)果，圖11為數(shù)據(jù)RMS誤差和模型參數(shù)誤差的變化曲線，表3為四種方案的網(wǎng)格變化。由表3可知最后四種網(wǎng)格優(yōu)化方案的網(wǎng)格數(shù)量接近，由圖11可知四種方案的反演參數(shù)變化曲線大致相同，RMS誤差平穩(wěn)下降，模型誤差單調(diào)上升，說明算法具有較好的穩(wěn)定性。四種網(wǎng)格優(yōu)化自適應(yīng)反演都能識別出橫向0～2.0km和4.5～6.5km的地下高阻地質(zhì)體，且左側(cè)的電阻率更大。基于“邊—角”檢測和模型梯度信息網(wǎng)格細(xì)化方案能夠清晰地刻畫出地質(zhì)體邊界； 基于模型靈敏度的網(wǎng)格加密方案，地質(zhì)體邊界網(wǎng)格較大，邊界分辨率略低； 基于模型改變量的網(wǎng)格加密方案在地質(zhì)體內(nèi)部明顯加密，但目標(biāo)體邊界網(wǎng)格明顯大于梯度信息和“邊—角”檢測方案，對邊界的描述與基于模型靈敏度的加密方案相當(dāng)。此外，基于“邊—角”檢測的反演方案能夠?qū)蓚€方向的模型梯度信息進行重新平衡，減小異常值之間的差異，避免網(wǎng)格在異常大的地質(zhì)體單元周圍過度加密。

圖10 實測數(shù)據(jù)四種方案自適應(yīng)反演結(jié)果(a)基于“邊—角”檢測； (b)基于靈敏度信息； (c)基于模型改變量； (d)基于梯度信息

圖11 實際數(shù)據(jù)兩種方案的反演參數(shù)隨迭代變化曲線(a)數(shù)據(jù)的RMS誤差； (b)模型參數(shù)誤差

表3 實際數(shù)據(jù)四種方案自適應(yīng)反演的網(wǎng)格單元數(shù)

綜上所述，實際的地質(zhì)情況往往較為復(fù)雜，且可能存在較大的地質(zhì)體，因此在實測數(shù)據(jù)反演中基于模型梯度信息和“邊—角”檢測網(wǎng)格優(yōu)化的自適應(yīng)反演在刻畫地質(zhì)體邊界方面更具優(yōu)勢。

4 結(jié)論

本文研究了自適應(yīng)漸進網(wǎng)格細(xì)化反演算法，并應(yīng)用于二維MT反演，得到以下結(jié)論。

(1)自適應(yīng)反演方案，在反演初期使用粗網(wǎng)格，通過減小反演單元數(shù)能夠有效降低反問題的不適定性； 隨反演迭代進行，基于網(wǎng)格優(yōu)化方案自適應(yīng)加密網(wǎng)格，能夠逐步改善反演效果。

(2)對基于模型靈敏度、模型改變量、模型梯度、“邊—角”檢測的四種網(wǎng)格細(xì)化方案形成的Hessian矩陣特征值分布研究表明：總體上，反問題不適定性會隨網(wǎng)格數(shù)量增加而提高，基于模型靈敏度的網(wǎng)格優(yōu)化方案，隨網(wǎng)格的增加對反演不適定性的影響相比其他三種方法小。

(3)四種網(wǎng)格優(yōu)化方案有三種主要特點，其中基于模型梯度方法與“邊—角”檢測的特征相似，集中在異常體邊界加密，能夠有效提高對異常體邊界刻畫的效果； 基于模型改變量方案主要對異常體分布區(qū)進行加密，對異常不明顯區(qū)域識別較差； 基于模型靈敏度方案具有一定的全局加密特點。

在反演中，可根據(jù)網(wǎng)格優(yōu)化特點綜合應(yīng)用多種方案進一步改善反演效果。基于靈敏度的方法可以用于反演的初期，以提高網(wǎng)格總體質(zhì)量甚至可用于初始網(wǎng)格生成，基于梯度與“邊—角”檢測的方案有利于得到精確的異常邊界，基于模型改變量的方案可以用于中期以提高異常體響應(yīng)精度。以后將進一步研究多種優(yōu)化方案的組合技術(shù)。此外，本文只采用了由粗到細(xì)的優(yōu)化方案，由細(xì)到粗的反向優(yōu)化策略是進一步的研究方向。