电测深曲线数字解释原理及自动反演算法

电测深曲线数字解释原理及自动反演算法 1. 拟合核函数法数字解释原理和方法 2. 视电阻率转换函数滤波计算及滤波系数 3. 电测深曲线自动反演解释的最优化算法 4. 拟合核函数自动反演方法及程序设计 5. 直接拟合视电阻率法反演原理及程序 6. 7. 电测深曲线的其他反演算法 17:37:16 1 法 （一）由视电阻率计算核函数的褶积表达式 理论对称四极电测深视电阻率表达式 上式可变为 利用汉克尔变换，可得 17:37:16 2 法 （一）由视电阻率计算核函数的褶积表达式 引入对数变量 代入上两式，可得 继续变换 17:37:16 3 法 （一）由视电阻率计算核函数的褶积表达式 则可简写为 )( )( )(1 反演 )(1 )( )( 正演 17:37:16 4 法 （二）拟合核函数法数字解释原理和方法 首先由实测 电测深曲线 计算出视电阻率转换 函数曲线 ，然后与根据层参数初值计算出理论核函数 曲线进行 对比，如果二者的拟合差大于给定拟合精度，则修改初始模型的层参数，重新计算理论核函数曲线，并重新进行比较。这一过程反复进行，直到拟合差小于给定精度为止。 基本思想 17:37:16 5 法 （二）拟合核函数法数字解释原理和方法 基本步骤 确定初始模型：（层参数初值 0P ） ； 由实测 ()s r 曲线 数 字 滤 法 视电阻率转换函数 T 曲线 由初值 0P 计算 T 理论曲线 按照一定法则进行曲线对比，判断 T 与 T 拟合程度是否达到要求。根据拟合差确定是否进行层参数调查 调整 0P ，以改善 T 与 T 的拟合程度 重复 步骤，直到 T 与 T 曲线达到拟合要求为止 17:37:16 6 法 （二）拟合核函数法数字解释原理和方法 迭代方法 （ 1）自动迭代方法 单纯形法 最小二乘法 阻尼最小二乘法 梯度法 牛顿法 变尺度法 （ 2）非自动迭代方法 17:37:16 7 （一）视电阻率转换函数滤波计算的数学模型 上节提到的褶积积分表达式 根据采样定理，视电阻率离散 代入化简 17:37:16 8 （一）视电阻率转换函数滤波计算的数学模型 令采样间隔 代换 或 为反演滤波系数 17:37:16 9 （二）反演滤波系数 反演滤波系数求取与正演相同 实测视电阻率曲线外延方法 17:37:16 10 法 （一）最优化算法的数学实质 寻找模型 电测深曲线定量解释的关键 建立法则 吻合 评定拟合程度的方法 如何自动修改模型参数 评定拟合程度的方法 采样值偏差的平方和 17:37:16 11 法 （一）最优化算法的数学实质 即 设离散实测数据 ),2,1( 实际曲线 建立数学模型 ),( 即理论曲线 评定理论与实测曲线的拟合程度 目标函数 改变 X 的值，使得 m i n)( X 如果给定模型参变量初值为 )0(X ， 并设修正量为 X , 则修正后的参变量为 )1(X 参变量修正使目标函数 m i n)( )1( X 17:37:16 12 法 （二）高斯法与广义最小二乘法 对于 理论模型函数为： 是层参数 向量 n=2元非线性函数 )( k 个采样点的理论值为 显然 17:37:16 13 法 （二）高斯法与广义最小二乘法 17:37:16 14 法 （二）高斯法与广义最小二乘法 具体做法 在 0x 点附近将 ),( 成泰勒级数， 略去 X 的高次项 ,2,1 离散点数 ,2,1 参数个数 代入目标式 17:37:16 15 法 （二）高斯法与广义最小二乘法 为了求出 )( X 的极小值，其各个参量 导数应为 0 化简后得到 17:37:16 16 法 （二）高斯法与广义最小二乘法 式中 写成矩阵形式 17:37:16 17 法 （二）高斯法与广义最小二乘法 实际计算步骤 设定 0x （初值），计算 ； 解 法 方程 得到 x ，有 判断 x 是否满足 )( X （给定精度），否 则继续下步； 把 x 当作 0x ，进行， 步，有 2 以上步骤反复进行，直到满足条件； 最终的 x 是 解释结果。 17:37:16 18 法 （二）高斯法与广义最小二乘法 由于泰勒展开式是近似的，所以要反复迭代修正，求得参数的近似最决于 的大小，逐次迭代所得的参数 逐步逼近真值， 也逐渐减小，当小到二次项可忽略不计时，便成精确等式。如果二次和二次以上的偏导数均为零，则不论 多大，式也总是精确成立，这就是线性模型的的情况。 17:37:16 19 法 （二）高斯法与广义最小二乘法 高斯法 要求系数矩阵是正定和非奇异，否则不能确定极小点，无法求解 迭代过程不稳定，当变量过大，泰勒展开式高次项不容忽略，不嫩呢保证平稳收敛 广义最小二乘法 根据实际应用经验，对高斯法提出了一个改进方案，这时不把 X 直接作为修正量，而是乘以一个常数 作为修正量，将问题变为求 的值 为极小值，这是单变量求极小值的问题，然后取： )()1( ，作为新的初值。通常 1,0 17:37:16 20 法 （三）最速下降法 最早提出来的一种最优化方法，其基本思想是： 在某点 附近以沿负梯度方向下降 最快 因此取修正步长： 新参数： 17:37:16 21 法 （三）最速下降法 如何计算 t ？ 将 的邻域内展成泰勒 级数，省略高次项 在 代入目标函数式 由于 的选择应使 趋于最小 故 17:37:16 22 法 （三）最速下降法 如何计算 t ？ 将 的邻域内展成泰勒 级数，省略高次项 在 代入目标函数式 由于 的选择应使 趋于最小 故 17:37:16 23 法 （三）最速下降法 最速下降法计算步骤 1 给定出始发点 0x 和精度 ， ； 2 计算方向 )( xg r a 3 判断 若成立，转步 4 求 t ，使 )( 极小 转 2 5输出极小值点 17:37:16 24 法 （三）最速下降法 一 维寻查法求 t 二分法 : 当 ,函数递减 0)( 当 , 函数递增 0)( 设 a,b 上有 0)(,0)( 则 a, b 内必有极小值 x 取20 0)(0 小点在 0x 右边，令 101 , 0)( 0 小点在 0x 左边，令 011 , 得到新的极小点区域 a 1 ,b 1 , )(2111 17:37:16 25 法 （三）最速下降法 一 维寻查法求 t 如上逐次分半 , 得到一系列极小点在域 , 11 kk 当 k ， 2 1 17:37:16 26 法 （四）阻尼最小二乘法 最优化问题的求解可归结为： 1. 确定修正量的方向（即搜索方向），使目标函数 收敛到极小点 的 下降方向 2. 确定修正向量的大小，即迭代步长 3. 选取修正后的参数作为初值，重新进行迭代计算 实际应用中，最小二乘法与负梯度法各有特点 17:37:16 27 法 （四）阻尼最小二乘法 高斯法： 迭代过程常常不稳定 。 当 x 过大， 如果初值选择和合适，则能很快的收敛，如果初值选择不合适，则易于发散， 使01 发散。 高斯法修正向量 x 的 方向 称为最小二乘方向 （记为t） ，对于二次目标函数 )( X 来说，直接 指向 的极小值点。 最速下降法： 最速度下降方向 垂直于过点 )( 等值面，是目标函数 在 )( 处的 梯度方向 的反方向（记为g） 。 能够 保证收敛01 ，如果 修正步长 太小，收敛速度过慢，求解十分费时。 17:37:16 28 法 （四）阻尼最小二乘法 马奎特提出阻尼最小二乘法 在迭代的每一步，最好尽量使用最小二乘法作为搜索方向，以便使修正向量的步长尽可能地增大，从而节省求解时间。但是当修正步长过大不能保证收敛时，则改用接近最速下降方向的搜索方向，缩小步长以保证收敛。 这种方法是在高斯法和最小二乘法之间取某种插值，它力图以最大步长前进，同时又能紧靠目标函数的负梯度方向，从而既保证收敛，又加快收敛 速度。 17:37:16 29 法 （四）阻尼最小二乘法 具体做法 法方程加入阻尼因子 单位矩阵 用来控制修正量方向与步长的任意正数，称为阻尼因子 （ 1 ） 0 时，阻尼最小二乘法退化为高斯法 （ 2 ） 时，阻尼最小二乘法修正微量的方向是最速下降方向，为最速下降法 （ 3 ） 当 由零开始逐渐增大时 ，搜索方向由最小二乘法方向转向最速下降方向； （ 4 ） 越大，阻尼最小二乘法修正向量的步长 x 较越小，可保证迭代稳定收敛， 速度较慢 （ 5 ） 越 小时，阻尼最小二乘法的搜索方向接近于最小二乘法方向，从而可取得足够大的步长 x ，收敛速度加快 17:37:16 30 （一）拟合核函数自动反演的数学模型 对数型误差形式的目标函数 视电阻率转换函数，第 k 个采样点的值 ,T k p 理论模型核函数 1 2 , 1 , 2 , 1, h h h 层函数，电性层数 N m 采样点数 17:37:16 31 （一）拟合核函数自动反演的数学模型 根据阻尼最小二乘法，给定初值 ，则 层参数修正量 由法方程决定 17:37:16 32 （一）拟合核函数自动反演的数学模型 实际计算中，为节省计算量，利用 中元素的对称性 定义 则 17:37:16 33 （一）拟合核函数自动反演的数学模型 若令 则有 法方程 17:37:16 34 （二）雅可比矩阵的形成 核函数递推公式 利用复合函数求导的连锁法则求 1当 jk 时， 不含 k ， 如 3T 中不含 2 ， 2h 如 :320T 17:37:16 35 （二）雅可比矩阵的形成 当 j=k 时，如：11T ， 11 ， 33T ， 33 17:37:16 36 （二）雅可比矩阵的形成 当 jk 时 , 如 :13T，13其中 : 221 111/j j 17:37:16 37 （二）雅可比矩阵的形成 具体计算时，利用递推公式，反向， ，1次计算 出 的值 1 ， 11 ， 12 2 ， 22 ， 23 （ 2T ， 11T 38 17:37:16 （二）雅可比矩阵的形成 17:37:16 39 （三）阻尼因子的选择 17:37:16 40 阻尼最小二乘法方程 增大： ， 一般为 2 10 减少： / 这种设置的缺点及改进 1 、 各电性层层参数s曲线形态变化幅度主要取决于/ 。 如： 1 10 ，2 1 ，3 100 1 1 , 2 , 1 0 1P 1 1 2 1 0 1 1 0 % 1 1 1 . 1 1 1 0 （三）阻尼因子的选择 17:37:16 41 改进：构造阻尼因子矩针 D ，主对角线上的元素为： 21 j j 1, 2 , , n 则法方程： B （三）阻尼因子的选择 17:37:16 42 2 、 根据线性化指标 R 来增大 ，或减少 式中： k k 1, ：以 k k 1p , p 为参数的目标函数值 k k 1 , ： T k , p k 作线性化处理后的目标函数值 （三）阻尼因子的选择 17:37:16 43 2 、 根据线性化指标 R 来增大 ，或减少 具体上：设常数 1 , 2： 0 1R 2R 1 2 时，目标函数线性化程度较好 , 迭代收敛 , 1 时，目标函数线性化程度较差 , 迭代不收敛 , 向最速下降方向靠拢 ,保证收敛 21 时， 线性化一般， 不变 （四） 拟合核函数法自动反演实用程序 17:37:16 44 （ 1 ） 给出实测 s 值，计算视电阻率转换函数 T 值 （ 2 ） 设定 初值，常数 1 2 1 2, , R , R , 其 它 指标 （ 3 ） 给出模型参数初值 0P ，计算理论 0T , , （ 4 ） 计算雅可比矩阵，形成法方程 B （ 5 ） 解法方程求 P ， 10P P P （ 6 ） 由 1P 计算目 标函数 1 ，计算 R (7) 比较 R 与 12R , R 大小： 1 , 发散 , 2 转； 2时 , 收敛 , 转 12R R R 时， 不变，令 01 ，转 (8) 令 01 ； 1 ， 转 步继续迭代，直到取得满意的结果为止 （四） 拟合核函数法自动反演 实用程序框图 17:37:16 45 5直接拟合视电阻率法反演原理及程序设计 （一）非自动迭代反演算法 17:37:16 46 特点 （ 1）最简单、最直接的方法 （ 2）层参数的修改人工完成 （ 3）根据工作区已知的地质、钻探、电测井和物性资料给出初始模型 （ 4）根据理论与实际曲线实时显示拟合情况，给出修正量 （ 5）适合于有经验的解释者 步骤 （ 1）给出实测采样值及采样点坐标 （ 2）建立地电模型，给出层参数初始值 （ 3）计算理论电测深曲线 （ 4）实测、理论曲线拟合对比 （ 4）修改、输出结果 5直接拟合视电阻率法反演原理及程序设计 （二）自动迭代反演算法 17:37:16 47 自动迭代反演的数学模型 对数型目标函数 首先对 ,s 线性化，即在初值 0P 的邻域内泰勒级数展开 根据多元函数求极值的条件 0 1 , 2 , , 可建立法方程： A P B 同样： J 5直接拟合视电阻率法反演原理及程序设计 （二）自动迭代反演算法 17:37:16 48 自动迭代反演的数学模型 5直接拟合视电阻率法反演原理及程序设计 （二）自动迭代反演算法 17:37:16 49 自动迭代反演的数学模型 阻尼最小二乘法 雅可比矩阵的形成 求 ,s 的偏导数，有二种方法 （ 1）数字滤波法 求导 5直接拟合视电阻率法反演原理及程序设计 （二）自动迭代反演算法 17:37:16 50 将上式离散化 计算步骤 计算采样点 1, ,的 1 的值 对 1 进行滤波 此方法特点 优点：精度高 缺点：计算量 大 多计算 的采样点数 5直接拟合视电阻率法反演原理及程序设计 （二）自动迭代反演算法 17:37:16 51 （ 2）差商法 s T 的 s 曲线首尾延长问题 以 21 y i b i 设 12: 2 , 3b i n n 2b 1b 0b 1b 2b 3b 1T 1s 0s 1s 2s 3s 4s 0T 2s 1s 0s 1s 2s 3s 1T 3s 2s 1s 0s 1s 2s 6T 8s 7s 6s 5s 4s 3s 5直接拟合视电阻率法反演原理及程序设计 （三）拟合核函数法与视电阻率直接拟合法比较 17:37:16 52 （ 1）拟合 核函数法 速度 快，迭代较稳定 拟合 精度低 须 对 实测视电阻率曲线延拓 （ 2）直接 拟合视电阻率法 计算 量大，迭代初值要求高（差商法） 拟合 精度高 不需对实测视电阻率曲线延拓 5 （一） 17:37:16 53 设 n 层层状介质模型 , 切出高 ，底面为 11 m 的柱体 总纵向电导 ： 1n in 总横向电阻 1nn i 5 （一） 17:37:16 54 则有 / ， / 各向异性系数 1n 根据：一个厚度为 H ，各向异性系数为 的电性层与一个各向同性的厚度为 H ，电阻率为m 的电性层在电性上是等效的 因此，上面给出的地电模型的等效各向同性电性层层参数为 m 拟真电阻率，又称 阻率 拟真深度，又称 度 5 （一） 17:37:16 55 当在 n 层增加 h 厚度时 即 m 与 在某种函数关系 5 （二） 计算 17:37:16 56 当 10 时 1 1 证明： 10, h 之间，是一条渐近线 当 12h h h 时 5 （二） 计算 17:37:16 57 当 2 时 同样，当 3 4 1, , , Nh h h h 时 有： 3 3 4 4( , ) , ( , ) ,m e q m e 这些点称为 ， N （不存在 ） 说明： N 层断面有 N 1 个 。 D Z 点之间的曲线为分只枝曲线 5 （二） 计算 17:37:16 58 第 j 分枝曲线 据 以 反推出各电性层的厚度和真电阻率 5 （三） 征 17:37:16 59 （ 1）由 （ 2）各分枝曲线的交点为折点（不光滑），此交点称为 应了电性层的分界面） （ 3） 与 视电阻率曲线对比： 二者首支重合，都以 1 为渐近线 n 为非零有限值时，二者都以 n 为尾支水平渐近线 n 时，二者都以 045 角直线上升 0n 时， 线以 045 下降， 曲线以 063 下降 二者曲线变化趋势基本一致 ()s r 圆滑， 线不圆滑 5 （四） 17:37:16 地电断面的层参数 ( , )变化时 引起的 线变化与电测深曲线变化规律是一致的 条件 依据 线中： 坐标与层参数 ( , )可方便地进行互算 线 i i 对 应 曲 线 60 5 （五）