殷長春, 任秀艷*, 劉云鶴,2, 蔡晶
1 吉林大學地球探測科學與技術學院, 長春 130026 2 中國地質大學(武漢)地球內部多尺度成像湖北省重點實驗室, 武漢 430074
?
航空瞬變電磁法對地下典型目標體的探測能力研究
殷長春1, 任秀艷1*, 劉云鶴1,2, 蔡晶1
1 吉林大學地球探測科學與技術學院, 長春 130026 2 中國地質大學(武漢)地球內部多尺度成像湖北省重點實驗室, 武漢 430074
航空電磁法的探測能力受飛行高度、發射波形、發射磁矩和發射基頻等因素的影響,致使不同分量間的勘探能力存在差異.航空電磁如對所有磁場和磁感應分量、on-和off-time數據進行觀測和解釋,不僅數據量大、耗時長,而且出現大量冗余數據.目前國內針對此問題尚無系統解決方法.本文針對吊艙式直升機航空電磁系統,采用積分方程法求解頻率域響應,經漢克爾變換轉換到時間域,計算了地下三維目標體的B和dB/dt時間域響應.利用異常體響應與背景場響應作比值,并通過設定響應閥值定義最大勘探深度,進而分析不同發射波形、不同分量以及on-和off-time期間的航空電磁系統的探測能力.基于本文分析手段,可根據實際勘探目標,確定一套探測能力較強的航空電磁最佳參數組合,為野外測量和數據處理提供技術指導,高效完成勘探任務.
航空電磁; 時間域; 探測能力; 系統設計; 磁場和磁感應; on-time和off-time
航空電磁法(Airborne Electromagnetic簡稱AEM) 是一種以飛機為運載工具進行地下目標體探測的方法.它通過測量大地的二次磁場來研究地層沿水平和垂直方向的電性差異,描述地電斷面的電性特征,進而了解地質構造情況(黃皓平,1991).AEM具有速度快、成本低、通行性好、可用于海域和大面積覆蓋區域勘測等優勢.時間域航空電磁系統用于地質填圖和礦產普查已超過50年,其中大部分系統主要使用的發射波形為半正弦波,發射基頻介于25 Hz和300 Hz之間,隨著地質條件不同而改變.
AEM系統誕生于1948年,在20世紀70年代得到迅速發展.90年代起,時間域航電系統開始向大磁矩發射、多分量測量方向發展,并通過選擇合適的脈沖頻率及寬度,以期在導電地層中獲得大的穿透深度.黃皓平和王維中(1990)對時間域航空電磁數據進行了反演研究,Smith和Annan (1998)分析了時間域航電系統磁場數據相對于磁感應數據的優點,Liu(1998)研究了脈沖、方波、半正弦波等幾種不同發射波形對航空瞬變電磁響應的影響, 而Balch等(2003)等對三角波激勵源下AeroTEM系統的on-time數據進行了反演解釋.羅延鐘等(2003)導出了層狀大地條件下時間域航空電磁法的正演計算公式,Yin 等(2008)對半正弦和梯形波激勵源下均勻半空間模型的on-和off-time電磁響應特征進行模擬和分析.陳曙東等(2012)利用自由空間回線作為目標體,推導出方波、梯形波、半正弦波和三角波的瞬變電磁響應.時間域航空電磁系統由于發射線圈和接收線圈感應的影響,理想的階躍波形不可實現,因此,數據的處理和成像主要在off-time進行,這常常導致地面淺部信息丟失(許洋鋮等,2012).半正弦波和梯形波在一定程度上可以彌補不足,對其進行全波形的分析和研究具有重要意義.裴易峰等(2014)介紹了由dB/dt數值計算磁場B的積分算法,并分析了磁場B和磁感應dB/dt對地下良導體的探測能力.
本文利用積分方程進行了三維頻率域航空電磁響應的正演計算,通過漢克爾變換及高斯積分得到航空電磁三維時間域響應.然后,探討了異常體電阻率對時間域響應的影響特征.在分析探測能力過程中,我們首先將異常響應最大值的位置設為記錄點位置,并通過不斷地增大異常體的埋深,使其對應的異常響應不斷減小,當異常響應值到達背景(均勻半空間)響應值的10%時,定義此時的異常體埋深為最大探測深度.進而,本文針對不同發射波形、磁場和磁感應分量、on-和off-time期間的航空電磁勘探深度進行系統研究,并對不同情況下的航空電磁探測能力進行了討論和分析.鑒于航空電磁系統對于高阻體探測能力有限的情況,本文對航空電磁系統用于高阻目標體探測能力分析不作詳細討論.本文的研究思路有益于針對不同勘探目標體設計最佳航空電磁系統參數組合,從而使得設計的儀器系統以最優的觀測方式實現地下目標體探測.另外,由于我國航空電磁勘查技術目前相對落后,本文航空電磁系統探測能力研究有益于我國航空電磁技術標準和規范的制定.
本文采用的頻率域三維正演算法是基于Raiche(2001)建立的.其基本思路是利用并矢格林函數理論建立二次感應場和一次場的關系,求解異常體內感應電流分布,并通過對異常體內感應電流的體積分計算航空電磁系統的頻率域響應.
由麥克斯韋旋度方程出發,可得到頻率域三維電磁散射問題的矢量場問題:
(1)
(2)
式中,波數k2=iωμ(σ-iωε),其中μ為磁導率,ε為介電常數,Jf代表源電流密度.
引入并矢格林函數求解磁場積分表達式,其滿足的方程為
(3)
結合(1)和(3)式,得到計算任意一點電場表達式為
(4)
式中,EP(r)為入射電場,σ*和ε*分別代表目標體與圍巖的電導率和介電常數之差.利用迭代方法求解(4)式可得目標體內電場,進而利用歐姆定律計算目標體內電流密度.
根據Raiche(1974)的思想,利用電磁場和并矢格林函數滿足的方程,可得到磁場積分:
(5)
利用傅里葉變換將頻率域航空電磁響應轉換為時間域響應,即
ω.
(6)
通過簡單的變量代換,我們可得到與(6)式相對應的階躍波響應:

