翁愛(ài)華，郭峻豪

(吉林大學(xué) 地球探測(cè)科學(xué)與技術(shù)學(xué)院，吉林 長(zhǎng)春 130026)

0 引 言

蘭州地區(qū)地處歐亞大陸中部，其西部可能存在塔里木地幔柱[1-2]，而南部可能存在峨眉地幔柱[3]。最新的地震成像發(fā)現(xiàn)歐亞大陸中北部下地幔中存在低速異常[4]，動(dòng)力學(xué)模擬推測(cè)其可能是太平洋深部下地幔中大低速體的一個(gè)分支[5]。地幔柱的存在將導(dǎo)致局部高溫異常，在轉(zhuǎn)換帶中引起電導(dǎo)率的增高。而蘭州地區(qū)位于上述可能的地幔柱的中間位置，其深部電導(dǎo)率可能對(duì)這些地幔柱產(chǎn)生響應(yīng)。因此，研究蘭州地區(qū)地幔轉(zhuǎn)換帶的電導(dǎo)率模型，能夠從電性結(jié)構(gòu)角度間接給出地幔柱可能存在的相關(guān)信息。

地球深部，尤其是在410～1 600 km深度的巖石，其電導(dǎo)率主要依靠地磁測(cè)深C-響應(yīng)的變化進(jìn)行有效表達(dá)[6]。地磁測(cè)深C-響應(yīng)的估計(jì)使用地磁臺(tái)站記錄的水平磁場(chǎng)和垂直磁場(chǎng)數(shù)據(jù)[7-8]，它們受到地核中的導(dǎo)電對(duì)流流體以及地幔、地殼、電離層、磁層和洋流體系的影響而發(fā)生磁場(chǎng)強(qiáng)度隨時(shí)間的改變[9-10]。因此，通過(guò)研究地表磁場(chǎng)的變化，可以獲得地球內(nèi)部電導(dǎo)率，進(jìn)而研究地球內(nèi)部與巖石電性特征密切相關(guān)的參數(shù)和地球內(nèi)部動(dòng)力學(xué)特征[11-12]。

Kelbert等采用非線(xiàn)性共軛梯度法[13]，反演獲得地幔的三維電性結(jié)構(gòu)[14]。Kuvshinov等采用有限內(nèi)存擬牛頓法進(jìn)行地磁測(cè)深反演[15]，發(fā)現(xiàn)歐洲、非洲地區(qū)地幔轉(zhuǎn)換帶具有高阻特征，而中國(guó)東北地區(qū)地幔轉(zhuǎn)換帶具有低阻特征[16]。Munch等采用隨機(jī)優(yōu)化和模型探索技術(shù)進(jìn)行地磁測(cè)深一維反演，結(jié)合溫度、壓力和含水量實(shí)驗(yàn)數(shù)據(jù)，獲得歐洲、非洲等地區(qū)地幔轉(zhuǎn)換帶的溫度信息[17]。張艷輝等基于全局光滑約束技術(shù)將得到的地磁數(shù)據(jù)進(jìn)行了光滑模型一維反演，獲得中國(guó)東部地區(qū)地幔轉(zhuǎn)換帶的含水量分布[18]。這些不同反演方法獲得的結(jié)果既有相似性，也有互異性。

