非結構化網格漸進自適應正則化反演算法

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

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

0 引言

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

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

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

1 漸近自適應正則化反演

1.1 二維MT正演

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

(1)

(2)

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

KU=P

(3)

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

(4)

(5)

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

1.2 正則化反演目標函數

建立二維MT正則化反演目標函數，其中模型穩定函數采用Zhdanov的一般表達式[25]，則有

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

(6)

式中：m為離散模型；Φd(m)為數據誤差函數；λ為正則化因子；Φm(m)為模型穩定函數，本文采用最小結構穩定函數[26]，非結構化網格模型粗糙度采用最小二乘算法計算[27]；Wm為模型權重矩陣；mapr為先驗模型；Wd為數據權重矩陣；dobs為觀測數據；A(m)為與m相關的正演算子。自適應反演過程中考慮到粗網格正演存在誤差，則Wd可表示為

(7)

式中：ξobs為觀測數據誤差；ξcal為正演計算誤差，通過網格細化后的正演數據評估得到；N為觀測數據個數。

1.3 目標函數優化

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

(8)

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

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

(9)

(10)

(11)

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

mn=mn-1+αΔm

(12)

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

1.4 自適應反演算法

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

圖1 自適應反演流程

(1)設置反演網格最大細化次數及自適應反演中每一重網格最大迭代次數，根據反演需求設置每重網格細化的單元比例和最小單元面積約束等；

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

(3)對當前網格，采用高斯—牛頓最優化法求解目標函數，直至達到當前重網格的最大反演迭代次數；

(4)根據網格優化方案計算單元優化參考信息，并根據優化比例選擇優化單元，然后對網格自適應加密；

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

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

1.5 網格優化策略

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

(13)

(14)

(15)

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

(16)

(17)

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

圖2 各重網格自適應反演結果(a)第一重網格； (b)第二重網格； (c)第三重網格； (d)第四重網格

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

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

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

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

2 模型試算

2.1 水平地形模型

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

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

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

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

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

圖6 固定網格反演結果

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

表1 水平地形模型四種自適應加密方案的單元數

2.2 起伏地形模型

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

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

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

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

表2 起伏地形模型兩種方案的自適應反演單元數

3 實測數據反演

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

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

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

表3 實際數據四種方案自適應反演的網格單元數

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

4 結論

本文研究了自適應漸進網格細化反演算法，并應用于二維MT反演，得到以下結論。

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

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

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

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