劉云鶴，殷長(zhǎng)春

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

1 引 言

三維電磁法正反演技術(shù)近年發(fā)展較快，已經(jīng)開始從單純的理論研究走向?qū)嶋H應(yīng)用.三維航空電磁模擬需要在每個(gè)源位置單獨(dú)做一次正演導(dǎo)致計(jì)算量過大，無(wú)法實(shí)現(xiàn)快速數(shù)值計(jì)算，所以其數(shù)據(jù)反演相對(duì)于其他固定源或者源位置較少的電磁法發(fā)展緩慢.最近，Cox等提出一種利用航空電磁系統(tǒng)“footprint”概念來減小反演模型規(guī)模的算法［1］.利用該算法，整個(gè)反演區(qū)域被離散成每個(gè)測(cè)點(diǎn)下方“footprint”大小的子區(qū)域，而該點(diǎn)航空電磁響應(yīng)和數(shù)據(jù)的靈敏度矩陣只在這個(gè)子區(qū)域內(nèi)計(jì)算.另外，Oldenburg等利用矩陣分解技術(shù)實(shí)現(xiàn)了多發(fā)射源電磁法的快速正反演［2］.這兩種技術(shù)加速了航空電磁正演的過程，進(jìn)而使得大規(guī)模三維航空電磁反演成為可能.然而，在正演速度提高的基礎(chǔ)上我們還需要尋找一種穩(wěn)定和精確的最優(yōu)化算法以實(shí)現(xiàn)三維航空電磁的高效反演.

過去的十多年中，各種優(yōu)化算法被引用到求解多維情況下的最小值問題.主要包括Gauss－Newton（GN）算法、非線性共軛梯度法（NLCG）和擬牛頓法（QN）等.這幾種方法在內(nèi)存需求和計(jì)算速度上可以滿足三維電磁反演的需要.Mackie等、Siripunvaraporn等和林昌洪等利用GN方法從事大地電磁的三維反演［3－5］，而 Sasaki和 Oldenburg 等 利 用 GN 方 法 實(shí)現(xiàn)了三維航空電磁數(shù)據(jù)反演［2，6］.基于梯度的直接最優(yōu)化算法NLCG和有限內(nèi)存L－BFGS（標(biāo)準(zhǔn)的QN方法）都是節(jié)省內(nèi)存的有效反演方法，目前被廣泛應(yīng)用于大地電磁、海洋電磁和地面可控源類電磁法的數(shù)據(jù)反演中［7－16］.對(duì)于大規(guī)模的三維數(shù)據(jù)反演，為了節(jié)省內(nèi)存和計(jì)算時(shí)間，上述三種方法都非顯式地計(jì)算靈敏度矩陣，而是通過伴隨方法計(jì)算靈敏度矩陣和向量的乘積.GN方法通常利用預(yù)處理的共軛梯度法求解反演方程（PCG），在每次的共軛梯度法迭代過程中需要兩次額外的伴隨正演來計(jì)算靈敏度矩陣與向量的乘積.為了減少計(jì)算量一般都采用截?cái)嗟姆绞絹斫K止共軛梯度的迭代過程.NLCG和L－BFGS方法在每次迭代過程中只需要計(jì)算一次梯度，即只需要一次正演和一次伴隨正演.與GN方法相比，NLCG和L－BFGS方法每次迭代的正演計(jì)算量較少，然而反演收斂速度較慢.據(jù)Schwarzbach和Haber等的研究，GN方法在進(jìn)行伴隨正演過程中需要存儲(chǔ)每個(gè)源的全域解向量.當(dāng)發(fā)射源較多時(shí)，有限的內(nèi)存空間無(wú)法存儲(chǔ)所有源的全域解向量，因此該方法不適合多發(fā)射源的電磁反演［17］.Newman和Boggs曾對(duì)NLCG和L－BFGS的電磁反演效果進(jìn)行對(duì)比［9］.然而，其步長(zhǎng)搜索存在缺陷，導(dǎo)致 LBFGS算法的迭代步長(zhǎng)總是不能達(dá)到最優(yōu)理論值1；由此，L－BFGS較之于NLCG在反演效果上的優(yōu)越性沒有得到體現(xiàn).本文中，我們采用Egbert和Kelbert提出的改進(jìn)NLCG步長(zhǎng)搜索技術(shù)及Nocedal和 Wright提出的L－BFGS步長(zhǎng)搜索技術(shù)進(jìn)行三維航空電磁數(shù)據(jù)反演［18］，并以垂直磁偶極子發(fā)射多分量觀測(cè)方式為例，比較NLCG和L－BFGS在三維航空電磁反演中的應(yīng)用效果.