結(jié)果上的差異有可能是反演方法或反演時(shí)初始模型不同造成的，而非啟發(fā)式尋優(yōu)算法應(yīng)用在地球物理反演工作中能很好地?cái)[脫初始模型對(duì)反演結(jié)果的影響[19-21]，另外其尋得的最優(yōu)解是在整個(gè)約束范圍搜索，可有效避免解陷入局部最小的問(wèn)題。因此，本文首次嘗試將模擬退火(Simulated Annealing,SA)非啟發(fā)式算法應(yīng)用于地磁測(cè)深反演工作中，為研究地球深部結(jié)構(gòu)提供新的參考依據(jù)[22-24]。

模擬退火全局優(yōu)化算法目前已經(jīng)被應(yīng)用于地球物理反演中[25-27]。Ingber等針對(duì)模擬退火算法的效率問(wèn)題，提出了非?？焖倌M算法(Very Fast Simulated Annealing，VFSA)[28-29]。師學(xué)明等應(yīng)用模擬退火算法進(jìn)行一維大地電磁測(cè)深反演，獲得了與地震剖面相一致的解釋結(jié)果[30]。Sharma用模擬退火算法對(duì)印度東部甘加盆地進(jìn)行一維直流電阻率測(cè)深反演，獲得準(zhǔn)確的地下含水層信息[31]。據(jù)此可以看到，模擬退火算法在電磁探測(cè)中都取得了較好的應(yīng)用效果。

模擬退火算法反演一般用最后一次迭代的模型作為反演結(jié)果，而本文的不同之處在于對(duì)所有有效的隨機(jī)反演模型進(jìn)行統(tǒng)計(jì)分析，將統(tǒng)計(jì)模型作為地磁測(cè)深數(shù)據(jù)一維隨機(jī)反演結(jié)果。為此，本文首先簡(jiǎn)要介紹地磁測(cè)深的基本原理；接著給出結(jié)合統(tǒng)計(jì)分析的隨機(jī)反演技術(shù)；在此基礎(chǔ)上，利用理論模型討論影響反演結(jié)果的因素和效果；最后，對(duì)甘肅蘭州臺(tái)站實(shí)測(cè)地磁C-響應(yīng)數(shù)據(jù)進(jìn)行反演，并討論該地區(qū)地幔轉(zhuǎn)換帶的性質(zhì)。

1 地磁測(cè)深原理

地球磁層中的電流可以激發(fā)產(chǎn)生磁場(chǎng)，在地磁測(cè)深中常被稱(chēng)為外源場(chǎng)，其能在地球內(nèi)部導(dǎo)電介質(zhì)中引起感應(yīng)場(chǎng)。因此，在地表觀(guān)測(cè)到的感應(yīng)場(chǎng)可攜帶地幔導(dǎo)電信息[32]。一般通過(guò)定義參數(shù)C-響應(yīng)建立地表觀(guān)測(cè)數(shù)據(jù)和地球內(nèi)部電性結(jié)構(gòu)的關(guān)系[33]。其表達(dá)式為

(1)

式中：tanθ為源空間的補(bǔ)償項(xiàng)；a0為地球半徑，取值6 370 km；在地磁測(cè)深研究中，采用球坐標(biāo)系，原點(diǎn)定義在地心處，則任意點(diǎn)的磁場(chǎng)H均包含指向地心的分量(Hr)、水平北向分量(Hθ)和水平東向分量(Hφ)，其中，φ為經(jīng)度，θ為緯度，r為指向地心的距離。

(2)

(3)

在球坐標(biāo)系中，對(duì)式(3)進(jìn)行變量分離，并取時(shí)諧因子eiωt，得到的磁場(chǎng)分量[35-36]為

(4)

(5)

而每層中Rn(r)與層導(dǎo)電性相關(guān)，滿(mǎn)足

[n(n+1)-(iωμσ)2r2]Rn(r)=0

(6)

考慮層界面上磁場(chǎng)連續(xù)性，相鄰的l-1和l層的r=rl邊界上Rn(r)滿(mǎn)足邊界條件

Rn,l-1(rl)=Rn,l(rl)

(7)

(8)

2 反演算法

2.1 算法流程

(9)

式中：λ為正則化因子；Nω為頻率個(gè)數(shù)；Cobs(ω)為觀(guān)測(cè)數(shù)據(jù)計(jì)算得到C-響應(yīng)；Cmod(ω)為模型理論C-響應(yīng)；fm為模型約束項(xiàng)，用以描述模型粗糙度；σ為電導(dǎo)率模型向量，是需要反演的參數(shù)。

需要說(shuō)明的是，組成地球的各薄球?qū)雍穸裙潭?，不參與反演。目標(biāo)函數(shù)最小化的過(guò)程與固體物質(zhì)冷卻退火達(dá)到能量最小穩(wěn)定狀態(tài)的過(guò)程相似。因此，模擬退火全局優(yōu)化算法應(yīng)該可以實(shí)現(xiàn)一維地磁測(cè)深反演。為此，將式(9)定義的目標(biāo)函數(shù)作為退火過(guò)程中的能量函數(shù)，通過(guò)合理降低溫度，當(dāng)目標(biāo)函數(shù)相應(yīng)的能量達(dá)到最小時(shí)各個(gè)分子對(duì)應(yīng)的狀態(tài)就是反演問(wèn)題的解。

