本發(fā)明屬于地球物理勘探領(lǐng)域，具體涉及一種適用于復(fù)雜地質(zhì)條件下的航空自然場源電磁法三維正演模擬方法。

背景技術(shù)：

1、高原、山地、冰川、沼澤等地質(zhì)環(huán)境復(fù)雜的區(qū)域，地面勘查工作難于開展。航空自然場電磁法通過利用天然的正交場源電磁探測地下電性結(jié)構(gòu)，是一種新興的地球物理勘查手段。航空自然場電磁法通常采用直升機(jī)或無人機(jī)平臺(tái)搭載磁傳感器，在空中觀測磁場的不同分量，同時(shí)在地面設(shè)置少量的基站，觀測電場或磁場分量。由于采用航空數(shù)據(jù)采集方式，該方法具有高效、便捷的優(yōu)勢，非常適合應(yīng)用于復(fù)雜地貌環(huán)境和大面積勘察，在礦產(chǎn)勘探以及工程勘察等領(lǐng)域具有巨大的發(fā)展?jié)摿Α：娇兆匀粓鲭姶欧ㄈS正演方法是三維反演的基礎(chǔ)，但目前航空自然場電磁法三維正演研究存在以下不足：(1)?航空自然場電磁法不同于傳統(tǒng)的地面電磁觀測方式，現(xiàn)有地面大地電磁三維正演模擬方法不能直接用于航空自然源電磁法模擬；(2)航空自然源電磁法的應(yīng)用區(qū)域通常具有復(fù)雜的地形和地質(zhì)結(jié)構(gòu)，但目前尚未有專用于復(fù)雜地質(zhì)結(jié)構(gòu)的航空自然場源電磁法三維正演模擬方法；?(3)?巖石地層中廣泛存在各向異性與高磁導(dǎo)率介質(zhì)，現(xiàn)有的航空自然源電磁三維正演法中未考慮這兩種現(xiàn)象，從而影響航空自然源電磁正演模擬的計(jì)算準(zhǔn)確性。

2、中國專利申請cn118884546a公開了一種基于物理信息神經(jīng)網(wǎng)絡(luò)的航空電磁正演模擬方法及介質(zhì)，包括：確定航空電磁勘探區(qū)域的三維立體模型；基于三維立體模型劃定多個(gè)航空電磁測點(diǎn)；構(gòu)建物理信息神經(jīng)網(wǎng)絡(luò)；以不同航空電磁測點(diǎn)的三維坐標(biāo)、邊界測點(diǎn)的二次電場值、航空電磁系統(tǒng)參數(shù)為輸入，以損失函數(shù)的均方誤差最小為目標(biāo)對物理信息神經(jīng)網(wǎng)絡(luò)進(jìn)行訓(xùn)練優(yōu)化；將目標(biāo)航空電磁測點(diǎn)的三維坐標(biāo)輸入訓(xùn)練優(yōu)化后的物理信息神經(jīng)網(wǎng)絡(luò)進(jìn)行頻率域航空電磁三維正演模擬，得到目標(biāo)航空電磁測點(diǎn)的二次電場值；利用法拉第電磁感應(yīng)定律計(jì)算目標(biāo)航空電磁測點(diǎn)的二次磁場值。專利cn118884546a所提出的方法用于可控源航空頻率域電磁法，不適用于航空自然源電磁法。

3、中國專利申請cn114880902a公開了一種磁性源瞬變電磁法時(shí)域矢量有限元正演方法及裝置，方法包括：建立模擬地質(zhì)模；對網(wǎng)格剖分后的模擬地質(zhì)模型施加激勵(lì)源，獲取模擬地質(zhì)模型的麥克斯韋方程組得到電場控制方程；采用加權(quán)余量法離散控制方程并選擇第一類whitney基函數(shù)建立有限元方程；選擇時(shí)間差分格式對有限元方程進(jìn)行時(shí)間離散；組裝總體剛度矩陣，并調(diào)用padiso求解器求解大型線性稀疏方程組；根據(jù)激勵(lì)源的電流關(guān)斷特征確定迭代的時(shí)間步長，再以可變的時(shí)間步長進(jìn)行迭代得到電場各個(gè)分量，并通過基函數(shù)得到各個(gè)節(jié)點(diǎn)電場值；根據(jù)法拉第電磁感應(yīng)定律，求得瞬變電磁觀測值和視電阻率參數(shù)。cn114880902a公開的方法適用于磁性源瞬變電磁法，而不適用于航空自然源電磁法。

技術(shù)實(shí)現(xiàn)思路