2 正演模擬

三維頻率域航空電磁正演計(jì)算基于下面的二次電場(chǎng)的偏微分方程：

其中Ep為發(fā)射偶極在背景電導(dǎo)率σb中產(chǎn)生的一次背景場(chǎng)，Es為二次電場(chǎng)，σ為模型電導(dǎo)率，ω為角頻率，μ為磁導(dǎo)率［19］.總的電場(chǎng)為E＝Ep＋Es，而磁場(chǎng)H由公式

得到.本文采用交錯(cuò)網(wǎng)格有限差分技術(shù)離散方程（1），最終形成的大型方程組采用擬最小殘差法（QMR）求解.圖1給出的正演模型采用Newman和Alumbaugh在1995年給出的模型［19］.發(fā)射和接收裝置離地面20m，發(fā)射－接收源距為10m.我們考慮水平共面HCP與垂直共軸VCX兩種裝置，發(fā)射頻率為900Hz.計(jì)算區(qū)域采用長(zhǎng)方體單元進(jìn)行離散，網(wǎng)格數(shù)為50×50×40.圖2中給出本文算法與Newman和Alumbaugh結(jié)果的比較.由圖可以看出，本文的結(jié)果與Newman和Alumbaugh的結(jié)果整體吻合很好.虛部在異常體正上方的誤差稍大，但仍小于1ppm.

圖1 三維航空電磁正演模型Fig.1 3Dearth for AEM modeling（Newman and Alumbaugh，1995）

圖2 三維頻率域航空電磁正演結(jié)果驗(yàn)證Fig.2 Model results for 3Dfrequency－domain AEM systems in comparison to those of Newman and Alumbaugh（1995）

3 反演理論

3.1 正則化

依據(jù)Egbert和Kelbert的闡述［15］，本文通過最小化目標(biāo)函數(shù)

來解決正則化反演問題.上式中，m是M維地球電導(dǎo)率模型參數(shù)向量，d是Nd維數(shù)據(jù)向量，Cd是數(shù)據(jù)方差矩陣，f（m）為正演算子，m0是先驗(yàn)或初始估計(jì)模型，λ是正則化因子，Cm是模型方差矩陣.通過變換（3）式可以進(jìn)一步簡(jiǎn)化.假設(shè)

目標(biāo)函數(shù)的梯度為

其中J是靈敏度矩陣，最終

的模型參數(shù)可以通過

來計(jì)算.為簡(jiǎn)便起見，我們?cè)诤竺娴挠懻撝泻雎阅Ｐ蛥?shù)上面的波浪線.

3.2 反演約束