本文提出的改進(jìn)模擬退火(Improved Simulated Annealing，ISA)算法是在傳統(tǒng)模擬退火算法反演的基礎(chǔ)上，結(jié)合統(tǒng)計(jì)分析改進(jìn)了反演的效果。其主要步驟包括：

步驟1，初始化T0、M、σmax、σmin、σ0、L、Rmax，j=0。

步驟2，置l=1，執(zhí)行鏈內(nèi)循環(huán)L次。

②當(dāng)R

④更新l=l+1。

步驟3，降低溫度(T)；j=j+1。

步驟4，當(dāng)j>M時(shí)，退火結(jié)束；否則，返回步驟2。

步驟5，統(tǒng)計(jì)最優(yōu)解。對(duì)擬合差小于Rmax的所有有效模型進(jìn)行統(tǒng)計(jì)分析用以確定最終解。

在上述迭代過(guò)程中，T0為初始溫度；j為當(dāng)前迭代次數(shù)；M為迭代次數(shù)上限；l為鏈內(nèi)當(dāng)前循環(huán)次數(shù)，受馬爾可夫鏈的鏈長(zhǎng)(L)控制，鏈長(zhǎng)越大則同一溫度下全局搜索次數(shù)越多，更有機(jī)會(huì)得到全局最優(yōu)解[37]，考慮時(shí)間成本，結(jié)合實(shí)驗(yàn)驗(yàn)證，本文采用L=10可以達(dá)到搜索效果；σmax和σmin分別為電導(dǎo)率參數(shù)取值上、下限，它們可以根據(jù)已知的巖石物理模型[38]或者前人經(jīng)驗(yàn)[13]進(jìn)行估計(jì)；σ0為模型初始電導(dǎo)率，因?yàn)槟M退火算法反演不需要固定的初始模型，所以本文采用有效范圍內(nèi)的隨機(jī)數(shù)作為初始值；Rmax表示擬合差閾值，其中數(shù)據(jù)擬合差(R)的定義為

(10)

擬合差閾值的合理選擇對(duì)反演至關(guān)重要。合理的閾值選擇根據(jù)C-響應(yīng)的數(shù)據(jù)質(zhì)量、迭代次數(shù)等綜合考慮。根本原則是要求該閾值能在退火過(guò)程中篩選出足夠大且符合對(duì)數(shù)正態(tài)分布的樣本量，否則需要在反演過(guò)程中調(diào)整擬合差閾值。

需要指出的是，在上述步驟中，傳統(tǒng)模擬退火算法執(zhí)行步驟1～4，而目標(biāo)函數(shù)的最優(yōu)解由最后一次退火結(jié)果確定。本文增加步驟5，使得反演的結(jié)果不單由最后一次退火結(jié)果確定，而是由所有符合約束條件解的統(tǒng)計(jì)期望值確定，從而更能反映隨機(jī)反演的特點(diǎn)。

2.2 模型參數(shù)隨機(jī)生成

隨機(jī)電導(dǎo)率參數(shù)的產(chǎn)生基于溫度的似柯西分布法，具體公式為

(11)

(12)

抽樣和反演的結(jié)果受模型參數(shù)的范圍約束。如果模型參數(shù)區(qū)間選擇不合理，則反演結(jié)果集中在區(qū)間邊界處(圖1)[31]，在區(qū)間不合理時(shí)應(yīng)適當(dāng)擴(kuò)大參數(shù)區(qū)間的上、下限。

2.3 Metropolis準(zhǔn)則

當(dāng)前被隨機(jī)修改后的模型需要依據(jù)Metropolis準(zhǔn)則去判斷是否被接受。該準(zhǔn)則需要計(jì)算新?tīng)顟B(tài)的接納概率(P)，其表達(dá)式為

(13)

Δ=(σj+1)-(σj)

(14)

陰影部分為預(yù)置參數(shù)變化區(qū)間；虛線(xiàn)代表理論真實(shí)值(黑色或者紅色)；實(shí)線(xiàn)代表反演結(jié)果圖1 模型約束范圍對(duì)反演結(jié)果的影響Fig.1 Influence of Constraint Range of Model on Inversion Results

