一種時(shí)域有限差分的三維感應(yīng)-極化雙場(chǎng)數(shù)值模擬方法
【專利摘要】本發(fā)明涉及一種時(shí)域有限差分的三維感應(yīng)?極化雙場(chǎng)數(shù)值模擬方法，目的在于快速計(jì)算三維模型的感應(yīng)?極化雙場(chǎng)的電磁響應(yīng)。主要包括采用逆拉普拉斯變換獲得迪拜模型電導(dǎo)率的時(shí)域表達(dá)式，構(gòu)建電導(dǎo)率參數(shù)的e指數(shù)輔助方程，通過(guò)梯形積分法獲得歐姆定律時(shí)域離散遞推表達(dá)式，由四維數(shù)值運(yùn)算降低為三維運(yùn)算，再將其代入無(wú)源Maxwell旋度方程中，基于三維時(shí)域有限差分方法推導(dǎo)電場(chǎng)和磁場(chǎng)的迭代方程，進(jìn)而完成三維模型的感應(yīng)?極化雙場(chǎng)電磁響應(yīng)數(shù)值計(jì)算。本發(fā)明目的在于解決歐姆定律的時(shí)域卷積運(yùn)算耗時(shí)長(zhǎng)、內(nèi)存占用大等問(wèn)題，最終實(shí)現(xiàn)三維模型的感應(yīng)?極化雙場(chǎng)電磁響應(yīng)快速數(shù)值計(jì)算。
【專利說(shuō)明】
一種時(shí)域有限差分的三維感應(yīng)-極化雙場(chǎng)數(shù)值模擬方法
技術(shù)領(lǐng)域
[0001] 本發(fā)明涉及一種地球物理勘探領(lǐng)域中的時(shí)域電磁場(chǎng)數(shù)值模擬方法，尤其適用快速 計(jì)算三維模型的感應(yīng)-極化雙場(chǎng)電磁響應(yīng)。
【背景技術(shù)】
[0002] 瞬變電磁法(Time-domain Electromagnetic，簡(jiǎn)稱TEM)是基于電磁感應(yīng)禍流效應(yīng) 為主的地下近地表探測(cè)方法。當(dāng)?shù)皖l時(shí)忽略位移電流和極化電荷產(chǎn)生的電磁場(chǎng)，僅利用地 下介質(zhì)的禍流效應(yīng)觀測(cè)二次瞬變場(chǎng)。時(shí)域激發(fā)極化法（Time-domain Induced Polarization，簡(jiǎn)稱IP)是基于地下介質(zhì)極化效應(yīng)為主的地下淺層探測(cè)方法，主要采用電性 源施加直流或交流電磁場(chǎng)激勵(lì)下，忽略傳導(dǎo)電流產(chǎn)生的電磁場(chǎng)以及高頻介電常數(shù)引起的極 化響應(yīng)，僅利用時(shí)域電流斷開(kāi)后低頻激勵(lì)下地下介質(zhì)產(chǎn)生的極化場(chǎng)，觀測(cè)極化電荷產(chǎn)生的 二次電位變化過(guò)程。在時(shí)域激發(fā)極化測(cè)量中總是去除感應(yīng)場(chǎng)的耦合。
[0003] 但是，無(wú)論是瞬變電磁還是時(shí)域激發(fā)極化方法，都是以電磁場(chǎng)麥克斯韋方程為理 論依據(jù)并進(jìn)行波場(chǎng)簡(jiǎn)化近似，實(shí)際探測(cè)中感應(yīng)場(chǎng)和極化場(chǎng)同時(shí)存在，互相伴生。無(wú)論瞬變電 磁還是時(shí)域激發(fā)極化方法，將兩種電流場(chǎng)割裂開(kāi)分別進(jìn)行研究，是不符合實(shí)際電磁場(chǎng)的傳 播規(guī)律的，可能導(dǎo)致實(shí)測(cè)數(shù)據(jù)與理論計(jì)算結(jié)果不符，甚至導(dǎo)致解釋結(jié)果出現(xiàn)錯(cuò)誤。為此，只 有準(zhǔn)確計(jì)算地下結(jié)構(gòu)的感應(yīng)和極化雙場(chǎng)電磁響應(yīng)，理論模擬結(jié)果與實(shí)際測(cè)量數(shù)據(jù)互相吻 合。
[0004] 中國(guó)專利CN104408021A公開(kāi)了一種電偶源三維時(shí)域有限差分正演成像方法，對(duì)海 洋空氣、海水和海底大地三分空間均建立麥克斯韋方程組和本構(gòu)方程，采用時(shí)域有限差分 法得到海水和海底大地中任意時(shí)刻電磁場(chǎng)的分布。
[0005] 中國(guó)專利CN105277980A公開(kāi)了一種高精度空間和時(shí)間任意倍數(shù)可變網(wǎng)格有限差 分正演方法，通過(guò)建立地下介質(zhì)的正演速度模型，對(duì)正演速度模型中的聲波波場(chǎng)進(jìn)行二維 網(wǎng)格離散化，同時(shí)對(duì)最佳匹配層邊界條件進(jìn)行網(wǎng)格離散化，進(jìn)而通過(guò)聲波波動(dòng)方程進(jìn)行時(shí) 域有限差分正演模擬。
[0006] 加拿大專利CA2388271公開(kāi)了一種基于時(shí)域有限差分方法的電磁場(chǎng)計(jì)算方法，并 提出了一種基于二維導(dǎo)體的網(wǎng)格劃分方法和數(shù)據(jù)儲(chǔ)存方式，但只涉及感應(yīng)場(chǎng)的數(shù)值計(jì)算。
[0007] 加拿大，D Marchant(2015，University of British Columbia)米用有限元方法， 在時(shí)間域直接求解并模擬計(jì)算了三維空間下存在激發(fā)極化效應(yīng)的電磁響應(yīng)。
[0008] 以上所述方法公布了基于時(shí)域有限差分(FDTD)方法的電磁響應(yīng)計(jì)算方法，國(guó)內(nèi)外 專利還未涉及時(shí)域有限差分方法計(jì)算三維模型的感應(yīng)和極化雙場(chǎng)的電磁響應(yīng)方面的，為 此，本發(fā)明基于H)TD方法直接計(jì)算三維模型的感應(yīng)和極化雙場(chǎng)電磁響應(yīng)，并解決歐姆定律 的時(shí)域卷積計(jì)算耗時(shí)長(zhǎng)、內(nèi)存占用大等問(wèn)題，最終實(shí)現(xiàn)三維模型的感應(yīng)-極化雙場(chǎng)電磁響應(yīng) 快速數(shù)值計(jì)算。