(7)
應用褶積理論和脈沖響應與階躍響應間的導數關系,我們可得到任意發射波形時間域響應(殷長春等,2013):
(8)
(9)
式中I(t)代表發射電流,*號代表褶積.(8)和(9)式的褶積可用高斯積分進行計算(Yin 等,2008;殷長春等,2013).
假設均勻半空間中存在一個100 m×100 m×400 m的地質體,圍巖電阻率為100 Ωm,網格剖分為10 m×10 m×40 m.發射和接收線圈中心的高度分別為30 m和50 m,發射與接收線圈水平距離為10 m.異常體中心位置在地面上的投影作為坐標原點,取Z軸向下為正,如圖1所示.發射波形考慮兩種情況:第一種為半正弦波,基頻30 Hz,脈沖寬度為4 ms;第二種為梯形波,基頻30 Hz,上升沿和下降沿時間均為0.2 ms,穩定電流時間為3.6 ms.發射偶極距為615000 Am2.
首先對基頻為30 Hz的半正弦波和梯形波的垂直磁場及磁感應的時間域響應進行了計算,得到其隨地質體埋藏深度變化的曲線,如圖2所示.
從圖2可以看出,磁場B和磁感應dB/dt保持良好的微分/積分關系,且無論半正弦波或梯形波,B和dB/dt在on-和off-time都存在響應極大值.地下存在異常體時的時間域響應與均勻半空間的響應

圖1 模型示意圖Fig.1 Sketch map of AEM system and the 3D earth model

圖2 B和dB/dt時間域響應隨深度變化(a,c)分別為發射半正弦波時的B和dB/dt響應曲線;(b,d)分別為發射梯形波時的B和dB/dt響應曲線.Fig.2 Time-domain AEM responses B and dB/dt versus depth(a) and (c) are half-sine B and dB/dt responses, respectively; (b) and (d) are trapezoid B and dB/dt responses, respectively.

圖3 響應比隨異常體電阻率變化曲線(a)、(c)和(e)分別代表發射半正弦波時Bx、Bz和dBz/dt與半空間響應的比值隨電阻率變化曲線;(b)、(d)和(f)分別代表發射梯形波時Bx、Bz和dBz/dt與半空間響應的比值隨電阻率變化曲線.Fig.3 Time-domain EM response ratios versus resistivity(a), (c) and (e) are ratios of Bx, Bz and dBz/dt for a half-sine wave to half-space responses versus resistivity; (b), (d) and (f) are ratios of Bx, Bz and dBz/dt for a trapezoid wave to half-space responses versus resistivity