1、本發(fā)明的目的是為克服上述現(xiàn)有技術(shù)的不足，提供一種基于非結(jié)構(gòu)化四面體網(wǎng)格建模技術(shù)與矢量有限元方法的航空自然源電磁法復(fù)雜介質(zhì)并行三維正演模擬方法，彌補(bǔ)了尚未有方法統(tǒng)一解決航空自然源電磁正演模擬中復(fù)雜地質(zhì)結(jié)構(gòu)、電導(dǎo)率各向異性與磁導(dǎo)率影響難題的現(xiàn)狀，從而可以獲取真實(shí)的復(fù)雜地球介質(zhì)的航空自然場源電磁響應(yīng)。

2、為實(shí)現(xiàn)上述目的，本發(fā)明采用下述技術(shù)方案：

3、一種航空自然源電磁法復(fù)雜介質(zhì)并行三維正演模擬方法，包括以下步驟：

4、s1：構(gòu)建研究區(qū)域的三維地質(zhì)幾何模型；基于cad隱式建模方法，利用地形、地層和地質(zhì)體數(shù)據(jù)構(gòu)建研究區(qū)域的三維地質(zhì)幾何模型；

5、s2：生成模型非結(jié)構(gòu)化四面體網(wǎng)格；使用非結(jié)構(gòu)化四面體網(wǎng)格對幾何模型進(jìn)行剖分，在地面基站和空中的測點(diǎn)處進(jìn)行局部加密，核心區(qū)域以外的采用粗網(wǎng)格進(jìn)行剖分，從而生成整個(gè)研究區(qū)域地質(zhì)模型的非結(jié)構(gòu)化四面體網(wǎng)格；

6、s3：模型非結(jié)構(gòu)化四面體網(wǎng)格各向異性電導(dǎo)率張量和磁導(dǎo)率設(shè)置；根據(jù)已知的地質(zhì)、鉆孔和地球物理等信息，為模型四面體網(wǎng)格單元分配任意各向異性電導(dǎo)率張量和磁導(dǎo)率。

7、構(gòu)建頻率域電場控制方程；頻率域麥克斯韋方程組消去磁場項(xiàng)后可得到以下矢量電場控制方程公式（1）：

8、公式（1）：，

9、其中,其中,i=，ω表示角頻率，?為電場，為磁導(dǎo)率，為電導(dǎo)率張量,，、、、、、、、和分別為各方向電導(dǎo)率分量；

10、s5：采用矢量有限元離散電場控制方程；利用galerkin有限元分析對頻域電場控制方程公式?(1)?進(jìn)行處理，得到微分方程的弱形式積分方程:

11、公式（2）：

12、其中，?ω?表示模擬的區(qū)域，??為矢量插值基函數(shù)，??為研究區(qū)域內(nèi)待求解的電場。

13、采用一階線性矢量基函數(shù)離散非結(jié)構(gòu)化四面體內(nèi)的電場，四面體單元內(nèi)的電場?可以由每條棱邊中心的電場?定義：

14、公式（3）：，

15、其中表示第e個(gè)四面體單元，j表示四面體內(nèi)棱邊的局部編號(hào)。

16、利用第一個(gè)矢量格林定理在每個(gè)四面體元素中得到公式（2）的離散形式，利用模型非結(jié)構(gòu)化四面體網(wǎng)格、電導(dǎo)率和磁導(dǎo)率參數(shù)組裝成稀疏的有限元系數(shù)方程組：

17、公式（4）：，

18、其中，?為施加邊界條件后的方程右端項(xiàng)，?和?分別為剛度矩陣和質(zhì)量矩陣，具體離散的形式為，

19、公式（5）：，

20、公式（6）：，

21、其中，?和?分別表示剛度矩陣和質(zhì)量矩陣，?為矢量基函數(shù)，?為求解的單元區(qū)域，為體積微元，??和??分別表示單元的磁導(dǎo)率和電阻率張量，上標(biāo)“e”表示第e個(gè)四面體單元，下標(biāo)“j”，“k”表示單元的棱邊編號(hào)，采用解析方法計(jì)算以上單元積分。

22、s6：設(shè)置邊界條件；設(shè)置狄利克雷邊界條件。模型的頂部施加激勵(lì)源和?，模型的側(cè)向和底部邊界的電場可以通過一維正演得到。設(shè)在三維模型邊界的棱邊的中心坐標(biāo)為，通過一維正演算法可以得到棱邊中心的電場，將一維解析解投影到三維模型邊界的棱邊中心：

23、公式（7）：，

24、其中，是三維模型邊界棱邊的切向單位矢量，?是棱邊中心電場，是一維正演得到的電場分量。

