非結(jié)構(gòu)化網(wǎng)格漸進(jìn)自適應(yīng)正則化反演算法

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

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

0 引言

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

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

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

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

1.1 二維MT正演

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

(1)

(2)

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

KU=P

(3)

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

(4)

(5)

式中：ρs為視電阻率；Z為阻抗；φ為相位； Im(·)表示取復(fù)數(shù)虛部； Re(·)表示取復(fù)數(shù)實(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)格模型粗糙度采用最小二乘算法計(jì)算[27]；Wm為模型權(quán)重矩陣；mapr為先驗(yàn)?zāi)Ｐ?；Wd為數(shù)據(jù)權(quán)重矩陣；dobs為觀測(cè)數(shù)據(jù)；A(m)為與m相關(guān)的正演算子。自適應(yīng)反演過(guò)程中考慮到粗網(wǎng)格正演存在誤差，則Wd可表示為

(7)

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

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

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

(8)

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

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

(9)

(10)

(11)

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

mn=mn-1+αΔm

(12)

式中α為改進(jìn)量的搜索步長(zhǎng)[25]。

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

自適應(yīng)反演算法在整個(gè)反演過(guò)程中采用可變網(wǎng)格。在反演初期使用粗網(wǎng)格，在反演迭代過(guò)程中基于網(wǎng)格優(yōu)化方案進(jìn)行自適應(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)化三角單元進(jìn)行網(wǎng)格離散生成網(wǎng)格，為確保正演精度，在測(cè)量點(diǎn)周?chē)捎镁植考用芗夹g(shù)；

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

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

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

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

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

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

(13)

(14)

(15)

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

(16)

(17)

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

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

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

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

2 模型試算

2.1 水平地形模型

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

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

圖7a為五種反演方案的數(shù)據(jù)均方根(RMS)誤差曲線，圖7b為模型參數(shù)誤差曲線； 表1為四種方案的網(wǎng)格變化。由圖7可知，四種優(yōu)化方案的自適應(yīng)反演，初始RMS均快速下降，整個(gè)自適應(yīng)反演過(guò)程中RMS誤差均呈穩(wěn)定下降趨勢(shì)并逐漸趨近于1，以及模型參數(shù)誤差的穩(wěn)定上升說(shuō)明反演過(guò)程穩(wěn)定[25]； 固定網(wǎng)格方案反演RMS誤差呈穩(wěn)定下降趨勢(shì)并逐漸趨近于1，說(shuō)明反演過(guò)程穩(wěn)定，但下降速率較自適應(yīng)反演略慢，說(shuō)明粗網(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)反演中的有效性?；凇斑叀恰睓z測(cè)技術(shù)與物性梯度信息的網(wǎng)格加密方式都以物性兩個(gè)方向的梯度變化為基礎(chǔ)，兩者網(wǎng)格加密方式較為相似，反演結(jié)果相近； 基于模型靈敏度的加密方式，反演區(qū)域整體由上至下逐漸加密，這是由于地表觀測(cè)的地球物理方法模型靈敏度由淺到深逐漸減??； 基于模型改變量的網(wǎng)格加密能夠有效在高、低阻異常區(qū)進(jìn)行網(wǎng)格加密，異常體周?chē)W(wǎng)格較大，能夠有效反演出異常體。

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

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

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

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

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

2.2 起伏地形模型

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

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

選擇A花崗巖區(qū)一條MT實(shí)測(cè)剖面驗(yàn)證自適應(yīng)反演方案能力。該測(cè)線長(zhǎng)約為7km，點(diǎn)距約為50m，采用四種網(wǎng)格細(xì)化方案進(jìn)行自適應(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)上升，說(shuō)明算法具有較好的穩(wěn)定性。四種網(wǎng)格優(yōu)化自適應(yīng)反演都能識(shí)別出橫向0～2.0km和4.5～6.5km的地下高阻地質(zhì)體，且左側(cè)的電阻率更大?；凇斑叀恰睓z測(cè)和模型梯度信息網(wǎng)格細(xì)化方案能夠清晰地刻畫(huà)出地質(zhì)體邊界； 基于模型靈敏度的網(wǎng)格加密方案，地質(zhì)體邊界網(wǎng)格較大，邊界分辨率略低； 基于模型改變量的網(wǎng)格加密方案在地質(zhì)體內(nèi)部明顯加密，但目標(biāo)體邊界網(wǎng)格明顯大于梯度信息和“邊—角”檢測(cè)方案，對(duì)邊界的描述與基于模型靈敏度的加密方案相當(dāng)。此外，基于“邊—角”檢測(cè)的反演方案能夠?qū)蓚€(gè)方向的模型梯度信息進(jìn)行重新平衡，減小異常值之間的差異，避免網(wǎng)格在異常大的地質(zhì)體單元周?chē)^(guò)度加密。

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

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

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

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

4 結(jié)論

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

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

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

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

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