圖5 半正弦波B和dB/dt的勘探深度對比(a) on- 和off-time期間的Bx和Bz響應比隨深度變化;(b) on- 和off-time期間的dBx/dt和dBz/dt響應比隨深度變化;(c) on-和off-time期間的Bx和dBx/dt響應比隨深度變化;(d) on-和off-time期間的Bz和dBz/dt響應比隨深度變化.Fig.5 Exploration depth of B and dB/dt for a half-sine wave(a) on- and off-time Bx and Bz response ratios varying with depth; (b) on- and off-time dBx/dt and dBz/dt response ratios varying with depth; (c) on- and off-time Bx and dBx/dt response ratios varying with depth; (d) on- and off-time Bz and dBz/dt response ratio varying with depth.

圖6 梯形波B和dB/dt的勘探深度對比(a) on- 和off-time期間的Bx和Bz響應比隨深度變化;(b) on- 和off-time期間的dBx/dt和dBz/dt響應比隨深度變化;(c) on- 和off-time期間的Bx和dBx/dt響應比隨深度變化;(d) on- 和off-time期間的Bz和dBz/dt響應比隨深度變化.Fig.6 Exploration depth of B and dB/dt for a trapezoid wave(a) on- and off-time Bx and Bz response ratios varying with depth; (b) on- and off-time dBx/dt and dBz/dt response ratios varying with depth; (c) on- and off-time Bx and dBx/dt response ratios varying with depth; (d) on- and off-time Bz and dBz/dt response ratio varying with depth.

圖7 不同波形B和dB/dt的勘探深度對比(a)—(d)發射波形為半正弦波; (e)—(h)發射波形為梯形波. (a)和(e)為Bz響應的on- 和off-time勘探深度對比圖;(b)和(f)為dBz/dt響應的on-和off-time勘探深度對比圖;(c)和(g)為on-time的Bx和dBx/dt勘探深度對比圖;(d)和(h)為on-time的Bz和dBz/dt勘探深度對比圖.Fig.7 Exploration depth of B and dB/dt for different transmitting waveforms(a)—(d) Half-sine transmitting waves; (e)—(h) Trapezoid transmitting waves. (a) and (e) are on- and off-time Bz exploration depth; (b) and (f) are on- and off-time dBz/dt exploration depth; (c) and (g) are on-time Bx and dBx/dt exploration depth; (d) and (h) are on-time Bz and dBz/dt exploration depth.
曲線趨勢十分一致,并隨著深度的不斷增加,響應值不斷接近半空間響應.由于探測分辨率限制,航空電磁系統無法識別所有高于半空間背景響應的異常.本文以異常響應值到達均勻半空間背景響應的10%的信號作為最小可識別信號,其對應的探測深度定義為最大探測深度.
地質體的電阻率變化會引起異常響應變化,使探測能力發生變化.當異常體電阻率從0.1 Ωm變化到10000 Ωm時,總響應與均勻半空間背景響應之比的變化規律如圖3所示.

