具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明实施方式作进一步地描述。

请参考图1，本发明提供了一种基于点扩散函数的Hessian矩阵的计算与存储方法，应用于地震勘探数据处理中的照明分析，最小二乘偏移等方面；

在地震勘探的成像问题中，地震数据与反射系数的模型的关系是非线性的。

d＝L(m) (1)

其中，为d为地震数据，m为反射系数的模型，L为传播算子。

这个非线性的过程在一阶Born近似下可以近似表示为线性，从而得到

d＝Lm (2)

给定地震观测资料估算反射系数模型是一个典型的地震反问题，可以用以下目标函数来求解；

式中，d_obs为地震观测资料，目标函数对模型的一阶偏导数是梯度，二阶偏导数是Hessian算子，让梯度为零，我们可以得到下面的方程：

L^Hd＝L^HLm (4)

其中L^Hd为传统的偏移成像结果将其记为I_mig，L^HL为Hessian算子记为H，上面的式子可以写成如下：

I_mig＝Hm (5)

Hessian矩阵可以被视为一个模糊算子，进而将理想图像m与偏移结果I_mig联系起来，见上式，为了求得理想图像m，上式可以化为：

m＝H^-1I_mig (6)

因此Hessian矩阵的计算就显得十分重要了，在考虑线性化成像问题时，Hessian矩阵可以表达为一个对频率ω和炮检对ξ的二重积分，其中积分的内核为“散射格林函数”，见公式(2)。

H(i,j)＝∫_ω∫ξ[-ω²G(ω,x_r；x_i)G(ω,x_i；x_s)f(ω；s_x)]^*×[-ω²G(ω,x_r；x_j)G(ω,x_j；x_s)f(ω；s_x)]dξdω

＝∫_ω∫ξω⁴|f(ω；s_x)|²G^* _scatter(ω,x_r；x_i；x_s)G_scatter(ω,x_r；x_j；x_s)dξdω (7)

其中，G代表格林函数，G_scatter代表散射波格林函数，T代表旅行时；ξ代表炮检对的索引，f(ω；s_x)代表震源子波，ω代表频率，x_r表示检波点坐标，x_s表示炮点坐标，x_j，x_i表示地下空间点坐标。

根据公式的可以恢复得到理想的图像m，但是包含Hessian矩阵的反演是非平凡的，面临着计算量巨大的挑战。从(2)式中可以看出，计算成本来源于于两方面：

Green函数的大量计算，频率域格林函数是一个频率、炮点和地下空间为独立变量的高维函数，见公式：

其中，表示炮点局部平面波的初始波矢量方向，u_GB(x,x_s,p^s,ω)表示笛卡尔坐标系中一个高斯束，其中，x表示成像点坐标，x_s表示炮点坐标，ω为地震波的频率，i表示为复数。

Hessian矩阵大规模的计算，因为Hessian矩阵为一个大规模的稀疏矩阵，大小为模型维度*模型维度，如地下空间模型的维度为500*500，参考图2；

本发明主要从两个方面来解决上述困难：

1)在局部窗中对Hessian进行计算，从公式(2)可以得到，Hessian的一个元素为模型中某点所有炮检对Green函数的乘积与某点所有炮检对Green函数乘积的内积，参考图3。在模型中，计算某一点的Hessian时，离该点较远的点与该点内积值很小，可以忽略，所以可以选取离该点较近的周围n个点与该点内积，计算出该行Hessian中较大的元素即可，并且忽略接近零的元素显著的降低了Hessian矩阵的规模，参考图4。

2)通过采用WKBJ近似下的渐近格林函数见公式(9)，可以对积分方程(2)进行化简，因此，近似的将Hessian的解析形式在公式(4)中表达为射线振幅与震源子波自相关的积分。由于自相关函数R_ss的衰减特性，其强度随x_i和x_j距离的增大而衰减，这也与PSF的局部支撑性相符合。经过WKBJ近似下的Hessian矩阵计算公式从二重积分转换为一重积分，并且积分内核变为射线振幅以及传播时间控制的自相关函数，且在本发明中的射线的传播时间以及振幅相对于格林函数来说都是低维函数，大大的降低了运算规模。图5中图(a)的箭头代表所需公式(10)中的射线路径(x_r←x_i←x_s)和(x_r←x_j←x_s)。

由此，本发明的具体步骤如下：

S1、基于高斯束并根据背景速度进行射线追踪计算地下空间旅行时；具体如下：

S11、申请结构体数组Cells(Nz，Nx，Npx)储存不同方向的旅行时，申请旅行时数组储存单炮time(Nz，Nx，Nshot)；对Cells(Nz，Nx，Npx)，time(Nz，Nx，Nshot)及涉及的循环相关变量(如：Nz，Nx,Npx,Nshot，Ix，Iz，Ipx等)进行定义；

S12、炮点循环，即对炮点shot进行Nshot次循环；

S13、炮方向循环，即对炮方向shot进行Npx次循环，对一个方向ipx的高斯束进行正演获得射线参数，获得该方向的Cells(Nz，Nx，Npx)；

S14、模型空间循环，循环1-Nx*NZ次；对模型空间内的Iz,Ix网格点，进行不同方向旅行时的累加，即得到单炮旅行时time(Iz，Ix，Ishot)，参考图6；

其中，z和x指空间点二维坐标，Nshot指炮点的循环次数，shot指炮点；其中I＝(1，N)。

S2、采用短波逼近法WKBJ近似下的渐近格林函数以及在局部窗中近似估计得到PSFs，包括：

S21、申请Hessian数组H(Ni，Nj)

S22、炮点循环，即对炮点Ishot进行Nshot次循环；

S23、检波点循环，对任意检波器Ireceiver进行处理，对检波点进行1-Ntrace次循环；

S24、模型空间循环，循环1-Ni*Nj次,Ni*Nj为Hessian矩阵模型空间大小，模型空间积分输出：

其中，R_ss为自相关函数，A_ray代表振幅可通过旅行时获得，T代表旅行时，ξ代表炮检对的索引，f(ω；s_x)代表震源子波，ω代表频率，x_r表示检波点坐标，x_s表示炮点坐标，x_j，x_i表示地下空间点坐标，参考图7；

S3、利用得到的PSFs对偏移图像进行校正，得到校正后的图像，公式如下：

m＝H^-1I_mig

其中，偏移成像结果记为I_mig，理想图像记为m。

为了测试本发明的有效性，选择了sigsbee2a模型进行测试；图8为sigabee2a模型真实反射系数，图9是sigabee2a模型通过背景速度得到的偏移图像，图10为sigabee2a模型通过PSFs对偏移图像进行校正之后得到的图像。可以看出本发明计算的Hessian矩阵对偏移图像进行校正后得到的像层位清楚，能量比较平衡。

对所公开的实施例的上述说明，使本领域专业技术人员能够实现或使用本发明。对这些实施例的多种修改对本领域的专业技术人员来说将是显而易见的，本文中所定义的一般原理可以在不脱离本发明的精神或范围的情况下，在其它实施例中实现。因此，本发明将不会被限制于本文所示的这些实施例，而是要符合与本文所公开的原理和新颖特点相一致的最宽的范围。