积分方程在大地电磁三维反演的应用.docVIP

积分方程在大地电磁三维反演的应用 目前, 大地电磁二维反演已趋于成熟, 达到实用化程度, 而大地电磁三维反演成为研究的前沿课题, 并取得了长足的发展。三维反演方法主要有共轭梯度法、拟线性近似反演、快速松弛反演、贝叶斯统计反演和人工神经网络反演等。作者在本文中采用积分方程法计算均匀半空间中三维异常体的大地电磁响应, 利用Eaton等提出的方法, 将灵敏度矩阵分为线性项和体现反演过程非线性的偏微分项来求解, 提高了计算精度; 然后利用阻尼最小二乘技术进行最优化反演, 得到地下三维异常体电阻率。由于数据点只分布在电阻率异常区域而不是整个三维空间, 因而与有限元法和有限差分法相比其计算量相对较小, 效果较好。 1 积分方程法基本理论 根据Maxwell方程、积分方程理论, 以及相应的电磁张量格林函数, 可以得到均匀大地中三维异常体电磁响应的积分方程。 其中E（r）和H（r）是r 处的总电场和总磁场；EP ( r) 和HP ( r ) 是以天然平面电磁波为场源, 在大地中产生的一次电场和磁场； 是三维异常体与大地电导率差值; 和为电和磁并矢格林函数；为异常体内电场。 为进行数值计算, 将大地中三维异常体剖分成N个小立方体单元, 并假设在每个剖分单元内部电阻率是均匀分布的, 则各单元内电场可用如下有限求和近似： 整理上式最终可得如下矩阵方程: 解方程( 4)求得异常体内各单元总电场后, 再利用式( 1)和式( 2)即可得到空间内任意一点电磁场值。在求得空间总的电磁场后, 可由式( 5)得到视电阻率 其中 i = x, y; j = x, y 。 2 积分方程非线性反演方法 在反演计算中, 我们采用阻尼最小二乘最优化技术, 利用正演模拟得到的电磁场数据拟合观测到的电磁数据, 并逐步修改地电模型参数初值, 最终达到最优拟合, 反演得到地下介质电阻率异常分布情况。 阻尼最小二乘法可表述为如下求解模型修改量的线性方程组: 其中为电导率模型修改量; P 为灵敏度矩阵的元素; 右端矢量S表示实测场值和理论正演场值之间的相对偏差; K为阻尼因子; D 为平滑矩阵。利用式( 6)我们可以求出异常体电导率模型参数 , 并以作为新的模型参数初值, 重新计算目标函数。这样反复迭代, 以达到目标函数小于预先给定的一个小量, 此时的X 即为所求的反演结果。 3 灵敏度矩阵的求取 灵敏度矩阵的求取是反演过程中非常重要的一步, 根据式( 5)可知: 其中x表示地下剖分单元的电导率值;下标k表示第k个单元; i表示第i个观测点。 可以看出, 若要求得地表视电阻率或阻抗对地下模型参数的灵敏度矩阵, 只要求出地表电场和磁场对地下模型参数的灵敏度矩阵。利用Eaton等提出的方法, 将灵敏度矩阵分为线性项和体现反演过程非线性的偏微分项。首先计算地下异常体各剖分单元电导率扰动, 对于各单元中心点处电磁场的影响; 然后结合相应的张量格林函数, 得到各单元电导率扰动对于观测点处电磁场的影响。 根据式( 1), 且由于背景电导率是已知的, 故。则： 其中为异常体第n个单元, 对地表观测点i处的电张量格林函数。进一步求导可得 在计算式( 10)右端后一项时, 首先要求解异常体各单元中心点处电场对各单元电导率的偏导数, 然后利用相应的张量格林函数叠加, 求得地表总场对各单元电导率的偏导数。式( 10)中Ei 也是有x、y、z 三个分量, 以x 方向分量为例: 是异常体第n个单元电场对第k 块单元电导率的偏导数。的求解可仿照正演求解中对式( 3)的求解 方法, 离散后的方程可表示为: 这样, 通过求解n 次n 维线性方程组, 即可得到 ( n = 1, 2, N; k = 1, 2, N ) ,然后利用式( 10)、式( 11)可求得地表观测点电场值对地下异常体各单元电导率的偏导数。同理, 可得地表观测点处, 磁场值对地下异常体各单元电导率的偏导数, 然后代入式( 7) , 可最终求得地表观测的视电阻率对地下异常体各单元电导率的偏导数。 4 模型反演结果 设计如图1 所示模型, 围岩电阻率为200,异常体模型大小为600m600m300m，埋深500 m, 被剖分为66 3个小单元,各小单元电阻率分布均匀。 图2为地表观测点分布示意图, 每个测点之间横向和纵向间距都为200m, 工作频率采用五个频点, 分别为0.01 H z、0.1 Hz、1 H z、10 H z、100H z。以地表观测的数据进行反演, 反演时初始模型电阻率都为200, 初始阻尼因子为10,阻尼因子的缩放系数为2。经过十二次迭代, 反演拟合误差小于10-6, 反演结束。 图1 异常体模型剖分示意图 图2 地表观测点分布示意图 图3、图4、图5 分别为异常体各层的反演结果。异常体剖分为三层: 图3为剖分异常体第一层各个小单元