基于光滑約束的磁梯度張量3D正則化反演方法

李金朋，張英堂，李志寧，范紅波，尹 剛

?

基于光滑約束的磁梯度張量3D正則化反演方法

李金朋，張英堂，李志寧，范紅波，尹 剛

(解放軍軍械工程學(xué)院七系，石家莊 050003)

針對(duì)鐵磁物質(zhì)反演方法中反演解的病態(tài)性問題，提出基于光滑約束的磁梯度張量3D正則化反演方法。首先，利用積分靈敏度矩陣和粗糙度矩陣對(duì)經(jīng)典Tikhonov正則化理論框架下的反演模型進(jìn)行約束，避免由于反演參數(shù)多于采集點(diǎn)數(shù)而導(dǎo)致的不穩(wěn)定；然后，通過定義正則化修正系數(shù)確定合理的初始正則化因子，從而減少因正則化因子引入而在反演結(jié)果中介入的誤差；最后，利用物性參數(shù)上下限約束函數(shù)，將反演過程中的解轉(zhuǎn)換到合理的值域范圍內(nèi)，得到與原始模型更加吻合的反演結(jié)果。仿真及實(shí)驗(yàn)結(jié)果表明：該反演方法能夠準(zhǔn)確還原磁性異常體的輪廓形態(tài)，具有較好的橫向和縱向分辨率。

光滑約束；磁梯度張量；積分靈敏度矩陣；粗糙度矩陣；正則化因子

磁性目標(biāo)反演技術(shù)是通過地面或航空等實(shí)測(cè)數(shù)據(jù)利用某種手段推算出磁性體(未爆彈、水雷、潛艇等)物性參數(shù)在地下的分布規(guī)律從而達(dá)到尋找磁性體的目的[1]。根據(jù)反演方法的不同將反演分為物性反演和形態(tài)反演兩類。其中基于磁場(chǎng)信息的物性反演方法因具有探測(cè)精度高、虛警率低、定位能力強(qiáng)等優(yōu)點(diǎn)得到了較為廣泛的關(guān)注，逐漸成為軍事偵察，水文及工程地質(zhì)勘探等領(lǐng)域的關(guān)鍵技術(shù)。磁梯度張量場(chǎng)是指磁場(chǎng)三分量在3個(gè)坐標(biāo)方向上的變化率，共有9個(gè)分量。相比磁總場(chǎng)、磁總場(chǎng)梯度和磁場(chǎng)三分量數(shù)據(jù)，磁梯度張量場(chǎng)具有更高的分辨率，能更好地描繪地下小尺度磁性體的空間形態(tài)及位置[2?6]。

在磁梯度張量數(shù)據(jù)三維反演方法中，反演結(jié)果的多解性主要是由于反問題的病態(tài)性引起的。TIKHONOV等[7?8]和ZHDANOV[9]指出在反演過程中通過加入正則化參數(shù)重新構(gòu)造目標(biāo)函數(shù)，從而達(dá)到獲得穩(wěn)定解的目的；WILSON等[10]利用獲得的航空重力梯度測(cè)量數(shù)據(jù)，采用共軛梯度法進(jìn)行反演；LI[11]提出了基于平滑約束的重力梯度張量數(shù)據(jù)反演方法，并應(yīng)用于實(shí)測(cè)資料，取得良好的效果；吳小平[12]提出將共軛梯度法應(yīng)用于電阻率三維反演中，并討論了不同正則化因子對(duì)反演結(jié)果的影響；朱自強(qiáng)等[13]利用混合正則化反演方法，減弱了反演的聚焦效應(yīng)，能有效獲得異常邊界；楊嬌嬌等[14]通過向目標(biāo)函數(shù)中引入深度加權(quán)函數(shù)來提高對(duì)目標(biāo)深度信息的識(shí)別能力。上述文獻(xiàn)研究方法為磁梯度張量反演提供了較好的研究思路，但是所得反演成像結(jié)果存在“趨膚效應(yīng)”，同時(shí)正則化因子的選取常采用經(jīng)驗(yàn)定值方法，導(dǎo)致計(jì)算精度較低，難以達(dá)到最佳的擬合效果。針對(duì)此問題，本文作者提出基于光滑約束的磁梯度張量3D正則化反演方法，首先，利用積分靈敏度矩陣和粗糙度矩陣解決反演過程中的病態(tài)解問題；然后，通過定義正則化修正系數(shù)減少介入誤差提高擬合精度；最后，引入物性參數(shù)上下限約束函數(shù)將反演結(jié)果轉(zhuǎn)換到合理的值域內(nèi)，得到與原始模型更加吻合的反演結(jié)果。

