j)的值越大，其对应的图像梯 度值也越大，但这还不足以确定边缘，因为这里仅将图像快速变换的问题转化成求幅值矩 阵M(i，j)的局部最大值问题；为确定边缘，必须细化幅值图像中的屋脊带，即只保留幅值 局部变化最大的点；这一过程叫非极大值抑制（NMS)，它会生成细化的边缘。
[0039] 上述方案中，所述步骤24包括：首先利用累计统计直方图得到一个高阀值T1，然 后再取一个低阀值T2 ;如果图像信号的响应大于高阀值，那么它一定是边缘；如果低于低 阀值，那么它一定不是边缘；如果在低阀值和高阀值之间，就看它的8个邻接像素有没有大 于高阀值的边缘，如果有，那么它是边缘，否则它不是边缘；观测到模型的边界信息后，用得 到的边界信息对数据的残差的梯度做一定的修正，使得模型梯度对界面更加敏感。
[0040] 上述方案中，步骤3中所述修正后的残差梯度Rgk为：
[0041]
[0042] 上述公式中梯度符合边界是指模型的自己的梯度方向和目标函数梯度模型的梯 度方向同向，使得边界不会消减和平滑。
[0043] 上述方案中，步骤4中所述双边滤波是一种非线性的滤波方法，是结合图像的空 间邻近度和像素值相似度的一种折衷处理，同时考虑空域信息和灰度相似性，双边滤波器 比高斯滤波多了一个高斯方差sigma-d，它是基于空间分布的高斯滤波函数，所以在边缘附 近，离的较远的像素不会太多影响到边缘上的像素值，这样就保证了边缘附近像素值的保 存，其中：
[0044] 高斯滤波器可以定义如下：
[0045]
(13)
[0046] 和似的，但源汲益疋乂卯卜：
[0047]
(14)
[0048] 结合几何域和值域过滤，从而同时考虑几何和值域两个方面，合并的滤波方法被 定义如下：
[0049] w(i, j) = ws(i, j)wr(i, j) (15)
[0050] 合并的几何域和值域的滤波将被表示为双边滤波，双边滤波器被定义为：
[0051]
(16)
[0052] 利用上述定义的双边滤波器对残差梯度Rgk进行滤波，去除噪声对模型的干扰，得 到保持边缘结构信息的残差梯度RFgk。
[0053] 上述方案中，步骤6中所述在步骤5得到最佳步长ak后，即为模型在共 轭梯度的方向lif的改变量，故能够计算模型的更新而基于该模型的 更新由低频到高频进行反演，具体包括：在反演过程中频率f的选择是由低频到高频依次 进行反演，即在公式2中选取对应的《 (?= 2 Jr /f)，再将较小频率的反演结果作为下一次 更高频率的初始模型，这样逐步由低频到高频依次进行反演。
[0054] (三）有益效果
[0055] 从上述技术方案可以看出，本发明具有以下有益效果：
[0056] 1、本发明提供的这种基于canny边缘检测算子和双边滤波的边缘引导和结构约 束的全波形反演方法，利用canny算子检测模型边缘，即在全波形反演中速度和界面同时 反演，再利用双边滤波对模型的残差的梯度进行保持结构的正则化和结构约束的正则化， 解决了常规全波形反演方法的边缘信息和结构信息模糊的问题，并且可以克服强噪声去噪 过程中对模型结果和边界的损坏，增强反演结果结构和边界信息。
[0057] 2、本发明提供的这种基于canny边缘检测算子和双边滤波的边缘引导和结构约 束的全波形反演方法，由于在第2步中省去了由偏移成像得到边界信息来引导全波形反 演，故本方法具有节约快速有效的特点；由于在第3步中采用双边滤波这种较常规滤波方 法（如高斯滤波）能更好的保护和增强模型的边界和结构的滤波方法，同时又能很好的去 除噪声的影响，故能使反演结果的结果和边界更加清晰，反演结果也会更加准确。同时对一 些含强噪声的数据（如陆地数据）会有更好的效果；由于在第5步中加入了结果约束的正 则化，使得全波形朝着结果更简单的方向反演，避免了反演的多解性；进而，显著的提高了 全波形反演的效率和反演结果的准确度。
[0058] 3、相对于传统的方法，本发明提供的这种基于canny边缘检测算子和双边滤波的 边缘引导和结构约束的全波形反演方法，通过将图像处理里面的边缘检测和自动识别模型 结构的双边滤波方法引入波形反演，显著提高了计算效率，反演结果的准确度。本发明是采 用canny算子检测模型边缘，即在全波形反演中速度和界面同时反演，再利用双边滤波对 模型的残差的梯度进行保持结构的正则化和结构约束的正则化，在频率域由低频到高频进 行反演，有效地解决了常规全波形反演方法的边缘信息和结构信息模糊的问题，并且在一 定程度上减小了全波形反演得周期跳跃的问题。另外，本发明采用的步骤少，加快了收敛速 度；并且速度和界面同时迭代，大幅度地提升了反演结果的准确度。
[0059] 4、本发明提供的这种基于canny边缘检测算子和双边滤波的边缘引导和结构约 束的全波形反演方法，在每次全波形反演迭代中，可以通过canny边缘检测算子得到反演 模型的边界信息，相对于全波形反演此过程仅仅需要很少的计算量；观测到模型的边界信 息后，用得到的边界信息对数据的残差的梯度做一定的修正，使得模型梯度对界面更加敏 感。此外本发明们利用得到边界信息对全波形反演做一定结构约束正则化，使反演结果的 结构更加简单。
【附图说明】
[0060] 图1是本发明提供的基于canny边缘检测算子和双边滤波的边缘引导和结构约束 的全波形反演快速方法的流程图；
[0061] 图2是本发明提供的用canny边缘检测算子检测模型边缘的方法流程图；
[0062] 图3是依照本发明实施例的基于canny边缘检测算子和双边滤波的边缘引导和结 构约束的全波形反演快速方法的流程图；
[0063] 图4是依照本发明实施例的第七步中抛物线拟合法求取优化步长的函数曲线示 意图。
[0064] 图5是依照本发明实施例的反演结果界面和常规全波形反演结果界面的对比，其 中（a)是真实速度模型边界；（c)是常规全波形反演界面结果；（c)是边缘引导全波形反演 界面结果。
【具体实施方式】
[0065] 为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照 附图，对本发明进一步详细说明。
[0066] 本发明利用canny算子检测模型边缘，即在全波形反演中速度和界面同时反演， 再利用双边滤波对模型的残差的梯度进行保持结构的正则化和结构约束的正则化，解决常 规全波形反演方法的边缘信息和结构信息模糊的问题，并且可以克服强噪声去噪过程中 对模型结果和边界的损坏，增强反演结果结构和边界信息。
[0067] 如图1所示，图1是本发明提供的基于canny边缘检测算子和双边滤波的边缘引 导和结构约束的全波形反演快速方法的流程图，该方法包括以下步骤：
[0068] 步骤1 :在频率域计算目标函数J(m)和目标函数梯度gk;
[0069] 频率域全波形反演的目标是使得下面数据残差最小，为了让反演结果的结构更加 简单，本发明利用得到边界信息对全波形反演做一定正则化，因此给目标函数增加一个正 则化项J e(m);
[0070] 现在，目标函数定义如下：
[0071] J (m) = min {JD (m) + A Je (m)} (1)
[0072] 其中JD(m)是数据的拟合差，Je(m)是正则化项；
[0073] 则反演目标函数为：
[0074]
[0075] 其中m是反演的速度模型，《是单个频率ns和nr代表震源和检波器的个数，d°bs 和cfal是观测和对应计算的地震波形，Lm是模型的一阶差分算子；
[0076] 通过计算目标函数J(m)，能够得知反演残差的下降趋势，目标函数J(m)对模型的 导数就是目标函数梯度g k，通过目标函数对模型求导计算目标函数梯度gk:
[0077]
[0078] 计算目标函数梯度gk能够得知数据残差在本次反演迭代过程中下降时，模型改 变的最优方向。
[0079] 在本步骤中计算的目标函数J(m)和目标函数梯度gk在后续步骤3中有使用。
[0080] 步骤2 :用canny边缘检测算子检测模型边缘，在全波形反演中速度和界面同时反 演，同时得到速度模型和速度边缘模型，利用得到的边缘信息对下一次迭代做结构约束；
[0081] Canny边缘检测算法是被业界公认的性能最为优良的边缘检测算法。严格来说，它 不是一个简单算子（如R〇berts、Prewitt、Sobel等）或锐化