线，并求取极小值点，以此作为优化步长；
[0136] 第八步，用优化步长和梯度场来更新速度场，直到满足迭代终止条件或者达到最 大迭代次数，然后进行下个频率（组）的反演，即重复第二步到第十步；
[0137] 第九步，当所有的频率均反演完毕，反演完毕；
[0138] 第十步，循环结束。
[0139] 如图5所示，图5是依照本发明实施例的反演结果界面和常规全波形反演结果界 面的对比，其中（a)是真实速度模型边界；(b)是常规全波形反演界面结果；（c)是边缘引 导全波形反演界面结果。本实施例用marmousi模型来验证本发明的有效性。
[0140] 所选用模型大小为8940mX 3000m，纵横采样间距均为30m，所选炮的数据为炮点 个数为200炮，每炮445道，记录长度为48秒，时间采样间隔为4毫秒，炮间距为110m。
[0141] 反演参数设置为反演频率从4Hz到25Hz选取6个频率组，反演按照频率组进行循 环，每个频率组选取两个频率。
[0142] 迭代终止条件为：第i次的目标函数减去第i_l次目标函数的模为第i次目标函 数的千分之一。
[0143] 以上所述的具体实施例，对本发明的目的、技术方案和有益效果进行了进一步详 细说明，所应理解的是，以上所述仅为本发明的具体实施例而已，并不用于限制本发明，凡 在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保 护范围之内。
【主权项】
1. 一种边缘引导和结构约束的全波形反演快速方法，该方法包括： 步骤1 :在频率域计算目标函数J(m)和目标函数梯度gk; 步骤2 :用canny边缘检测算子检测模型边缘，在全波形反演中速度和界面同时反演， 同时得到速度模型和速度边缘模型，利用得到的边缘信息对下一次迭代做结构约束； 步骤3 :利用步骤2中描述的canny算子及其探测边界的方法探测模型边界，利用探测 得到的模型边界信息修正步骤1中得到的目标函数梯度gk，得到修正后的残差梯度Rgk; 步骤4 :利用双边滤波对残差梯度Rgk进行滤波，去除噪声对模型的干扰，得到双边滤 波后的残差梯度RFgk; 步骤5 :在步骤4中得到双边滤波后的残差梯度RFgk，结合上一次迭代的方向llji+4 得到本次迭代的方向ftf，再利用抛物线拟合法在共轭梯度的方向hT搜索一个最佳步长ak，其中本次迭代方向和上一次迭代方向存在共轭的数学关系； 步骤6 :在步骤5得到最佳步长ak后，计算模型的更新叫+尸基于该模型 的更新由低频到高频进行反演。2. 根据权利要求1所述的边缘引导和结构约束的全波形反演快速方法，其中，所述步 骤1包括： 目标函数定义如下：J(m) =min{JD(m)+AJe(m)} (1) 其中JD(m)是数据的拟合差，J?是正则化项； 则反演目标函数为：其中m是反演的速度模型，《是单个频率ns和nr代表震源和检波器的个数，d°b5P(T1是观测和对应计算的地震波形，Lm是模型的一阶差分算子； 通过计算目标函数J(m)，能够得知反演残差的下降趋势，目标函数J(m)对模型的导数 就是目标函数梯度gk，通过目标函数对模型求导计算目标函数梯度gk:计算目标函数梯度gk能够得知数据残差在本次反演迭代过程中下降时，模型改变的最 优方向。3. 根据权利要求1所述的边缘引导和结构约束的全波形反演快速方法，其中，所述步 骤2包括： 步骤21 :引入高斯滤波器，对原始图像进行高斯平滑滤波，以提高算法的抗噪性； 步骤22 :计算平滑滤波后图像的梯度强度和方向信息； 步骤23 :利用得到的梯度强度和方向信息划分进行梯度强度的非极大抑制，获取单像 素边缘点； 步骤24:利用获得的单像素边缘点，采用双阈值进行边缘的二值化，其中进行边缘的 二值化时仅选择在双阈值内的边缘点。4. 根据权利要求3所述的边缘引导和结构约束的全波形反演快速方法，其中，所述步 骤21包括： 一个二维高斯函数为在某一方向n上是G(x，y)的一阶方向导数为：式中m式方向矢量，▽G是梯度矢量； 用f(i，j)表示输入图像，使用可分离滤波方法求图像与高斯平滑滤波器卷积，得到的 结果是一个已平滑的数据矩阵： S(i,j) =G(i,j；〇)*f(i,j) (7) 其中〇是高斯函数的标准差，用于控制着平滑程度。5. 根据权利要求3所述的边缘引导和结构约束的全波形反演快速方法，其中，所述步 骤22包括： 已平滑数据矩阵S(i，j)的梯度是使用2X2 -阶有限差分近似式来计算x与y偏导数 的两个矩阵P(i，j)与Q(i，j): P(i,j) ^(S(i+1,j)-S(i,j)+S(i+1,j+1)-S(i,j+1)) /2 (8) Q(i,j) ^(S(i,j+1)-S(i,j)+S(i+1,j+1)-S(i+1,j)) /2 (9) 在这个2X2方形内求有限差分的均值，以便在图像中的同一点计算x和y的偏导数梯 度；幅值和方向角采用直角坐标到极坐标的坐标转化公式来计算：M(i，j)反映了图像的边缘强度；0 (i，j)反映了边缘的方向，使得M(i，j)取得局部最 大值的方向角9 (i，j)就反映了边缘的方向。6. 根据权利要求3所述的边缘引导和结构约束的全波形反演快速方法，其中，所述步 骤23包括： 幅值图像阵列M(i，j)的值越大，其对应的图像梯度值也越大，但这还不足以确定边 缘，因为这里仅将图像快速变换的问题转化成求幅值矩阵M(i，j)的局部最大值问题；为确 定边缘，必须细化幅值图像中的屋脊带，即只保留幅值局部变化最大的点；这一过程叫非极 大值抑制（匪S)，它会生成细化的边缘。7. 根据权利要求3所述的边缘引导和结构约束的全波形反演快速方法，其中，所述步 骤24包括： 首先利用累计统计直方图得到一个高阀值T1，然后再取一个低阀值T2;如果图像信号 的响应大于高阀值，那么它一定是边缘；如果低于低阀值，那么它一定不是边缘；如果在低 阀值和高阀值之间，就看它的8个邻接像素有没有大于高阀值的边缘，如果有，那么它是 边缘，否则它不是边缘；观测到模型的边界信息后，用得到的边界信息对数据的残差的梯度 做一定的修正，使得模型梯度对界面更加敏感。8. 根据权利要求1所述的边缘引导和结构约束的全波形反演快速方法，其中，步骤3中 所述修正后的残差梯度Rgk为：上述公式中梯度符合边界是指模型的自己的梯度方向和目标函数梯度模型的梯度方 向同向，使得边界不会消减和平滑，t是设定的参数，将影响反演结果边界的清晰度。9. 根据权利要求1所述的边缘引导和结构约束的全波形反演快速方法，其中，步骤4 中所述双边滤波是一种非线性的滤波方法，是结合图像的空间邻近度和像素值相似度的一 种折衷处理，同时考虑空域信息和灰度相似性，双边滤波器比高斯滤波多了一个高斯方差 sigma-d，它是基于空间分布的高斯滤波函数，所以在边缘附近，离的较远的像素不会太多 影响到边缘上的像素值，这样就保证了边缘附近像素值的保存，其中： 高斯滤波器可以定义如下：相似的，值域滤波器定义如下：结合几何域和值域过滤，从而同时考虑几何和值域两个方面，合并的滤波方法被定义 如下： w(i,j) =ws(i,j)wr(i,j) (15) 合并的几何域和值域的滤波将被表示为双边滤波，双边滤波器被定义为：利用上述定义的双边滤波器对残差梯度Rgk进行滤波，去除噪声对模型的干扰，得到保 持边缘结构信息的残差梯度RFgk。10. 根据权利要求1所述的边缘引导和结构约束的全波形反演快速方法，其中，步骤6 中所述在步骤5得到最佳步长ak后，Offcrt+I1即为模型在共轭梯度的方向的改变量，故 能够计算模型的更新狗^尸_ - 而基于该模型的更新由低频到高频进行反演，具 体包括： 在反演过程中频率f?的选择是由低频到高频依次进行反演，即在公式2中选取对应的u(u = 2JT/f)，再将较小频率的反演结果作为下一次更高频率的初始模型，这样逐步由 低频到高频依次进行反演。
【专利摘要】本发明公开了一种边缘引导和结构约束的全波形反演快速方法，利用canny算子检测模型边缘，即在全波形反演中速度和界面同时反演，再利用双边滤波对模型的残差的梯度进行保持结构的正则化和结构约束的正则化，解决了常规全波形反演方法的边缘信息和结构信息模糊的问题，并且克服了强噪声去噪过程中对模型结果和边界的损坏，增强了反演结果结构和边界信息。
【IPC分类】G01V1/28
【公开号】CN104991269
【申请号】CN201510309007
【发明人】向世明, 张海江
【申请人】中国科学技术大学
【公开日】2015年10月21日
【申请日】2015年6月4日