[发明专利]一种基于交替方向法和全变分理论的地震数据处理方法

权利要求书

1.一种基于交替方向法和全变分理论的地震数据处理方法，其基于建立全变分的地震数据模型，利用经典的二次罚函数方法将此模型转换为满足能够迭代与交替方向最小化的求解问题，通过交替方向法来迭代获得恢复后的地震数据；

所述基于交替方向法和全变分理论的地震数据处理方法主要包括以下步骤：

首先在野外勘探目标区中在地表以人工方法激发地震波，利用地震采集设备，获得地震数据；将野外采集到的地震数据进行格式转换，满足数据处理的格式需求；

然后采用下述步骤进行数据处理：

步骤1：建立基于全变分的地震数据模型；将全变分理论融入到带有噪声的原始地震反褶积模型中，建立成基于全变分的地震数据模型，能免于受噪声的影响，其目标是：在给定合适的褶积算子情况下，可从含有噪声的观测地震数据中恢复原始地震数据；

步骤2：对建立的地震数据恢复目标函数进行重写和转化，即利用二次罚函数方法将“有约束优化”转化为“无约束优化”问题，使其能够迭代与交替方向最小化求解；

步骤3：利用交替方向法(ADM)对目标函数公式(8)进行迭代与交替方向最小化求解；即在步骤2的基础上，利用交替方向法(ADM)对转化后的模型进行迭代与交替方向最小化求解，这种求解方法可以得到每一步迭代过程中的近似最优解；

步骤4：利用第三步获得数据，通过交替方向法来迭代获得恢复后的地震数据，能够在合理的时间内得到全局最优解；

其中，所述建立基于全变分的地震数据模型具体为：设有n个地震道记录的地震剖面，每个地震道记录有m个采样点，则用矩阵r来表示此n·m个数据，其元素为x_ij，i＝1,2,…,n；j＝1,2,…,m，则地震数据r被表示为：

令x^*＝vec(r)为矩阵r的拉伸向量，即x^*∈R^nm；K∈R^nm×nm为褶积算子，w∈R^nm为添加的随机噪声，f∈R^nm为观测的地震数据，满足如下关系：

f＝Kx^*+w(2)

对于给定的褶积算子K，其目标是从含有噪声的观测地震数据f中恢复原始地震数据x^*，即视为去噪和反褶积处理；由于||Kx-f||₂≈||w||₂，恢复原始地震数据x^*为最小化公式(3)：

${\mathrm{\Φ}}_{fid}(x,f)=\left|\right|Kx-f|{|}_2^2---\left(3\right)$

其中，||·||₂表示L₂范数，Φ_fid(x,f)称为随机噪声的数据保真项；对于公式(2)，直接从f中恢复x^*为不适定问题；为了使恢复x^*具有稳定性，引入TV正则项，即全变分(TV)的离散形式：

$TV\left(x\right)=\underset{i=1}{\overset{n\·m}{\Σ}}\left|\right|D_{i}x|{|}_2---\left(4\right)$

其中，D_i为局部有限差分算子，D_ix∈R²为x在第i个采样点处包含水平和垂直两个方向的一阶有限差分；TV正则化也有利于保留数据锋利的边缘和边界；地震数据的恢复为解决如下最小化问题：

$\underset{x}{min}TV\left(x\right)+{\mathrm{\μ\Φ}}_{fid}(x,f)---\left(5\right)$

即地震数据的重建模型为：

$\underset{x}{\min }\underset{i=1}{\overset{n\·m}{\Σ}}\left|\right|D_{i}x|{|}_2+\frac{\mathrm{\μ}}{2}||Kx-f|{|}_2^2---\left(6\right)$

其中，μ＞0为平衡二者的参数；x^*为原始地震数据；x为被恢复的地震数据；

对建立的地震数据恢复目标函数进行重写和转化具体为，将目标函数公式(6)重写为公式(7)形式，则有：

