本發(fā)明屬于地球物理正演，具體涉及一種球坐標系任意各向異性電導率張量大地電磁仿真方法。

背景技術(shù)：

1、地球的活動和演化不僅給人類帶來資源同時也伴隨自然災害，研究地球相關的科學問題對人類可持續(xù)發(fā)展具有重要意義。地球物理是一種通過物理場探測地球物質(zhì)結(jié)構(gòu)分布的有效方法，一直被廣泛應用于災害監(jiān)測、資源勘探、地球科學問題研究。地球物理領域針對地幔尺度的探測方法主要包括地震觀測和大地電磁測深。其中大地電磁測深法利用天然場源進行對地探測，具有勘探深度大、成本低等諸多優(yōu)勢，是地球深部探測的一把利刃。如今，為了充分發(fā)揮大地電磁法的探測優(yōu)勢，實現(xiàn)地球深部的電性結(jié)構(gòu)精準成像，國際上已部署了多項洲際尺度大地電磁測深計劃。洲際尺度的大地電磁觀測可在大尺度上實現(xiàn)地下介質(zhì)的整體成像，減少多解性，但同時也為數(shù)據(jù)的反演帶來了新的挑戰(zhàn)。

2、正演模擬是反演的基礎，隨著大地電磁觀測尺度越來越大，球坐標系下的大地電磁模擬成為地球物理研究領域的熱點與難點。首先是大尺度模擬往往需要離散化求解，基于交錯網(wǎng)格的有限差分法的球坐標正演發(fā)展較快。眾多學者分析了球坐標系下大地電磁響應和笛卡爾坐標系的下響應差距，一致認為地球曲率對于大尺度長周期的大地電磁數(shù)據(jù)影響不容忽視?；诎瞬鏄渚W(wǎng)格的有限元法實現(xiàn)了球坐標系大地電磁正演。然而，考慮到復雜的海陸邊界、地核以及電離層對電磁響應會產(chǎn)生顯著影響，全球尺度的大地電磁正演需要多尺度的模型離散技術(shù)。

3、非結(jié)構(gòu)四面體可以靈活地擬合多種尺度的結(jié)構(gòu)，是全球大地電磁三維正演的優(yōu)選方法。其次隨著模擬尺度的增加，場源的形式發(fā)生了顯著變化。傳統(tǒng)的大地電磁場源被假設來自遠離地球表面的電離層系統(tǒng)（4-10倍地球半徑處），對于小尺度測區(qū)，可以認為是均勻的平面極化波垂直入射到地球介質(zhì)。然而，當測區(qū)為洲際尺度時，這一假設將不再適用，此時應建立符合電離層結(jié)構(gòu)的球坐標系場源?；谟邢薏罘炙惴▽崿F(xiàn)了球坐標系各向同性介質(zhì)大地電磁模擬，并對比了平面波場源和非均一電場場源的響應差異，認為場源會影響長周期數(shù)據(jù)。最后，地幔存在不均勻性導致電各向異性會顯著影響大地電磁觀測信號，而現(xiàn)有的球坐標系大地電磁方法基于各向同性假設容易帶來錯誤的結(jié)果，在球坐標系下進行電各向異性介質(zhì)模擬的研究較少，其一在球坐標系下沿用笛卡爾坐標系下的電導率各向異性張量模擬地磁測深響應，但這顯然具有嚴重偏差，其各向異性張量元素沿著經(jīng)緯度方向是不連續(xù)的。

4、因此，為了準確模擬球坐標系大地電磁場以適應深地觀測中洲際大地電磁陣列數(shù)據(jù)處理需求，亟需研發(fā)球坐標系任意各向異性介質(zhì)大地電磁三維仿真方法，服務于未來大尺度陣列數(shù)據(jù)反演與深地科學研究。

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

1、針對現(xiàn)有大地電磁仿真技術(shù)無法靈活模擬洲際尺度任意各向異性介質(zhì)大地電磁響應的技術(shù)問題，本發(fā)明提出一種球坐標系任意各向異性電導率張量大地電磁仿真方法，其采用非結(jié)構(gòu)四面體網(wǎng)格靈活刻畫真實地球環(huán)境并采用球諧函數(shù)模擬大地電磁正交場源，提出了新的球坐標系任意各向異性介質(zhì)電導率張量及其笛卡爾坐標系等效形式，最后基于矢量有限元方法給出仿真方程以及求解方法。為洲際尺度陣列大地電磁數(shù)據(jù)反演研究提供正演技術(shù)，未來有望服務于深地科學研究。

2、為實現(xiàn)上述發(fā)明目的，實施例提供的一種球坐標系任意各向異性電導率張量大地電磁仿真方法，包括以下步驟：

