本发明属于地震波形反演领域，具体涉及一种基于模型平滑算法的全波形反演梯度预处理方法。
背景技术：
上世纪八十年代，tarantola等人首先基于最小平方理论提出了时间域全波形反演方法，该方法充分利用实测地震记录的走时、振幅以及相位等信息来重建地下介质的速度结构，理论上对于理想观测系统其精度可达到波长数量级，因此全波形反演被认为是能够进一步提升地震勘探的油气勘探能力的重要方法，有望在未来的油气勘探中发挥重要作用。
全波形反演是以理论地震记录与实际地震记录误差的最小平方泛函作为反演的目标函数，应用双程波动方程有限差分数值模拟方法分别模拟炮点震源子波的正向传播波场与接收点理论地震记录与实际地震记录残差的逆时传播波场，然后应用正向及逆时波场的时间导数值生成梯度方向并计算共轭梯度方向，最后通过梯度寻优实现速度模型参数迭代更新。具体实现流程如图1所示。
常规地震全波形反演计算时，通常情况下模型浅层地震波能量要强于深层，并且浅层的地震射线密度也要高于深层，这就使得在迭代反演的过程中浅层的速度更新一直占据主导地位，从而导致深层速度反演收敛慢且精度不高。当前主要有两种策略解决深层反演精度低下的问题，一种是基于“层剥法”的思想，由浅层到深层逐层反演，其存在着难以做到精确层位划分的问题；另一种策略则是在反演的过程中对梯度进行预处理，以校正几何扩散对梯度计算的影响。因此梯度预处理方法在实际全波形反演中得到广泛应用。
pratt与1998年提出了应用hessian矩阵或近似hessian矩阵对梯度进行预处理的方法，随后诸多学者在hessian矩阵类梯度预处理方向上做了大量的研究工作。hessian矩阵或近似hessian的梯度预处理方法通常都需要存储hessian矩阵并需进行求逆运算，而hessian矩阵或近似hessian维数通常都非常大(为模型网格点个数的平方)，因此若应用该类方法基于较大模型进行反演，将产生巨大的内存消耗且由此带来的大型矩阵的求逆也极易产生不稳定现象，因此该类方法仅适合于小模型数据的反演。
zhang等于2012年seg年会上提出了一种应用地震波能量进行梯度预处理的方法，其应用地下每个网格点正时、逆时波场的地震波能量对梯度进行加权，其本质是以波场能量来作为每个点的几何扩散校正因子，以使得浅中深部的梯度值达到均衡。zhang的方法是基于当前迭代得到的模型进行波动方程模拟，获得各点的能量值，并以其作为几何扩散及照明度不均衡的校正因子。理论上，应用透射波场的能量作为几何扩散及照明度不均衡的校正因子更为精确，但由于当前迭代得到的模型往往已具备成层性，因此基于当前迭代得到的模型模拟得到的波场中除透射波能量外还包含了大量的反射波/多次波，反射波/多次波的能量会干扰梯度预处理的精度。
本发明基于模型平滑算法对全波形反演中的梯度进行预处理(具体流程见图2所示)，可有效降低波场中的反射波/多次波的影响，提高与处理的精度，从而显著提升反演的效果。
技术实现要素：
本发明要解决的技术问题在于提供一种基于模型平滑算法的全波形反演梯度预处理方法。在全波形反演的每次迭代时，首先将迭代得到的速度模型进行平滑，然后基于此平滑速度模型，应用声波方程有限差分模拟技术，模拟得到炮点震源子波的正时波场能量和以实际地震记录为扰动的逆时波场能量，最后应用此正时波场能量和逆时波场能量对梯度进行预处理。该方法可有效降低波场能量计算中反射波及多次波的影响，提高梯度预处理的精度，从而显著提升全波形反演的效果
本发明采取以下技术方案：
一种基于模型平滑算法的全波形反演梯度预处理方法，其特征在于它具体包括以下步骤：
(1)基于模型，模拟得到炮地震记录作为反演的观测地震记录；在全波形反演的第i(i>1)次迭代计算时，输入上次的梯度值gi-1(x)，速度模型mn、地震子波w及观测地震记录pobs(xr,xs)，其中x表示空间位置坐标，r表示检波点，s表示炮点，xr和xs分别表示检波点和炮点位置；
(2)应用五点三次平滑方法，将速度模型mn平滑n次(n>＝50)，得到新的平滑速度模型
(3)基于平滑速度模型
(4)基于平滑速度模型
(5)利用公式(1)计算基于模型
(6)利用公式(2)计算基于模型
(7)应用公式(3)对上次的梯度gi-1(x)进行预处理，得到新的梯度γi(x)；所述的公式(3)为
全波形反演中每次梯度预处理时均需重复以上7个步骤。
本发明与现有技术相比的有益效果：
本发明提出的一种基于模型平滑算法的全波形反演梯度预处理方法，每次迭代反演时，首先将迭代得到的速度模型进行平滑，然后基于此平滑速度模型，应用声波方程有限差分模拟技术，模拟得到炮点震源子波的正时波场能量和以实际地震记录为扰动的逆时波场能量，最后应用此正时波场能量和逆时波场能量对梯度进行预处理。模型实验结果表明本方法可有效降低波场能量计算中反射波及多次波的影响，提高梯度预处理的精度，从而显著提升全波形反演的效果，有望在实际数据的全波形反演中发挥重要作用。
附图说明
图1基于波场能量的全波形反演流程图；
图2基于模型平滑算法的梯度预处理的全波形反演流程图；
图3overthrust速度模型；
图4初始速度模型；
图5基于波场能量梯度预处理的反演结果；
图6基于模型平滑算法梯度预处理的反演结果；
图7基于波场能量梯度预处理的反演结果(局部)；
图8基于模型平滑算法梯度预处理的反演结果(局部)。
具体实施方式
这里采用overthrust模型(overthrust速度模型见图3)进行具体实施方式阐述。overthrust模型是一个测试地震数据处理、反演方法准确性的著名模型，模型横向长10000m、纵向深2325m。
基于overthrust模型，采用全排列的观测方式，模拟得到共521炮地震记录作为反演的观测地震记录。该炮集记录每炮为521道接收，炮间隔和道间隔均为25m，炮点深度和接收点深度也均为25m。
全波形反演的初始模型如图4所示，下面详细阐述本发明的具体实施过程：
(1)输入上次迭代的梯度值gi-1(x)(i为当前迭代次数)，速度模型mn、地震子波w及观测地震记录pobs(xr,xs)，其中x表示空间位置坐标，r表示检波点，s表示炮点，xr和xs分别表示检波点位置和炮点位置；
(2)应用五点三次平滑方法，将速度模型mn平滑50次，得到新的平滑速度模型
(3)基于平滑
(4)基于平滑速度模型
(5)利用公式(1)计算基于模型
(6)利用公式(2)计算基于模型
(7)应用公式(3)对上次的梯度gi-1(x)进行预处理，得到新的梯度γi(x)。
全波形反演过程中每次迭代时，均需重复以上7个步骤进行梯度预处理。需要特别说明的是在首次迭代(即i＝1)时，不需要进行梯度预处理(即不需要以上步骤处理)。
为了说明本发明所述方法的有效性，这里与基于常规的梯度预处理方法进行对比。如图5为基于常规梯度预处理的全波形反演方法得到的结果，图6为基于模型平滑算法的梯度预处理全波形反演得到的结果。为了更清晰的对比两个结果，将这两个速度剖面选取同样的区域进行局部放大显示，分别如图7、8所示。对比图7和图8可知，基于模型平滑算法的梯度预处理算法，其全波形反演结果的反演精度明显高于基于常规梯度预处理算法的全波形反演。