2.4 溫度降低機(jī)制

合理的退火溫度降低機(jī)制可以高效控制迭代過(guò)程的收斂，從而獲得近似最優(yōu)解。反演過(guò)程中控制溫度的公式為

T(j)=T0exp(-cj1/m)

(15)

式中：T為當(dāng)前溫度；c為常數(shù)。

過(guò)高的初始溫度會(huì)使解在迭代初期收斂很慢，因此，選擇過(guò)高的初始溫度是沒(méi)有必要的。Sharma在直流電測(cè)深反演實(shí)驗(yàn)中發(fā)現(xiàn)，c取值為1，T0取值為0.01，m取值為2.5，可以保證反演的效率，并且在反演過(guò)程中它們不依賴(lài)于模型參數(shù)[31]。本文通過(guò)實(shí)驗(yàn)驗(yàn)證了這組參數(shù)對(duì)地磁測(cè)深也較為適用。

2.5 解分布統(tǒng)計(jì)分析

由于電導(dǎo)率參數(shù)不能為負(fù)值且為隨機(jī)數(shù)，所以迭代過(guò)程中所有的有效模型參數(shù)可以用對(duì)數(shù)正態(tài)分布來(lái)描述，且有l(wèi)nσ服從期望值為μ、方差為ε2的正態(tài)分布。就某一層而言，反演最終的電導(dǎo)率等于該層有效電導(dǎo)率參數(shù)集合的期望值E(σ)，其表達(dá)式為

E(σ)=eμ+ε2/2

(16)

其中，ε作為標(biāo)準(zhǔn)差可以用來(lái)衡量結(jié)果的不確定性。當(dāng)ε越小，表示估計(jì)值越可靠。距期望值一個(gè)標(biāo)準(zhǔn)差以?xún)?nèi)的范圍所包含的數(shù)據(jù)占總體數(shù)據(jù)的68%，本文將此范圍作為可靠區(qū)間。

估計(jì)最終解的可靠性也要通過(guò)檢驗(yàn)樣本數(shù)據(jù)是否符合對(duì)數(shù)正態(tài)分布加以驗(yàn)證，并采用對(duì)數(shù)正態(tài)概率分布P-P圖來(lái)判斷[39]。P-P圖是根據(jù)變量的累積概率對(duì)應(yīng)于自然對(duì)數(shù)正態(tài)理論分布累積概率所繪制的。其橫坐標(biāo)表示電導(dǎo)率樣本累積概率密度；縱坐標(biāo)為理論累積概率密度p(σ)，由概率密度函數(shù)計(jì)算得到。其表達(dá)式為

(17)

當(dāng)被檢測(cè)的數(shù)據(jù)符合對(duì)數(shù)正態(tài)分布時(shí)，在P-P圖中樣本累積概率和理論累積概率近似相等，從而在圖中表現(xiàn)為45°的似直線(xiàn)分布。從圖2可見(jiàn)，反演的電導(dǎo)率σ4樣本分布基本符合對(duì)數(shù)正態(tài)分布規(guī)律。

圖2 σ4樣本分布直方圖和對(duì)數(shù)正態(tài)概率分布P-P圖Fig.2 Histogram and P-P Graph of Lognormal Probability Distribution of Sample σ4

3 結(jié)果分析與討論

3.1 模型和理論C-響應(yīng)

為驗(yàn)證模擬退火算法應(yīng)用于地磁測(cè)深一維反演中的可行性及其準(zhǔn)確性，對(duì)合成模型進(jìn)行反演測(cè)試并對(duì)結(jié)果進(jìn)行討論。理論模型(表1)以Kelbert等的一維正演層狀模型為基礎(chǔ)[13]，其C-響應(yīng)(圖3)作為觀(guān)測(cè)數(shù)據(jù)進(jìn)行反演。反演過(guò)程中模型變化范圍根據(jù)Püthe等應(yīng)用蒙特卡羅抽樣方法獲得的400個(gè)獨(dú)立一維全球平均模型的最大變化區(qū)間確定[21]。

