曾昭發(fā)，霍祉君，李文奔，李 靜，趙雪宇，何榮欽

1.吉林大學地球探測科學與技術(shù)學院，長春 130026 2.國土資源部應用地球物理重點實驗室，長春 130026 3.河北地質(zhì)大學信息工程學院，石家莊 050031

0 引言

航空電磁法(airborne electromagnetic method)因其不受地形影響、測量面積大等優(yōu)點而廣泛應用于復雜地區(qū)的礦產(chǎn)勘查、地質(zhì)填圖、地下水勘查及環(huán)境檢測等領(lǐng)域。1984年Stanmac和McPhar公司在加拿大成功實現(xiàn)了固定翼飛機AEM系統(tǒng)的首次飛行,標志著AEM的誕生，之后航空電磁法測量平臺和系統(tǒng)配置不斷創(chuàng)新和發(fā)展[1]。我國于20世紀60年代初開始研制長導線式航空電磁系統(tǒng)，之后又陸續(xù)研制了雙頻和三頻硬架式航電儀[2]。隨著航空電磁技術(shù)的發(fā)展，航空電磁正演計算方法一直是理論研究的重要內(nèi)容。Newman等[3]應用交錯有限差分(FD)分別進行了二維和三維帶地形的三維航空電磁響應。朱凱光等[4]開展了一維層狀模型航空電磁法正演模擬，同時研究了航空電磁響應的影響因素和一定條件下有效探測深度。高亮等[5]運用Impulse系統(tǒng)(頻率域吊艙式直升機航空電磁測量系統(tǒng))對頻率域航空電磁系統(tǒng)進行了磁性條件下頻率域航空電磁正演研究。李小康[6]基于MPI(message passing interface，用于并行計算的通信協(xié)議)運用有限元法并行計算了2.5維頻率域航空響應。曲昕馨[7]基于有限差分法實現(xiàn)了三維頻率域航空電磁法正演模擬，同時分析了線圈姿態(tài)對結(jié)果的影響并給出了校正方案。王衛(wèi)平等[8]在實現(xiàn)二維和三維航空電磁正演模擬的基礎(chǔ)上研究了地形影響,并使用地形校正法對典型地形地電模型進行校正。殷長春等[9]運用非結(jié)構(gòu)化網(wǎng)格的有限元法進行了2.5維航空電磁正演模擬，分析總結(jié)了典型地形對航空響應的影響及其特征。黃威等[10]利用三維矢量有限元正演模擬研究了覆蓋層和垂直接觸帶等典型地電構(gòu)造對航空電磁響應的影響特征。Li等[11]運用等參有限元法進行了三維源/二維地電模型(2.5維)的頻率域航空電磁法正反演算法研究。

由于航空電磁方法探測環(huán)境復雜，地下介質(zhì)的各向異性特征越來越受到重視。Yin 等[12-13]開發(fā)了用于計算層狀大地任意各向異性中直流電場和磁場的算法;處理了任意各向異性介質(zhì)中正演問題。Avdeev等[14]研究了各向異性對航空電磁測量的影響，但其研究結(jié)果只適用于三維簡單各向異性地質(zhì)模型。Yin等[15-16]研究各向異性介質(zhì)對于一維介質(zhì)模型的影響；考慮了各向異性對機載電磁系統(tǒng)響應的影響,并分析了實測航空電磁數(shù)據(jù)的各向異性主軸方向。Liu等[17]運用有限差分法進行了任意三維各向異性模型的航空電磁響應計算，取得了較好的應用效果。Yin等[18]進行了起伏地形任意各向異性模型下的三維時域航空電磁響應正演計算。Huang等[19]采用譜元法對航空電磁各向異性響應進行三維正演模擬，得出了不同各向異性異常體的電磁響應的收斂條件及其識別方式。上述研究中，各向異性介質(zhì)電磁模擬研究主要以有限差分為主，對于地下復雜結(jié)構(gòu)模擬難以達到滿意的模擬精度。本文采用矢量有限元法開展三維任意各向異性介質(zhì)頻率域航空電磁響應模擬算法，利用矢量有限元法求解散射電場耦合方程，通過將總場分解為一次場(背景場)和二次場(散射場)來進行各向異性介質(zhì)中的頻率域三維航空電磁正演模擬。針對三維模型計算量大的問題，本文采用共享內(nèi)存直接求解器PARDISO對大規(guī)模稀疏矩陣并行計算，大大提高了計算效率。

