楊克思，朱培民，趙 娜，徐 陽

（1.中國石油吉林油田分公司a.地球物理勘探研究院，b.新民采油廠工藝所，松原 138000；2.中國地質(zhì)大學(xué) 地球物理與空間信息學(xué)院，武漢 430074）

YANG Ke－si1a，ZHU Pei－min2，ZHAO Na2，XU Yang1b

（1.Jilin Oilfield Branch，PetroChina a.Geophysical Exploration Research Institute，b.Xinmin Oil Production Factory，Songyuan 138000，China；2.China University of Geosciences，Institute of Geophysics and Geomatics，Wuhan 430074，China）

一種提高稀疏觀測數(shù)據(jù)區(qū)反演精度的策略

楊克思1a，朱培民2，趙 娜2，徐 陽1b

（1.中國石油吉林油田分公司a.地球物理勘探研究院，b.新民采油廠工藝所，松原 138000；2.中國地質(zhì)大學(xué) 地球物理與空間信息學(xué)院，武漢 430074）

在一個(gè)觀測數(shù)據(jù)稀疏的地區(qū)，僅用觀測數(shù)據(jù)來獲得接近于真實(shí)的模型解是非常困難的，因此在反演時(shí)加入先驗(yàn)信息，考慮基于參考模型的反演方法是合理的。在一個(gè)觀測數(shù)據(jù)有限，先驗(yàn)信息也很少的區(qū)域，獲取合適的參考模型非常困難。在無法選取合適參考模型的情況下，提出了一種在參考模型遠(yuǎn)離真實(shí)模型的情況下，通過反演迭代過程中對參考模型的不斷修改，使參考模型逐步逼近真實(shí)模型，從而提高了反演精度的策略。這里設(shè)計(jì)了均勻?qū)訝钅Ｐ图昂惓ｓw的模型對算法進(jìn)行了驗(yàn)證，證實(shí)了此方法對反演大型地質(zhì)異常體具有較好的適用性。

反演；天然地震層析成像；稀疏觀測數(shù)據(jù)；模型修改方法

0 引言

在一個(gè)稀疏觀測數(shù)據(jù)區(qū)，其反演問題往往是欠定問題，因此在反演過程中，不能獲得唯一的反演結(jié)果。為了解決這一問題，可以考慮基于模型的反演，對模型施加合理的約束條件（例如讓模型長度‖m‖2最?。５Ｐ烷L度最小并不是一個(gè)最好的度量（例如解有關(guān)海洋中密度擾動的反演問題，就不能找到一個(gè)模型長度最小的解），而選擇一個(gè)參考模型m0，使要求解的模型與參考模型的差距‖mm0‖2最小，會更合理一些，這就是基于參考模型的反演［1］。

基于參考模型的反演需要給定一個(gè)合適的參考模型，而參考模型選擇的適當(dāng)與否，關(guān)系到反演結(jié)果的好壞及迭代的效率。如果參考模型選擇的不適當(dāng)，可能會使迭代次數(shù)增多，甚至收斂到一個(gè)與真實(shí)模型相差很大的反演結(jié)果。

給定參考模型的方法很多，常見的方法有兩種：①人工給出一個(gè)均勻模型［2］；②根據(jù)目標(biāo)區(qū)現(xiàn)有的相關(guān)資料，通過適當(dāng)?shù)奶幚斫o出一個(gè)比較粗略的參考模型。第一種方法簡單易行，但給出的參考模型隨意性很大，有可能造成算法收斂緩慢甚至不收斂；第二種方法合理可靠，但在數(shù)據(jù)稀疏、噪聲較大的情況下或調(diào)查的初期階段，仍無法得到合適的參考模型。

在無法選取合適參考模型的情況下，作者提出了一種在參考模型遠(yuǎn)離真實(shí)模型的情況下，通過反演迭代過程中對參考模型的不斷修改，使參考模型逐步逼近真實(shí)模型，從而提高反演精度的策略。這種方法的主要思想，是將每次迭代反演后的結(jié)果進(jìn)行處理，作為下次迭代的參考模型，并將最后輸出的參考模型作為反演結(jié)果。

