-1)/2)，；[1 1={0，1，2，，"，1111-1}，11=0，1， 2,3,函数floor()表示向下取整，kikXix表示四维向量kik与Xix的内积。并且ik与i0、ii、i2、i3 保持关系式ik=((i0 · m+ii) · n2+i2) · n3+i3,其中，11()、111、112、113的定义见步骤(3)。本例中， c0 - 5, ci - 5，C2 - 32，C3 - 8 0
[0194] (6)、对每一个频率切片ico = {〇，1，2,…，icomax}，执行如下循环，并存储yiWD
[0195] biw = FHdiw
[0196] yiW = CG(biw ,Ψ?ω)
[0197] UiGj = conj(yiGj) · *yi〇j〇
[0198] Wi〇j+i = (sqrt(UiGj/sum(UiGj) )+Wi〇j)/2
[0199] 步骤(6)所述的icomax是根据实际需要选取的频率上限索引，选取原则是使得地震 数据d(Xlx，t lt)的有效信号能量刚好处于的频带内，其中=
[0200] 步骤(6)所述的wo为一个包含n〇mn2n3个元素的列向量，且wo的每一个元素都等于
其中，WQ是Wiu中下角标i ω =〇时的情况。
[0201 ] 步骤(6)所述的F=F(Xix，kik)，为一个Νχ行、η〇η?Π 2η3列的矩阵。
[0202] 步骤(6)所述的 = [(6/(X。,竭^)，火xp辦x2,qJ,. .,?/(χ νι,^ Nx个元素的列向量。
[0203] 步骤(6)所述的函数conj()表示对输入向量的每一个元素求复共辄，并输出一个 向量;运算符·*表示两个向量的对应元素相乘，并得到一个向量；函数sqrt()表示对输入 向量的每一个元素开平方，并输出一个向量；函数sum()表示对输入向量的所有元素求和， 并输出一个标量;运算符η表示求矩阵的共辄转置，下同。
[0204] 步骤(6)所述的函数y = CG(b，w)表示如下预条件共辄梯度算法，①设定常量。包括 迭代终止相对容差tol = l(T4,最大迭代次数1]^1；^ = 30，纟=?乍;@计算迭代初始变量。包括
' r〇 = b，ρ-ι = 1;③计算迭代终止绝对容
i = {0，l，2，'"，maxit-l}，执行如下循环
[0205]
，或者i=maxit，跳出循环，返回y = yi，否则继续循环
[0206] zi=w ·
[0207] Pi = riHZi
[0208] qi = Api
[0209] ai = Pi/piHqi
[0210] yi+i = yi+aipi
[0211 ] ri+i = ri-〇iqi
[0212] pi+i = Zi+(Pi/Pi-i)pi
[0213] 步骤(6)所述的函数y = CG(b,w)中，每次迭代需要计算一次矩阵与向量乘积运算q =Ap，其中A = FHF具有多级分块Toeplitz矩阵结构，可以采用如下快速算法计算q = Ap。
[0214] ①设 Ν〇 = ηιη2η3、Νι = η2η3、Ν2 = η3、Ν3 = 1= ，Π = 〇，1，2，3。按如下方式生成 长度为mn的向量g和略，其中η = 〇，1，2,3。本例中，~=12870、沁=1170、吣=18、吣=1。_ = 21，mi = 21，m2 = 129，m3 = 35。
[0217] ②设向量a= {a(i) I i = ia(i。，ii，i2, i3)，in=0,1，2,…，mn-1，η = 〇, 1，2,3}是一个 四维数组的向量化表示，其中，：[3(：[()4142 43) = ((：[().1111+；[1).1]12+12).1113+13 4为向量索 引，io，U，i2，i3为四维数组索引。按如下方式对a的元素赋值。初始化变量 ζ=0., t\=nun,n-,n. -1,
[0218] 对于i〇= {〇，l，2,···，m『l}，执行如下循环⑴至（12)
[0219] (l):d+i;;(，0)
[0220] (2) = /；,+!〇〇；)
[0221] (3)对于ii= {0,1，2，···，执行如下循环(4)至（12)
[0222] (4) < =(+i;'(z·!)
[0223] (5) ^=^+15(/,)
[0224] (6)对于i2= {0,1，2，···，m2-l}，执行如下循环(7)至（12)
[0225] (7.): i3r =十i;(/2).
[0226] W 4 =ξ+\\(?2)
[0227] (9)对于i3= {0,1，2，···，m3-l}，执行如下循环（10)至（12)
[0228] (10)=广；十Γ,(/,)
[0229] (11) t\ =t\ +13(/,)
[0230] (12) a(4(44,44)>= A(i:，〇;
[0231] 其中，i(i)表示向量i的第i个元素，a(i)表示向量a的第i个元素，A(ir，i。)表示矩 阵A的第i r行、i。列的元素。
[0232] ③对向量a所表示的四维数组做四维快速傅里叶变换，并返回向量i。其中，向量 3 =丨>'>|/_二4(/;|,/。/2,/_.;),/:7=_体1又~.，^-丨,/广 (〕，丨二3|是一个四维数组的向量化表示， , /p /2, /;)二((/:, ' /?:i + / j· ?ζ2 + /2) _ "2; + /;，i 为向量索引，i。，i i，i 2，i 3 为四维数组索引。 ? = FfT4 (a)，函数FFT4()表示对四维数组做四维快速傅里叶变换，并返回一个四维数组。
[0233] ④设向量p={p(i) | i = ip(i〇，ii，i2，i3)，in=0，l，2，···，ηη_1，η = 〇，1，2,3}是一个 四维数组的向量化表示，其中，ip(i()，il，i2，i3) = ((i()·nl+il)·n2+i2)·n3+i3，i为向量索 引，io, ii，i2，i3为四维数组索弓丨。设向量ρ7 = {p7 (i) I i = iP，（io，ii，i2，i3)，in = 0,1，2，···， 11111-1，11 = 〇，1，2,3}是一个四维数组的向量化表示，其中，；[1/(；[。414243) = ((：[。*1]11+；[1)· 1112+12)*1113+13，：[为向量索引，：[()4142 43为四维数组索引。按照如下方式对向量口/的各个元 素赋值
[0235] ⑤计算//7厂/#·〒/^4(P'))，其中，函数FFT4()和IFFT 4()分别表示对四维数组 做四维快速傅里叶变换和逆变换，并返回一个四维数组;运算符·*表示两个四维数组的对 应元素相乘，并得到一个四维数组；向量q' = {q' (i) I ? = (io，ii，i2，i3)，in = 〇, 1，2,…， 11111-1，11 = 〇，1，2,3}是一个四维数组的向量化表示，其中山/(；[。414243) = ((：[。*1]11+；[1)· m2+i2) · m3+i3，i为向量索引，:1()414243为四维数组索弓|。
[0236] ⑥设向量q={q(i) | i = iq(i〇，ii，i2，i3)，in=0，l，2，···，ηη_1，η = 〇，1，2,3}是一个 四维数组的向量化表示，其中山(：1()4142 43) = ((：[()*111+；[1)*112+12)*113+13，：[为向量索 引，1041424 3为四维数组索引。按照如下方式对向量9的各个元素赋值，9(1<1(1()414 2， i3)) = q' (iq，（io+no, ii+m，i2+n2, i3+n3))，其中0 仝 in 仝 Πη-1，η = 0，1，2,3。
[0237] ⑦整个方法流程中，以上第①~③步只需要计算一次，运算量基本可以忽略，因 此，如第⑤步所述，矩阵与向量乘积运算q = Ap的计算量主要包括一次四维快速傅里叶变换 和一次四维快速傅里叶逆变换。
[0238] (7)、对yiu执行函数=，并存储L，其中祕={0,口，…,
[0239] 步骤（7)所述的函数5' = AT(y)表示对四维数组y的元素位置进行调整，并返回一 个四维数组f，其中，
[0240] y = {y(i) | i = iy(i〇，ii，i2，i3)，in = 0,1，2，···，Πη-1，n = 〇, 1，2,3}是一个四维数组 的向量化表示，1丫(1〇414243) = ((：[。*111+；[1)*112+12)*113+13，：[为向量索引，：[。4142 43 为四维数组索引；
[0241] $叫? ".Ψ. =Λ,丨,2,…,，7;/ -丨 w 表示，& (I。,I;, 4，匕）=((4 _ ?1 + 0 _~+ 4 ^ % Η，i为向量索引，i0，i 1，i2，i3为四维数组索引。
[0242] 具体调整方法如下，(4 44)) = y(&(Μι，?·2,4))，其中 I,圈d(/,, - floor((nn-l)/2)，in= {〇，1,2,···，nn-l}，n = 〇, 1,2,3,函数mod(vi，V2)表示vi对V2求0~V2-1 之间的余数，函数floorO表示向下取整。
[0243] (8)、对|祕执行函数％ =liF；T4(5^),并存储Siu，其中 i ω = {〇，1，2，…，i ω max} 〇
[0244] 步骤(8)所述的函数s = /FF7；(5〇农示对四维数组#做四维快速傅里叶逆变换，并 返回一个四维数组s。其中，
[0245] 夕=以0|卜4私44:，4),(=〇，1,2,"_，^；-切=〇,1,2,3}是一个四维数组的向量化 表示，& (/W, ) = /v+乂）_/V-6 WvH ^
[0246] s = {s(i) | i = is(io，ii，i2，i3)，in = 0,1，2，···，Πη-1，n = 〇, 1,2,3}是一个四维数组 的向量化表示，15(：[。414243) = ((：[。*111+；[1)*112+12)*113+13，：[为向量索引，：[。4142 43 为四维数组索引。
[0247] (9)、利用Slu生成插值后的频率域五维数据其中，1为插值后的空间 采样坐标，ix = {〇, 1，2,…，η〇ηιη2η3_1}，iw = {〇，l，2，···，Nt_l}。
[0248] 步骤（9)所述的 t = (.r,,(/,)),5(/1),1 办),4(/凡，其中，夂以 {0，1，2,…，1^-1}，11=〇，1，2,3，觸,的定义见步骤(3)。并且11与1()、11、12、13保持关系式11 = ((io · m+ii) · n2+i2) · n3+i3。
[0249] 步骤(9)所述的?^η.,%)生成方式如下，为了方便表述，用向量t代表, 的一部分数据，具体地，l = 以仏,l)3(仁，l),.··ΜΚΛηΛ-ν ω>?〇))^\
[0251 ] 其中，函数con j ()表示对输入向量的每一个元素求复共辄。0表示由11()111112113个0元 素构成的向量。
[0252] (10)、利用快速傅里叶逆变换将Ai,、，％.,)变换的时间域，得到插值后的五维数据
把,、At )。
[0253] 步骤（10)所述的乂,）与, &)满足如下关系式，
,tit = it· At，it={0，l，2，···，Nf 1丨，ix= {0,1，2,…，η〇ηιη2Π3-1}。
[0254] 图2是本发明实施例中插值前的地震数据观测系统示意图；图3是本发明实施例中 插值前的同一个CMP面元的地震数据观测系统示意图；图4是本发明实施例中与图3对应的 插值后的的同一个CMP面元的地震数据观测系统示意图；图5是本发明实施例中采用不同快 速算法的运算时间对比图表示意图；图6是本发明实施例中与图3对应的插值前的CMP道集 地震数据示意图；图7是本发明实施例中与图4对应的插值后的CMP道集地震数据。
[0255] 与现有技术相比较，通过图2至图7所示和上述实施例的介绍可知，本发明实施例 提供的不规则采集地震数据的五维插值快速算法，核心是提出了一种新的迭代方案，大大 提高了计算效率。整个计算过程仅涉及少量NFFT运算，绝大部分的矩阵与向量乘积运算都 可以利用FF