1 磁梯度張量三維反演理論

磁梯度張量正演公式表示為

式中：矩陣的元素為第(=1, 2, 3, …,)個(gè)觀測(cè)點(diǎn)觀測(cè)的第(=1, 2, 3, …,)個(gè)單元的響應(yīng)；維的行向量為任意張量分量數(shù)據(jù)；為反演物性參數(shù)。

根據(jù)Tikhonov正則化理論，構(gòu)建三維目標(biāo)函數(shù)：

(2)

2 光滑約束正則化理論

在反演過程中，由于重構(gòu)模型的核函數(shù)隨深度快速衰減，因此導(dǎo)致反演結(jié)果產(chǎn)生“趨膚效應(yīng)”。為了解決這一問題，利用積分靈敏度矩陣對(duì)模型進(jìn)行加權(quán)。觀測(cè)數(shù)據(jù)與模型參數(shù)的變化可以表示為

數(shù)據(jù)靈敏度對(duì)模型m的積分可以表示為

故積分靈敏度矩陣為對(duì)角矩陣：

(5)

(6)

對(duì)應(yīng)的Tiknonov正則化方程(2)變?yōu)?/p>

(7)

式中：為粗糙度矩陣的權(quán)重系數(shù)，數(shù)量級(jí)范圍在10?8~10?9之間。

利用共軛梯度法對(duì)目標(biāo)函數(shù)進(jìn)行求解，式(7)可以改寫為方程組：

(9)

具體算法實(shí)施步驟如下：

，，

2.1 自適應(yīng)正則化參數(shù)的選擇

本文在傳統(tǒng)的自適應(yīng)的正則化算法基礎(chǔ)上引入了正則化修正系數(shù)提高反演精度：

式中：正則化修正系數(shù)；為衰減因子；0為物性參數(shù)初值。ZHDANOV[9]指出衰減因子的取值范圍為[0.5, 0.9]，但是并未指出一種明確的衰減機(jī)制。

定義一個(gè)長(zhǎng)方體模型，尺寸為5 m×4 m×4 m,長(zhǎng)方體中心坐標(biāo)為(11, 11, 5)。測(cè)區(qū)范圍是21 m×21 m。在僅考慮感應(yīng)磁化條件下，異常體的磁化強(qiáng)度為40 A/m，磁傾角為70°，磁偏角為20°。設(shè)置最高迭代次數(shù)為60次。為確定自適應(yīng)正則化的各部分參數(shù)：

定義逼近誤差函數(shù)：

圖1(a)所示為各個(gè)分量在不同迭代次數(shù)下的逼近誤差。由圖1(a)可知，隨迭代次數(shù)增加，各分量的擬合誤差逐漸減小,其中和分量迭代過程相比其他分量較為穩(wěn)定，而且收斂速度較快。由圖1(b)可以看出，各分量值在10到50之間的模型逼近情況較好。圖1(c)可以看出，當(dāng)在[0.8, 0.9]區(qū)間時(shí)，模型逼近誤差小，擬合效率更好。綜合以上因素，設(shè)置為0.9，為10，選取進(jìn)行單分量反演。利用光滑約束進(jìn)行處理，處理結(jié)果如圖2(a)和圖2(c)所示。根據(jù)圖像可以看出，反演結(jié)果具有較好的橫向和縱向分辨率，與之相比如果不進(jìn)行光滑約束處理，異常體分布于地表附近，其縱向分辨率很差(見圖2(b)和2(d))。

2.2 物性參數(shù)分布的上下限約束函數(shù)

KIM等[15]通過對(duì)數(shù)約束實(shí)現(xiàn)了將電導(dǎo)率限制在合理的值域目標(biāo)中。CARDARELLI等[16]將對(duì)數(shù)轉(zhuǎn)換方法成功應(yīng)用于電阻率層析成像中[16]。本研究使用的上下限約束函數(shù)[17]如下所示：

