杜興忠，朱海東

（中國電建集團(tuán)貴陽勘測設(shè)計研究院有限公司，貴陽 550081）

航空瞬變電磁正演中頻時轉(zhuǎn)換方法分析及系數(shù)的選取

杜興忠，朱海東

（中國電建集團(tuán)貴陽勘測設(shè)計研究院有限公司，貴陽 550081）

在瞬變電磁法正演計算中，常常需要頻時轉(zhuǎn)換來獲得垂直磁場隨時間的變化率。目前已經(jīng)提出了多種頻時轉(zhuǎn)換方法，其中G－S變換法只涉及實數(shù)計算，具有較快的速度，因此已經(jīng)成為航空瞬變電磁正演模擬的候選方法。由于航空瞬變電磁數(shù)據(jù)巨大，時間數(shù)較多，因此在保證精度的條件下怎樣提高計算速度成為當(dāng)前亟待解決的問題。鑒于大量頻點的正弦變換和余弦變換具有較高的精度，這里使用三層和四層典型地電模型，通過比較余弦變換、正弦變換和不同點數(shù)G－S變換的正演計算結(jié)果來確定最佳的G－S變換系數(shù)。實驗結(jié)果比較和分析認(rèn)為，采用8個系數(shù)的G－S變換法不僅有較快的計算速度，而且也有較高的計算精度。

航空瞬變電磁法；G－S變換；正弦變換；余弦變換

0 引言

在瞬變電磁法中，常常需要通過頻時轉(zhuǎn)換來實現(xiàn)時間域的電磁響應(yīng)計算。通常的頻時轉(zhuǎn)換法有G －S變換法［1－2］、余弦變換多項式近似法［3］、余弦變換法［4］、正弦變換法［4－5］、數(shù)值線性濾波法［6］、滯后褶積法［7］、延遲譜法［8］、混合算法［9－10］等。盡管在上述方法的基礎(chǔ)上，許多研究者也提出了相應(yīng)的改進(jìn)方法，羅延鐘［11］基于拉氏變換的延遲定理，建立了一種新的快速Gaver－Stehfest拉氏逆變換；阮百堯［12］提出了改進(jìn)的Guptasarma算法并用于瞬變電磁的正演計算中；為了提高計算效率，蔣淑芬［13］等提出了一種基于復(fù)積分的正余弦變換算法；徐振平等［14］給出了改進(jìn)的余弦變換方法，但是都未能在計算效率上有顯著改善。昌彥君［15］等通過比較余弦變換多項式近似法、G－S變換法和Guptasarma濾波算法，認(rèn)為Guptasarma濾波算法計算速度快，但是晚期誤差范圍大，適用的時間范圍最小；余弦變換法精度高，適用的時間范圍最寬，但是耗時最多；G －S變換法的計算精度和耗時介于兩者之間［15］。此外，孫業(yè)發(fā)［5］等通過比較正弦變換、余弦變換和余弦變換的折線逼近法，最后認(rèn)為正弦變換在瞬變電磁響應(yīng)晚期的轉(zhuǎn)換精度明顯高于其他兩種方法。

由于航空瞬變電磁法數(shù)據(jù)龐大，涉及該方法的正演計算耗時較大，為此作者以大量頻點（300個）的正弦變換法和余弦變換法為基準(zhǔn)，通過與不同系數(shù)的G－S變換法的正演計算結(jié)果進(jìn)行比較，確定出G－S變換法的最佳系數(shù)，以保證航空瞬變電磁正演計算在擁有較高計算精度的條件下，具有較高的計算速度。

1 航空瞬變電磁法正演計算

航空瞬變電磁法是以飛機為運載工具，發(fā)射脈沖電流，在發(fā)射機關(guān)斷電流時，觀測二次磁場隨時間的衰減特性，進(jìn)而推斷地下不同介質(zhì)的分布，從而解決各類地質(zhì)問題的方法［16］。

圖1 航空瞬變電磁正演模型Fig.1 The forward model of airborne transient electromagnetic

對于圖1所示的航空瞬變電磁正演模型，中心回線裝置的航空瞬變電磁響應(yīng)計算公式表示為［17］式（1）：

其中：J1表示一階第一類貝塞爾函數(shù)；a表示發(fā)射線圈半徑；λ是積分變量；RTE表示反射系數(shù)，具體為：

并且

