陳 強 葉險峰

（湘潭大學 土木工程與力學學院，湖南 湘潭 411105）

反演具有多解性，單一的反演方法很難反演出準確且清晰的地下模型，Vozoff 和Jupp（1975）[1]首次提出聯(lián)合反演的概念，通過聯(lián)合反演來耦合不同的物性信息。在此之后聯(lián)合反演經(jīng)過不斷的發(fā)展，Gallardo 和Meju（2003、2004）[2-3]基于不同的物性參數(shù)具有同一地質(zhì)結構的特性，提出了交叉梯度聯(lián)合反演方法。在Bezdek（1974）[4]提出的模糊c-均值聚類算法的基礎上，研究者們提出了模糊c-均值（FCM）聚類聯(lián)合反演算法，后來Sun和Li（2015，2016）[5-6]又做出了一些改進，加入了先驗地質(zhì)物性信息，隨后Sun 和Li（2016）[7]又提出了模糊c-回歸（FCRM）聚類方法。Sun 和Li 在FCRM 聚類反演中使用高斯牛頓法解方程組，計算效率低，所以本文采用共軛梯度法實現(xiàn)FCRM 聚類反演。共軛梯度法的概念是Hestenes[8]與Stiefel[9]于1950 年首次提出的，由于收斂速度快、內(nèi)存需求小，在地球物理反演中被廣泛采用[10]。近年來，預條件共軛梯度法[11]和重加權共軛梯度法[12]也被廣泛應用于反演中。

本文在模型目標函數(shù)的構建中采用了光滑約束，使用有限差分法計算模型加權矩陣，所以，在求解大型線性方程組的計算中，核函數(shù)和大型稀疏矩陣的存儲與計算是一個不容忽視的問題。本文在涉及核函數(shù)矩陣的運算中采用了矩陣分解算法[13]，用三元組的形式表示大型二維稀疏矩陣。除此之外，本文在編程時用OpenMP 指令實現(xiàn)CPU 多核環(huán)境下的多線程計算，用于提高計算效率。

1 反演理論

重磁反演的基礎是正演基本原理，本文使用的是一般直立六面體的正演計算公式[14]。

1.1 光滑約束

1.2 模糊c-回歸聚類聯(lián)合反演

聯(lián)合反演可以有效的耦合多種物性信息，降低反演多解性。模糊聚類聯(lián)合反演的目標函數(shù)如下：預期的線性聚類的距離，引導反演結果的線性聚類靠近預期的線性聚類，達到預期的地質(zhì)分化效果。

1.3 重加權共軛梯度法

1.4 核函數(shù)矩陣的分解運算

1.5 大型稀疏矩陣的優(yōu)化存儲與運算

2 模型試驗

2.1 簡單模型

加入所述的約束條件后，反演結果如圖1（c）（g）所示，根據(jù)它們可以看出傳統(tǒng)約束反演的結果雖然顯示出了中心直立長方體的位置，但是它的邊界是發(fā)散的，模型和圍巖的邊界模糊不清。相比于磁力反演的結果，重力反演的結果很明顯地下沉了，這是由于假設的地下空間的底層圍巖仍有剩余密度，這就相當于存在一個擁有小剩余密度的模型，深度分辨率較差的重力反演的結果就出現(xiàn)了下沉的情況。

圖1（d）（h）是模糊c-回歸聚類聯(lián)合反演得到的地下空間的密度分布圖和磁化強度分布圖，相比于單獨反演的結果，聯(lián)合反演加強了兩種反演結果的共同部分，削弱了它們之間的差異部分，在圖中呈現(xiàn)出的現(xiàn)象就是直立長方體所在的位置的密度值和磁化強度值更為密集，且高于圍巖。而密度模型中底層的物性聚類是傳統(tǒng)的約束反演得到的密度模型的向下發(fā)散導致，在底層位置，傳統(tǒng)約束反演得到的密度模型和磁化強度模型在底層位置的物性產(chǎn)生了線性關系，這被模糊c-回歸聚類聯(lián)合反演誤識別為模型間物性的線性關系，所以使密度模型的下沉部分加強了。

圖1

2.2 復雜模型

圖2（a）（b）是根據(jù)湖南省有色地質(zhì)勘查研究院提供的坪寶-黃沙坪鉛鋅銀多金屬礦某處的地下剖面建模建立的復雜剖面模型。為了進行三維的反演試驗，對該剖面模型進行了理想化的假設：XZ 剖面模型在Y 方向進行復制，獲得多層完全相同的XZ 剖面模型。為了減少數(shù)據(jù)量，減少三維的反演試驗的時間，我們在Y 方向只設置了5 層。

圖2（c）（g）是傳統(tǒng)約束反演的結果，圖2（d）（h）是模糊c-回歸聚類聯(lián)合反演得到的結果，可以看出聯(lián)合反演加強了密度模型和磁化強度模型之間共同的部分，削弱了它們之間差異的部分，改善了磁化強度模型的物性分布。雖然三維反演的多解性問題仍然嚴重，約束條件遠遠少于未知條件，但是總體來說，F(xiàn)CRM 聚類聯(lián)合反演方法相比于傳統(tǒng)約束反演方法是有所改進的。

圖2

3 結論

針對重、磁反演問題的趨膚現(xiàn)象，本文采用了深度加權函數(shù)抵消核函數(shù)隨深度幾何衰減的影響；針對重、磁反演的多解性，本文采用了光滑約束、正定約束等多種約束；針對大型線性方程組的計算難度問題，本文對核函數(shù)矩陣進行矩陣分解計算，采用三元組實現(xiàn)大型稀疏矩陣的存儲和計算，采用基于PRP 公式的重加權共軛梯度最優(yōu)化方法，同時采用OpenMP 并行計算技術實現(xiàn)CPU 多線程加速，最后，采用模糊c-回歸聚類算法實現(xiàn)重、磁數(shù)據(jù)的聯(lián)合反演，并先后采用簡單直立長方體模型和復雜模型進行驗證，得出以下結論：

3.1 矩陣分解方法和三元組數(shù)組的引入改善了大型二維數(shù)組的存儲內(nèi)存大和計算時間長的問題，特別是在對兩個二維數(shù)組間的計算中，這種改善尤為明顯，達到了減少算法運行時間的目的；

3.2 光滑約束、正定約束等方法的引入改善了重、磁物性反演的多解性問題，深度約束和重加權反演方法很好的解決了反演的趨膚效應；

3.3 根據(jù)直立長方體模型和復雜模型的試驗，對比重、磁單獨反演和聯(lián)合反演的結果，重、磁數(shù)據(jù)的聯(lián)合反演加強了兩種物性數(shù)據(jù)中一致的特征，同時又對差異的部分進行了互補和改善，聯(lián)合反演得到的結果相比單獨反演的結果更加真實準確。

附錄A

糊聚類聯(lián)合反演的目標函數(shù)如下：

為了使目標函數(shù)在線性聚類上最小化，我們求出（A-1）式對線性聚類參數(shù)vk1的導數(shù)，然后令導數(shù)方程等于0，求解得到線性聚類參數(shù)（斜率）的計算公式：

同理，線性聚類參數(shù)vk2（截距）的計算公式如下：

求出（A-1）式對隸屬度函數(shù)ujk的導數(shù)，然后令導數(shù)方程等于0，求解得到隸屬度函數(shù)ujk的計算公式：