1 正演理論

E=Es+Ep；

(1)

H=Hs+HP。

(2)

式中：E為電場強度；H為磁場強度；下標p和s分別代表一次場和二次場。為了在保證精度的同時能提高計算效率，本文以空氣為介質(zhì)的均勻全空間為背景對一次場采用解析式計算，采用的源為垂直磁偶極源(對應水平共面裝置，簡稱HCP)，引入謝坤諾夫矢量勢F，通過解偏微分方程可得矢量勢F表達式：

(3)

由此得到一次場電場和磁場的表達式：

Ep=-

(4)

(-k2r2+3ikr+3)+(k2r2-ikr-1)uz]。

(5)

對于二次場的計算，下式為正演模擬中二次場形式下的麥克斯韋方程組：

×Es=-iμmHs-iω(μ-μ0)Hp，

(6)

(7)

(6)式和(7)式中：Js=[(σ-σp)+iω(ε-εp)]Ep,為等效場源項，σp為背景電導率。在各向異性介質(zhì)中，σ及σp用一個3×3的張量矩陣表示[20]。為實現(xiàn)任意各向異性電導率介質(zhì)模型，利用歐拉旋轉(zhuǎn)從主軸坐標系轉(zhuǎn)換到模擬坐標系，主軸電導率張量用如下形式表示：

(8)

其下標x、y和z為笛卡爾坐標系的方向。歐拉旋轉(zhuǎn)計算過程如下：

σ=DTσ0D，D=DxDyDz。

(9)

式中：D為主旋轉(zhuǎn)矩陣；Dx、Dy、Dz分別是繞著x軸、y軸和z軸旋的對應矩陣(圖1)，φ、ψ、χ分別為繞著x軸、y軸、z軸旋轉(zhuǎn)的角度，即

(10)

最后得到如下形式的任意各向異性介質(zhì)中的電導率張量，該張量具有對稱和正定性質(zhì)[21]：

(11)

假定磁導率為真空中的磁導率μ0，不考慮存在位移電流的情況，對式(6)兩邊同時取旋度，將(7)式

圖1 坐標旋轉(zhuǎn)示意圖 Fig.1 Coordinate rotation diagram

代入得到二次電場的雙旋度方程：

××Es+iωμ0σEs=

-iωμ0ΔσEp，Δσ=σ-σp。

(12)

采用矢量有限元法，用六面網(wǎng)格將模型剖分，并將模擬區(qū)域離散成Ne個單元，對式(12)采用伽遼金加權(quán)余量法可得：

iωμ0σEs+iωμ0ΔσEp]dΩe。

(13)

(14)

式中，ns為面外法線單位向量。

使用狄利克雷邊界條件，n為邊界上的法向量，即假設(shè)在遠離異常體的邊界處二次異常場衰減為0，棱邊元賦值為0，該條件使得電磁場呈指數(shù)衰減，從而邊界處不產(chǎn)生反射并遠離異常區(qū)[22]：

n×E|?Ω=0。

(15)

根據(jù)物性變化界面處電磁場切向連續(xù)性及所選用邊界條件，式(14)的第二項可以忽略，改寫成：

(16)

根據(jù)矢量有限元法，即將自由度賦于棱邊而不是單元節(jié)點，可得到任意單元內(nèi)場的展開式：

(17)

(18)

(19)

(20)

圖2 xyz坐標系與ξηζ坐標系轉(zhuǎn)換[23]Fig.2 Conversion of xyz coordinate system and ξηζ coordinate system

Es=[Ee]，將式(17)代入式(16)得：

(21)

Ke和Me即為單元剛度矩陣，具體如下：

(22)

(23)