與傳統(tǒng)的模型約束方式不同，本文中的Cm不是僅考慮每個(gè)離散單元和其相鄰的單元之間的影響，而是將整個(gè)模型所有元素之間的相互影響關(guān)聯(lián)起來.Cm的計(jì)算采用Egbert和Siripunvaraporn等提出的一種遞歸自回歸算法［20－21］.該算法的實(shí)質(zhì)為求解一個(gè)擴(kuò)散方程，實(shí)現(xiàn)過程為，先在x、y、z三個(gè)方向每個(gè)參數(shù)都在與前一個(gè)參數(shù)和光滑因子的乘積求和后，再與光滑因子相乘加到下一個(gè)參數(shù)，這樣便將所有參數(shù)關(guān)聯(lián)起來，然后對(duì)已經(jīng)光滑的所有參數(shù)按z、y、x方向再反向光滑一遍.（6）式中，C－1／2m除了通過正則化因子λ對(duì)模型電導(dǎo)率進(jìn)行約束外，還直接約束目標(biāo)函數(shù)梯度表達(dá)式中的數(shù)據(jù)梯度部分.類似的直接修正梯度方法已得到廣泛應(yīng)用.Plessix和Mulder提出用變尺度的趨膚深度冪指數(shù)對(duì)梯度進(jìn)行修正［22］；Mackie和Newman等將近似的Hessian矩 陣 作 為 梯 度 的 預(yù) 處 理 算 子［9，23－24］；Avdeev 和Avdeeva提出了一種附加正則化方法來修正梯度［14］.這些修正梯度的方法可以使淺層異常的反演變得穩(wěn)定，加速了反演過程.本文中，我們進(jìn)一步深入研究中光滑因子對(duì)反演結(jié)果的影響.進(jìn)而，為改變反演的穩(wěn)定性和速度，我們?cè)诜囱葜胁捎昧俗兓墓饣蜃?在反演初始階段采用大的光滑因子以確定準(zhǔn)確的異常體空間位置；而當(dāng)數(shù)據(jù)擬合差下降緩慢時(shí)，則通過減小光滑因子獲得準(zhǔn)確的異常體大小和電阻率值.

3.3 反演流程

NLCG方法反演流程：

1）令i＝1，選擇初始模型mi及模型光滑因子，并計(jì)算ri＝－Δφ（mi，d）；

2）如果‖ri‖很小則停止，否則令M為預(yù)處理矩陣；

3）搜索αi使φ（mi＋αiui，d）最小化.當(dāng)步長(zhǎng)αi太小時(shí)減小正則化因子λ，如果λ達(dá)到閥值則反演終止；

4）令mi＋1＝mi＋αiui和ri＋1＝－Δφ（mi＋1，d）；

5）如果‖ri＋1‖很小則停止，否則向下繼續(xù)執(zhí)行；

L－BFGS方法反演流程：

1）令i＝1，選擇初始模型mi、模型光滑因子及初始對(duì)稱正定矩陣Hi，設(shè)定導(dǎo)數(shù)控制精度ε＞0；

2）計(jì)算ri＝－Δφ（mi，d）.如果 ‖ri‖ ≤ε，則迭代終止，輸出最終解mi；

3）否則，令ui＝Hiri，搜索αi.當(dāng)步長(zhǎng)αi太小時(shí)減小正則化因子λ.如果λ達(dá)到閥值則反演終止；

4）令mi＋1＝mi＋αiui，如果 ‖ri＋1‖ ≤ε，則得到最終解mi＋1；

6）令i＝i＋1，然后跳到2）.

NLCG和L－BFGS反演過程中均采用自適應(yīng)方式選擇正則化因子λ.具體實(shí)現(xiàn)方式為：當(dāng)數(shù)據(jù)擬合差下降小于0.002時(shí)，正則化因子λ變?yōu)樵瓉淼氖种?正則化因子λ在反演中約束步長(zhǎng)的大小，每減小一次λ導(dǎo)致正則化部分在目標(biāo)函數(shù)的比重變小，模型修正空間變大，可以搜索較大的步長(zhǎng).當(dāng)λ減小到閥值時(shí)，反演基本只擬合數(shù)據(jù)部分，如果數(shù)據(jù)擬合差仍緩慢下降，說明目標(biāo)函數(shù)值達(dá)到極小.另外，如果目標(biāo)函數(shù)梯度過小，目標(biāo)函數(shù)值也達(dá)到極小.因?yàn)榧词故抢碚摂?shù)據(jù)也不能使數(shù)據(jù)擬合差總降低到0，因此本文的理論數(shù)據(jù)反演并沒有給定數(shù)據(jù)擬合差閥值，而是采用λ和梯度是否達(dá)到閥值確定終止反演.

NLCG和L－BFGS反演過程中需要計(jì)算目標(biāo)函數(shù)的梯度和搜索步長(zhǎng).下面分別介紹梯度的計(jì)算方法和步長(zhǎng)的搜索方法.

3.4 梯度計(jì)算

有限差分方法的最終方程組可以寫成

其中K為系數(shù)矩陣，s為方程（1）的右端源項(xiàng).將（8）