4.1 記錄點位置
非對稱航空電磁系統的記錄點由模型試驗確定.取異常體電阻率為1 Ωm,異常體埋深為100 m,發射基頻為30 Hz.圖4為沿飛行剖面的各磁場和磁感應分量的時間域響應.圖中橫坐標為發射線圈中心對應地面測線上的位置,縱坐標為存在目標體的時間域響應與均勻半空間響應的比值.
根據剖面結果,以時間域二次場響應最大值處作為記錄點位置,可以確定各場分量Bx、Bz及磁感應分量dBx/dt和dBz/dt在探測異常體最大埋深時的發射線圈的位置.經計算,對于半正弦波,可得Bx、Bz、dBx/dt、dBz/dt記錄點位置對應x、y坐標分別為(40.94 m,0),(-4.45 m,0),(33.49 m,0),(-4.45 m,0);對于梯形波,其Bz和dBz/dt響應值最大時記錄點位置對應坐標均為(-4.45 m,0).
4.2 探測能力分析
在確定異常記錄點位置的基礎上,本文對不同發射波形、不同時間道情況下的Bx、Bz和dBx/dt、dBz/dt進行了勘探深度的比較和分析.圖5為發射半正弦波時,磁場分量和磁感應各分量響應比隨異常體埋深的變化;圖6為發射梯形波時,響應比隨異常體埋深的變化.其中橫軸表示深度,縱軸為航空電磁響應與均勻半空間響應之比.
由圖5可見,隨著地質體埋藏深度的不斷增加,on-和off-time的最大響應不斷減小.圖5a顯示,同線磁場分量Bx在on-和off-time期間的探測深度均大于垂向磁場分量Bz;同時,off-time的磁場分量的探測深度大于on-time;圖5b顯示off-time的垂向磁感應分量,從異常體埋深較淺至最大深度變化過程中,其二次場響應都相對較小;圖5(c,d)顯示,無論on-time或off-time,磁感應同線分量dBx/dt的探測深度大于磁場Bx的探測深度;垂向分量dBz/dt的探測深度大于Bz的探測深度.
如圖6所示,當發射波形為梯形波時,對于同一磁場或磁感應分量,其on-time和off-time期間的探測深度十分接近,但不同分量間存在一定差別.通過對比圖6(a—d)發現,Bx的探測深度大于Bz,且垂向磁感應分量dBz/dt探測能力弱于dBx/dt;同時,dBx/dt的探測能力強于Bx,dBz/dt的探測能力強于Bz,這與發射半正弦波時結論相似.
通過計算,得到半正弦波和梯形波各磁場分量和磁感應分量的on-和off-time所對應的勘探深度情況如表1所示.
表1數據顯示,對于柱狀良導體,在on-time和off-time過程中,不同發射波形、不同磁場和磁感應分量,對應不同的勘探深度.發射半正弦波時的探測能力分析結論與發射梯形波時結論相似.航空測量可根據實際地質條件選擇合適的參量組合進行最優觀測和數據解釋.

表1 磁場B和磁感應dB/dt的探測深度Table 1 Exploration depth for magnetic field B and magnetic induction dB/dt
取異常體電阻率為10 Ωm,其他參數選擇與前述柱狀良導體一致.利用類似于低阻區勘探能力分析方法,可以獲得該地質條件下不同磁場及磁感應分量間的探測能力.
如圖7所示,當發射半正弦波時,垂向磁場分量Bz及磁感應分量dBz/dt在on-time的勘探深度大于off-time的勘探深度,且磁感應分量的探測深度在on-time還是off-time都大于對應的磁場分量;當發射波形為梯形波時,垂向磁場和磁感應分量各自的探測能力在on-和off-time期間相當,且on-和off-time時的同線磁感應分量的探測深度比磁場分量的探測深度稍大.經過計算,不同波形時的各分量的探測深度對比見表2.

表2 磁場B和磁感應dB/dt的探測深度(模型電阻率10 Ωm)Table 2 Exploration depth for B and dB/dt for a model of 10 Ωm
表2顯示的是相同地質體條件下(中阻區),不同發射波形時,相應磁場和磁感應分量的探測深度值.可以看出,大部分同一分量的on-和off-time的相對探測能力結論與低阻區的結論相反,其他各交差分量間的探測能力分析結論與低阻區一致.
當地質體形狀、尺寸和電阻率發生變化時,航空電磁響應及探測能力隨之變化.這里僅介紹垂直薄板和水平良導板的探測情況,所使用參數與低阻區探測能力分析時采用的模型一致.模型尺寸與勘探深度對比如表3所示.

