[發明專利]一種圓錐型場源瞬變電磁反演方法

有效

• 申請號: 202111157516.3
• 申請日: 2021-09-30
• 公開（公告）號: CN113933905B
• 公開（公告）日: 2023-09-19
• 發明（設計）人: 楊海燕;劉志新;汪凌;儲丹華;顧艷艷;李哲;楊夫杰
• 申請（專利權）人: 中國礦業大學
• 主分類號: G01V3/28
• 分類號: G01V3/28;G01V3/38;G06F30/20;G06F17/16
• 代理公司: 南京經緯專利商標代理有限公司 32200
• 代理人: 姜慧勤

1.一種圓錐型場源瞬變電磁反演方法，其特征在于，該方法包括如下步驟：

步驟1，讀入初始模型m₁即初始電阻率，根據初始模型m₁構造粗糙度矩陣，設置初始拉格朗日乘子

步驟2，對于第k次迭代，計算當前模型m_k對應的正演數據F[m_k]、雅克比矩陣J_k以及數據擬合差X²；

具體過程如下：

對于第k次迭代，向圓錐型場源通入電流，對于圓錐型場源的第i匝圓形發射線圈，其中心處垂直方向磁場的頻率域響應為：

式中，H(ω)表示頻率域響應，r_i為第i匝圓形發射線圈的半徑，I(ω)為發射電流的傅式譜，λ為積分變量，h_i為第i匝圓形發射線圈離地高度，u₀為真空中的磁導率，Z₀為表層介質的輸入波阻抗，Z⁽¹⁾為第一層的輸入阻抗，J₁為第一類零階貝塞爾函數；k₁為第一層介質的波數，σ₁為第一層介質的電導率，ω為角頻率，μ₀為真空中磁導率；

對頻率域響應公式進行140點漢克爾變換得到：

式中，H_z,i(ω)表示140點漢克爾變換得到的頻率域，m為抽樣點的位置，b為回線半徑，a＝-7.91001919000e+00，s＝8.79671439570e-02，為漢克爾變換濾波系數；

設定激勵電流I(t)為單位階躍電流：

將140點漢克爾變換得到的頻率域，通過傅里葉變換求得時間域中響應：

式中，t為時間，f(t)為時間域電磁場分布，F(ω)為頻率域的電磁場分布，j為虛數單位；

將頻率域響應公式代入上式即得時間域的磁場響應：

式中，Im表示虛部，Re表示實部，為角頻率，H_z(t)為時間域電磁場，I₀為通入回線的電流，B_z(t)為時間域磁感應強度；

根據法拉第電磁感應定律，得到二次場的感應電位為v_i，將圓錐型場源產生的瞬變電磁感應電位視為n個線圈感應電位的疊加，則圓錐型場源下層狀模型的瞬變電磁感應電位V(t)為：

式中，n為圓錐型場源纏繞的圓形發射線圈匝數，Δt_i是第i匝線圈與第1匝線圈距離產生的延時，對于第1匝線圈，Δt₁＝0；采用全區視電阻率算法將感應電位V(t)轉換成反映電性變化的視電阻率，即當前模型對應的正演數據；

雅克比矩陣J_k為：

其中，表示求F[m_k]的梯度，F[m_k]表示當前模型m_k對應的正演數據；

數據擬合差X²為：

X²＝‖Wd-WF[m_k]‖²

其中，W＝diag{1/δ₁,1/δ₂,…,1/δ_M}，δ_c為第c個數據的標準差，c＝1,…,M，M為場源裝置測得的數據個數，d為場源裝置測得的數據集合；

步驟3，搜索使數據擬合差X²最小時的模型m_k+1，同時計算模型m_k+1對應的粗糙度矩陣R_k+1以及步長a_k；

所述模型m_k+1定義為：

其中，為第k+1次迭代的拉格朗日乘子，θ為N×N的矩陣，N為地層模型分段得到的總層數，J_k為雅克比矩陣，d_k＝d-F[m_k]+J_km_k，F[m_k]表示當前模型m_k對應的正演數據；

所述步長a_k的計算公式為：

其中，表示拉格朗日乘子函數的模型集，m_k、m_k+1分別第k、k+1次迭代的模型，為第k次迭代的拉格朗日乘子；

步驟4，判斷步驟3得到的最小數據擬合差是否小于反演所要求達到的擬合差若是，則進入步驟5，否則，判斷步驟3得到的最小數據擬合差是否小于第k-1次迭代得到的最小數據擬合差；若是，則進入步驟6，否則使步長a_k減半并返回步驟3；

步驟5，確定當數據擬合差X²等于反演所要求達到的擬合差時的拉格朗日乘子及對應的模型m′_k+1，根據模型m′_k+1計算粗糙度矩陣R′_k+1，判斷R′_k+1是否大于R_k+1，若是，則使步長a_k減半并返回步驟3，否則進入步驟6；

步驟6，保存步長a_k和模型m_k+1；

步驟7，判斷是否達到最大迭代次數或是否達到反演所要求達到的擬合差若達到最大迭代次數或達到反演所要求達到的擬合差則迭代結束，輸出反演結果，否則，令k＝k+1并返回步驟2。

2.根據權利要求1所述圓錐型場源瞬變電磁反演方法，其特征在于，步驟1所述根據初始模型m₁構造粗糙度矩陣為：

R₁＝‖θm₁‖²

其中，R₁表示根據初始模型m₁構造的粗糙度矩陣，θ為N×N的矩陣，定義如下：

N為地層模型分段得到的總層數。