式兩端對(duì)第k個(gè)模型參數(shù)mk求導(dǎo)得到［9］

將方程（1）中源項(xiàng)表達(dá)式帶入（9）式并進(jìn)行化簡(jiǎn)，得到

將方程（2）簡(jiǎn)寫為Hs＝LEs，其中L為空間插值算子，則磁場(chǎng)的導(dǎo)數(shù)為

假設(shè)如下Nd×M維矩陣

則有

將（13）式帶入（6）式，則（6）式右端第一項(xiàng)變?yōu)?/p>

令

則可以通過求解伴隨正演

得到v.將v帶入到（14）式中就可以得到數(shù)據(jù)部分的梯度值，而模型部分的梯度值是解析的，這樣只通過一次正演和一次伴隨正演就可計(jì)算出目標(biāo)函數(shù)的梯度.

3.5 步長(zhǎng)搜索

與其他可以通過解析方法得到步長(zhǎng)的方法（比如共軛梯度法等）不同，NLCG和L－BFGS兩種方法需要通過線性搜索得到反演迭代步長(zhǎng).NLCG方法需要沿負(fù)梯度及其共軛方向搜索步長(zhǎng)；而L－BFGS方法則沿經(jīng)近似Hessian矩陣的逆修正后的負(fù)梯度方向搜索步長(zhǎng).最優(yōu)化過程的步長(zhǎng)一般需要滿足如下充分下降條件

和曲率條件

其中，f為正演算子，mk為第k次迭代的模型參數(shù)向量，αk為迭代步長(zhǎng)，p為搜索方向，c1和c2為常數(shù)，滿足0＜c1＜c2＜1，k為迭代次數(shù).這兩種條件統(tǒng)稱為 Wolfe條件.在實(shí)際應(yīng)用中，一般取c1＝0.0001；對(duì)于GN或者 QN方法，c2＝0.9，而對(duì)于非線性共軛梯度法c2＝0.1.另外，還存在一種逆向追蹤方法（back tracing）用來約束反演步長(zhǎng)，其作用和（18）式類似.逆向追蹤法的基本步驟為：選取一個(gè)在0到1之間變化的值，將其與未滿足充分下降條件的步長(zhǎng)做乘積，直到步長(zhǎng)滿足充分下降條件（17）為止.實(shí)際反演中，通常選取多次函數(shù)插值方式來進(jìn)行逆向追蹤［18］.目前，NLCG方法廣泛采用充分下降條件和逆向追蹤方法組合約束反演步長(zhǎng)，而LBFGS利用充分下降條件和曲率條件來約束步長(zhǎng).

本文采用的NLCG步長(zhǎng)搜索方式為：

1）第一次初始迭代步長(zhǎng)設(shè)為＝0.001.以后每次迭代中，初始步長(zhǎng)首先設(shè)為

2）利用二次方函數(shù)插值技術(shù)得到步長(zhǎng).如果滿足充分下降條件（17），則判斷f（mk＋）是否成立.如果成立則取，否則取；如果不滿足充分下降條件（17），則利用三次方函數(shù)插值技術(shù)進(jìn)行逆向追蹤.

L－BFGS步長(zhǎng)搜索方式為：

1）每次迭代反演初始步長(zhǎng)均為＝1，判斷此步長(zhǎng)是否滿足Wolfe條件，如果滿足則取初始步長(zhǎng)為該次迭代步長(zhǎng)；

2）否則，如果初始步長(zhǎng)不滿足Wolfe條件，則利用線性搜索算法求取迭代步長(zhǎng).具體算法包括兩層循環(huán)：外層循環(huán)為將整個(gè)步長(zhǎng)求取區(qū)間［0，αmax］利用插值技術(shù)分成數(shù)個(gè)子區(qū)間，然后依次在每個(gè)子區(qū)間中搜索滿足Wolfe條件的步長(zhǎng)，當(dāng)找到滿足條件的步長(zhǎng)后循環(huán)自動(dòng)結(jié)束；內(nèi)層循環(huán)為在每個(gè)子區(qū)間搜索滿足Wolfe條件的步長(zhǎng)時(shí)，利用二次方函數(shù)插值以及三次方函數(shù)插值技術(shù)計(jì)算迭代步長(zhǎng)αj，如果得到的步長(zhǎng)不合適，則利用該步長(zhǎng)更新子區(qū)間的邊界，然后在更新后的子區(qū)間，重復(fù)上面計(jì)算步長(zhǎng)的步驟.內(nèi)層循環(huán)直到找到合適步長(zhǎng)或者在子空間內(nèi)已搜索完畢后終止.如果在區(qū)間［0，αmax］中得不到滿足Wolfe條件的步長(zhǎng)αk，則重新計(jì)算新的目標(biāo)函數(shù)值和導(dǎo)數(shù)，然后重新搜索.pk）＜f（mk＋