式中：和分別為強(qiáng)制約束邊界的上下限。通過圖3可以看出，值反映出物性參數(shù)轉(zhuǎn)換的速度。值取0.5~1之間可以使反演結(jié)果較為合理的轉(zhuǎn)換到符合地球物理意義的數(shù)值范圍。處理實(shí)測(cè)數(shù)據(jù)時(shí)，將反演結(jié)果視為式(12)中的，通過該式可以實(shí)現(xiàn)反演數(shù)據(jù)的上下限約束。

圖1 各分量不同參數(shù)選取方案的數(shù)據(jù)統(tǒng)計(jì)

Fig. 1 Statistics of different parameters program for each component: (a) Data approximation error under different iterations; (b) Model approximation error under different regularization correction factors; (c) Model approximation error under different attenuation factors

圖2 單直立長(zhǎng)方體單分量反演結(jié)果切片圖

3 仿真分析

為驗(yàn)證本文方法的有效性，設(shè)計(jì)包含了兩個(gè)異常體的組合模型進(jìn)行三維測(cè)試。將地下待測(cè)空間劃分為21×21×10=4410個(gè)單元格，每個(gè)單元格均為邊長(zhǎng)為1 m的正方體，地面的觀測(cè)點(diǎn)為21×21=441個(gè)網(wǎng)格。在僅考慮感應(yīng)磁化的條件下異常體由兩個(gè)長(zhǎng)方體組成，其中長(zhǎng)方體1的中心坐標(biāo)為(6, 10.5, 4)，尺寸為5 m×4 m×4 m；長(zhǎng)方體2的中心坐標(biāo)為(13, 10.5, 4)，尺寸為9 m×8 m×6 m，假設(shè)磁化傾角為70°，磁化偏角為20°，磁化強(qiáng)度為40 A/m。

地下正演模型如圖4所示，圖4(a)所示黑色部分表示待測(cè)異常體，箭頭表示磁化方向，切片圖為模型在待測(cè)平面上產(chǎn)生的磁總場(chǎng)異常數(shù)據(jù)平面圖；圖4(b)表示待測(cè)異常體產(chǎn)生的理論磁梯度張量數(shù)據(jù),從磁梯度張量各分量云圖可以看出，分量正極值與待測(cè)異常體水平位置具有良好的對(duì)應(yīng)關(guān)系，而且通過前面討論可知分量在計(jì)算過程中具有良好的穩(wěn)定性，由此，本文作者利用磁梯度張量數(shù)據(jù)并加入5%的高斯噪聲作為實(shí)測(cè)數(shù)據(jù)進(jìn)行反演。反演的物性參數(shù)分布在不同切面處的分布特點(diǎn)如圖5所示，反演過程中設(shè)置最高迭代次數(shù)為55次，為50，取0.9。黑色框所圈定的范圍是正演組合模型體所在位置及分布范圍，從圖5中結(jié)果可以看出：反演方法能夠較好地反映長(zhǎng)方體組合模型的輪廓形態(tài)，具有較好的橫向和縱向分辨率。

圖3 上下限約束特征中參數(shù)變化對(duì)轉(zhuǎn)換空間值的影響

圖4 組合模型三維空間形態(tài)及磁梯度張量數(shù)據(jù)平面圖

圖5 組合異常體反演切片圖及三維姿態(tài)

4 實(shí)測(cè)資料驗(yàn)證

由于地理環(huán)境和測(cè)量誤差的影響，會(huì)導(dǎo)致實(shí)際測(cè)量的磁梯度張量數(shù)據(jù)精度下降。選擇垂向和水平分辨率相對(duì)較高的分量進(jìn)行方法的有效性驗(yàn)證。測(cè)區(qū)位于石家莊某地，側(cè)區(qū)內(nèi)放置一個(gè)南北走向的未爆彈，其直徑約20 cm，沿軸坐標(biāo)為0.5~1.3 m，未爆彈的軸線距離觀測(cè)面0.5 m。圖6展示了測(cè)量臺(tái)架及傳感器陣列圖，將待測(cè)空間劃分為22×22×10=4840個(gè)單元格，每個(gè)單元格均為邊長(zhǎng)為0.1 m的正方體，地面的觀測(cè)點(diǎn)為22×22=484個(gè)網(wǎng)格。圖7所示為未爆彈的實(shí)測(cè)磁梯度張量數(shù)據(jù)，反演過程中設(shè)置最高迭代次數(shù)為55次，為50，取0.9。