式中：ω是角頻率；kn是第n層的波數(shù)；σn是第n層的電導(dǎo)率；通常取大地的導(dǎo)磁率μn和自由空間的導(dǎo)磁率μ0相等，在準(zhǔn)靜態(tài)情況下，通常取u0＝λ。

根據(jù)上述公式可以實現(xiàn)航空瞬變電磁法一維正演模擬，其計算步驟如下：

1）設(shè)置系統(tǒng)參數(shù)和地電模型參數(shù)。

2）計算反射系數(shù)。

3）漢克爾變換計算。

4）頻時轉(zhuǎn)換計算。

5）計算電磁響應(yīng)。

在上述計算步驟，影響電磁響應(yīng)精度和計算效率的關(guān)鍵是漢克爾變換和頻時轉(zhuǎn)換。

漢克爾變換數(shù)值濾波計算方法有很多，這里采用Guptasarma和Singh［18］所提出的數(shù)字濾波方法。

漢克爾積分變換的一般形式可以表示為式（6）：

其中：K（λ）表示待變換的函數(shù)；Ji表示i階第一類貝塞爾函數(shù)；r是一個已給出的定值；λ表示被積函數(shù)。數(shù)字濾波器有一組濾波系數(shù)Wi和兩個常數(shù)（位移a和采樣間隔s）。位移a決定了輸入函數(shù)采樣的起始點，間隔s決定抽樣時的時間間隔。文獻(xiàn)［18］給出了漢克爾變換的計算公式：

一般常用的漢克爾變換濾波系數(shù)有安德森［19］引入的283個系數(shù)和801個系數(shù)，盡管其精度較高，但是計算速度較慢，所以這里主要采用文獻(xiàn)［18］中給出的濾波系數(shù)，即零階有61或120個系數(shù)，一階有47或120個系數(shù)。

2 頻時轉(zhuǎn)換

2.1 正弦、余弦變換

根據(jù)航空瞬變電磁響應(yīng)計算公式（1）可以求得Hz。令激發(fā)波為階躍波：

對式（8）進(jìn)行傅氏變換，可得式（9）。

將式（8）代入式（9）再經(jīng)過傅氏反變換，得到式（10）與式（11）。

根據(jù)正余弦變換數(shù)值濾波算法［4］，式（10）和式（11）可以分別化為：

2.2 G－S變換

Gaver－Stehfest逆拉氏變換（簡稱G－S變換）是純實數(shù)運算，只需對較少的拉氏變換量s值做計算，因而計算速度較快。根據(jù)文獻(xiàn)［1］，可以獲得用G－S變換作逆拉氏變換的算法為：對給定的時間t，時間域的響應(yīng)值f（t），由拉氏變換域中的變量（其中m＝1，2，…，n）的拉氏域響應(yīng)值F（sm）算出：

表1 G－S變換系數(shù)Tab.1 G－S transformation coefficients

式中：n是決定于計算機位數(shù)的正偶整數(shù)；Km是G －S變換系數(shù)，且可由公式（15）計算得出。

式中：m＝1，2，…，n；求和下限i1是的整數(shù)部分。

在航空瞬變電磁響應(yīng)計算中，為了在保證計算精度的條件下獲得較高的計算效率，作者考查濾波系數(shù)個數(shù)分別取n＝6、8、10、12、14、16時所對應(yīng)的響應(yīng)值。根據(jù)式（15）求得不同濾波系數(shù)個數(shù)所對應(yīng)的G－S變換系數(shù)（表1）。

對于G－S變換系數(shù)個數(shù)n的取值問題，目前存在多種不同的觀點，Wooden［20］等認(rèn)為n取8或者10比較合適，當(dāng)n取得更大時，會使得變換變得不穩(wěn)定，數(shù)據(jù)結(jié)果會發(fā)生振蕩；羅延鐘［21］等認(rèn)為n取16合適；昌彥君和張桂青等［15］則認(rèn)為n取12合適。作者結(jié)合航空瞬變電磁數(shù)據(jù)特性，通過正演模擬計算確定G－S變換系數(shù)個數(shù)n的取值，同時和使用正、余弦變換的響應(yīng)計算結(jié)果進(jìn)行對比分析，驗證n的取值的正確性。

3 電磁響應(yīng)計算結(jié)果及分析