25、采用“01”賦值法施加邊界條件，得到待求解的稀疏的有限元線性方程組：

26、公式（8）：，

27、其中，和??兩種極化方式的有限元系數(shù)矩陣??相同，而方程右端的源向?不同。

28、s7：求解有限元線性方程；對不同計(jì)算頻率采用并行計(jì)算，采用直接求解器pardiaso求解每個(gè)計(jì)算頻率的有限元大型線性方程組，先對系數(shù)矩陣進(jìn)行矩陣分解，進(jìn)行兩次回代后得到兩種極化模式下的解?和?。

29、s8：計(jì)算測點(diǎn)、基站處電場和磁場三分量；求解測點(diǎn)電場分量和磁場分量，測點(diǎn)在第 e個(gè)四面體單元中時(shí)，測點(diǎn)電場?可以由矢量形函數(shù)表示為：

30、公式（9）：

31、根據(jù)法拉第定律，測點(diǎn)磁場分量?由電場和矢量形函數(shù)表示為：

32、公式（10）：。

33、s9：計(jì)算測點(diǎn)處阻抗、視電阻率、相位、傾子數(shù)據(jù)；航空自然場源電磁法的響應(yīng)數(shù)據(jù)包括阻抗張量、視電阻率、相位、傾子矢量。航空自然場源電磁法的阻抗張量各個(gè)分量的計(jì)算表達(dá)式為：

34、公式（11）：，

35、，

36、，

37、，

38、其中，，，，為空中觀測的觀測磁場水平分量，，，，為地面固定基站觀測電場水平分量，下標(biāo)“x”，“y”代表兩個(gè)水平方向，下標(biāo)“1”，“2”代表兩種極化方式。

39、在計(jì)算阻抗張量后可以獲得視電阻率和相位，即：

40、公式（12）：

41、公式（13）：

42、傾子矢量的計(jì)算表達(dá)式為：

43、公式（14）：,?，

44、當(dāng)在空中只觀測磁場垂直分量時(shí)，，為空中磁場垂直分量，，，，為地面固定基站觀測磁場水平分量；當(dāng)在空中觀測磁場三分量時(shí)，，，，，，均為空中觀測的磁場分量。下標(biāo)“x”，“y”，“z”代表三分量方向，下標(biāo)“1”，“2”代表兩種極化方式。

45、優(yōu)選的，步驟s4所述的矢量插值基函數(shù)采用一階線性基函數(shù)表示：?其中?l?為節(jié)點(diǎn)基函數(shù)，l為棱邊的長度，下標(biāo)“”和“”代表棱邊的起點(diǎn)和終點(diǎn)。

46、優(yōu)選的，步驟s6所述的并行計(jì)算方式包括兩層并行，第一層是為每個(gè)頻率的有限元系數(shù)矩陣計(jì)算分配一個(gè)計(jì)算進(jìn)程，第二層是使用pardiaso求解器內(nèi)部openmp線程并行。

47、本發(fā)明提出了涉及一種適用于復(fù)雜地質(zhì)結(jié)構(gòu)的航空自然場源電磁法三維正演方法，克服了現(xiàn)有航空自然源電磁法正演技術(shù)無法綜合考慮復(fù)雜地質(zhì)結(jié)構(gòu)、各向異性介質(zhì)與磁導(dǎo)率影響的局限性，其創(chuàng)新點(diǎn)及優(yōu)勢如下：

48、1.?本發(fā)明可以準(zhǔn)確模擬各向異性電導(dǎo)率和磁導(dǎo)率雙參數(shù)對航空自然源電磁數(shù)據(jù)的影響，適用于真實(shí)地球介質(zhì)中廣泛存在各向異性與磁性介質(zhì)影響。

49、2.本發(fā)明采用三維地質(zhì)建模技術(shù)與基于非結(jié)構(gòu)化四面體網(wǎng)格矢量有限元方法，可以精確模擬具有復(fù)雜形態(tài)的真實(shí)地質(zhì)結(jié)構(gòu)的航空自然源電磁響應(yīng)。

50、3.本發(fā)明采用并行技術(shù)計(jì)算不同頻率觀測數(shù)據(jù)，具有較高的計(jì)算速度。

51、本發(fā)明為航空自然場源電磁法提供了三維正演模擬算法，為航空自然場源電磁法理論的研究以及實(shí)際觀測數(shù)據(jù)的反演解釋提供了基礎(chǔ)，進(jìn)而可以讓航空自然場源電磁法更好的服務(wù)于礦產(chǎn)勘察、工程勘察等應(yīng)用領(lǐng)域。