以天然地震層析成像為例，作者設(shè)計(jì)了三個(gè)理論模型對算法進(jìn)行了驗(yàn)證，證實(shí)這種反演策略在參考模型與真實(shí)模型相差較遠(yuǎn)的情況下，對于反演較大型的地質(zhì)構(gòu)造有較好的適用性。

1 原理

在數(shù)據(jù)稀疏的情況下，獲得合理的反演結(jié)果，選用了如下的目標(biāo)函數(shù)［3］：

其中：m代表需要反演求解的模型向量：Ψ（m）為數(shù)據(jù)擬合項(xiàng)；Φ（m）為基于參考模型的模型約束項(xiàng)；Ω（m）為對模型進(jìn)行光滑約束的項(xiàng)，其中Φ（m）和Ω（m）的作用主要是正則化；ε和η分別代表阻尼因子和光滑因子，其作用是平衡Ψ（m）、Φ（m）和Ω（m）對目標(biāo)函數(shù)S（m）的貢獻(xiàn)，以使得反演過程能夠穩(wěn)定進(jìn)行，并獲得較為理想的模型解。

采用上述的目標(biāo)函數(shù)，其作用可以總結(jié)為兩點(diǎn)：①在觀測數(shù)據(jù)少、觀測系統(tǒng)無法完全覆蓋或觀測系統(tǒng)布置不合理的情況下，解決反演中的多解問題，獲得唯一的確定解；②加入了正則化項(xiàng)，可以解決反演過程中的不穩(wěn)定性。

若數(shù)據(jù)d與模型之間的關(guān)系為d＝g（m），式（1）中的第一項(xiàng)Ψ（m）具體形式為式（2）：

其中：dobs為觀測數(shù)據(jù)；第二項(xiàng)Φ（m）為模型約束項(xiàng)，形式為式（3）。

Φ（m）的作用是促使模型解m逼近于參考模型m0，在反演問題是欠定或混定問題的情況下，獲得唯一的確定解。由于各模型參數(shù)對數(shù)據(jù)的敏感程度不同，需要用一個(gè)權(quán)矩陣（模型協(xié)方差）Cm平衡每個(gè)模型參數(shù)誤差對預(yù)測總誤差Φ（m）的貢獻(xiàn)，起到對每個(gè)模型參數(shù)加權(quán)的作用。

目標(biāo)函數(shù)中最后一項(xiàng)的作用是定義模型參數(shù)之間的關(guān)系，以獲得最平滑模型解。此項(xiàng)的形式為式（4）。

其中：Dm是一個(gè)特定空間導(dǎo)數(shù)的有限差分估計(jì)；D是二階導(dǎo)數(shù)算子，因此模型越圓滑，Ω（m）會越小。

作者在使用阻尼最小二乘算法進(jìn)行理論模型反演時(shí)發(fā)現(xiàn)，當(dāng)參考模型選取合理時(shí)，迭代后的反演結(jié)果會逐漸趨近于真實(shí)模型?；谶@樣的考慮，在無法獲取合適的參考模型的情況下，將每次迭代后的結(jié)果經(jīng)過處理作為下次迭代的參考模型，經(jīng)過多次迭代后，可以獲取比較合適的參考模型，從而提高了反演的可靠性。圖1為參考模型逐步逼近真實(shí)模型的反演策略流程圖。在圖1中，“反演結(jié)果處理”項(xiàng)的目的是對當(dāng)前的反演結(jié)果進(jìn)行適當(dāng)?shù)奶幚?，從而獲取新的參考模型。

圖1 參考模型逐步逼近真實(shí)模型的反演策略流程Fig.1 Illustration of the inversion strategy to encourage reference model near to real model step by step

