重磁三維數(shù)據(jù)的模糊聚類聯(lián)合快速反演

陳 強(qiáng) 葉險(xiǎn)峰

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

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

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

1 反演理論

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

1.1 光滑約束

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

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

1.3 重加權(quán)共軛梯度法

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

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

2 模型試驗(yàn)

2.1 簡(jiǎn)單模型

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

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

圖1

2.2 復(fù)雜模型

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

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

圖2

3 結(jié)論

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

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

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

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

附錄A

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

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

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

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