表1 理論模型電導(dǎo)率參數(shù)Tab.1 Electric Parameters of Theoretical Model

圖3 理論模型的C-響應(yīng)Fig.3 C-response of Theoretical Model

圖(b)中層電導(dǎo)率是反演結(jié)束時(shí)的模型參數(shù)圖4 理論模型的地磁響應(yīng)模擬退火算法反演結(jié)果Fig.4 Results of SA Algorithm Inversion of Geomagnetic Response of Theoretical Model

圖5 反演結(jié)果穩(wěn)定性對(duì)比Fig.5 Stability Comparisons of Inversion Result

3.2 模擬響應(yīng)數(shù)據(jù)反演

首先，對(duì)圖3中的純理論數(shù)據(jù)進(jìn)行反演測(cè)試。模擬退火時(shí)，迭代次數(shù)上限設(shè)置為20 000，鏈長(zhǎng)為1，其他參數(shù)同Sharma在直流電測(cè)深反演實(shí)驗(yàn)中的參數(shù)設(shè)置[31]。圖4給出反演過(guò)程中擬合差和模型電導(dǎo)率變化及迭代結(jié)束時(shí)各層電導(dǎo)率。由圖4可以看出，擬合差在反演過(guò)程中整體為降低趨勢(shì)。由于理論C-響應(yīng)不含噪聲且足夠精確，所以根據(jù)擬合差計(jì)算公式，其值可以降到接近0。當(dāng)?shù)螖?shù)超過(guò)3 000后，各個(gè)層的電導(dǎo)率變化集中在理論真實(shí)值的附近；當(dāng)?shù)螖?shù)大于10 000，反演得到的各層電導(dǎo)率與真實(shí)值相差極小，但其中依然存在經(jīng)過(guò)Metropolis準(zhǔn)則的判斷，較差模型被接受，而擬合差升高的情況出現(xiàn)。

由于噪聲在實(shí)際地磁數(shù)據(jù)中普遍存在[18]，一般認(rèn)為，受過(guò)大噪聲影響的C-響應(yīng)不能用于反演研究，所以本文對(duì)圖3中理論C-響應(yīng)添加可接受的8%高斯噪聲進(jìn)行反演方法測(cè)試。反演采用與理論數(shù)據(jù)反演相同的初始條件，分別進(jìn)行了3次獨(dú)立模擬退火算法反演，反演結(jié)束時(shí)的電導(dǎo)率模型如圖5所示。從圖5(a)可見(jiàn)，在反演結(jié)束時(shí)，3次反演數(shù)據(jù)擬合效果基本一樣；但從圖5(b)可見(jiàn)，受噪聲影響，反演結(jié)束時(shí)，傳統(tǒng)模擬退火算法獲得的最終解(3條藍(lán)色虛線(xiàn))存在較大差異。例如，σ3分別為0.160、0.168和0.367 S·m-1，都與0.200 S·m-1的理論電導(dǎo)率差異較大。這表明反演結(jié)果除了受方法本身的模型靈敏度差異影響外，即使采用全局反演方法，受噪聲等因素的影響，依然存在明顯多解性，不容易確定最合理的模型。

3.3 改進(jìn)模擬退火算法反演結(jié)果

根據(jù)圖5可知，傳統(tǒng)模擬退火算法反演只取最后一次迭代值作為結(jié)果，顯然忽略電導(dǎo)率參數(shù)在迭代過(guò)程中變化的分布特征，因此，如果結(jié)合統(tǒng)計(jì)分析方法可能會(huì)更準(zhǔn)確地計(jì)算模型電導(dǎo)率的估計(jì)值。本文采用上述解分布特征分析的方法后，再對(duì)3次反演結(jié)果進(jìn)行統(tǒng)計(jì)分析，得到的結(jié)果模型如圖5(b)中3條紅色實(shí)線(xiàn)所示。圖5(b)中，紅色虛線(xiàn)表示估計(jì)值上、下限，可以作為可靠區(qū)間。