4 理論數(shù)據(jù)反演

為了進(jìn)行理論數(shù)據(jù)的反演測(cè)試，本文首先設(shè)計(jì)一個(gè)電阻率為100Ωm的背景模型，并將其離散成28×28×25個(gè)網(wǎng)格，其中在x和y方向包括8個(gè)擴(kuò)邊，在z方向包括4個(gè)向下的擴(kuò)邊和8個(gè)空氣層.模型中心區(qū)域三個(gè)方向的網(wǎng)格邊長(zhǎng)均為10m，擴(kuò)邊的大小依次為20m，40m，80m，160m.前兩個(gè)空氣層在z方向的長(zhǎng)度均為10m，上方的空氣層厚度的遞增系數(shù)為3，最頂層厚度為30000m.觀測(cè)方式選擇水平共面裝置（HCP），發(fā)射與接收裝置的高度設(shè)為20m，發(fā)射頻率為900Hz和5000Hz，觀測(cè)數(shù)據(jù)包含三個(gè)方向的磁場(chǎng)分量.測(cè)點(diǎn)間距在x和y方向均為20m，總計(jì)100個(gè)觀測(cè)點(diǎn).根據(jù)“footprint”概念，某測(cè)點(diǎn)的響應(yīng)大小只與其下方一定范圍的地電模型相關(guān)［1］.在反演大區(qū)域數(shù)據(jù)時(shí)可以將整個(gè)大區(qū)域離散成每個(gè)測(cè)點(diǎn)下方的“footprint”大小的模型進(jìn)行靈敏度矩陣或靈敏度矩陣與向量乘積的計(jì)算.由于本文選取的模型大小與所選的航空電磁系統(tǒng)的“footprint”尺度接近，因此本文算例分析實(shí)質(zhì)上反應(yīng)了航空電磁基于“footprint”反演技術(shù)的效率.

反演測(cè)試的第一個(gè)模型為將上述背景模型中嵌入一個(gè)30Ωm的異常體.異常體頂部埋深為20m，三個(gè)方向的長(zhǎng)度均為40m（見圖1）.反演開始時(shí)，式（1）中的正則化因子λ設(shè)為10，且在反演過程中采用自適應(yīng)方式遞減.當(dāng)正則化因子減小到0.0001時(shí)，數(shù)據(jù)擬合誤差基本不再變化，反演過程自動(dòng)終止.圖3給出理論數(shù)據(jù)的反演結(jié)果.從圖可以看出，當(dāng)選擇光滑因子為0.3時(shí)，L－BFGS和NLCG的反演結(jié)果中，異常體的電阻率比較接近真實(shí)電阻率值，但是位置均較真實(shí)位置更靠近地表；當(dāng)選擇光滑因子為0.5時(shí)，L－BFGS和NLCG的反演結(jié)果中異常體的位置接近正演模型中的異常體真實(shí)位置，但是反演電阻率值與真實(shí)值偏差較大，且異常分布范圍較廣；當(dāng)光滑因子初始值選為0.5，而在數(shù)據(jù)擬合差下降小于0.002時(shí)，光滑因子調(diào)整到0.3（判斷條件與正則化因子自適應(yīng)算法相同，但光滑因子先于正則化因子變化），發(fā)現(xiàn)L－BFGS和NLCG的反演結(jié)果中，異常體的位置和電阻率較光滑因子分別為0.3和0.5時(shí)更接近模型真實(shí)值.對(duì)于圖3中所選的三種光滑因子，L－BFGS均比NLCG方法得到的結(jié)果準(zhǔn)確.另外，從圖4給出的數(shù)據(jù)擬合誤差隨迭代次數(shù)的變化可以看出，L－BFGS方法的數(shù)據(jù)擬合差下降速度比NLCG快.如果從正演次數(shù)來考慮反演效率，L－BFGS方法較NLCG方法體現(xiàn)出明顯優(yōu)勢(shì).事實(shí)上，當(dāng)反演迭代步長(zhǎng)為1時(shí)，L－BFGS方法每次迭代需要的正演次數(shù)為2次；而對(duì)于圖4中迭代步長(zhǎng)不為1的情況，L－BFGS每次迭代需要進(jìn)行4次正演.綜合起來，L－BFGS每次反演迭代大約需要作2.2次正演計(jì)算.根據(jù)Newman等的結(jié)論，NLCG每次反演迭代至少需要3次正演［7］.