3.1 G－S變換系數(shù)對應(yīng)的電磁響應(yīng)對比實驗

為了比較不同數(shù)量的G－S變換系數(shù)對應(yīng)的電磁響應(yīng)，設(shè)置系統(tǒng)參數(shù)如下：發(fā)射電流為階躍電流，其歸一化值為1A，發(fā)射線圈半徑為1m；接收線圈有效面積歸一化為1m2；飛行高度為30m。設(shè)置三層K型地電模型的參數(shù)為：第一層電阻率為10Ω· m，層厚為50m，第二層電阻率為100Ω·m，層厚為50m，第三層的電阻率為50Ω·m。在上述系統(tǒng)參數(shù)和地電模型的參數(shù)設(shè)置下，使用表1中不同的G－S變換系數(shù)，通過正演計算，獲得的電磁響應(yīng)如圖2所示。

從圖2中可以看出，不同系數(shù)的G－S變換所對應(yīng)的電磁響應(yīng)差異主要在晚期，說明G－S變換系數(shù)對航空瞬變電磁階躍響應(yīng)的晚期影響較大，而早期響應(yīng)幾乎不受影響。此外，6點G－S變換和16點G－S變換的電磁響應(yīng)在晚期發(fā)生振蕩，可以認(rèn)為6點G－S變換和16點G－S變換不適用于航空瞬變電磁響應(yīng)計算，實際上，我們在實驗中也發(fā)現(xiàn)，小于6點或大于16點的G－S變換系數(shù)在電磁響應(yīng)晚期都發(fā)生振蕩，因此對于航空瞬變電磁正演模擬，只余8點、10點、12點和14點G－S變換可以考慮作為候選頻時轉(zhuǎn)換方法。同時從圖2中也發(fā)現(xiàn)，10點、12點、14點G－S變換的階躍響應(yīng)值幾乎一致，而8點G－S變換的階躍響應(yīng)值和這3種G －S變換在晚期存在一定差異。

圖2 不同系數(shù)個數(shù)的G－S變換所對應(yīng)的階躍響應(yīng)Fig.2 The step response values corresponding to the G－S transform of different coefficients

對于航空瞬變電磁響應(yīng)計算，我們在8點、10點、12點和14點G－S變換進(jìn)行選擇，分別將8點、10點、12點、14點G－S變換和余弦變換的K型地電模型階躍響應(yīng)值進(jìn)行對比（圖3）。從圖3中可以看出，在早期，8～14點G－S變換和余弦變換所分別對應(yīng)的航空瞬變電磁響應(yīng)的形態(tài)和趨勢基本一致，只是幅值上稍有差異；但是在晚期，10～14點G －S變換和余弦變換所對應(yīng)的電磁響應(yīng)則有很大差異，而8點G－S變換和余弦變換的電磁響應(yīng)比較接近，而當(dāng)頻點足夠多時，余弦變換的晚期響應(yīng)值更準(zhǔn)確和穩(wěn)定并且具有更寬的時間范圍［15］，這里假設(shè)8點G－S變換更能適用于航空瞬變電磁法。

3.2 G－S變換和余弦變換實驗對比

為了驗證8點G－S變換選取的正確性和適用性，這里考慮三層地電模型A型、H型、Q型和四層地電模型AA型、AK型、KH型、KQ型、QQ型、QH型、HA型和HK型，并對比8點G－S變換、10點G－S變換和余弦變換得到的階躍響應(yīng)值。系統(tǒng)參數(shù)的設(shè)置同上，地電模型的設(shè)置為：A型地電模型的電阻率從上到下依次為1 0Ω·m、5 0Ω·m、1 0 0 Ω·m；H型地電模型的電阻率從上到下依次為50 Ω·m、10Ω·m、100Ω·m；K型地電模型的電阻率從上到下依次為10Ω·m、100Ω·m、50Ω·m；Q型地電模型的電阻率從上到下依次為100 Ω·m、50Ω·m、10Ω·m；AA型地電模型的電阻率從上到下依次為10Ω·m、50Ω·m、100Ω·m、200Ω·m；AK型地電模型的電阻率從上到下依次為10Ω·m、50Ω·m、100Ω·m、50Ω·m；KH型地電模型的電阻率從上到下依次為10Ω·m、50 Ω·m、10Ω·m、50Ω·m；KQ型地電模型的電阻率從上到下依次為10Ω·m、100Ω·m、50Ω·m、10Ω·m；QQ型地電模型的電阻率從上到下依次為200Ω·m、100Ω·m、50Ω·m、10Ω·m；QH型地電模型的電阻率從上到下依次為100Ω·m、50Ω·m、10Ω·m、50Ω·m；HA型地電模型的電阻率從上到下依次為50Ω·m、10Ω·m、50 Ω·m、100Ω·m；HK型地電模型的電阻率從上到下依次為50Ω·m、10Ω·m、50Ω·m、10Ω·m。每種地電模型各層之間的層厚都為50m。在典型的三層地電模型和四層地電模型中所獲得的電磁響應(yīng)值分別如圖4和圖5所示。