式中:k=1，…，12;l=1，…，12。之后，引入雅克行列式將函數(shù)轉(zhuǎn)換為以ξ、η和ζ表示的函數(shù)，則式(22)和(23)轉(zhuǎn)化為

(24)

(25)

將單元剛度矩陣Me和Ke按式(22)進行全局剛度矩陣合成，可以得到大型線性方程組：

Ax=b

(26)

式中:A為稀疏復對稱系數(shù)矩陣；x為帶求解的場值；b為源端項。引入PARDISO求解器，運用其內(nèi)嵌套分割并行算法對方程共享內(nèi)存并行直接求解。求解出電場值后，磁場值可根據(jù)法拉第定律求得：

(27)

2 正演模擬結(jié)果

2.1 精度驗證

本文借鑒Liu等[17]的三維各向異性模型進行精度驗證。三層介質(zhì)模型如圖3所示，其中各向異性低阻層厚度為100 m，取主軸電導率元素為σxx=1 S/m、σyy=10 S/m、σzz=1 S/m，并令主電導率張量繞y軸旋轉(zhuǎn)45°，上方為空氣層，下方為σ=0.01 S/m的各向同性介質(zhì)高阻均勻半空間。采用水平磁偶極子源，接收線圈和發(fā)射線圈均離地表20 m，發(fā)射頻率為900 Hz。三維有限元計算結(jié)果與一維擬解析解結(jié)果對比如圖4所示。由圖4可見，一維解析解和本文三維矢量有限元模擬結(jié)果能很好地吻合，說明了本文算法的可靠性。

圖3 精度驗證模型Fig.3 Precision verification model

圖4 驗證結(jié)果Fig.4 Validation results

2.2 各向異性介質(zhì)對航空電磁響應影響特征

為了進一步分析不同電阻率各向異性特征對航空電磁響應的影響，建立如圖5所示的三維地電模型。該模型為一個高導異常體埋藏在均勻半空間中，模型大小為1 540 m×1 540 m×1 540 m，異常體大小為50 m×50 m×50 m，異常體中心位于坐標系原點，埋深為30 m，三維模型剖分為42×42×46個單元，包含擴邊區(qū)域和目標區(qū)域。其中，擴邊區(qū)單元為呈2倍擴展的漸變粗網(wǎng)格，目標區(qū)域為均勻細網(wǎng)格。同時，針對空中場的(空間)衰減較地中有耗衰減得慢的問題，對于空氣層的擴邊多，地下區(qū)域擴邊相對少，盡可能減少截斷邊界對模擬區(qū)域的影響，從而在保證了計算精度的同時提高計算速度。長源發(fā)射頻率為900 Hz，記錄磁場Hz分量，發(fā)射線圈和接收線圈距離為10 m，離地面20 m，點距為10 m，每條測線共計21個測點。

圖5 三維地電模型與測線分布Fig.5 Three dimensional model and survey line distribution

首先，當各向異性異常體的主電導率張量繞不同主軸旋轉(zhuǎn)后，得到傾斜各向異性電導率張量，設(shè)各向異性異常體主電導率參量為

(28)

然后分別討論以下不同各向異性情況對航空電磁響應影響規(guī)律：1)將(28)式的主電導率參量分別繞著x、z軸旋轉(zhuǎn)0°、45°和90°，分別研究磁場垂直分量的實部和虛部電磁響應特征；2)各向異性圍巖電阻率變化響應特征。

2.2.1 繞z軸旋轉(zhuǎn)各向異性磁場響應特征

當異常體主電導率繞z主軸旋轉(zhuǎn)0°、45°、90°時，相當于目標體從各向異性介質(zhì)變?yōu)閮A斜各向異性介質(zhì)的過程。圖6為上述電導率張量繞z軸旋轉(zhuǎn)不同角度后得到的磁場Hz分量在xy面的響應特征。由于垂直磁偶極子所產(chǎn)生的電流在地層內(nèi)主要沿xy面流動，因此磁場響應Hz分量沿著電導率的旋轉(zhuǎn)角度旋轉(zhuǎn)，場的分布形態(tài)和主軸的方位對應良好，可用場值的形態(tài)判斷出主軸方位。根據(jù)磁場特征，可以確定異常體的各向異性主方向以及主電導率張量旋轉(zhuǎn)角度。