圖3 單個(gè)異常體模型及不同光滑因子L－BFGS和NLCG反演結(jié)果Fig.3 Single anomalous target model and L－BFGS and NLCG inversions for different smoothing parameters

圖4 單個(gè)異常體反演的數(shù)據(jù)擬合誤差及反演步長(zhǎng)隨迭代次數(shù)變化Fig.4 Data misfit and model modification vs.iterations for single anomalous target model

為了進(jìn)一步在復(fù)雜模型中比較這兩種反演方法的有效性，我們將圖3中的單個(gè)異常體替換為8個(gè)異常體.如圖5a所示，靠近地表的4個(gè)異常體大小為50m×50m×10m，電阻率分別為30Ωm和300Ωm；另外4個(gè)異常體埋深均為30m，大小為50m×50m×30m，電阻率分別為20Ωm和2000Ωm.對(duì)于這個(gè)模型，正則化因子λ在反演開始時(shí)設(shè)為10000，最小值仍為0.0001.反演過程中采用變光滑因子0.5～0.3.比較圖5b和5c中的反演結(jié)果可以看出，LBFGS較NLCG方法更好地恢復(fù)了真實(shí)模型.特別是對(duì)地下高阻異常的反映，L－BFGS和NLCG方法存在較大差異.NLCG方法沒有得到高阻異常的準(zhǔn)確消息，而L－BFGS方法基本上確定了高阻異常的位置和形狀.從圖5中還可以看出，即使是利用LBFGS方法，地下高阻異常的電阻率也只反演到300Ωm左右，遠(yuǎn)低于真實(shí)的2000Ωm；而NLCG方法只反演到了170Ωm左右，距離真實(shí)值差距更大.Oldenburg等認(rèn)為這種現(xiàn)象的出現(xiàn)是正常的，因?yàn)槔酶袘?yīng)源和磁場(chǎng)數(shù)據(jù)本質(zhì)上很難反演出高阻異常體［2］.特別是對(duì)于航空電磁，由于發(fā)射與接收機(jī)距地面均有一定距離，電磁響應(yīng)的體積效應(yīng)會(huì)掩蓋高阻異常體信息.圖6中數(shù)據(jù)擬合差隨反演迭代次數(shù)的變化特征表明NLCG方法過早地陷入了局部極小，導(dǎo)致反演結(jié)果沒有L－BFGS方法的準(zhǔn)確.

圖5 8個(gè)異常體模型與L－BFGS和NLCG變光滑因子反演結(jié)果Fig.5 Model with 8anomalous targets and L－BFGS and NLCG inversion results for variable smoothing parameters

圖6 8個(gè)異常體反演的數(shù)據(jù)擬合誤差及反演步長(zhǎng)隨迭代次數(shù)變化Fig.6 Data misfit and model modification vs.iterations for the 8anomalous targets model