【發(fā)明內(nèi)容】

[0009] 本發(fā)明所要解決的關(guān)鍵問(wèn)題是提供一種時(shí)域有限差分的三維感應(yīng)-極化雙場(chǎng)數(shù)值 模擬方法，先獲得迪拜模型的電導(dǎo)率時(shí)域形式，通過(guò)構(gòu)建電導(dǎo)率參數(shù)e指數(shù)輔助方程，并采 用梯形積分法獲得歐姆定律的時(shí)域離散遞推形式，最后基于H)TD方法，實(shí)現(xiàn)了時(shí)域內(nèi)快速 計(jì)算三維模型的感應(yīng)-極化雙場(chǎng)電磁響應(yīng)。
[0010] 本發(fā)明是這樣實(shí)現(xiàn)的，一種時(shí)域有限差分的三維感應(yīng)-極化雙場(chǎng)數(shù)值模擬方法包 括：
[0011] 1)、基于迪拜模型(頻率相關(guān)系數(shù)c = l時(shí)的Cole-Cole復(fù)電阻率模型），先獲得迪拜 模型電導(dǎo)率的頻域形式，再通過(guò)逆拉普拉斯變換，得到電導(dǎo)率的時(shí)域表達(dá)式 〇(t);
[0012] 2)、將電導(dǎo)率的時(shí)域表達(dá)式〇(t)代入歐姆定律中，獲得歐姆定律時(shí)域卷積積分表 達(dá)式，結(jié)合電導(dǎo)率〇(t)近似于aeT et形式的特點(diǎn)，構(gòu)建電導(dǎo)率參數(shù)的e指數(shù)輔助方程r(t)，代 入歐姆定律時(shí)域卷積積分表達(dá)式，將其轉(zhuǎn)換成線性積分形式；
[0013] 3)、將計(jì)算時(shí)間進(jìn)行密集剖分，利用梯形積分法，獲得歐姆定律時(shí)域離散遞推表達(dá) 式；
[0014] 4)、將歐姆定律時(shí)域離散遞推表達(dá)式代入Maxwell方程中，基于三維時(shí)間域有限差 分方法，在時(shí)間和空間上進(jìn)行中心差分離散，推導(dǎo)電場(chǎng)E(t)、磁場(chǎng)H(t)的迭代方程；
[0015] 5)、采用非均勻三維Yee氏網(wǎng)格對(duì)計(jì)算區(qū)域進(jìn)行剖分，設(shè)置電性參數(shù)，時(shí)間步長(zhǎng)，并 對(duì)各個(gè)網(wǎng)格進(jìn)行電導(dǎo)率參數(shù)賦值；
[0016] 6)、計(jì)算初始場(chǎng)，對(duì)三維模型進(jìn)行電場(chǎng)E(t)、磁場(chǎng)H(t)的迭代，加載狄利克雷邊界 條件，完成電磁場(chǎng)數(shù)值計(jì)算；
[0017] 7)、計(jì)算結(jié)束后，提取磁場(chǎng)或電場(chǎng)的各分量響應(yīng)進(jìn)行成圖。
[0018] 進(jìn)一步地，步驟2中，根據(jù)歐姆定律的頻域微分形式，先通過(guò)頻時(shí)變換得到歐姆定 律時(shí)域卷積積分形式;再將電導(dǎo)率時(shí)域表達(dá)式代入，獲得歐姆定律時(shí)域卷積積分表達(dá)式，如 式（1)所示;根據(jù)電導(dǎo)率〇(t)近似于的形式的特點(diǎn)，構(gòu)建電導(dǎo)率參數(shù)的e指數(shù)輔助方程r (t)，如式(2)所示，將r(t)的表達(dá)式代入歐姆定律的卷積積分表達(dá)式中，得到歐姆定律的線 性積分形式，如式(3)所示：
(1) ') (3)
[0022] 式中J(t)是電流密度，E(t)為電場(chǎng)值，t為時(shí)間；電導(dǎo)率絕對(duì)值表達(dá)式為
A為時(shí)間常數(shù)、n為極化率，為頻率趨近無(wú)窮大時(shí)的電導(dǎo)率值。
[0023] 進(jìn)一步地，步驟3中，將計(jì)算時(shí)間進(jìn)行密集剖分，利用梯形積分法，得到積分求和| (tn)，提取Ut n)最后一項(xiàng)構(gòu)建遞推公式，如式(4)所示，進(jìn)而實(shí)現(xiàn)降維計(jì)算;最后將式⑷代 入式(3)，并對(duì)時(shí)間進(jìn)行離散，最終得到n時(shí)刻歐姆定律時(shí)域離散遞推表達(dá)式，如式(5)所示：
[0026] 式中n表示第n個(gè)時(shí)刻，r (tn)、利:tn>、E (tn)、U tn)為第n個(gè)時(shí)刻所對(duì)應(yīng)的數(shù)值，A t (n)為第n個(gè)時(shí)間步長(zhǎng)；
[0027] 接著在
時(shí)刻對(duì)式(4)進(jìn)行離散，再代入(5)式中，進(jìn)行整理，得到
時(shí)刻歐 姆定律時(shí)域離散遞推表達(dá)式，如式所示：
[0030] 進(jìn)一步地，步驟4中，利用步驟3中的推導(dǎo)結(jié)果，將式(6)中的
I項(xiàng)，采用平均 值近似法進(jìn)行計(jì)算，有
，代入無(wú)源Maxwell方程，采用時(shí)域有限差分 方法，構(gòu)建出電場(chǎng)E(tn+1)與E(tn)的迭代方程;對(duì)于
1中的Utn)項(xiàng)，先按照步驟3中式 ⑷的遞推公式進(jìn)行計(jì)算，獲得n時(shí)刻的|(tn)值，再將其代入電場(chǎng)迭代公式，最終獲得n+1時(shí) 刻電場(chǎng)的迭代方程，如式(8)所示：
[0033] (i，j，k)為Yee式網(wǎng)格的坐標(biāo)點(diǎn)，￡f+i)(m)為n+1時(shí)亥丨jx軸方向電場(chǎng)分量，E x(n)(m)為n 時(shí)刻x軸方向電場(chǎng)分量
：時(shí)刻對(duì)應(yīng)坐標(biāo)點(diǎn)上的y軸方向的磁場(chǎng)分 量
?時(shí)刻對(duì)應(yīng)坐標(biāo)點(diǎn)上的Z軸方向的磁場(chǎng)分量，A yj、A Zk為Yee氏 網(wǎng)格的y、z方向的步長(zhǎng)。
[0034]進(jìn)一步地，步驟5中，根據(jù)空間、時(shí)間的先后順序進(jìn)行離散，先對(duì)三維Yee網(wǎng)格的& (x，y，z)賦值，將每個(gè)網(wǎng)格〇~(x，y，z)代入到電導(dǎo)率時(shí)域表達(dá)式中，得到〇(x，y，z，t);再結(jié)合 電流密度7 (^+^)的計(jì)算時(shí)刻特點(diǎn)，設(shè)定時(shí)間步長(zhǎng)A t(n)，將時(shí)間密集剖分成N+1個(gè)時(shí)刻，計(jì) 算時(shí)間半步長(zhǎng)
，再對(duì)時(shí)間進(jìn)行肘刻剖分，對(duì)應(yīng)為t
共2N+1個(gè)
時(shí)刻，在2N+1個(gè)時(shí)刻分別對(duì)r(t)、6(t)進(jìn)行賦值后，存儲(chǔ)r(t)的2N+1個(gè)值。對(duì)于$(t)的數(shù)據(jù)， 先舍棄n時(shí)刻的值，再保留
時(shí)刻的值，這樣只需存儲(chǔ)N個(gè)數(shù)據(jù)。
[0035]進(jìn)一步地，步驟6中，在三維模型的基礎(chǔ)上，建立地面中心回線發(fā)射與接收線圈的 系統(tǒng)，計(jì)算出初始場(chǎng)后，代入到控制方程中，計(jì)算x、y、z方向的電場(chǎng)與x、y方向的磁場(chǎng)，為了 確保計(jì)算結(jié)果的唯一性，采用散度方程17必=0,計(jì)算z方向的磁場(chǎng)，應(yīng)用狄利克雷邊界條件， 完成三維模型的感應(yīng)-極化雙場(chǎng)電磁響應(yīng)數(shù)值計(jì)算。
[0036] 本發(fā)明與現(xiàn)有技術(shù)相比，有益效果在于，將時(shí)域電磁探測(cè)中的感應(yīng)和極化雙場(chǎng)同 時(shí)進(jìn)行理論計(jì)算，通過(guò)構(gòu)建輔助方程，推導(dǎo)歐姆定律的時(shí)域離散遞推表達(dá)式，對(duì)計(jì)算場(chǎng)量進(jìn) 行降維處理，避免了卷積和的求取，大大減少了內(nèi)存的使用，可以更高效、準(zhǔn)確地模擬三維 模型的感應(yīng)-極化雙場(chǎng)電磁響應(yīng)。
【附圖說(shuō)明】
[0037] 圖1是三維有限差分?jǐn)?shù)值算法的整體示意圖；
[0038] 圖2是歐姆定律時(shí)域卷積形式的離散算法示意圖；
[0039]圖3是三維異常體模型的示意圖。
[0040]圖4是迪拜模型均勻半空間的磁場(chǎng)響應(yīng)圖；
[0041 ]圖5是迪拜模型層狀大地的磁場(chǎng)響應(yīng)圖；
[0042]圖6是迪拜模型三維異常體的磁場(chǎng)響應(yīng)圖；
[0043] 圖7是歐姆定律時(shí)域卷積運(yùn)算降維前后的計(jì)算時(shí)間對(duì)比圖。
[0044] 圖8是歐姆定律時(shí)域卷積運(yùn)算降維前后的占用內(nèi)存對(duì)比圖。
【具體實(shí)施方式】
[0045] 為了使本發(fā)明的目的、技術(shù)方案及優(yōu)點(diǎn)更加清楚明白，以下結(jié)合實(shí)施例，對(duì)本發(fā)明 進(jìn)行進(jìn)一步詳細(xì)說(shuō)明。以均勻半空間模型、層狀的大地模型、三維異常體模型為例，進(jìn)行時(shí) 域有限差分的三維感應(yīng)-極化雙場(chǎng)電磁響應(yīng)計(jì)算。應(yīng)當(dāng)理解，此處所描述的具體實(shí)施例僅僅 用以解釋本發(fā)明，并不用于限定本發(fā)明。參見(jiàn)圖1結(jié)合圖2所示，一種時(shí)域有限差分的三維感 應(yīng)-極化雙場(chǎng)數(shù)值模擬方法，包括：
[0046] 1、基于迪拜模型(頻率相關(guān)系數(shù)c = l時(shí)的Cole-Cole復(fù)電阻率模型），先獲得迪拜 模型電導(dǎo)率的頻域形式，再通過(guò)逆拉普拉斯變換，得到時(shí)域電導(dǎo)率表達(dá)式 〇(t);迪拜模型的 頻域表達(dá)式如式(1)所示
(1)
[0048]式（1)中極化率為n，時(shí)間常數(shù)為T(mén)1;經(jīng)拉普拉斯逆變換后得到迪拜模型時(shí)域電導(dǎo) 率表達(dá)式為
[0050] 式(2)又可表示為 (2)
[0051 ] = (7^5(1) - a(t) ? u(t) (3)
[0052] 式(3)中
J(t)、u(t)分別是單位脈沖函數(shù)和單位階躍函數(shù)。
[0053] 2、將電導(dǎo)率的時(shí)域表達(dá)式〇(t)代入歐姆定律中，獲得歐姆定律時(shí)域卷積積分表達(dá) 式，結(jié)合電導(dǎo)率〇(t)近似于aeTet形式的特點(diǎn)，構(gòu)建電導(dǎo)率參數(shù)的e指數(shù)輔助方程r(t)，代入 歐姆定律時(shí)域卷積積分表達(dá)式，將其轉(zhuǎn)換成線性積分形式；
[0054] 首先，將歐姆定律的頻域微分形式進(jìn)行頻時(shí)變換后得到時(shí)域卷積形式，表達(dá)式為： J(t) = f_\o(t-r)E(r)dT (4)
[0055] 將〇(t)代入到式(4)中得歐姆定律時(shí)域卷積積分表達(dá)式為：
[0056] /(〇 = a^-EU) -- /〇C a(t: - T)￡(i)dr (5)
[0057] 根據(jù)電導(dǎo)率〇(t)近似于的形式的特點(diǎn)，構(gòu)造輔助方
代入歐姆定 律時(shí)域卷積積分表達(dá)式將其轉(zhuǎn)換成線性積分形式，如式(6)所示：
[0058] ./(t) = ' . -孩(t) ?. /(|>(7：)后(T)dr (6)
[0059] 3、對(duì)計(jì)算時(shí)間進(jìn)行密集剖分，利用梯形積分法，獲得歐姆定律時(shí)域離散遞推表達(dá) 式：
[0060] 1)、將時(shí)間軸劃分成N+1個(gè)時(shí)刻，則有⑴、t⑵~t(n)這樣的時(shí)間點(diǎn)，At( n) = t At(n)是可變步長(zhǎng)，隨時(shí)間的增加逐漸加長(zhǎng)；
[0061] 采用梯形積分法
(?)
[0062] 將各個(gè)時(shí)間步長(zhǎng)的積分用梯形積分法近似后的表達(dá)式為
(B) (9)
[0065]對(duì)式(9)進(jìn)行拆分，提取最后一項(xiàng)，構(gòu)建可遞推公式，實(shí)現(xiàn)數(shù)據(jù)的降維計(jì)算，表達(dá)式 為
(10)
[0067] 2)、特殊時(shí)間點(diǎn)遞推公式的求解方法，當(dāng)t = ti時(shí)，〖(t1)表達(dá)式為
(11)
[0069] 由于電磁場(chǎng)計(jì)算中￡(#)=0,所以忽略掉Ht13) ?￡(#)的乘積項(xiàng)，之后仍按照梯形 積分法來(lái)近似各個(gè)時(shí)間步長(zhǎng)上的線性積分。
[0070] 3)、時(shí)域有限差分方法中以"
時(shí)刻為電場(chǎng)的觀察點(diǎn)，為了與時(shí)域有限差分方法 更好的結(jié)合，繼續(xù)推導(dǎo)n $時(shí)刻/ (Vn+9方程為 (12)
[0073] 4、將離散后的歐姆定律代入Maxwell方程中，利用三維時(shí)間域有限差分方法，在時(shí) 間和空間上進(jìn)行中心差分離散，推導(dǎo)得到電場(chǎng)E(t)、磁場(chǎng)H(t)的迭代方程：
[0074] 1)、將步驟3的
|中拆分成兩項(xiàng)，兩者形式相同，表達(dá)式為
_ (14) (卩）
[0077]對(duì)（14)式單獨(dú)進(jìn)行遞推計(jì)算，對(duì)（15)式采用平均值近似法把
替換成
項(xiàng)，從而構(gòu)建出E(tn+1)與E(t n)的迭代方程。
[0078] 2)、將離散歐姆定律的卷積遞推公式代入無(wú)源Maxwell方程，推導(dǎo)E、H迭代公式。無(wú) 源Maxwell旋度方程組表達(dá)式及其本構(gòu)關(guān)系表達(dá)式為
(16)
[0081 ]將旋度方程在x、y、z方向展開(kāi)表達(dá)式如下：
(17) (IB)
[0084] 進(jìn)一步地，采用三維時(shí)域有限差分方法，將Maxwell方程組在空間和時(shí)間進(jìn)行差分 離散，將微商轉(zhuǎn)換為差商。令f (x，y，z，t)代表直角坐標(biāo)系下的E、H的某一分量，在時(shí)間和空 間上離散有式(19)
[0085] f(x,y,z,t)=f(i Ax,j Ay,kAz,nAt)=f(n)(i,j,k) (19)
[0086] 對(duì)函數(shù)的時(shí)間和空間的一階偏導(dǎo)數(shù)取中心差分近似，推導(dǎo)迭代公式。
[0087] 例如:選取觀察點(diǎn)（x，y，z)為Ex的節(jié)點(diǎn)，即.
點(diǎn)在時(shí)間
?的時(shí)間節(jié)點(diǎn) 觀察，將式(12 )、（ 14 )、（ 15)代入后推導(dǎo)表達(dá)式為
[0090] 同理可以推導(dǎo)出其他方向的分量。
[0091] 同樣，選取觀察點(diǎn)(x，y，z)為Hx的節(jié)點(diǎn)，即
）點(diǎn)在n時(shí)刻進(jìn)行觀察，獲得
表達(dá)式為
同理，可以繼續(xù)推導(dǎo)磁場(chǎng)其他方向的分量。
[0094] 5、采用非均勻三維Yee氏網(wǎng)格對(duì)計(jì)算區(qū)域進(jìn)行剖分，設(shè)置電性參數(shù)，時(shí)間步長(zhǎng)，并 對(duì)各個(gè)網(wǎng)格進(jìn)行電導(dǎo)率參數(shù)賦值：
[0095] 1)、劃分非均勻Yee式網(wǎng)格，根據(jù)時(shí)域電磁響應(yīng)在近源處幅值大且變化快，離源較 遠(yuǎn)時(shí)幅值小且變化慢的特點(diǎn)，滿足網(wǎng)格步長(zhǎng)近源處密集，離源遠(yuǎn)處稀疏的條件，設(shè)置網(wǎng)格數(shù) 目為101 X 101 X50,其中x、y方向上的網(wǎng)格數(shù)目均為101個(gè)，z方向上網(wǎng)格數(shù)50個(gè)，三個(gè)方向 上的最小和最大網(wǎng)格步長(zhǎng)均為l〇m和120m;在計(jì)算區(qū)域內(nèi)磁導(dǎo)率均設(shè)置為真空磁導(dǎo)率，均勻 半空間模型 0~設(shè)置為0.01 s/m;層狀大地模型上層50m，上層設(shè)置〇~為0.01 S/m，下層2300m， 設(shè)置為〇. lS/m;三維異常體模型，如圖3所示，背景場(chǎng)尺寸2860 X 2860 X 2350(m3)，設(shè)置 為0.01S/m，異常體尺寸1600 X 1600 X 1700(m3)，〇~設(shè)置為0 ? lS/m。
[0096] 2)、初始時(shí)刻t(Q)按照公式進(jìn)行取值，其中，^為最上層的磁導(dǎo) 率，這里取真空磁導(dǎo)率，01為最頂層的電導(dǎo)率〇~，A 1為Yee氏網(wǎng)格中最小的空間步長(zhǎng);在時(shí)間 域有限差分方法中，一般采用
(22)
[0098] 式(22)中，a的取值范圍是0.1~0.2。
[0099]迭代前就對(duì)A t(n)進(jìn)行賦值并完成迭代所需的所有中間參數(shù)的計(jì)算，將時(shí)間t預(yù)先 分成幾個(gè)較大的時(shí)間段，近似地依照式(22)，為每一段的A t賦值。這樣的處理方式加快了 迭代速度。算例中，初始時(shí)刻為1微秒，總時(shí)長(zhǎng)為10毫秒，一共分成了 15段，前14段每段有100 個(gè)時(shí)間步長(zhǎng)，最后一段有500個(gè)時(shí)間步長(zhǎng)，其中，時(shí)間步長(zhǎng)At最小為0.05微秒，最大為10微 秒。
[0100] 3)、根據(jù)空間、時(shí)間的先后順序進(jìn)行離散，先對(duì)三維Yee網(wǎng)格的〇~(X， y，Z)賦值，將 每個(gè)網(wǎng)格〇c"(x，y，z)代入到電導(dǎo)率時(shí)域表達(dá)式中，得到〇(x，y，z，t);再結(jié)合電流密度
丨的計(jì)算時(shí)刻特點(diǎn)，設(shè)定時(shí)間步長(zhǎng)A t(n)，將時(shí)間密集剖分成N+1個(gè)時(shí)刻，計(jì)算時(shí)間半 步長(zhǎng)^，再對(duì)時(shí)間進(jìn)行$時(shí)刻剖分，對(duì)應(yīng)為-
，共2N+1個(gè)時(shí)刻，在 2N+1個(gè)時(shí)刻分別對(duì)r (t)、#(:￡)進(jìn)行賦值后，存儲(chǔ)r (t)的2N+1個(gè)值。對(duì)于4⑴的數(shù)據(jù)，先舍棄n 時(shí)刻的值，再保留
時(shí)刻的值，這樣只需存儲(chǔ)N個(gè)#(t)數(shù)據(jù)。
[0101] 6、計(jì)算初始場(chǎng)，對(duì)三維模型進(jìn)行電場(chǎng)E(t)、磁場(chǎng)H(t)的迭代，加載狄利克雷邊界 條件，完成電磁場(chǎng)數(shù)值計(jì)算：
[0102] 在三維模型的基礎(chǔ)上，建立地面中心回線發(fā)射與接收線圈的系統(tǒng)，計(jì)算出初始場(chǎng) 后，代入到控制方程中，計(jì)算x、y、z方向的電場(chǎng)與x、y方向的磁場(chǎng)，為了確保計(jì)算結(jié)果的唯一 性，采用散度方程17 * 6 = 0,計(jì)算z方向的磁場(chǎng)，應(yīng)用狄利克雷邊界條件，完成三維感應(yīng)-極 化雙場(chǎng)電磁響應(yīng)計(jì)算。最后提取z分量的磁場(chǎng)響應(yīng)。
[0103] 7、計(jì)算結(jié)束后，提取磁場(chǎng)或電場(chǎng)的各分量響應(yīng)進(jìn)行成圖。
[0104] 計(jì)算成圖結(jié)果詳見(jiàn)圖4-圖6。
[0105] 以上所述僅為本發(fā)明的較佳實(shí)施案例而已，并不用以限制本發(fā)明，凡在本發(fā)明的 精神和原則之內(nèi)所作的任何修改、等同替換和改進(jìn)等，均應(yīng)包含在本發(fā)明的保護(hù)范圍之內(nèi)。
【主權(quán)項(xiàng)】
1. 一種時(shí)域有限差分的三維感應(yīng)-極化雙場(chǎng)數(shù)值模擬方法其特征在于，包括如下步驟： 1) 、基于迪拜模型(頻率相關(guān)系數(shù)C = I時(shí)的Cole-Cole復(fù)電阻率模型），先獲得迪拜模型 電導(dǎo)率的頻域形式，再通過(guò)逆拉普拉斯變換，得到電導(dǎo)率的時(shí)域表達(dá)式 σ (t); 2) 、將電導(dǎo)率的時(shí)域表達(dá)式〇(1)代入歐姆定律中，獲得歐姆定律時(shí)域卷積積分表達(dá)式， 結(jié)合電導(dǎo)率〇(t)近似于ae+形式的特點(diǎn)，構(gòu)建電導(dǎo)率參數(shù)的e指數(shù)輔助方程r(t)，代入歐姆 定律時(shí)域卷積積分表達(dá)式，將其轉(zhuǎn)換成線性積分形式； 3 )、將計(jì)算時(shí)間進(jìn)行密集剖分，利用梯形積分法，獲得歐姆定律時(shí)域離散遞推表達(dá)式； 4) 、將歐姆定律時(shí)域離散遞推表達(dá)式代入Maxwell方程中，基于三維時(shí)間域有限差分方 法，在時(shí)間和空間上進(jìn)行中心差分離散，推導(dǎo)電場(chǎng)E(t)、磁場(chǎng)H(t)的迭代方程； 5) 、采用非均勻三維Yee氏網(wǎng)格對(duì)計(jì)算區(qū)域進(jìn)行剖分，設(shè)置電性參數(shù)，時(shí)間步長(zhǎng)，并對(duì)各 個(gè)網(wǎng)格進(jìn)行電導(dǎo)率參數(shù)賦值； 6) 、計(jì)算初始場(chǎng)，對(duì)三維模型進(jìn)行電場(chǎng)E(t)、磁場(chǎng)H(t)的迭代，加載狄利克雷邊界條件， 完成電磁場(chǎng)數(shù)值計(jì)算； 7) 、計(jì)算結(jié)束后，提取磁場(chǎng)或電場(chǎng)的各分量響應(yīng)進(jìn)行成圖。2. 按照權(quán)利要求1所述的一種時(shí)域有限差分的三維感應(yīng)-極化雙場(chǎng)數(shù)值模擬方法，其特 征在于：步驟2中，根據(jù)歐姆定律的頻域微分形式，先通過(guò)頻時(shí)變換得到歐姆定律時(shí)域卷積積分 形式;再將電導(dǎo)率時(shí)域表達(dá)式代入，獲得歐姆定律時(shí)域卷積積分表達(dá)式，如式（1)所示;根據(jù) 電導(dǎo)率〇(t)近似于CteT iit的形式的特點(diǎn)，構(gòu)建電導(dǎo)率參數(shù)的e指數(shù)輔助方程r(t)，如式(2)所 示，將r(t)的表達(dá)式代入歐姆定律的卷積積分表達(dá)式中，得到歐姆定律的線性積分形式，如 式(3)戶，> (1) (2) (3) 式中J(t)是電流密度，E(t)為電場(chǎng)值，t為時(shí)間；電導(dǎo)率絕對(duì)值表達(dá)式為噸)= 為時(shí)間常數(shù)、η為極化率，為頻率趨近無(wú)窮大時(shí)的電導(dǎo)率值。3. 按照權(quán)利要求1所述的一種時(shí)域有限差分的三維感應(yīng)-極化雙場(chǎng)數(shù)值模擬方法，其特 征在于： 步驟3中，將計(jì)算時(shí)間進(jìn)行密集剖分，利用梯形積分法，得到積分求和ξ(〇,提取ξ(〇 最后一項(xiàng)構(gòu)建遞推公式，如式(4)所示，進(jìn)而實(shí)現(xiàn)降維計(jì)算;最后將式(4)代入式(3)，并對(duì)時(shí) 間進(jìn)行離散，最終得到η時(shí)刻的歐姆定律時(shí)域離散遞推表達(dá)式，如式(5)所示：式中11表示第11個(gè)時(shí)刻4(〇、利^)4(〇、|(〇為第11個(gè)時(shí)刻所對(duì)應(yīng)的數(shù)值，/^(11)為 第η個(gè)時(shí)間步長(zhǎng)； 接著在η + ^時(shí)刻對(duì)式⑷進(jìn)行離散，再代入(5)式中，進(jìn)行整理，得到η + #時(shí)刻的歐姆 Z Z. 定律時(shí)域離散遞推表達(dá)式，如式所示： (6)(7)。4. 按照權(quán)利要求1所述的一種時(shí)域有限差分的三維感應(yīng)-極化雙場(chǎng)數(shù)值模擬方法，其特 征在于： 步驟4中，利用步驟3中的推導(dǎo)結(jié)果，將式(6)中的|:項(xiàng)，采用平均值近似法進(jìn)行計(jì) 算，笮?戈入無(wú)源Maxwell方程，采用三維時(shí)域有限差分方法，構(gòu)建出 電場(chǎng)E(tn+1)與E(tn)的迭代方程;對(duì)于式(6)中的|(tn)項(xiàng)，先按照步驟3中式(4)的遞推公式 進(jìn)行計(jì)算，獲得η時(shí)刻的|(tn)值，再將其代入電場(chǎng)迭代公式，最終獲得n+1時(shí)刻電場(chǎng)的迭代 方程，如式(8)所示：格的y、z方向的步長(zhǎng)。5. 按照權(quán)利要求1所述的一種時(shí)域有限差分的三維感應(yīng)-極化雙場(chǎng)數(shù)值模擬方法，其特 征在于： 步驟5中，根據(jù)空間、時(shí)間的先后順序進(jìn)行離散，先對(duì)三維Yee網(wǎng)格的〇~(X，y，Z)賦值，將 每個(gè)網(wǎng)格〇~(X，y，Z)代入到電導(dǎo)率時(shí)域表達(dá)式中，得到0(x，y ,Z，t);再結(jié)合電流密度的計(jì)算時(shí)刻特點(diǎn)，設(shè)定時(shí)間步長(zhǎng)A t(n)，將時(shí)間密集剖分成N+1個(gè)時(shí)刻，計(jì)算時(shí)間半步七 再對(duì)時(shí)間進(jìn)行：;時(shí)刻剖分，對(duì)應(yīng)： 共2N+1個(gè)時(shí)刻，在 2 2N+1個(gè)時(shí)刻分別對(duì)r (t)、外〇進(jìn)行賦值后，存儲(chǔ)r (t)的2N+1個(gè)值，對(duì)于(6((6)的數(shù)據(jù)，先舍棄η 時(shí)刻的值，再保留Tl + ^時(shí)刻的值，這樣只需存儲(chǔ)N個(gè)#(?)的數(shù)據(jù)。 2
【文檔編號(hào)】G06F17/50GK105893678SQ201610202882
【公開(kāi)日】2016年8月24日
【申請(qǐng)日】2016年4月1日
【發(fā)明人】嵇艷鞠, 吳燕琪, 關(guān)珊珊, 黃廷哲, 吳瓊, 王遠(yuǎn)
【申請(qǐng)人】吉林大學(xué)