圖7所示為繞z軸旋轉(zhuǎn)HCP裝置Hz分量觀測記錄結(jié)果。對比圖中實虛分量各向同性(黑實線)與不同旋轉(zhuǎn)角度各向異性(黑虛線)在信號響應振幅上有較大的差異，可見各向異性參數(shù)對航空電磁觀測結(jié)果影響很大。由各向同性變?yōu)楦飨虍愋詶l件下，實分量響應信號由單峰值信號逐漸變?yōu)殡p峰值信號，并且信號振幅有較大幅度的衰減(圖7a)。Hz虛分量相對影響較小，信號相位沒有明顯變化(圖7b)。

2.2.2 繞x軸旋轉(zhuǎn)各向異性磁場響應特征

當各向異性電性參數(shù)沿著x軸分別旋轉(zhuǎn) 0°、45°和 90°。HCP 裝置磁場Hz分量響應如圖 8所示。隨著旋轉(zhuǎn)角度的變化， HCP 裝置Hz分量的實部隨著旋轉(zhuǎn)角度達到90°時，磁場分量響應從雙峰異常變成了關(guān)于x軸對稱的單峰異常，而虛部變化異常較小。同理，當沿y軸旋轉(zhuǎn)時，HCP裝置Hz分量的實部隨著旋轉(zhuǎn)角度的變化從單峰異常變成關(guān)于y軸對稱的單峰異常。

圖9所示為繞x軸旋轉(zhuǎn)HCP裝置Hz分量觀測記錄結(jié)果。與上述繞z軸旋轉(zhuǎn)結(jié)果類似，實虛分量各向同性(黑實線)與不同旋轉(zhuǎn)角度各向異性(黑虛線)在信號響應振幅上有較大的差異，實分量在不同旋轉(zhuǎn)角度情況下，響應信號振幅影響更大，在目標位置響應信號振幅隨著旋轉(zhuǎn)角度的變化，振幅峰值由-400×10-9A/m(旋轉(zhuǎn)90°)轉(zhuǎn)變?yōu)?100×10-9A/m(旋轉(zhuǎn)0°)。相比而言，虛分量幅值變化率相對較小。

a. 各向異性異常體繞z軸旋轉(zhuǎn)0°Hz分量實部；b. 各向異性異常體繞z軸旋轉(zhuǎn)45°Hz分量實部；c. 各向異性異常體繞z軸旋轉(zhuǎn)90°Hz分量實部；d. 各向異性異常體繞z軸旋轉(zhuǎn)0°Hz分量虛部；e. 各向異性異常體繞z軸旋轉(zhuǎn)45°Hz分量虛部；f. 各向異性異常體繞z軸旋轉(zhuǎn)90°Hz分量虛部。圖6 各向異性異常體繞z軸旋轉(zhuǎn)HCP裝置磁場Hz分量響應Fig.6 Hz component of the magnetic response of the HCP device when the anisotropic body rotates around the z axis

a.各向同性異常體及各向異性異常體繞z軸旋轉(zhuǎn)不同角度Hz分量實部對比；b. 各向同性異常體及各向異性異常體繞z軸旋轉(zhuǎn)不同角度Hz分量虛部對比。圖7 各向異性異常體繞z軸旋轉(zhuǎn)HCP裝置觀測記錄Hz分量 Fig.7 Hz component that is recorded by an observation of the HCP device when an anisotropic anomalous body rotates around the z axis

a. 各向異性異常體繞x軸旋轉(zhuǎn)0°Hz分量實部；b. 各向異性異常體繞x軸旋轉(zhuǎn)45°Hz分量實部；c. 各向異性異常體繞x軸旋轉(zhuǎn)90°Hz分量實部；d. 各向異性異常體繞x軸旋轉(zhuǎn)0°Hz分量虛部；e. 各向異性異常體繞x軸旋轉(zhuǎn)45°Hz分量虛部；f. 各向異性異常體繞x軸旋轉(zhuǎn)90°Hz分量虛部。圖8 各向異性異常體繞x軸旋轉(zhuǎn)HCP裝置磁場Hz分量響應Fig.8 Hz component of the magnetic response of the HCP device when the anisotropic body rotates around the x axis