如3.3節(jié)所述，正則化因子λ的變化會(huì)增大模型修正空間，對(duì)于NLCG算法，導(dǎo)致反演步長(zhǎng)變大，數(shù)據(jù)擬合差下降變快（擬合差下降曲線突變的地方）；而對(duì)于L－BFGS算法，其步長(zhǎng)選擇不受正則化因子變化的影響，數(shù)據(jù)擬合差曲線較為平緩.另外，從上述兩個(gè)反演測(cè)試?yán)涌梢钥闯觯琇－BFGS方法基本采用1作為迭代步長(zhǎng)；相比之下，NLCG方法的迭代步長(zhǎng)特別小.其原因是NLCG方法的反演過程類似于最速下降法的線性優(yōu)化算法，它需要很小的步長(zhǎng)使反演收斂過程體現(xiàn)出超線性（super linear）反演算法特征.

5 反演效率分析

本文中所有的數(shù)值測(cè)試都是在PC上完成，其配備有Intel?CoreTMi3－2100CPU ＠3.10GHz和4GB RAM.兩組模型反演結(jié)果表明，L－BFGS方法每次迭代大約需要2.2次正演，而NLCG大約需要3次正演計(jì)算.對(duì)于本文設(shè)計(jì)的航空電磁采集方式，每次正演包括100個(gè)位置的2個(gè)頻率的子正演，共200次子正演.統(tǒng)計(jì)結(jié)果表明，每次正演大約耗時(shí)8min，即單個(gè)頻率每個(gè)位置的子正演耗時(shí)大約需2.4s.以單個(gè)測(cè)點(diǎn)來考慮，每個(gè)頻率100次迭代NLCG耗時(shí)約720s，而L－BFGS耗時(shí)大約需528s.本文中討論的兩種算法其內(nèi)存占用量接近，均為100MB左右.本文中的算例進(jìn)一步表明，如果設(shè)定相同的數(shù)據(jù)擬合差，L－BFGS較NLCG方法節(jié)約一半以上計(jì)算時(shí)間.

6 結(jié) 論

本文通過引入變化的光滑因子，改善了三維航空電磁數(shù)據(jù)反演的效果.反演結(jié)果中異常體的位置和電性較之于固定光滑因子條件下的反演結(jié)果更接近真實(shí)模型.以此為基礎(chǔ)，通過對(duì)L－BFGS和NLCG方法的對(duì)比分析，發(fā)現(xiàn)L－BFGS方法不論是在反演速度上還是反演結(jié)果的準(zhǔn)確性上都明顯優(yōu)于NLCG，因此更適合于處理大規(guī)模三維航空電磁數(shù)據(jù)的反演問題.由于采用的反演模型尺寸接近于航空電磁系統(tǒng)的“footprint”，所以本文單個(gè)測(cè)點(diǎn)的反演效率與目前廣泛使用的“moving footprint”方法的反演效率相當(dāng)［1］.然而，分析表明，對(duì)于本文設(shè)計(jì)的理論模型，即使是利用L－BFGS方法，每個(gè)測(cè)點(diǎn)單個(gè)頻率仍需要10min以上的反演耗時(shí)，尚不能滿足航空電磁海量數(shù)據(jù)的三維反演需求.本文的主要任務(wù)是尋找一種高效的最優(yōu)化算法以實(shí)現(xiàn)精確、穩(wěn)定的三維航空電磁反演，對(duì)于在此基礎(chǔ)上如何進(jìn)一步提高三維航空電磁正演速度及反演效率還需要深入研究.未來可通過結(jié)合模型分區(qū)和矩陣分解技術(shù)改善航空電磁正演速度，進(jìn)而提高航空電磁數(shù)據(jù)的三維反演效率.這將是本研究團(tuán)隊(duì)的下一步研究目標(biāo).

（

）

［1］ Cox L H，Wilson G A，Zhdanov M S.3Dinversion of airborne electromagnetic data using a moving footprint.ExplorationGeophysics，2010，41（4）：250－259.

［2］ Oldenburg D W，Haber E，Shekhtman R.3Dinversion of multi－source time domain electromagnetic data. http：／／www.math.ubc.ca／～haber／Eldad＿Haber＿Web＿Page／Publications.html，2012［2013－03－21］.

［3］ Mackie R L，Madden T R.Three－dimensional magnetotelluric inversion using conjugate gradients.Geophys.J.Int.，1993，115（1）：215－229.