3、根據(jù)大地電磁仿真目標確定測點陣列的坐標和觀測頻率；

4、構(gòu)建地球背景模型并采用非結(jié)構(gòu)四面體網(wǎng)格進行離散得到非結(jié)構(gòu)網(wǎng)格；

5、應用球坐標任意各向異性電導率張量給非結(jié)構(gòu)網(wǎng)格的電性屬性賦值；

6、對非結(jié)構(gòu)網(wǎng)格的外層網(wǎng)格施加正交的非平面波場源，并計算非平面波場源的源項積分；

7、計算非結(jié)構(gòu)網(wǎng)格中單元的質(zhì)量矩陣和剛度矩陣，并根據(jù)編號關系形成總體單元矩陣；

8、根據(jù)總體單元矩陣和源項積分并按照當前觀測頻率構(gòu)建對應的仿真方程組；

9、求解仿真方程組得到所有單元棱邊的電場，并基于電場通過插值函數(shù)計算大地電磁響應數(shù)據(jù)，實現(xiàn)大地電磁仿真。

10、優(yōu)選地，應用球坐標任意各向異性電導率張量給非結(jié)構(gòu)網(wǎng)格的電性屬性賦值，包括：

11、給出其球坐標下的表達式為：

12、；

13、其中上標 s表示物理量在球坐標下，即表示球坐標下的各向異性電導率張量，公式右端為的展開形式，其中每個分量的下標指示沿著不同球坐標系方向提供其對應的標量電導率，將利用空間歐拉旋轉(zhuǎn)進一步把張量轉(zhuǎn)換成用三個球坐標主軸電導率和三個球坐標主軸旋轉(zhuǎn)角來描述：

14、；

15、其中上標t表示矩陣的轉(zhuǎn)置，而，，為球坐標下三個主軸、以及對應的空間歐拉旋轉(zhuǎn)矩陣，表示為：

16、，，；

17、球坐標系任意各向異性電導率張量滿足歐姆定律，此處上標 s表示物理量在球坐標下，表示球坐標下的電流密度，表示球坐標下的電場；

18、給出了球坐標系任意各向異性電導率張量的笛卡爾坐標系等效形式，此處上標 c表示物理量在笛卡爾坐標系下，表示笛卡爾坐標下的電流密度，表示笛卡爾坐標下的電場，得到球坐標系任意各向異性電導率張量的笛卡爾等效形式為：

19、；

20、其中表示笛卡爾坐標系下的各向異性電導率張量，和為兩個方向的旋轉(zhuǎn)矩陣其表達式如下：

21、，；

22、基于此，賦值的非結(jié)構(gòu)網(wǎng)格的電性屬性包括。

23、優(yōu)選地，對非結(jié)構(gòu)網(wǎng)格的外層網(wǎng)格施加正交的非平面波場源，包括：

24、給定兩類正交的磁場 h作為外部邊界的場源，假定空氣絕緣，磁場用磁標勢 u的負梯度表示為：

25、；

26、上標 ext表示物理量在網(wǎng)格的外邊界上，即表示外邊界的磁場，表示外邊界的磁標勢，表示對物理量求梯度，為真空中的磁導率，在球坐標系下的磁標勢 u表示成球諧函數(shù)：

27、；

28、其中表示虛數(shù)單位，為地球半徑，表示外源系數(shù)，exp表示指數(shù)函數(shù)計算， p為連帶勒讓德多項式，m和n分別為多項式的次數(shù)和階數(shù)，本發(fā)明取靜日狀態(tài)下的場分布即m=0,n=1，將r取5倍地球半徑得到外邊界的第一類場源為：

29、；

30、其中，和以及為在球坐標系下三個方向的磁場，根據(jù)球坐標和笛卡爾坐標的轉(zhuǎn)換關系，笛卡爾坐標系下的第一類場源（，，）為：

31、；

32、第二類場源（，，）取第一類場源的正交場為：

33、；

34、其中，和以及為在笛卡爾坐標系下第一類場源三個方向的磁場，和以及為在笛卡爾坐標系下第二類場源三個方向的磁場。

35、優(yōu)選地，計算非平面波場源的源項積分，包括：

36、采用高斯積分來計算右端源項的積分，積分項內(nèi)部看成只與空間位置有關的函數(shù)表示成f，應用空間三角形的積分節(jié)點得源項積分為：

37、；

38、其中，上標1，2表示第幾類場源，即和分別表示第1類場源和第2類場源對應的源項積分，a為加權(quán)系數(shù)，下標 k表示積分節(jié)點索引。