a.各向同性異常體及各向異性異常體繞x軸旋轉(zhuǎn)不同角度Hz分量實部對比；b. 各向同性異常體及各向異性異常體繞x軸旋轉(zhuǎn)不同角度Hz分量虛部對比。圖9 各向異性異常體繞x軸旋轉(zhuǎn)HCP裝置觀測記錄Hz分量Fig.9 Hz component that is recorded by an observation of the HCP device when an anisotropic anomalous body rotates around the x axis

2.2.3 各向異性圍巖與各向異性異常體的電磁磁場分布特征

當?shù)叵陆橘|(zhì)呈交替分布時，圍巖電性參數(shù)呈各向異性特征。在上述模型(圖5)的基礎(chǔ)上改變圍巖電性參數(shù)各向異性特征，使圍巖沿y方向電導率增大，圍巖主電導率張量為

(29)

同時，各向異性異常體的主電導率張量分別繞z、x軸旋轉(zhuǎn)45°，計算得到異常體在各向異性圍巖內(nèi)的電磁響應 (圖10)。對比圖6和圖8圍巖為各向同性旋轉(zhuǎn)相同角度的磁場Hz分量響應圖可見，圍巖各向異性對磁場響應影響較大，特別是在目標信號能量分布上產(chǎn)生較大變化，由對稱的雙峰或單峰結(jié)構(gòu)變?yōu)椴灰?guī)則的能量團分布。當圍巖沿y軸電導率增大則在該方向呈良導性質(zhì)，而相應的x軸方向呈高阻性質(zhì)，由于電磁波在良導和高阻的情況下有不同的衰減速度，因此可以看出磁場Hz分量平面分布沿著y方向延展。

圖11所示為HCP裝置記錄所得Hz響應曲線。圖中黑實線為常規(guī)各向同性記錄結(jié)果，黑虛線為異常體各向異性旋轉(zhuǎn)45°結(jié)果，黑點劃線為背景介質(zhì)同時為各向異性記錄結(jié)果。由圖11實虛分量繞x、z軸的旋轉(zhuǎn)HCP裝置Hz分量記錄結(jié)果可見，

a．各向異性圍巖-繞z軸旋轉(zhuǎn)各向異性異常體Hz分量實部；b. 各向異性圍巖-繞z軸旋轉(zhuǎn)各向異性異常體Hz分量虛部；c. 各向異性圍巖-繞x軸旋轉(zhuǎn)各向異性異常體Hz分量實部；d. 各向異性圍巖-繞x軸旋轉(zhuǎn)各向異異常體性Hz分量虛部。圖10 在主軸電導率不同的圍巖下各向異性異常體繞z、x軸旋轉(zhuǎn)磁場響應Hz分量Fig.10 Hz component of the magnetic field response of an anisotropic anomalous body rotating around z and x axes under different conductivities of surrounding rock

a.全各向同性模型至全各向異性模型(異常體繞z軸旋轉(zhuǎn))Hz分量實部曲線；b. 全各向同性模型至全各向異性模型(異常體繞z軸旋轉(zhuǎn))Hz分量虛部曲線；c. 全各向同性模型至全各向異性模型(異常體繞x軸旋轉(zhuǎn))Hz分量實部曲線；d. 全各向同性模型至全各向異性模型(異常體繞x軸旋轉(zhuǎn))Hz分量虛部曲線。圖11 在主軸電導率不同的圍巖下各向異性異常體HCP裝置Hz分量記錄結(jié)果Fig.11 Hz component of an anisotropic abnormal body HCP device under different conductivity of surrounding rock

當背景介質(zhì)存在各向異性特征時，記錄信號在形狀上與背景各向同性-目標各向異性類似，但在信號振幅值上有較大的差異，與完全各向同性相比，目標響應信號較弱，說明各向異性造成目標體的識別能力降低。