圖3 K型地電模型中G－S變換和余弦變換的階躍響應(yīng)值Fig.3 The step response values of G－S transform and cosine transform in K Geoelectric model

圖4 8、10點G－S變換和余弦變換在三層地電模型中的階躍響應(yīng)值Fig.4 The step response values of 8，10points G－S transform and cosine transform in three－layer model

圖5 8、10點G－S變換和余弦變換在四層地電模型中的階躍響應(yīng)值Fig.5 The step response values of 8，10points G－S transform and cosine transform in four－layer model

圖6 8、10點G－S變換與余弦變換和正弦變換在K型和H型上的階躍響應(yīng)Fig.6 The step response values of 8，10points G－S transform，cosine transform and sine transform in K－type and H－type geoelectric model respectively （a）K型；（b）H型

在圖4和圖5中，藍(lán)色實線為余弦變換得到的階躍響應(yīng)曲線，綠色實線為8點G－S變換得到的階躍響應(yīng)曲線，紅色虛線為10點G－S變換得到的階躍響應(yīng)曲線。從A型、K型、AA型、AK型地電模型圖中可以看出，在晚期的時候綠色實線較紅色虛線明顯更靠近藍(lán)色實線，而在其他幾種地電模型圖中綠色實線和紅色虛線幾乎重合，這說明8點和10點G－S變換在一些模型中，其階躍響應(yīng)基本一致；而在另一些模型中，8點G－S變換所得到的階躍響應(yīng)在晚期更接近于余弦變換所得到的階躍響應(yīng)，這就證明了8點G－S變換適用于航空瞬變電磁法，因為其不僅擁有較高的計算精度，同時能獲得更快的計算速度。

3.3 G－S變換、余弦和正弦變換實驗對比

孫業(yè)發(fā)［5］等認(rèn)為正弦變換在瞬變電磁響應(yīng)晚期的轉(zhuǎn)換精度明顯高于余弦變換和余弦變換的折線逼近法。為了考查8點G－S變換的轉(zhuǎn)換精度，這里以三層K型和H型地電模型為例，對8點G－S變換、10點G－S變換、余弦變換和正弦變換的電磁響應(yīng)進(jìn)行比較。在地電模型的參數(shù)設(shè)置中，K型模型的電阻率由上至下分別為：10Ω·m、100Ω·m、50 Ω·m，H型模型的電阻率由上至下分別為：100 Ω·m、10Ω·m、100Ω·m，兩種地點模型各層之間的層厚均為50m。在K型和H型地電模型上，8 點G－S變換、10點G－S變換、余弦變換和正弦變換的電磁響應(yīng)如圖6所示。

從圖6中可以看出，正弦變換得到的響應(yīng)結(jié)果和余弦變換、G－S變換得到的響應(yīng)結(jié)果在早期走勢基本一致，而在晚期它們之間存在一定差異?？梢钥吹皆谕砥诘臅r候8點G－S變換（綠色實線）比10點G－S變換（粉紅色虛線）更靠近正弦變換（紅色實線）或者余弦變換（藍(lán)色實線），驗證了8點GS變換具有較高的計算精度。

從運算時間上來說，在同等計算機平臺上（CPU：i3，內(nèi)存：3G；操作系統(tǒng)：Windows 7；語言：Matlab 2010b），G－S變換耗時為0.063s，而正弦變換和余弦變換耗時分別為4.768s和4.766s。這說明正、余弦變換盡管有較寬的時間范圍，但是其計算量太大，而G－S變換具有較高的計算效率，這是由于G－S變換系數(shù)較少，且不涉及復(fù)數(shù)運算，而正、余弦變換需要進(jìn)行復(fù)數(shù)計算。