39、優(yōu)選地，計算非結(jié)構(gòu)網(wǎng)格中單元的質(zhì)量矩陣和剛度矩陣，并根據(jù)編號關系形成總體單元矩陣，包括：

40、應用矢量有限元法對大地電磁場滿足的電場雙旋度方程進行分析，得到單元ele的質(zhì)量矩陣和剛度矩陣為：

41、；

42、；

43、其中，下標 i,j表示棱邊在單元ele的局部編號，v表示單元ele的體積，表示第 j條棱邊相對的第 i條棱邊對應的矢量插值基函數(shù)，表示第 j條棱邊的矢量插值基函數(shù)，×表示對矢量插值基函數(shù)求旋度，表示球坐標系任意各向異性電導率張量，每個單元的棱邊都有全局編號，因此形成總體單元矩陣s和m。

44、優(yōu)選地，根據(jù)總體單元矩陣和源項積分并按照當前觀測頻率構(gòu)建對應的仿真方程組，包括：

45、根據(jù)當前觀測頻率 f計算角頻率，為圓周率，對應的復頻率系數(shù)為，其中表示虛數(shù)單位，令矩陣，則當前頻率下的仿真方程組按下式計算：

46、；

47、其中，s為總體單元矩陣中的總體剛度矩陣，m為總體質(zhì)量矩陣，右端項中的向量b中的元素為源項積分，向量b中只有外邊界棱邊對應位置才有值其余為0，上標1，2表示第幾類場源，e表示電場。

48、優(yōu)選地，求解仿真方程組得到所有單元棱邊的電場，包括：

49、當仿真單元的個數(shù)少且計算資源足夠時采用直接求解器進行求解，其原理為對仿真方程組的k進行l(wèi)u分解，對于不同的方程右端項進行回代即得到所有棱邊的電場e1、e2。

50、優(yōu)選地，求解仿真方程組得到所有單元棱邊的電場，包括：

51、當仿真單元的個數(shù)多且計算資源有限時采用迭代法并行求解，采用基于附屬空間預條件子及塊狀對角右預條件的靈活廣義最小殘差迭代法，首先將原來仿真方程組呈現(xiàn)的復數(shù)方程轉(zhuǎn)化成等效實數(shù)形式，去掉不同場源的上標，給出一般形式：

52、；

53、其中上標 re表示物理量的實部，上標 im表示物理量的虛部，記一般形式的方程為kr·er=br，塊狀對角化矩陣p為：

54、；

55、將p的逆矩陣乘于一般形式方程的系數(shù)矩陣右側(cè)，則原方程為krp-1?·u=br，向量u為一個中間向量，求解過程分為內(nèi)外兩層循環(huán)，內(nèi)層用附屬空間預條件子及共軛梯度法求解向量，q是中間向量u的迭代更新步長向量，外層循環(huán)利用塊狀對角右預條件對應的塊狀對角化矩陣p及廣義最小殘差法求解電場er。

56、優(yōu)選地，基于電場通過插值函數(shù)計算大地電磁響應數(shù)據(jù)，實現(xiàn)大地電磁仿真，包括：

57、基于電場通過插值函數(shù)計算大地電磁響應數(shù)據(jù)包括當前觀測頻率的大地電磁張量阻抗響應、視電阻率、以及相位。

58、與現(xiàn)有技術(shù)相比，本發(fā)明具有的有益效果至少包括：

59、本發(fā)明采用非結(jié)構(gòu)四面體網(wǎng)格能夠靈活刻畫真實地球環(huán)境和精細刻畫海岸線，同時其具有斜邊的特性能夠模擬非規(guī)則異常體（如俯沖板塊，褶皺造山帶，地幔柱等），本發(fā)明采用球諧函數(shù)構(gòu)建非均一的正交場源打破傳統(tǒng)平面波場源無法適應地球曲率問題，同時本發(fā)明提出一種新的球坐標系任意各向異性介質(zhì)電導率張量，發(fā)明張量沿著球坐標系是連續(xù)分布的，解決了傳統(tǒng)笛卡爾坐標系電導率張量不連續(xù)的問題，同時發(fā)明基于歐拉旋轉(zhuǎn)給出了張量在笛卡爾坐標下的等效形式，使得基于非結(jié)構(gòu)矢量有限元方法能夠應用于構(gòu)建球坐標系任意各向介質(zhì)大地電磁仿真方程組，最后發(fā)明給出了電磁仿真方程組的求解方法，根據(jù)本發(fā)明提供的一個具體實施例可以看出發(fā)明提出的仿真技術(shù)能夠準確有效模擬球坐標任意各向異性介質(zhì)的大地電磁響應，這給未來進行洲際尺度大地電磁數(shù)據(jù)反演及深地研究提供了正演技術(shù)。