［4］ Sirepunvarpporn W，Egbert G D.Data space conjugate gradient inversion for 2－D magnetotelluric data.Geophys.J.Int.，2007，170（3）：986－994.

［5］ 林昌洪，譚捍東，佟拓.傾子資料三維共軛梯度反演研究.地球物理學(xué)報(bào)，2011，54（4）：1106－1113.Lin C H，Tan H D，Tong T.Three－dimensional conjugate gradient inversion of tipper data.ChineseJ.Geophys.（in Chinese），2011，54（4）：1106－1113.

［6］ Sasaki Y.Full 3－D inversion of electromagnetic data on PC.JournalofAppliedGeophysics，2001，46（1）：45－54.

［7］ Newman G A，Alumbaugh D L.Three－dimensional magnetotelluric inversion using non－linear conjugate gradients.Geophys.J.Int.，2000，140（2）：410－424.

［8］ Rodi W L，Mackie R L.Nonlinear conjugate gradients algorithm for 2－D magnetotelluric inversion.Geophysics，2001，66（1）：174－187.

［9］ Newman G A，Boggs P T.Solution accelerators for largescale three－dimensional electromagnetic inverse problems.InverseProblems，2004，20（6）：151－170.

［10］ Avdeev D B.Three－dimensional electromagnetic modelling and inversion from theory to application.Surv.Geophys.，2005，26（6）：767－799.

［11］ Haber E.Quasi－Newton methods for large－scale electromagnetic inverse problems.InverseProblems，2005，21（1）：305－323.

［12］ Plessix R E，Peter S.3DCSEM modeling and inversion in complex geologic settings.77thAnnual International Meeting，SEG，Expanded Abstracts，2007：589－593.

［13］ Kelbert A，Egbert G D，Schultz A.Non－linear conjugate gradient inversion for global EM induction：resolution studies.Geophys.J.Int.，2008，173（2）：365－381.

［14］ Avdeev D，Avdeeva A.3Dmagnetotelluric inversion using a limited－memory quasi－Newton optimization.Geophysics，2009，74（3）：F45－F57.

［15］ Egbert G D， Kelbert A. Computational recipes for electromagnetic inverse problems.Geophys.J.Int.，2012，189（1）：251－267.

［16］ 翁愛華，劉云鶴，賈定宇等.地面可控源頻率測(cè)深三維非線性共軛梯度反演.地球物理學(xué)報(bào)，2012，55（10）：3506－3515.Weng A H，Liu Y H，Jia D Y，et al.Three－dimensional controlled source electromagnetic inversion using non－linear conjugate gradients.ChineseJ.Geophys.（in Chinese），2012，55（10）：3506－3515.

［17］ Schwarzbach C，Haber E.Finite element based inversion for time－h(huán)armonic electromagnetic problems.http：／／www.math.ubc.ca／～haber／Eldad＿Haber＿Web＿Page／Publications.html，2012［2013－03－21］.

［18］ Nocedal J，Wright S J.Numerical Optimization.New York：Springer－Verlag，1999.

［19］ Newman G A，Alumbaugh D L.Frequency－domain modelling of airborne electromagnetic responses using staggered finite differences.GeophysicalProspecting，1995，43（8）：1021－1042.

［20］ Egbert G D.A new stochastic process on the sphere：application to characterization of long－period global scale external sources.14thWorkshop on Electromagnetic Induction in the Earth and Moon，Brest，F(xiàn)rance，1994.

［21］ Siripunvaraporn W，Egbert G D.An efficient data－subspace inversion method for 2－D magnetotelluric data.Geophysics，65（3）：791－803.

［22］ Plessix R E， Mulder W A. Resistivity imaging with controlled－source electromagnetic data： Depth and data weighting.InverseProblems，2008，24：1－22.

［23］ Mackie R L，Rodi W，Watts M D，et al.3Dmagnetotelluric inversion for resource exploration.71stAnnual International Meeting，SEG，Expanded Abstracts，2001：1501－1504.

［24］ Mackie R L，Watts M D，Rodi W，et al.Joint 3Dinversion of marine CSEM and MT data.77thAnnual International Meeting，SEG，Expanded Abstracts，2007：574－578.