[发明专利]应力值变化方向约束下的二维矿震全波形反演方法

权利要求书

1.一种应力值变化方向约束下的二维矿震全波形反演方法，其特征在于：

步骤一：确定目标区域网格剖分，目标区域被剖分为m×m的矩形网格，m为奇数；应力变化前的波速分布矩阵为VSB，应力变化后的波速分布矩阵为VSA，波速分布矩阵中的元素与网格剖分中的每格波速一一对应，波速均为大于0的实数；

步骤二：将VSB进行二维FFT(Fast Fourier Transform)后所得m×m矩阵为VFB，记此过程为VFB＝FFT2(VSB)，逆过程为IFFT2，有VSA＝IFFT2(VFA)；根据二维FFT性质，VFB中元素均为复数，且幅度谱关于坐标轴对称，相位谱关于原点对称，其中心点坐标为同样，得应力变化后波速分布矩阵的二维FFT矩阵VFA＝FFT2(VSA)；

步骤三：确定频域拟保留坐标长度D，即对于VFB、VFA，位于横坐标范围或纵坐标范围或的元素均为0，记频域截断后的VFB为VFBC，频域截断后的VFA为VFAC；同样的，VFBC、VFAC的幅度谱也关于坐标轴对称，相位谱也关于原点对称，因此VFBC、VFAC中的自由变量数目为4D²+4D+1；

将VFBC、VFAC中的自由变量用向量VFBCV、VFACV表示，并记此过程为M_TO_V，则有VFBCV＝M_TO_V(VFBC)和VFACV＝M_TO_V(VFAC)，记逆过程为V_TO_M，则有VFBC＝V_TO_M(VFBCV)和VFAC＝V_TO_M(VFACV)；

步骤四：确定应力值变化方向和变化位置矩阵SV

其中(svx_isv,svy_isv)为VSB、VSA中发生应力变化的元素坐标，满足1≤svx_isv≤m和1≤svy_isv≤m且为整数，下标满足1≤nsv≤m²且为整数，svd_isv为应力值变化方向，若应力值变大svd_isv＝1，应力值变小svd_isv＝-1，应力值不变svd_isv＝0；

步骤五、计算目标函数

入参：TruSeiB(应力变化前的真实地震记录)、TruSeiA(应力变化后的真实地震记录)、VFBCV(应力变化前的波速分布向量)、VFACV(应力变化前的波速分布向量)、SV(应力值变化方向和变化位置矩阵)、SSP(震源、检波器位置，地震子波等参数)

返回值：TFV(目标函数值)

计算过程：

获得应力变化前、后的理论波速分布矩阵VSB、VSA

VSB＝IFFT2(V_TO_M(VFBCV))

VSA＝IFFT2(V_TO_M(VFACV))

获得应力变化前、后的理论地震记录：

记常密度声波方程正演为FW，应力变化前的理论地震记录为TheoSeiB，应力变化后的理论地震记录为TheoSeiA，有：

TheoSeiB＝FW(VSB,SSP)

TheoSeiA＝FW(VSA,SSP)

计算目标函数值：

TFV＝CHECK(VSB,VSA,SV)×(||TruSeiB-TheoSeiB||₂+||TruSeiA-TheoSeiA||₂)

CHECK为应力分布检查函数，当VSB、VSA中应力值变化方向与SV中完全一致时，函数值为1，否则函数值为2。

2.根据权利要求1所述的应力值变化方向约束下的二维矿震全波形反演方法，其特征在于：所述步骤一中，若目标区域无法剖分为奇数行或奇数列，应补齐一行或一列以满足m为奇数。