用反演結(jié)果估計(jì)參考模型的方法可以很多，其基本原則是在反演的最初階段，新的參考模型應(yīng)該是整個(gè)模型的背景場，而在反演的最后階段，參考模型應(yīng)該逐漸接近真實(shí)模型?；谏鲜鲈瓌t，可以對反演結(jié)果進(jìn)行多尺度平滑處理來估計(jì)出參考模型（式（5））。

式（5）表示以第j個(gè)模型參數(shù)mj為中心，取周圍n個(gè)模型參數(shù)做平均（下文稱n為平滑的階數(shù)），估計(jì)出參考模型的第j個(gè)參數(shù)，其中N為模型參數(shù)的總數(shù)。

在反演的開始階段，n的取值可以較大，超過反演的最大異常體尺度；隨著反演迭代的進(jìn)行，n的取值逐漸減小；在反演的最后階段，減小到與最小異常體尺度相當(dāng)或更小。另外，當(dāng)多種不同類型的模型參數(shù)（例如，地層界面和速度）同時(shí)反演時(shí)，在“反演結(jié)果處理”項(xiàng)中就需要分別對不同的模型參數(shù)進(jìn)行多尺度平滑處理。

2 算法驗(yàn)證

為了對提出的反演策略進(jìn)行驗(yàn)證，作者設(shè)計(jì)了一個(gè)均勻?qū)訝钅Ｐ秃蛢蓚€(gè)含異常體的地質(zhì)模型，進(jìn)行地震波速度層析成像實(shí)驗(yàn)。首先采用球坐標(biāo)系下波前快速推進(jìn)（Fast Marching Method，F(xiàn)MM）方法正演計(jì)算出理論模型的走時(shí)作為觀測數(shù)據(jù)，然后給定與真實(shí)模型相差較遠(yuǎn)的參考模型，最后采用FMM正演算法及子空間反演法進(jìn)行迭代反演［4］，在迭代過程中通過對參考模型進(jìn)行修正，使其逐步逼近真實(shí)模型。

2.1 均勻?qū)訝钅Ｐ?/p>

圖2為均勻?qū)訝罾碚撃Ｐ图坝^測系統(tǒng)（為了方便，剖面圖在直角坐標(biāo)系下顯示），模型范圍為地球東經(jīng)15°～45°E，北緯15°～45°N，深度0km～ 1 000km之間，分為4層，各層具體參數(shù)見表1。觀測系統(tǒng)的分布范圍為地球東經(jīng)20°～40°E和北緯20°～40°N之間，震源點(diǎn)及接收點(diǎn)的具體分布位置見表2。模型及觀測系統(tǒng)的覆蓋范圍均在2 000 km2以上，可以構(gòu)建大小不同的地質(zhì)異常體，從而驗(yàn)證算法對不同規(guī)模異常體的分辨率。

圖2 均勻?qū)訝罾碚撃Ｐ图坝^測系統(tǒng)Fig.2 Homogeneously stratified model and acquisition geometry

表1 理論模型參數(shù)Tab.1 Parameters of theoretical model

反演初始模型如圖3所示，四層的速度均為6.75km／s，初始參考模型（這里均將第一次反演迭代中輸入的參考模型命名為初始參考模型）與初始模型相同。圖4和圖5分別是參考模型沿北緯30°和東經(jīng)30°的剖面，兩圖分別顯示了六次反演迭代過程中參考模型的變化過程。對比分析發(fā)現(xiàn)，隨著迭代次數(shù)的增加，參考模型與真實(shí)模型越來越接近，最后一次的參考模型與真實(shí)模型相差不遠(yuǎn)。在反演過程中，為了簡化計(jì)算，對每次反演迭代后的結(jié)果均進(jìn)行三維五點(diǎn)平滑（5×5×5＝125個(gè)點(diǎn)）。

表2 觀測系統(tǒng)設(shè)計(jì)Tab.2 Design of acquisition geometry