反演結(jié)果如圖8所示，根據(jù)反演結(jié)果可以看出，本方法所得反演結(jié)果與未爆彈的位置及輪廓相吻合，進(jìn)一步證明本方法的有效性。

圖6 實(shí)測(cè)區(qū)域及傳感器陣列

圖7 實(shí)測(cè)未爆彈磁梯度張量數(shù)據(jù)

圖8 反演結(jié)果圖

5 結(jié)論

1) 提出的光滑約束正則化反演方法能很好的壓制反演結(jié)果的“趨膚效應(yīng)”，提高3個(gè)方向的分辨率，使反演結(jié)果呈現(xiàn)出平滑結(jié)構(gòu)。

2) 提出改進(jìn)的自適應(yīng)正則化方法，通過確定合理的正則化因子，提高反演的精度及穩(wěn)定性。

3) 引入上下限約束函數(shù)，深入討論了參數(shù)的選擇及設(shè)定的原則，該函數(shù)可以將不具有實(shí)際物理意義的反演結(jié)果轉(zhuǎn)化到合理的范圍，給進(jìn)一步解釋工作帶來方便。

4) 通過仿真及實(shí)驗(yàn)驗(yàn)證，證明利用共軛梯度法對(duì)目標(biāo)函數(shù)進(jìn)行求解，能夠準(zhǔn)確地反演出異常體的空間姿態(tài)及物性參數(shù)特征。

REFERENCES

[1] 陳召曦, 孟小紅, 郭良輝. 重磁數(shù)據(jù)三維物性反演方法進(jìn)展[J]. 地球物理學(xué)進(jìn)展, 2012, 27(2): 503?511. CHEN Zhao-xi, MENG Xiao-hong, GUO Liang-hui. Review of 3D property inversion of gravity and magnetic data[J] . Progress in Geophysics, 2012, 27(2): 503?511.

[2] 卞光浪, 翟國(guó)軍, 黃謨濤. 顧及地磁背景場(chǎng)的多目標(biāo)磁異常分量換算方法[J]. 武漢大學(xué)學(xué)報(bào), 2011, 36(8): 914?918. BIAN Guang-lang, ZHUAI Guo-jun, HUANG Mo-tao. Transformation of multiobjects magnetic anomaly components with geomagnetic effect[J]. Geomatics and Information Science of Wuhan University, 2011, 36(8): 914?918.

[3] 張昌達(dá). 航空磁力梯度張量測(cè)量—航空磁測(cè)技術(shù)的最新進(jìn)展[J]. 工程地球物理學(xué)報(bào), 2007, 3(5): 354?361. ZHANG Chang-da. Airborne tensor magnetic gradiometry—the latest progress of airborne magnetomeric technology[J]. Chinese Journal of Engineering Geophysics, 2007, 3(5): 354?361.

[4] 石 磊, 郭良輝, 孟小紅. 磁總場(chǎng)異常垂直梯度三維相關(guān)成像[J]. 地球物理學(xué)進(jìn)展, 2012, 27(4): 1069?1614. SHI Lei, GUO Liang-hui, MENG Xiao-hong. 3D correlation imaging of the vertical gradient of magnetic total field anomaly[J]. Progress in Geophysics, 2012, 27(4): 1069?1614.

[5] 吳招才. 磁力梯度張量技術(shù)及其應(yīng)用研究[D]. 武漢: 中國(guó)地質(zhì)大學(xué), 2008. WU Zhao-cai. Magnetic gradient tensor technology and its application[D]. Wuhan: China University of Geosciences, 2008.

[6] 孟 慧. 磁梯度張量正演、延拓、數(shù)據(jù)解釋方法研究[D]. 吉林: 吉林大學(xué), 2012. MENG Hui. Forward modeling, continuation and data interpretation of magnetic gradient tensor[D]. Jilin: Jilin University, 2012.

[7] TIKHONOV A N, ARSENIN V Y. Solution of Ill-posed problems[M]. New York: V.H. Winston & Sons, 1997.