3 結(jié)論與討論

本文基于矢量有限元法并通過分離總場提出了任意各向異性情況下三維頻率域航空電磁正演模擬算法，通過典型的理論模型與一維解析解對比，驗證了該算法的精度和有效性。針對典型的各向異性電阻率模型，分別計算得到了不同參數(shù)條件下的磁場響應特征并總結(jié)得出以下規(guī)律：

1)磁場分布和主電導率張量的旋轉(zhuǎn)方向之間有很強的相關(guān)性，為識別各向異性異常體的主軸方向和主軸電導率提供依據(jù)，根據(jù)這些性質(zhì)可以進一步識別異常體的層理性質(zhì)，即其走向和傾角。

2)當圍巖呈各向異性時，由于各向異性介質(zhì)中的通道效應，即電流會沿呈導電性質(zhì)的方向聚焦，磁場響應會隨之改變，記錄信號相比于各向同性介質(zhì)在響應振幅和相位上都有較大的變化，且異常體響應信號變?nèi)酢?/p>

3)航空電磁響應受到各向異性介質(zhì)影響較大，但其也有明顯的特征，根據(jù)其分布可以判斷各向異性異常體的電性情況與層理方向，同時也能判斷圍巖的電性情況。本文的計算結(jié)果和各向異性影響特征為實際解釋航空電磁響應結(jié)果提供了依據(jù)。同時各向異性正演算法也為進一步開展航空電磁各向異性反演提供了基礎(chǔ)。

[1] 雷棟，胡祥云，張素芳. 航空電磁法的發(fā)展現(xiàn)狀[J]. 地質(zhì)找礦論叢，2006，21(1)：40-44,53.

Lei Dong, Hu Xiangyun, Zhang Sufang. Development of Airborne Electromagnetic Method[J].Contributions to Geology and Mineral Resources Research, 2006, 21(1):40-44,53.

[2] 王衛(wèi)平. 頻率域航空電磁法發(fā)展研究現(xiàn)狀與應用研究[C]//中國地球物理學會.中國地球物理學會第二十七屆年會論文集.北京：中國地球物理學會，2011：2.

Wang Weiping. Application and Prospect of the Research About Frequeney Domain Airborne Electromagneite Method[C]//Chinese Geophysical Society. Proceedings of the Twenty-Seventh Annual Meeting of China Geophysical Society. Beijing:China Geophysical Society, 2011:2.

[3]Newman G A, Alumbaugh D L.Frequency Domain Modeling of Airborne Electromagnetic Responses Using Staggered Finite Differences[J]. Geophysical Prospecting, 1995, 43:1021-1042.

[4] 朱凱光，林君，劉長勝，等. 頻率域航空電磁法一維正演與探測深度[J]. 地球物理學進展，2008，23(6)：1943-1946.

Zhu Kaiguang, Lin Jun, Liu Changsheng, et al. One-dimensional Forwand and Prospecting Depth for Airborne Frequency Domain Electromagnetic Method[J]. Progress in Geophysics, 2008, 23 (6):1943-1946.

[5] 高亮，胡祥云，王衛(wèi)平，等. 磁性條件下頻率域航空電磁法正演研究[J]. 工程地球物理學報，2009，6(4)：399-403.

Gao Liang, Hu Xiangyun, Wang Weiping, et al. A Study of Frequency Airborne Electromagnetic Forward Under Magnetic Condition[J]. Chinese Journal of Engineering Geophysics, 2009, 6(4):399-403.

[6] 李小康. 基于MPI的頻率域航空電磁法有限元二維正演并行計算研究[D]. 北京：中國地質(zhì)大學，2011.

Li Xiaokang. A MPI Based Parallel Calculation Investigation on Two Dimensional Finite Element Modelling of AEM[D]. Beijing：China University of Geosciences, 2011.

[7] 曲昕馨. 三維頻率域航空電磁法的數(shù)值模擬及姿態(tài)影響和校正研究[D]. 長春：吉林大學，2014.