綜合上述的計算精度和計算時間分析，可以看到，在計算量較大的航空瞬變電磁法中，8點G－S變換比較適用，因為其不僅具有較高的計算精度，而且還有很高的計算效率。

4 結(jié)論

在典型的三層和四層地電模型中，對8點、10點、12點、14點、16點G－S變換以及正弦變換和余弦變換的階躍響應(yīng)進(jìn)行對比分析，獲得了如下結(jié)論：

1）當(dāng)需求較高的轉(zhuǎn)換精度時，正弦變換和余弦變換所需的頻點數(shù)較多，因此導(dǎo)致計算量增大，是G －S變換耗時的幾十倍。

2）在不同點數(shù)的G－S變換中，8點G－S變換具有較高的轉(zhuǎn)換精度，僅次于高頻點的正弦變換和余弦變換。

3）在計算效率上，和其他的G－S變換以及正、余弦變換相比，8點G－S變換具有最高的計算效率。

4）由于8點G－S變換在正演計算中不僅具有較快的計算速度，而且也具有較高的計算精度，所以在數(shù)據(jù)量巨大的航空瞬變電磁正演計算中，應(yīng)該選擇其作為頻時轉(zhuǎn)換方案。

［1］ KNIGHT J H，RAICHE A P.Transient electromagnetic calculations using the Gaver－Stehfest inverse Laplace transform method［J］.GEOPHYSICS，1982，47（1）：47－50.

［2］ 樸華榮.電磁測深法原理［M］.北京：地質(zhì)出版社，1990.

PIAO H R.Principle of electromagnetic sounding method［M］.Beijing：Geological Publishing，1990.（In Chinese）

［3］ 考夫曼A A，凱勒G V.頻率域和時間域電磁測深［M］.地質(zhì)出版社，1987.

KAUFMAN A A，KELLER G V.Frequency and transient soundings［M］.Geology Press，1987.（In Chinese）

［4］ 王華軍.正余弦變換的數(shù)值濾波算法［J］.工程地球物理學(xué)報，2004，1（4）：329－335.

WANG HJ.Digital filter algorithm of the sine and cosine transform［J］.Chinese Journal of Engineering Geophysics，2004，1（4）：329－335.（In Chinese）

［5］ 孫業(yè)發(fā)，李桐林，范翠松，等.中心回線瞬變電磁響應(yīng)計算的頻－時域轉(zhuǎn)換方法研究［J］.地球物理學(xué)進(jìn)展，2014，29（3）：1243－1247.

SUN Y F，LI T L，F(xiàn)AN C S，et，al.Research on frequency－time domain transform method in centralloop transient electromagnetic field response computation［J］.Progress in Geophysics，2014，29（3）：1243 －1247.（In Chinese）

［6］ GUPTASARMA D.Computation of the time－domain response of a polarizable ground［J］.GEOPHYSICS，1982，47（11）：1574－1576.

［7］ ANDERSON W L.Fast Hankel transforms using related and lagged convolutions［J］.ACM TRANS.ON MATH.SOFTWARE，1982（8）：344－368.

［8］ NEWMAN G A，HOHMANN G W，ADNERSON W L.Transient electromagnetic response of a three－dimensional body in a layered earth［J］.GEOPHYSICS，1986，51（8）：1608－1626.

［9］ GOLDMAN M M.The integral－finite－difference method for calculating transient electromagnetic fields in a horizontally stratified medium［J］.GEOPHYSICAL PROSPECTING，1983，31：664－686.

［10］李汝傳，王書民，孫鴻雁.用混合算法計算電偶源瞬變電磁響應(yīng)［J］.物探與化探，2002，26（3）：215－217.

LI R CH，WANG SH M，SUN H Y.The application of mixed algorithm to the Calculation of transient electromagnetic response of couple source.Geophysical ＆Geochemical Exploration，2002，26（3）：215－217.（In Chinese）

［11］羅延鐘，昌彥君.G－S變換的快速算法［J］.地球物理學(xué)報，2000，43（5）：684－690.

LUO Y ZH，CHANG Y J.A rapid algorithm for GS transform.Ch inese Journal of Geophysics，2000，43 （5）：684－690.（In Chinese）