[8] TIKHONOV A N. Regularization of Ill-posed problems[J]. Doklady Akad, 1963, 153: 49?52.

[9] ZHDANOV M S. Geophysical inverse theory and regularization problems[M]. New York: Elsevier Science, 2002.

[10] WILSON G, CUMA M, ZHDANOV M S. Massively parallel 3D inversion of gravity and gravity gradiometer data[J]. Preview, 2011, 152(6): 29?34.

[11] LI Y. 3-D inversion of gravity gradiometer data[C]// 2001 SEG Annual Meeting, SEG Technical Program Expanded Abstracts. San Antonio: Society of Exploration Geophysicists, 2001: 1470?1473.

[12] 吳小平. 利用共軛梯度法的電阻率三維反演若干問題研究[J]. 地震地質(zhì), 2001, 23(12): 321?327. WU Xiao-ping. Study on some problems for 3d resistivity inversion using conjugate gradient[J]. Seismology and Geology, 2001, 23(12): 321?327.

[13] 朱自強(qiáng), 曹書錦, 魯光銀. 基于混合正則化的重力場(chǎng)約束反演[J]. 中國(guó)有色金屬學(xué)報(bào), 2014, 24(10): 2601?2608. ZHU Zi-qiang, CAO Shu-jin, LU Guang-yin. 3D gravity inversion with bound constraint based on hyper-parameter regularization[J]. The Chinese Journal of Nonferrous Metals, 2014, 24(10): 2601?2608.

[14] 楊嬌嬌, 劉 展, 陳曉紅, 徐凱軍. 基于深度加權(quán)的重力梯度張量數(shù)據(jù)的3D聚焦反演[J]. 世界地質(zhì), 2015, 34(1): 210?218. YANG Jiao-jiao, LIU Zhan, CHEN Xiao-hong, XU Kai-jun. Three-dimensional focusing inversion of gravity gradient tensor data based on depth weighting[J]. Global Geology, 2015, 34(1): 210?218.

[15] KIM H J, SONG Y, LEE K H. Inequality constraint in least-squared inversion of geophysical data[J]. Earth Planets Space, 1999, 51: 255?259.

[16] CARDARELLI E, FISCHAGER F. 2D data modelling by electrical resistivity tomography for complex subsurface geology[J]. Geophysical Prospecting, 2006, 54(2): 121?133.

[17] COMMER M, NEWMAN G A. New advances in three- dimensional controlled-source electromagnetic inversion[J]. Geophysical Journal International, 2008, 172(2): 513?535.

(編輯 王 超)

3D magnetic gradient tensor regularization inversion method based on smooth constrained

LI Jin-peng, ZHANG Ying-tang, LI Zhi-ning, FAN Hong-bo, YIN Gang

(Seventh Department, Ordnance Engineering College, Shijiazhuang 050003, China)

For the ill-posed problem of ferromagnetic materials inversion method, the method based on 3D magnetic gradient tensor regularization inversion method was proposed based on smooth constraint condition. First, the unstable defect caused by the number of inverse parameter far than observational points by using the integral sensitivity matrix and the roughness matrix in the classical Tikhonov regularization under the framework of the theory to constrain the model is avoided. Then, errors are reduced due to regularization factor introduced in the inversion results, by defining regularization correction factor to determine a reasonable initial regularization factor. Finally, by using the upper and lower constraint function of physical parameters, the solution to a reasonable value range in the inversion process is converted and the inversion results which are more consistent with the original model are obtained. Simulation and experimental results show that this inversion method can reflect the outline shade of the magnetic anomaly, and has good lateral and vertical resolution.

smooth constraint; magnetic gradient tensor; integral sensitivity matrix; roughness matrix; regularization factor

Project(51305454) supported by the National Natural Science Foundation of China

2015-11-23; Accepted date:2016-04-17

ZHANG Ying-tang; Tel: +86-311-87994748; E-mail: zyt01@mails.tsinghua.edu.cn

1004-0609(2016)-12-2640-07

P631

A

國(guó)家自然科學(xué)基金資助項(xiàng)目(51305454)

2015-11-23；

2016-04-17

張英堂，教授，博士；電話：0311-87994748；E-mail: zyt01@mails.tsinghua.edu.cn