圖6為改進(jìn)前算法與改進(jìn)后算法的反演結(jié)果與真實(shí)模型的速度差。對比圖6中算法改進(jìn)前后的結(jié)果可以看出，在初始參考模型與真實(shí)模型相差較遠(yuǎn)的情況下，其結(jié)果均向真實(shí)模型靠近，但在模型邊界區(qū)域兩種反演結(jié)果的分辨率都較差，其主要原因是由于邊界處穿過的射線較少，甚至沒有射線穿過。另外，采用未改進(jìn)算法反演出的結(jié)果中含有多個(gè)不真實(shí)異常體（如A、B、C、D、E等區(qū)域），有效區(qū)的最大誤差約為1.3km／s，最大相對誤差26%；改進(jìn)后算法反演出的結(jié)果消除了不真實(shí)異常體，有效區(qū)的最大誤差約為1km／s，最大相對誤差20%，而且較大的誤差主要集中于模型邊界區(qū)域。在有效反演區(qū)，改進(jìn)算法反演出的4層平均誤差，從上到下分別約為13%、4.84%、2.63%和6.1%，反演的整體速度結(jié)構(gòu)與真實(shí)模型基本一致。值得注意的是，中間兩層的誤差明顯小于第一層與第四層，分析其原因是由于給定的初始模型速度與中間兩層更加接近，而與其他兩層相差較遠(yuǎn)。雖然該反演方法準(zhǔn)確度較依賴于初始模型，但相對于改進(jìn)前的算法，改進(jìn)后算法的平均誤差大大降低，反演精度得到有效提高。

總的來說，在初始參考模型與真實(shí)模型相差較遠(yuǎn)的情況下，采用作者提出的反演策略及相應(yīng)的改進(jìn)算法，參考模型在反演過程中越來越逼近真實(shí)模型。相比常規(guī)算法，在一定程度上提高了反演結(jié)果的精度，這說明了這樣的反演策略及相應(yīng)的改進(jìn)算法是正確可行的。

圖3 初始模型及初始參考模型Fig.3 Initial model and initial reference model

2.2 含速度異常體模型

前述實(shí)驗(yàn)中的模型是比較簡單的均勻?qū)訝钅Ｐ?，而?shí)際的地球物理模型是十分復(fù)雜的。為了驗(yàn)證在有地球物理異常存在的情況下，我們所設(shè)計(jì)的反演策略與算法也能很好地工作，作者模擬月球質(zhì)量瘤的地球物理模型，進(jìn)行進(jìn)一步的算法實(shí)驗(yàn)。

月球質(zhì)量瘤是當(dāng)前月球探測中的熱點(diǎn)問題之一，其大小、分布及成因有助于研究月球的起源和演化等問題。質(zhì)量瘤是月球上高密度物質(zhì)聚集區(qū)，可能是由于月球形成早期月幔物質(zhì)上涌造成的［5－6］，其形狀近于圓盤狀。根據(jù)Konopliv［7］公布的質(zhì)量瘤分布圖，最大的質(zhì)量瘤直徑可達(dá)上千公里。相比于其他探測方法，月震層析成像是研究質(zhì)量瘤形態(tài)、規(guī)模及分布的最佳手段。但由于月球上僅有四個(gè)月震臺站的數(shù)據(jù)可用，臺站覆蓋密度小，而且相比于地震，月震發(fā)生次數(shù)非常少，射線覆蓋密度低。很多學(xué)者利用月震數(shù)據(jù)對月球進(jìn)行了大量研究，但至今仍未給出月球內(nèi)部可靠的結(jié)構(gòu)模型。基于這樣的考慮，作者設(shè)計(jì)了兩個(gè)大小不同的圓柱形速度異常體模型（本文分別稱為大型圓柱體及小型圓柱體）模擬質(zhì)量瘤，來進(jìn)行反演實(shí)驗(yàn)，驗(yàn)證改進(jìn)算法對不同規(guī)模異常體的分辨能力。

圖4 反演迭代過程中參考模型的變化（30°N剖面）Fig.4 Variation of reference model during inversion（cross section along 30°N）