表3 磁場B和磁感應dB/dt的探測深度Table 3 Exploration depth for B and dB/dt
由表3可知,對于直立良導薄板,無論發射半正弦波還是梯形波,垂向磁場分量探測能力明顯弱于同線分量的探測能力,這與電磁基本理論相吻合.在off-time期間,磁感應分量的探測深度均大于其對應的磁場分量的探測深度.對于水平良導板狀體,發射梯形波時,on- 和off-time的探測深度依然接近.實際航空電磁觀測中,需要針對特定地質情況進行模型設計和正演計算分析,以獲得特定探測目標的航空電磁參量組合,指導野外探測和數據解釋.
本文在時間域航空電磁法勘探原理和算法的基礎上,以均勻半空間中三維良導體為例,進行了半正弦波和梯形波激勵下的各磁場及磁感應分量的on-和off-time時間域響應的計算.通過對比分析,得出對深部異常體勘探能力的結論如下:
1) 對于非對稱航空電磁系統而言,記錄點的位置可以通過數值試驗來確定;
2) 對于柱狀良導體,無論on-還是off-time,同線磁場分量Bx的探測深度均大于垂向磁場分量Bz,磁感應分量的探測深度均大于其對應的磁場分量的探測深度;
3) 發射波形不同,航空電磁對地質體的探測能力不同,梯形波的同一個場分量在on-time和off-time期間的探測深度差別較小;
4) 地質體的形狀和尺寸對磁場和磁感應分量在on-和off-time的相對探測能力產生影響.實際航空電磁的觀測設計中,我們應根據勘探目標的特征,選擇具有最佳偶合的觀測系統,分析其對目標體的探測能力,針對設定的模型進行計算和分析,以獲得最佳參數.
根據以上結論及本文研究思想,可以針對具體的勘探目標,設計一套最佳航空電磁觀測系統參數組合,以保證在有效完成勘探目標探測的前提下,減少航空觀測的工作量、提高工作效率.
致謝 賁放和黃威博士在文章準備過程中提供的幫助及審稿人對本文提出的修改意見一并衷心感謝.
Balch S J, Boyko W P, Paterson N R. 2003. The AeroTEM airborne electromagnetic system.TheleadingEdge, 22(6): 562-566.
Chen S D, Lin J, Zhang S. 2012. Effect of transmitter current waveform on TEM response.ChineseJ.Geophys. (in Chinese), 55(2): 709-716, doi: 10.6038/ j.issn.0001-5733. 2012.02.035.Huang H P, Wang W Z. 1990. Inversion of time-domain airborne electromagnetic data.ChineseJ.Geophys. (in Chinese), 33(1): 88-97.Huang H P. 1991. The theoretical problems of the application of airborne electromagnetic method in engineering geophysical exploration (in Chinese). ∥ The Academic Essays of the 7th Chinese Geophysical Society.
Liu G. 1998. Effect of transmitter current waveform on airborne TEM response.ExplorationGeophysics, 29(2): 35-41.
Luo Y Z, Zhang S Y, Wang W P. 2003. A research on one-dimension forward for aerial electromagnetic method in time domain.ChineseJ.Geophys. (in Chinese), 46(5): 719-724, doi: 10.3321/j.issn:0001-5733.2003.05.021.
Pei Y F, Yin C C, Liu Y H, et al. 2014. Calculation and application of B-field for time-domain airborne EM.ProgressinGeophysics(in Chinese), 29(5): 2191-2196, doi: 10.6038/pg20140530.
Raiche A P. 1974. An integral equation approach to Three-Dimensional modelling.GeophysicalJournaloftheRoyalAstronomicalSociety, 36(2): 363-376.Raiche A. 2001. 3D EM modeling using integral equation algorithm. AMIRA Project Report P223E.
Smith R, Annan P. 1998. The use of B-field measurements in an airborne time-domain system: Part Ⅰ: Benefits of B-field versus dB/dt data.ExplorationGeophysics, 29(2): 24-29.
Xu Y C, Lin J, Li S Y, et al. 2012. Calculation of full-waveform airborne electromagnetic response with three-dimension finite-difference solution in time-domain.ChineseJ.Geophys. (in Chinese), 55(6): 2105-2114, doi: 10.6038/j.issn.00015733.2012.06.032.Yin C, Huang W, Ben F. 2013. The full-time electromagnetic modeling for time-domain airborne electromagnetic systems.ChineseJ.Geophys. (in Chinese), 56(9): 3153-3162, doi: 10.6038/cjg20130928.
Yin C, Smith R S, Hodges G, et al. 2008. Modeling results of on-and off-time B and dB/dt for time-domain airborne EM systems. ∥ 70th Annual EAGE Conference and Exhibition, Extended Abstract, Rome, 1-4.
附中文參考文獻
陳曙東, 林君, 張爽. 2012. 發射電流波形對瞬變電磁響應的影響. 地球物理學報, 55(2): 709-716, doi: 10.6038/j.issn.0001-5733.2012.02.035.
黃皓平, 王維中. 1990. 時間域航空電磁數據的反演. 地球物理學報, 33(1): 88-97.
黃皓平. 1991. 航空電磁法在工程物探中應用的理論問題. ∥ 中國地球物理學會第七屆學術年會論文集.
羅延鐘, 張勝業, 王衛平. 2003. 時間域航空電磁法一維正演研究. 地球物理學報, 46(5): 719-724, doi: 10.3321/j.issn:0001-5733.2003.05.021. 裴易峰, 殷長春, 劉云鶴等. 2014. 時間域航空電磁磁場計算與應用. 地球物理學進展, 29(5): 2191-2196, doi: 10.6038/pg20140530. 許洋鋮, 林君, 李肅義等. 2012. 全波形時間域航空電磁響應三維有限差分數值計算. 地球物理學報, 55(6): 2105-2114, doi: 10.6038/j.issn.00015733.2012.06.032.
殷長春, 黃威, 賁放. 2013. 時間域航空電磁系統瞬變全時響應正演模擬. 地球物理學報, 56(9): 3153-3162, doi: 10.6038/cjg20130928.
(本文編輯 汪海英)
Exploration capability of airborne TEM systems for typical targets in the subsurface
YIN Chang-Chun1, REN Xiu-Yan1*, LIU Yun-He1,2, CAI Jing1
1CollegeofGeo-explorationScienceandTechnology,JilinUniversity,Changchun130026,China2HubeiSubsurfaceMulti-scaleImagingLab(SMIL),ChinaUniversityofGeosciences(Wuhan),Wuhan430074,China
For airborne EM (AEM) systems, flight altitude, transmitting waveforms, transmitting dipole moment and base frequency can strongly affect the depth of exploration. In an AEM survey, if all the magnetic field and magnetic induction components are measured, the EM dataset will be huge,resulting in costly data processing. In this paper, we investigate the exploration capability of an AEM system to different targets in the subsurface. We try to optimize parameter combinations of the airborne system (e.g. the waveform, field components, on- and off-time etc.) to balance the survey cost and resolution of targets.For the towed-bird AEM system, we use an integral equation algorithm to calculate the frequency-domain EM field responses and convert them into the time domain via a Hankel′s transform. The on- and off-time magnetic fieldsBand magnetic induction dB/dtare calculated for 3D targets of a plate or a prism embedded in a conductive earth. We propose a response ratio to define the largest exploration depth, based on which we calculate the exploration depth to different underground targets for different transmitting waves, field components and on- and off-time signals.We study and analyze the influence of transmitting waveforms, field components, the on- and off-time signals on the target exploration capability. We find that the off-timeBx, dBx/dt, and dBz/dtfor a trapezoid transmitting wave is the best parameter combination for resolving an underground conductive prism target. For a vertical thin plate, however, the off-timeBx, dBx/dtand on-time dBz/dtfor a half-sine transmitting wave perform better. The combination of off-timeBx, dBx/dtand on-timeBz, dBz/dtfor a trapezoid waveform can make deeper exploration for a horizontal thick plate.From the research in this paper, we draw the conclusion that depending on the exploration target, there exists an optimal parameter combination of the AEM system that can achieve the maximum exploration capability. This will guide AEM survey design and data processing for an effective and efficient exploration. This research further lays the foundation for establishment of AEM standards and regulations.
Airborne EM; Time-domain; Exploration capability; System design; Magnetic field and magnetic induction; On- and off-time
殷長春,任秀艷, 劉云鶴.2015.航空瞬變電磁法對地下典型目標體的探測能力研究.地球物理學報,58(9):3370-3379,
10.6038/cjg20150929.
Yin C C, Ren X Y, Liu Y H,et al. 2015. Exploration capability of airborne TEM systems for typical targets in the subsurface.ChineseJ.Geophys. (in Chinese),58(9):3370-3379,doi:10.6038/cjg20150929.
10.6038/cjg20150929
P631
2015-02-05,2015-06-27收修定稿
國家自然科學基金項目(41274121),國家重大科研裝備研究項目(ZDYZ2012-1-03和20130523MTEM05),中國地質大學(武漢)地球內部多尺度成像湖北省重點實驗室開放基金項目(SMIL-2014-03)聯合資助.
殷長春,男,1965年生,教授,主要從事電磁勘探理論,特別是航空和海洋電磁方面的研究.E-mail:yinchangchun@jlu.edu.cn
*通訊作者 任秀艷,女,1989年生,碩士,主要從事航空電磁探測能力、時間域有限差分方法研究.E-mail:jdrxy@hotmail.com