對(duì)比圖5(a)、(b)可知，雖然傳統(tǒng)模擬退火算法反演結(jié)果相差較大，但改進(jìn)模擬退火算法反演結(jié)果不受隨機(jī)性影響，統(tǒng)計(jì)估計(jì)的解基本一樣，反映了解收斂的方向(圖中紅色實(shí)線(xiàn)基本重疊)。因此，改進(jìn)模擬退火算法反演具有更強(qiáng)的穩(wěn)定性和準(zhǔn)確性，受多解性影響小。3次反演過(guò)程中保留的數(shù)據(jù)集不同，但改進(jìn)模擬退火算法反演得到3個(gè)最終模型幾乎一致，且與理論模型高度近似，如σ3分別為0.214、0.218和0.220 S·m-1，接近0.200 S·m-1的理論電導(dǎo)率。雖然估計(jì)值與理論值有一定的差異，但每層電導(dǎo)率真實(shí)值都在可靠區(qū)間內(nèi)。需要指出的是，經(jīng)過(guò)實(shí)驗(yàn)，更多次的改進(jìn)模擬退火算法反演結(jié)果依舊保持高度近似。因此，在實(shí)際反演中，2、3次改進(jìn)模擬退火算法反演就可以實(shí)現(xiàn)電導(dǎo)率的可靠估計(jì)。

3.4 抗噪能力測(cè)試

為了進(jìn)一步考察改進(jìn)模擬退火算法反演技術(shù)抗噪能力，對(duì)圖3中理論C-響應(yīng)分別加上5%、10%、15%的高斯噪聲后進(jìn)行單次模擬退火算法反演，之后再進(jìn)行統(tǒng)計(jì)分析。除擬合差閾值外(因?yàn)樵搮?shù)的選擇需要考慮噪聲水平)，反演的其他參數(shù)都與圖5一致。隨著噪聲影響的增加，擬合差閾值受噪聲影響分別設(shè)置為：R5%,max=0.070、R10%,max=0.145、R15%,max=0.190。圖6展示了不同噪聲影響下C-響應(yīng)反演擬合情況。最終反演得到的擬合差依次為0.068、0.141、0.187，指示噪聲越大，擬合差越大。但從圖6可見(jiàn)，不同噪聲影響都得到與理論數(shù)據(jù)較好的擬合曲線(xiàn)。

圖6 不同噪聲影響下的C-響應(yīng)反演擬合Fig.6 Fitting Diagrams of C-response Inversion at Different Noises

圖(b)中不同顏色對(duì)應(yīng)的虛線(xiàn)表示可靠區(qū)間圖7 不同噪聲影響下的理論模型C-響應(yīng)反演結(jié)果Fig.7 Inversion Results of C-response of Theoretical Model at Different Noises

圖7(a)為4個(gè)樣本在3種噪聲影響下反演結(jié)果的P-P圖。從圖7(a)可見(jiàn)，反演得到的有效模型參數(shù)基本符合對(duì)數(shù)正態(tài)分布，且在當(dāng)前的噪聲水平條件下，基本不受噪聲的影響。圖7(b)給出了最終電導(dǎo)率模型和可靠區(qū)間。從圖7(b)可以看到，雖然受到不同噪聲影響，但真實(shí)值都在可靠區(qū)間內(nèi)。另外，反演電導(dǎo)率模型在410 km以淺和900 km以深基本重合。一方面可能是這些區(qū)間處于地磁測(cè)深的弱敏感區(qū)域，另一方面也可能反映本文方法壓制噪聲影響的能力。在地幔轉(zhuǎn)換帶深度附近，與理論模型比較，反演結(jié)果較好地反映理論模型參數(shù)。值得指出的是，根據(jù)噪聲水平對(duì)反演結(jié)果的影響，噪聲的存在導(dǎo)致反演的電導(dǎo)率稍微偏小，且存在噪聲越大，反演的導(dǎo)電性越差，可靠區(qū)間范圍越大。因此，當(dāng)數(shù)據(jù)存在小于15%的噪聲時(shí)，可以認(rèn)為本文的反演方法依然有效，能保證該方法對(duì)受噪聲影響的實(shí)測(cè)C-響應(yīng)數(shù)據(jù)反演的準(zhǔn)確性。

3.5 蘭州地區(qū)地磁測(cè)深數(shù)據(jù)反演