這里使用的觀測系統(tǒng)與2.1節(jié)中的相同，理論模型除了在第一層中加入了一個(gè)大型圓柱形異常體以外，其他參數(shù)均與2.1節(jié)中的理論模型相同（圖7）。大型圓柱異常體的頂?shù)酌姘霃綖槌嗟捞?°（經(jīng)緯度），其中心位置位于30°N／30°E，頂部埋深為50km，底部位于第一層與第二層的界面處，異常體的速度與第二層的速度相同。這樣設(shè)計(jì)的目的是對月球質(zhì)量瘤進(jìn)行模擬，模型的第二層作為月幔，圓柱形異常體為一個(gè)大型質(zhì)量瘤，質(zhì)量瘤與月幔物質(zhì)速度相同，相當(dāng)于月幔物質(zhì)上涌形成的質(zhì)量瘤。

圖5 反演迭代過程中參考模型的變化（30°E剖面）Fig.5 Variation of reference model during inversion（cross section along 30°E）

圖6 未改進(jìn)及改進(jìn)后算法的成像結(jié)果在30°N的剖面（用相對速度表示，即成像速度與理論模型速度之差）Fig.6 Cross sections of imaging results along 30°N before and after improved algorithm （difference between imaging result and theoretical model）

反演初始模型及初始參考模型與2.1節(jié)中的相同，各層均為6.75km／s，采用改進(jìn)算法經(jīng)過6次迭代計(jì)算，得到圖8和圖9所示的結(jié)果。由圖8和圖9可以看出，在含有地質(zhì)異常體的情況下，每次反演過程中的參考模型逐漸向真實(shí)模型靠近，從最后一次輸出的參考模型，即本算法的反演結(jié)果來看（圖8 （f），圖9（f）），模型的形態(tài)及大小可以比較清晰地分辨出來，并且異常的邊界較為清晰，且在有效反演區(qū)，改進(jìn)算法反演出的4層平均誤差，從上到下分別約為15.1%、4.98%、2.61%和6.12%，這都說明改進(jìn)算法對于反演大型異常體，具有較好的反演效果。

圖10為小型圓柱異常體模型。小型圓柱形速度異常體的頂、底面半徑為赤道處3°，其他各參數(shù)均與大型圓柱異常體的相同。采用改進(jìn)算法經(jīng)過6次反演迭代的結(jié)果如圖11所示。由圖11可以看出，異常體下部的延伸情況比較清晰，但在切片圖上的邊界較模糊。

圖7 大型圓柱形異常體模型Fig.7 Model with large size of cylindrical anomaly

圖8 反演迭代過程中參考模型的變化（100km深度處的切片）Fig.8 Variations of reference model during inversion（slices at depth of 100km）

由兩種不同規(guī)模異常體的反演結(jié)果可以看出，改進(jìn)后的算法對較大規(guī)模的地質(zhì)構(gòu)造有較好的適用性；對于小型地質(zhì)構(gòu)造，其分辨率明顯降低，異常體邊界相對模糊，但仍能分辨出來。值得指出的是，在反演過程中，為了簡化計(jì)算，在模型2和模型3的反演計(jì)算中，每次迭代后獲取參考模型的處理均為三維五點(diǎn)平滑。

圖9 反演迭代過程中參考模型的變化（30°N剖面圖）Fig.9 Variations of reference model during inversion（cross section along 30°N）

圖10 小型圓柱形異常體理論模型Fig.10 Model with small size of cylindrical anomaly

3 結(jié)論與討論

1）在觀測數(shù)據(jù)稀疏、觀測系統(tǒng)不能完全覆蓋或觀測系統(tǒng)布置不合理的情況下，且先驗(yàn)信息也很少的區(qū)域，作者提出了通過在反演過程中對參考模型的修正，使參考模型逐步逼近真實(shí)模型，從而提高反演精度的新策略。

2）三個(gè)理想模型的數(shù)值實(shí)驗(yàn)，說明了在參考模型遠(yuǎn)離真實(shí)模型的情況下，相比于常規(guī)方法，采用這種方法的反演結(jié)果與真實(shí)模型更加接近，證明了新反演策略及相應(yīng)的算法是正確可行的。