$\left\{\begin{array}{c}\underset{x,y}{\min }{\mathrm{\Σ}}_{i=1}^{n\·m}\left|\right|y_i|{|}_2+\frac{\mathrm{\μ}}{2}||Kx-f|{|}_2^2\\ \begin{array}{cc}s.t.& y_i=D_{i}x,i=1,...,n\·m\end{array}\end{array}\right.---\left(7\right)$

其中，y＝(y¹；y²)∈R^2nm，y¹和y²是长度为n·m维向量，满足((y¹)_i；(y²)_i)＝y_i∈R²；利用二次罚函数方法将公式(7)转化为无约束目标函数：

$\underset{x,y}{\min }\underset{i=1}{\overset{n\·m}{\Σ}}\left(\right|\left|y_i\right|{|}_2+\frac{\mathrm{\β}}{2}\left|\right|y_i-D_{i}x\left|{|}_2^2\right)+\frac{\mathrm{\μ}}{2}\left|\right|Kx-f|{|}_2^2---\left(8\right)$

其中，β＞＞0为罚参数；

利用交替方向法(ADM)对目标函数公式(8)进行迭代与交替方向最小化求解，初始x＝x^k，λ＝λ^k：

$\left\{\begin{array}{c}{y}^{k+1}\←\arg {\min }_y{L}_A(x^k,y,{\mathrm{\λ}}^k)\\ x^{k+1}\←\arg {\min }_x{L}_A(x,y^{k+1},{\mathrm{\λ}}^k)\\ {\mathrm{\λ}}^{k+1}\←{\mathrm{\λ}}^k-\mathrm{\β}(y^{k+1}-{\mathrm{Dx}}^{k+1})\end{array}\right.---\left(9\right)$

其中，L_A(x,y,λ)为公式(8)的增广Lagrangian函数：

$\begin{array}{c}{L}_A(x,y,\mathrm{\λ})=\underset{i=1}{\overset{n\·m}{\Σ}}\left(\left|\right|y_i|{|}_2-{\mathrm{\λ}}_i^T\left(y_i-D_{i}x\right)+\frac{\mathrm{\β}}{2}||y_i-D_{i}x|{|}_2^2\right)\\ +\frac{\mathrm{\μ}}{2}\left|\right|Kx-f|{|}_2^2\end{array}---\left(10\right)$

λ_i为R²空间中的向量；公式(9)中argmin_yL_A(x^k,y,λ^k)等价于如下形式：

$\underset{y_i\∈R^2}{\min }\left|\right|y_i|{|}_2+\frac{\mathrm{\β}}{2}||y_i-(D_i{x}^k+\frac{1}{\mathrm{\β}}{\left({\mathrm{\λ}}^k\right)}_i)|{|}_2^2,i=1,...,n\·m---\left(11\right)$

公式(11)通过二维收缩来解决：

$y_i^{k+1}=\max \left\{\left|\right|D_i{x}^k+\frac{1}{\mathrm{\β}}{\left({\mathrm{\λ}}^k\right)}_i|{|}_2-\frac{1}{\mathrm{\β}},0\right\}\frac{D_i{x}^k+\frac{1}{\mathrm{\β}}{\left({\mathrm{\λ}}^k\right)}_i}{\left|\right|D_i{x}^k+\frac{1}{\mathrm{\β}}{\left({\mathrm{\λ}}^k\right)}_i|{|}_2}---\left(12\right)$

在固定y^k+1和λ^k基础上，求解公式(9)中arg min_xL_A(x,y^k+1,λ^k)为一最小二乘问题，其对应的标准方程为：

$(D^{T}D+\frac{\mathrm{\μ}}{\mathrm{\β}}{K}^{T}K)x=D^T(y^{k+1}-\frac{1}{\mathrm{\β}}{\mathrm{\λ}}^k)+\frac{\mathrm{\μ}}{\mathrm{\β}}{K}^{T}f---\left(13\right)$

其中，D＝(D¹；D²)∈R^2nm×nm；

以上各个数学符号的物理意义

m：地震道记录的采样点个数；

n：地震道记录个数；

r：各地震道与每个地震道采样点组成的n×m矩阵；

x_ij：r中的元素，指第i个地震道、第j个采样点值，其中，i＝1,2,…,n；j＝1,2,…,m；

vec(r)：矩阵r的拉伸向量；

x：被恢复的地震数据；

K：褶积算子；

f：观测的地震数据；

R^nm：n·m维度空间；

w：添加的随机噪声；

||·||₂：指L₂范数；

Φ_fid(x,f)：随机噪声的数据保真项；

D：一阶全局有限差分算子；

D¹,D²：分别为水平和垂直方向的一阶全局有限差分算子；

D_i：局部有限差分算子，以D¹的第i行为第一行、以D²的第i行为第二行的两行矩阵；

TV(x)：全变分的离散形式；

μ：平衡数据保真项Φ_fid(x,f)和全变分形式TV(x)的参数；

y_i＝D_ix＝((y¹)_i；(y²)_i)：x在第i个采样点处包含水平和垂直两个方向的一阶有限差分，是一2维向量；

y＝(y¹；y²)：2nm维向量；

y¹,y²：长度为n·m维的向量，分别指水平和垂直两个方向的一阶有限差分；

β：二次罚函数方法中的罚参数；

λ：Lagrangian函数参数；

λ_i：R²空间中的向量；

x^k,y^k,λ^k：分别是x,y,λ的第k次迭代更新值；

L_A(x,y,λ)：增广Lagrangian函数。