Qu Xinxin. Numerical Simulation Attitude Effect and Correction of 3-D Frequency Domain Helicopterborne Electromagnetic Method[D]. Changchun: Jilin University, 2014.

[8] 王衛(wèi)平，曾昭發(fā)，李靜，等. 頻率域航空電磁法地形影響和校正方法[J]. 吉林大學學報(地球科學版)，2015，45(3)：941-951.

Wang Weiping, Zeng Zhaofa, Li Jing, et al. Topographic Effects and Correction for Frequency Airborne Electromagnetic Method[J]. Journal of Jilin University(Earth Science Edition), 2015, 45(3):941-951.

[9] 殷長春，張博，劉云鶴，等. 2.5維起伏地表條件下時間域航空電磁正演模擬[J]. 地球物理學報，2015，58(4)：1411-1424.

Yin Changchun, Zhang Bo, Liu Yunhe, et al. 2.5-D Forward Modeling of the Time-Domain Airborne EM Aystem in Areas with Topographic Relief[J]. Chinese Journal of Geophysics, 2015, 58(4):1411-1424.

[10] 黃威，殷長春，賁放，等. 頻率域航空電磁三維矢量有限元正演模擬[J]. 地球科學，2016，41(2)：331-342.

Huang Wei, Yin Changchun, Ben Fang, et al. 3D Forward Modelling for Frequency AEM by Vector Finite-Element[J].Earth Science, 2016, 41(2):341-342.

[11] Li Wenben, Zeng Zhaofa, Li Jing,et al. 2.5D Fo-rward Modeling and In-Version of Frequency-Domain Airborne Electromagnetic Data[J]. Applied Geophysics, 2016, 13(1):37-47,218.

[12] Yin C, Weidelt P. Geoelectrical Fields in a Layered Earth with Arbitrary Anisotropy[J]. Geophysics, 1999, 64(2): 426-434.

[13] Yin C, Maurer, H.M. Electromagnetic Induction in a Layered Earth with Arbitrary Anisotropy[J]. Geophysics, 2001, 66(5):1405-1416.

[14] Avdeev D B, Kuvshinov A V, Pankratov O V, et al. Three-Dimensional Frequency-Domain Modeling of Airborne Electromagnetic Responses[J]. Exploration Geophysics, 1998, 29(2):111-119.

[15] Yin C, Fraser D C. The Effect of the Electrical Ani-sotropy on the Response of Helicopter-Borne Frequency-Domain Electromagnetic Systems[J]. Geophysical Prospecting, 2004, 52(5):399-416.

[16] Yin C, Hodges G. Simulated Annealing for Airborne EM Inversion[J]. Geophysics, 2007, 72(4):F189-F196.

[17] Liu Y, Yin C.3D Anisotropic Modeling for Airborne EM Systems Using Finite-Difference Method[J]. Journal of Applied Geophysics, 2014, 109: 186-194.

[18] Yin C, Yan F, Liu Y. 3D Time-Domain Airborne EM Modeling for an Arbitrarily Anisotropic Earth[J]. Journal of Applied Geophysics, 2016, 131: 163-178.

[19]Huang X, Yin C,Cao X Y, et al. 3D Anisotropic Modeling and Identification for Airborne EM Systems Based on the Spectral-Element Method[J]. Applied Geophysics, 2017, 14(3):419-430.

[20] Yin C. Geoelectrical Inversion for a One-Dimensional Anisotropic Model and Inherent Non-Uniqueness[J]. Geophys J Int, 2000, 140(1):11-23.

[21] Onsager L. Reciprocal Relations in Irreversible Pro-cesses[J]. Physical Review, 1931, 37(37):405-426.

[22] Mitsuhata Y. 2-D Electromagnetic Modeling by Finite Element Method with a DipoleSource and Top--Ography[J]. Geophysics, 2000, 65(2):465-475.

[23] 金建銘. 電磁場有限元方法[M]. 西安：西安電子科技大學出版社，1998.

Jin Jianming. Electromagnetic Finite Element Method[M]. Xi’an:Xidian University Press, 1998.