3）新策略反演算法對那些數(shù)據(jù)稀疏，且只關(guān)心較大型地質(zhì)構(gòu)造的反演特別適用。其主要原因在于作者在改進(jìn)算法中加入了平滑處理以估計(jì)新的參考模型，當(dāng)異常體規(guī)模較小時(shí)，有可能會被平滑掉。在反演過程中，逐步地減少平滑的階數(shù)是必要的，可使小異常的邊界比較清晰。

圖11 小型圓柱形異常體反演成像結(jié)果Fig.11 Inverted result of the model with small size of cylindrical anomaly

［1］ 姚姚.地球物理反演基本理論與應(yīng)用方法［M］.北京：中國地質(zhì)大學(xué)出版社，2002.

YAO Y.Basic theory and method of geophysical inversion［M］.Beijing：China University of Geosciences Press，2002.（In Chinese）

［2］ 陳國金.井間地震層析成像中自動生成初始速度模型的方法研究［J］.地球物理學(xué)進(jìn)展，2007，22（6）：1831 －1835.

CHEN G J.Research of automatic generation of initial velocity model in crosswell seismic tomography［J］.Geophysical Prospecting For Petroleum，2007，22（6）：1831－1835.（In Chinese）

［3］ KENNETT B L N，SAMBRIDGE M S，WILLIAMSON PR.Subspace methods for large inverse problem with multiple parameter classes［J］.Geophysical Journal，1988，94：237－247.

［4］ JUPP D L B，VOZOFF K.Stable interative methods for inversion of geophysical data［J］.Geophysical Journal of the Royal Astronomical Society，1975，24：957 －976

［5］ RAWLINSON N.Fast marching tomography package ［EB／OL］.http：／／rses.anu.edu.au／seismology／fmmcode.

［6］ WISE D U，YATES M Y.Mascons as structural relief on a lunar“Moho”［J］.J.Geophys.Res.，1970，75：261－268.

［7］ 金丹.基于層析成像的阿波羅登陸地區(qū)月球質(zhì)量瘤空間分布研究［D］.武漢：中國地質(zhì)大學(xué)，2011.

JIN D.Spatial distribution of lunar mascons in apollo landing area using seismic tomography［D］.Wuhan：China University of Geosciences，2011.（In Chinese）

［8］ KONOPLIV A S，BINDER A B，HOOD L L，et al.Improved gravity field of the moon from lunar prospector.Science［J］.1998，281：1476－1480.

A strategy to improve inversion precision in target area with sparsely distributed observation

Since the observed data of the target area are very sparse，it is difficult to obtain the model which is close to real model only using the observed data.Hence，it is reasonable to add some priori information into inversion and apply the inversion method based on reference model.With the lack of observed data but also and priori information，it is also very difficult to select an appropriate reference model.If an appropriate reference model can＇t be selected，how do we invert for a reliable solution？This paper suggests a strategy to improve inversion precision，which encourages reference model near to real model by stepwise modifying reference model during inversion when the reference model is far away from real model.The authors designed a homogeneously stratified model and two models with the anomaly for testing algorithm.The inversion results illustrate that this strategy is appropriate for inversion of large size of anomaly.

inversion；seismic tomography；sparse observation；model modification method

YANG Ke－si1a，ZHU Pei－min2，ZHAO Na2，XU Yang1b

（1.Jilin Oilfield Branch，PetroChina a.Geophysical Exploration Research Institute，b.Xinmin Oil Production Factory，Songyuan 138000，China；2.China University of Geosciences，Institute of Geophysics and Geomatics，Wuhan 430074，China）

P 631.4

：A

10.3969／j.issn.1001－1749.2015.06.11

1001－1749（2015）06－0735－08

2014－11－28改回日期：2015－03－22

楊克思（1990－），男，碩士，主要研究方向?yàn)榈卣鹳Y料解釋與反演，E－mail：yangkesi＠foxmail.com。