一种边缘引导和结构约束的全波形反演快速方法
【技术领域】
[0001] 本发明涉及地震成像及全波形反演技术领域，尤其是一种基于canny边缘检测算 子和双边滤波的边缘引导和结构约束的全波形反演快速方法，其中canny边缘检测是John F. Canny于1986年开发出来的一个多级边缘检测算法。
【背景技术】
[0002] 地震反演问题是利用观测到的地震数据确定地下的结构，是和地震正问题紧密相 关的。根据反演时使用信息的不同，反演方法可以分为走时反演、振幅反演、全波形反演三 大类。利用波场运动学信息反演出来的结果反映的是较大尺度地下速度模型的宏观特征， 对于理想观测系统，它的分辨率可达到菲涅耳带（Fresnel zone)数量级。利用完整波场信 息反演的方法则具有精确刻画模型细节的潜力，对于理想观测系统，反演结果的分辨率可 达到波长数量级。
[0003] 全波形反演方法利用叠前地震波场的运动学和动力学信息重建地下速度结构，具 有揭示复杂地质背景下构造与岩性细节信息的潜力，其分辨率比传统基于射线的反演高。 与走时反演相比，全波形反演由于能拟合波形记录的振幅/相位信息，因此能得到一个更 准确更精细的速度模型。通常全波形反演包含两部分：波形的正演和反演。正演通过有限 差分，有限元的数学方法解声波或弹性波方程得到波形记录。对于全波形反演，通过共轭梯 度，高斯牛顿等非线性反演算法逐步迭代更新速度模型。
[0004] 根据研宄需要，全波形反演既可在时间域也可在频率域实现。频率域相对于时间 域反演具有计算高效、数据选择灵活等优势。国内频率域全波形反演理论研宄起步相对较 晚，反演理论研宄程度较低，尚未充分挖掘其应用潜力。频率域全波形反演方法具有计算高 效、反演频率选择灵活等优点。目前频率域波形反演理论在波场正演方法、反演频率选择 策略、目标函数设置方式、震源子波处理方式、梯度预处理方法等方面已取得了可喜的进 展，提高了全波形反演方法的实用性。理论模型试验和实际井间地震、广角地震以及垂直地 震剖面（VSP)资料的反演结果初步展示了该方法的应用潜力。多波多分量数据和三维短 偏移距反射地震资料的应用研宄将是频率域全波形反演的下一个重要突破方向（卞爱飞， 2010)。相比于走时成像，全波形反演具有更强的非线性，因此需要一个较为准确的初始速 度模型。因此在进行全波形反演之前，通常需要利用走时成像获得初始的速度模型。另外 对于波形数据，也需要进行必要的预处理。频率域波形反演的预处理主要包括选择波形反 演的时间窗，通过F-K滤波去除面波，去除S波和去掉噪音，对数据进行重采样，振幅矫正和 振幅归一化。
[0005] 现有技术中的全波形反演主要存在如下缺陷；（1)反演解非唯一性、易陷入周波 跃迀和局部极值；（2)强烈依赖初始模型、需要低频分量和初始速度模型很好的耦合；（3) 对噪声敏感以及计算量大等问题；而且，在数学上它是一个高度病态的非线性问题，需要解 决其多解性及收敛性问题。
[0006] 为了克服这些问题的影响，国内外发展出了基于偏移图像引导的全波形算法（Ma and Hale，2011)和基于偏移的全波形反演算法。偏移是提高成像质量的有效手段，其挑战 主要表现在成像算法、高精度的速度模型、计算效率等方面。全波形反演本身是一个高花费 的反演过程，偏移过程也同样如此。因此，基于偏移的全波形花费会更加巨大。界面信息是 全波形反演研宄中非常重要的信息，界面中包含着地下介质结构等构造信息，这决定着后 续的结果如地震偏移等是否准确。界面的信息主要包含在数据的高频成分里面，其对全波 形反演结构准确与否有着重要的指导意义。现有的常规全波形反演是仅仅反演速度模型， 没有考虑结构和边界的影响，并且在每次迭代中对速度模型进行了大量的平滑，这对去除 噪声有一定的抑制，却丧失了地下模型的界面和边界信息。
[0007] 因此，需要探索一种新方法解决上述问题，以提高全波形反演的速度结构和界面 的准确度。
【发明内容】
[0008] (一）要解决的技术问题
[0009] 有鉴于此，本发明的主要目的在于提供一种快速和有效的基于canny边缘检测算 子和双边滤波的边缘引导和结构约束的全波形反演方法，以解决常规全波形反演方法存在 的边缘信息和结构信息模糊的问题，克服去噪过程对模型结果和边界的损坏，增强反演结 果结构和边界信息。
[0010] (二）技术方案
[0011] 为达到上述目的，本发明提供了一种边缘引导和结构约束的全波形反演快速方 法，该方法包括：步骤1 :在频率域计算目标函数J (m)和目标函数梯度gk;步骤2 :用canny 边缘检测算子检测模型边缘，在全波形反演中速度和界面同时反演，同时得到速度模型和 速度边缘模型，利用得到的边缘信息对下一次迭代做结构约束；步骤3 :利用步骤2中描述 的canny算子及其探测边界的方法探测模型边界，利用探测得到的模型边界信息修正步骤 1中得到的目标函数梯度g k，得到修正后的残差梯度Rgk;步骤4 :利用双边滤波对残差梯 度Rgk进行滤波，去除噪声对模型的干扰，得到双边滤波后的残差梯度RFg k;步骤5 :在步骤 4中得到双边滤波后的残差梯度RFgk，结合上一次迭代的方向得到本次迭代的方向 ，再利用抛物线拟合法在共轭梯度的方向搜索一个最佳步长a k，其中本次迭代方 向和上一次迭代方向存在共轭的数学关系；步骤6 :在步骤5得到最佳步长a k后，计算模 型的更新7叫，基于该模型的更新由低频到高频进行反演。
[0012] 上述方案中，所述步骤1包括：
[0013] 目标函数定义如下：J(m) =min{JD(m) +XJe(m)} (1)
[0014] 其中JD(m)是数据的拟合差，上(m)是正则化项；
[0015] 则反演目标函数为：
[0016]
[0017] 其中m是反演的速度模型，《是单个频率，ns和nr代表震源和检波器的个数，d°bs 和cT1是观测和对应计算的地震波形，L是模型的一阶差分算子；
[0018] 通过计算目标函数J(m)，能够得知反演残差的下降趋势，目标函数J(m)对模型的 导数就是目标函数梯度g k，通过目标函数对模型求导计算目标函数梯度gk:
[0019]
[0020] 计算目标函数梯度gk能够得知数据残差在本次反演迭代过程中下降时，模型改变 的最优方向。
[0021] 上述方案中，所述步骤2包括：步骤21 :引入高斯滤波器，对原始图像进行高斯平 滑滤波，以提高算法的抗噪性；步骤22 :计算平滑滤波后图像的梯度强度和方向信息；步骤 23:利用得到的梯度强度和方向信息划分进行梯度强度的非极大抑制，获取单像素边缘点； 步骤24:利用获得的单像素边缘点，采用双阈值进行边缘的二值化，其中进行边缘的二值 化时仅选择在双阈值内的边缘点。
[0022] 上述方案中，所述步骤21包括：
[0023] 一个二维高斯函数为：
[0024] 在某一方向n上是G(x，y)的一阶方向导数为：
[0025]
[0026]
[0027] 式中：n式方向矢量，▽G是梯度矢量；
[0028] 用f(i，j)表示输入图像，使用可分离滤波方法求图像与高斯平滑滤波器卷积， 得到的结果是一个已平滑的数据矩阵：
[0029] S(i, j) =G(i, j ；〇)*f(i, j) (7)
[0030] 其中〇是高斯函数的标准差，用于控制着平滑程度。
[0031] 上述方案中，所述步骤22包括：已平滑数据矩阵S(i，j)的梯度是使用2X2 -阶 有限差分近似式来计算x与y偏导数的两个矩阵P(i，j)与Q(i，j):
[0032] P(i, j) ^ (S (i+1, j) -S (i, j) +S (i+1, j+1) -S (i, j+1))/2 (8)
[0033] Q(i, j) ^ (S (i, j+1) -S (i, j) +S (i+1, j+1) -S (i+1, j))/2 (9)
[0034] 在这个2X2方形内求有限差分的均值，以便在图像中的同一点计算x和y的偏导 数梯度；幅值和方向角采用直角坐标到极坐标的坐标转化公式来计算：
[0035]
[0036]
[0037] M(i，j)反映了图像的边缘强度；0 (i，j)反映了边缘的方向，使得M(i，j)取得局 部最大值的方向角0 (i，j)就反映了边缘的方向。
[0038] 上述方案中，所述步骤23包括：幅值图像阵列M(i，