本次研究將改進(jìn)模擬退火算法應(yīng)用在蘭州臺(tái)站的實(shí)測(cè)地磁測(cè)深數(shù)據(jù)的一維反演中。蘭州臺(tái)站地處內(nèi)陸，受海洋效應(yīng)小[40]，因此，直接利用基于BIRRP方法[41]處理得到的C-響應(yīng)[圖8(a)]進(jìn)行反演。反演時(shí)根據(jù)深部地幔導(dǎo)電性隨著主要礦物相變而躍變的規(guī)律，將地核以上劃分為9層[13]。各層電導(dǎo)率區(qū)間的取值依據(jù)Semenov等反演結(jié)果[16]估計(jì)。模擬退火算法反演控制參數(shù)同圖5。參考圖6，估計(jì)圖8(a)的噪聲水平約為10%，設(shè)置Rmax為0.145。從圖8(a)可以看出，最終反演得到與原始數(shù)據(jù)相同的變化趨勢(shì)且光滑的C-響應(yīng)曲線(xiàn)，計(jì)算得到的最終模型對(duì)應(yīng)的擬合差為0.135。

圖8(b)顯示模型的各個(gè)電導(dǎo)率概率分布P-P曲線(xiàn)呈對(duì)角線(xiàn)，表明反演中擬合差小于閾值的模型符合對(duì)數(shù)正態(tài)分布特征，從而可以計(jì)算得到可靠的近似最優(yōu)解。圖8(c)給出了本文估計(jì)的最終近似最優(yōu)解。在0～250 km深度內(nèi)，蘭州臺(tái)站下方0～250 km電導(dǎo)率高于全球平均值[13]，且可靠區(qū)間與張艷輝等研究結(jié)果的置信區(qū)間[18]有較好的重合。在地磁測(cè)深的敏感深度范圍內(nèi)(410～1 200 km)，反演得到了與Semenov等相似的蘭州臺(tái)站下方電導(dǎo)率結(jié)果[16,18]。在410 km的不連續(xù)面上，沒(méi)有明顯的電導(dǎo)率躍遷，這與Yoshino等得到的大陸轉(zhuǎn)換帶導(dǎo)電性特征[42]一致；670～900 km深度的電導(dǎo)率略高于中國(guó)地區(qū)平均水平，與李世文的研究結(jié)論[43]相一致。需要指出的是：地幔轉(zhuǎn)換帶上層，本文反演得到的電導(dǎo)率大致為0.074 S·m-1，接近全球平均電導(dǎo)率；但轉(zhuǎn)換帶下層電導(dǎo)率為0.147 S·m-1，相對(duì)前人的結(jié)果，本文反演得到的電導(dǎo)率偏低。這個(gè)特點(diǎn)和歐洲中部地區(qū)轉(zhuǎn)換帶具有相同特征[13]。另外，雖然反演結(jié)果在轉(zhuǎn)換帶附近偏低，但是其可靠區(qū)間很好地包含了Semenov等在該深度的反演電導(dǎo)率[16,18]；且從對(duì)理論數(shù)據(jù)的噪聲測(cè)試來(lái)看，存在噪聲時(shí)，反演得到的電導(dǎo)率會(huì)有偏小的可能，反演結(jié)果在圖7(b)中顯示會(huì)左偏。因此，可以認(rèn)為本文的反演結(jié)果是合理的，而實(shí)際導(dǎo)電性可能更大，但不會(huì)超過(guò)全球平均模型。

為了討論蘭州地區(qū)地幔轉(zhuǎn)換帶下層低的導(dǎo)電性，本文采用Yoshino等提出的高溫高壓導(dǎo)電模型[42]對(duì)轉(zhuǎn)換帶下層做巖石物理計(jì)算。模型表達(dá)式為

(18)

式中：σV為體積電導(dǎo)率；Cw為含水量；α為幾何因子；k為玻爾茲曼常數(shù)；Ttr為熱力學(xué)溫度；σOH、σOP為指數(shù)前因子；HH、HP分別為電子跳躍傳導(dǎo)和質(zhì)子傳導(dǎo)的活化焓。

Yoshino等的模型參數(shù)[42]如表2所示，其中，k=1.38×10-23J·K-1。假設(shè)地幔轉(zhuǎn)換帶是絕熱的，暫時(shí)不考慮溫度異常，以全球平局進(jìn)行分析，在410、520、670 km深度處對(duì)應(yīng)的溫度分別為1 780 K、1 850 K和1 925 K[42]。