［12］阮百堯.Guptasarma算法在瞬變電磁正演計算中的應(yīng)用［J］.桂林工學(xué)院學(xué)報，1996，16（2）：167－170.RUAN BY.Applied of Guptasarma algorithm in TEM ＇s forward［J］.Journal of Guilin Institute of Technology，1996，16（2）：167－170.（In Chinese）

［13］蔣淑芬，向淑晃.一種正余弦變換的高效算法［J］.工程地球物理學(xué)報，2007，4（5）：512－515.

JIANG SH F，XIANG SH H.An Eifficient Method of Calculating the Sine and Cosine Transform.Chinese Journal of Engineering Geophysics，2007，4（5）：512 －515.（In Chinese）

［14］徐振平，胡文寶.基于快速余弦變換的低頻電磁場數(shù)值濾波法正演［J］.石油天然氣學(xué)報（江漢石油學(xué)院學(xué)報），2011，33（1）：63－67.

XU ZH P，HU W B.Low frequency electromagnetic field numerical filtering method based on fast cosine transform［J］.Journal of Oil and Gas Technology，2011，33（1）：63－67.（In Chinese）

［15］昌彥君，張桂青.電磁場從頻率域轉(zhuǎn)換到時間域的幾種算法比較［J］.物探化探計算技術(shù)，1995，17（3）：25 －29.

CHANG Y J，ZHANG G Q.Comparison among three transformation algorithms of electromagnetic field from frequency domain to time domain.Computing Techniques for Geophysical and Geochemical Exploration，1995，17（3）：25－29.（In Chinese）

［16］NABIGHIAN M N.Electromagnetic Methods in Applied Geophysics：Volume 2，Application，Parts A and B［M］.Tulsa：Society of Exploration Geophysics，1991.

［17］NABIGHIAN M N.Electromagnetic Methods in Applied Geophysics：Volume 1，Theory［M］.Tulsa：Society of Exploration Geophysics.1988.

［18］GUPTASARMA D，SINGH B.New digital linear filter for Hankel J0and J1transforms［J］.GEOPHYSICAL PROSPECTING，1997，45（5）：745－762.

［19］ANDERSON W L.Improved digital filters for evaluating Fourier and Hankel transform integrals［R］.U.S.G.S.REPORT，GD－75－012，PAGES，223，1975.

［20］WOODEN B，AZARI M，SOLIMAN M.Well test analysis benefits from new method of Laplace space inversion［J］.OIL ＆GAS JOURNAL，1992，90（29）：108－110.

［21］羅延鐘，張勝業(yè)，王衛(wèi)平.時間域航空電磁一維正演研究［J］.地球物理學(xué)報，2003，46（5）：719－724.

LUO Y ZH，ZHANG SH Y，WANG W P.A research on one－dimension forward for aerial electromagnetic method in time domain［J］.Chinese Journal of Geophysics，2003，46（5）：719－724.（In Chinese）

Analysis and coefficients selection of frequency－time conversion methods for airborne transient electromagnetic forward modeling

DU Xing－zhong，ZHU Hai－dong
（Guiyang Hydropower investigation design and Research Institute，CHECC，Guiyang 550081，China）

The change rates of vertical magnetic field with time are usually obtained by the frequency－time conversion method in transient electromagnetic modeling.At present，despite many kinds of frequency－time conversion methods have been proposed，the G－S transform with faster computation speed has become an important candidate for airborne transient electromagnetic method.Due to the huge data and the long time of airborne electromagnetic method，how to improve the computation speed of modeling under the certain accuracy becomes a key problem to be solved.Due to the higher frequency－time conversion precision of sine and cosine transform with the larger number of frequency coefficients，this paper may determine the best coefficients of G－S transform by the comparison of the forward modeling results of sine transform，cosine transform and G－S transform with different frequency based on the three－layer and four－layer model.By comparison and analysis of the experimental results，this paper shows that the G－S transform with 8coefficients not only has faster speed，but also a high degree of accuracy.

airborne transient electromagnetic；G－S transform；sine transform；cosine transform

P 631.3

A

10.3969／j.issn.1001－1749.2015.06.01

1001－1749（2015）06－0671－09

2015－08－30改回日期：2015－09－25

貴州省科學(xué)技術(shù)基金（［2013］2297）

杜興忠（1973－），男，博士，高級工程師，主要從事工程物探技術(shù)研究和應(yīng)用，E－mail：xingzhongdu＠163.com。