表2 巖石物理模型參數(shù)

圖(b)中電導(dǎo)率是最終模型參數(shù)，對(duì)應(yīng)圖(c)中的黑色實(shí)線(xiàn)。圖(c)中黑色實(shí)線(xiàn)代表改進(jìn)模擬退火算法反演結(jié)果；黑色虛線(xiàn)代表可靠區(qū)間；黃色實(shí)線(xiàn)從左至右分別代表含水量(質(zhì)量分?jǐn)?shù))為0%、0.1%、0.5%、1.0%的電導(dǎo)率剖面圖8 蘭州臺(tái)站地磁測(cè)深反演結(jié)果Fig.8 Inversion Results of Geomagnetic Depth Sounding at Lanzhou Station

圖8(c)中黃色實(shí)線(xiàn)為計(jì)算的轉(zhuǎn)換帶電導(dǎo)率剖面。從圖8(c)可以發(fā)現(xiàn)，即使采用全球平均溫度，為了重現(xiàn)反演結(jié)果，在轉(zhuǎn)換帶上層也需要要求0.5%的含水量，這與全球觀(guān)測(cè)結(jié)果一致。而在轉(zhuǎn)換帶下層，由于反演的導(dǎo)電性接近全球平均電導(dǎo)率的下限，所以要求該層基本不含水。導(dǎo)電性同時(shí)取決于溫度和含水量，因此，如果適當(dāng)降低溫度，該地區(qū)的轉(zhuǎn)換帶下層也可以適當(dāng)含水。由于蘭州地區(qū)轉(zhuǎn)換帶具體的溫度參數(shù)難以獲取，不排除蘭州地區(qū)轉(zhuǎn)換帶濕冷的可能性。在這種情況下，蘭州地區(qū)地磁觀(guān)測(cè)不能提供地幔柱的證據(jù)。但需要指出的是，由于C-響應(yīng)噪聲導(dǎo)致轉(zhuǎn)換帶下層導(dǎo)電性可能偏低，保守估計(jì)，其最多可以接近全球平均電導(dǎo)率，顯然重現(xiàn)這個(gè)電導(dǎo)率也不需要用高溫進(jìn)行解釋。

4 結(jié) 語(yǔ)

(1)基于模擬退火算法反演，結(jié)合統(tǒng)計(jì)分析，提出了地磁測(cè)深一維反演的新思路。理論數(shù)據(jù)反演測(cè)試結(jié)果表明，結(jié)合統(tǒng)計(jì)分析的模擬退火算法反演得到的結(jié)果穩(wěn)定，且具有良好的抗噪能力，但較大噪聲影響會(huì)使得反演的導(dǎo)電率偏低。

(2)對(duì)蘭州臺(tái)站的實(shí)測(cè)數(shù)據(jù)進(jìn)行反演得到的地幔轉(zhuǎn)換帶導(dǎo)電模型，其上層與前人反演結(jié)果相一致,但轉(zhuǎn)換帶下層反演得到的電導(dǎo)率偏低，推測(cè)含水量非常低，接近干枯狀態(tài)。蘭州地區(qū)的地磁測(cè)深不容易給出歐亞大陸中部存在地幔柱的證據(jù)。

謹(jǐn)以此文慶祝母校七十周年華誕，祝母校桃李遍天下，綠野追唐裴！1989年金秋，我從安徽農(nóng)村帶著木制的行李箱，坐汽車(chē)倒火車(chē)來(lái)到西安地質(zhì)學(xué)院，圖書(shū)館大廳里的新生報(bào)道開(kāi)啟了我四年的求學(xué)生活。圖書(shū)館厚大的棗紅色閱讀書(shū)桌，物探樓階梯教室明亮的燈光，422教室老師們時(shí)而慷慨激昂、時(shí)而循循善誘的課堂，混合著蔬菜和肥肉香味的特有鋼絲面食堂，還有89211班情同手足的兄弟姐妹，是我人生中最寶貴的回憶。母校是燈塔，指引著我們前行，更是我心中的牽掛!