CN104035128A - 可控震源伪随机扫描信号生成方法 - Google Patents

Abstract

本发明是地球物理勘探中地震采集可控震源伪随机扫描信号生成方法,根据所要生成的Num组可控震源伪随机扫描信号长度生成Gold序列码，得到目标振幅谱,在时间域中对Num组振幅谱整形信号的互相关和自相关进行加窗处理得到压制后的Num组重构信号,使用相同时间点位置上各Num组高切滤波信号的之间的互相关最大值除以每组信号的零时刻自相关最小值得到Num组高切滤波信号的最大相对互相关值；迭代修改生成最终的Num组可控震源伪随机扫描信号,本发明能够使得各组伪随机信号间的互相关效弱，从而抑制高效采集过程中的邻炮干扰，改进炮记录的分离质量，瞬时频率能量并不随时间呈现单调变化，从而能够大大降低可控震源施工过程中引起周围建筑共振的风险。

Description

可控震源伪随机扫描信号生成方法

技术领域

本发明涉及地球物理勘探中地震采集的可控震源技术，是一种可控震源伪随机扫描信号生成方法。

背景技术

可控震源技术就是在陆地地震勘探中利用检波器接收由可控震源所激发出的地震波来获得地下构造信息。鉴于可控震源与炸药震源相比，具有安全环保、可重复使用、激发信号已知且可以控制等诸多优点，使其在陆上地震勘探中得到越来越多的应用。

常规的线性扫描信号目前仍然是可控震源勘探中使用最多的扫描方式。线性扫描信号有很多优点，如能量平均分配、单位频率所占有的能量相同、对扫描频率范围内的信号成分没有频率滤波作用等。然而，这种线性扫描信号也有其不足之处，由于这种信号的频率随时间曾线性单调变化趋势，使得信号中对应于每个时刻的单一频率成分的能量都较大。从而在进行震源采集时，会存在引起周围建筑共振的风险。并且常规震源系统在激发扫描信号时，通常会对信号低频成分的驱动幅度有一定程度地限制，这使得实际激发出的线性扫描信号往往不具备足够的低频能量。对此，目前已经出现了能够满足震源系统限制的低频扫描信号设计技术。

如针对线性扫描信号的自相关旁瓣严重问题，国内外有学者通过利用m序列码来控制一个单位周期长度、频率固定的正弦函数信号相位，进而生成伪随机可控震源扫描信号的方法。但这种伪随机信号的自相关旁瓣仅在主瓣峰值附近位置得到压制，对于远处位置上，其旁瓣比常规扫描信号要明显得多，这对于层位较多、反射时间较长的实际地下介质而言，反而会更加影响相关记录的质量。此外，这种伪随机信号的振幅谱与常规线性扫描信号存在较大差异，振幅谱的低频成分有较大缺失。

随着人们对高精度地震勘探采集技术的日益重视，许多高效可控震源采集技术也不断涌现出来。其中，独立同步扫描技术（ISS）作为一种十分高效的可控震源采集技术正被日益推广开来。然而，该方法直接采集到的是多炮混叠的未相关记录。当使用常规线性扫描信号作为震源信号进行激发时，使用扫描信号与振动信号进行互相关所得到的相关炮集记录会存在较强的邻炮干扰，因而无法对其进行现场质控，并对后续地震资料处理中的邻炮干扰压制工作提出了较高的要求。而使用上述直接利用m序列生成伪随机可控震源扫描信号的方法，同样无法生成若干组具有弱相关特性的扫描序列，因而仍然无法直接抑制高效采集所带来的邻炮干扰。

发明内容

本发明目的在于提供一种共振风险小、低频信息丰富、扫描信号之间的互相关能量弱且振幅谱与传统线性扫描信号相近的可控震源伪随机扫描信号生成方法。

本发明通过以下技术步骤实现：

1）根据所要生成的Num组可控震源伪随机扫描信号长度，在二元线性反馈移位寄存器上生成Gold序列码，任意选出Num组Gold序列码供生成后续的可控震源伪随机扫描信号；

所述的Num组可控震源伪随机扫描信号各组信号长度相等。

所述的Gold序列码生成是：

首先获取阶数为n的全部本原多项式，再将本原多项式的各项系数设置为二元线性反馈移位寄存器开关，利用此二元线性反馈移位寄存器即可生成与各组本原多项式对应的m序列；此后，从众多m序列中，确定一组互相关值最小的m序列优选对，并对此m序列优选对进行模2和运算，可生成一组Gold序列；每当改变两组m序列相对位置时，即生成一组新Gold序列。

当本原多项式阶数为n时，生成伪随机扫描信号的采样点长度为2ⁿ-1，此时采用相对位移方式可以得到（2ⁿ-1）个Gold序列。再加上两组原有的m序列，共可以生成（2ⁿ+1）个Gold序列码。根据所需要生成伪随机扫描信号的数量，即可从所生成的众多Gold序列码中，任意挑选出Num组，生成后续的可控震源伪随机扫描信号。

2）根据所要生成的Num组可控震源伪随机扫描信号频率域的采样间隔Δf，以及最大频率f_max，利用下式得到目标振幅谱：

Target(f)=g₁(f)g₂(f) （1）

此处，

$g_1\left(f\right)={\left(\frac{f}{F_1}\right)}^2[{(1-{\left(\frac{f}{F_1}\right)}^2)}^2+{\left(\sqrt{2}\left(\frac{f}{F_1}\right)\right)}^2{]}^{-0.5}---\left(2\right)$

并且

当F₂<f≤F₃时， $g_2\left(f\right)={\left(\frac{1+\cos \left(\left(\frac{f-F_2}{F_3-F_2}\right)\mathrm{\&pi;}\right)}{2}\right)}^2---\left(3\right)$

当f≤F₂时，g₂(f)=1 （4）

当f>F₃时，g₂(f)=0 （5）

公式（1）至（5）中，f表示频率采样位置，其间隔为Δf，单位为Hz，g₁(f)为低切滤波器，g₂(f)为高切滤波器，F₁表示低切频率值的1/2，F₂表示高切频率值，F₃表示高频振幅为零值的频率值，且F₃<f_max；

3）应用目标振幅谱与匹配滤波器对Num组Gold序列码进行振幅谱整形，生成Num组振幅谱整形信号；

步骤3）振幅谱整形处理是一个循环迭代的过程，当第一次执行步骤3）时，处理的信号为Gold序列码，当第二次及以后执行步骤3）时，此步处理的信号为步骤10）生成的Num组高切滤波信号。

4）对得到的Num组振幅谱整形信号进行相关运算；

步骤4）中的相关运算包括各组扫描信号间的互相关运算和每组信号的自相关运算，运算在频率域中实现。

步骤4）中的相关运算是首先对两组数据进行快速傅里叶变换，得到其频域的复数值，之后在频域中，使用第一组数据的复数样点值乘以另一组数据复数样点值的共轭；此后，再对得到的频域相关值进行快速傅里叶反变换，得到时域的相关值。

5）使用相同时间点位置上各Num组振幅谱整形信号之间的互相关最大值除以每组信号的零时刻自相关最小值得到Num组振幅谱整形信号的最大相对互相关值Rmax_i ^a；

此步骤得到的最大相对互相关值大小能够说明振幅谱整形信号之间的互相关能量强弱，进而为最终生成弱互相关能量的可控震源伪随机扫描信号提供定量评价依据。

6）在时间域中，对Num组振幅谱整形信号的互相关和自相关进行加窗处理，得到加窗后的相关序列；

步骤6）的加窗处理使用的时窗函数如下：

当RL<t≤RL+T时，Window(t)=1/2*(1-cos(π*（RL+T-t)/T)) （6）

当-(RL+T)≤t<-RL时，Window(t)=1/2*(1-cos(π*(RL+T+t)/T)) （7）

当t取其它值时，Window(t)=1 （8）

公式中，T为过渡时间长度，RL为听时间长度，t为时窗函数的时间值；

将公式（6）至(8)所示的时窗函数分别与各组扫描信号的互相关和自相关进行相乘，完成加窗处理。

7）利用公式（9）和（10）对Num组振幅谱整形信号进行重构，得到互相关压制后的Num组重构信号；

对于重构获得的任一组扫描信号可用如下公式获得其近似解：

$F\left\{\stackrel{\&OverBar;}{S_i}\right\}\&ap;F\left\{S_i\right\}-\mathrm{\&beta;}*\left({\mathrm{\&delta;}}_i\right)---\left(9\right)$

公式（9）中，

${\mathrm{\&delta;}}_i\&ap;F\left\{R(S_i,\stackrel{\&OverBar;}{S_1})\right\}/F\left\{R(\stackrel{\&OverBar;}{S_1},\stackrel{\&OverBar;}{S_1})\right\}*F\left\{\stackrel{\&OverBar;}{S_1}\right\}+$

$F\left\{R(S_i,\stackrel{\&OverBar;}{S_2})\right\}/F\left\{R(\stackrel{\&OverBar;}{S_2},\stackrel{\&OverBar;}{S_2})\right\}*F\left\{\stackrel{\&OverBar;}{S_2}\right\}+$

$............\left(10\right)$

$+F\left\{R(S_i,\stackrel{\&OverBar;}{S_{\mathrm{Num}-1}})\right\}/F\left\{R(\stackrel{\&OverBar;}{S_{\mathrm{Num}-1}},\stackrel{\&OverBar;}{S_{\mathrm{Num}-1}})\right\}*F\left\{\stackrel{\&OverBar;}{S_{\mathrm{Num}-1}}\right\}$

$+F\left\{R(S_i,\stackrel{\&OverBar;}{S_{\mathrm{Num}}})\right\}/F\left\{R(\stackrel{\&OverBar;}{S_{\mathrm{Num}}},\stackrel{\&OverBar;}{S_{\mathrm{Num}}})\right\}*F\left\{\stackrel{\&OverBar;}{S_{\mathrm{Num}}}\right\}$

式中，F{}表示对括号{}内的变量进行频率域运算；函数R(A,B)表示对信号A与信号B进行相关运算，且对相关运算结果进行加窗处理；S_i(i=1,2,……,Num)表示重构前的振幅谱整形信号；表示重构所获得的信号；β表示一组常量，根据其扫描间隔Δβ，对最小值β_min与最大值β_max之间的取值范围进行扫描搜索；

所述的扫描间隔Δβ取值为0.1或0.2。

所述的β取值为0.1～2.0之间。

8）应用目标振幅谱对Num组重构信号进行振幅谱整形，过程与步骤3）相同；

9）采用以下计算对Num组重构整形信号进行能级加强，得到Num组能级加强信号；

当-1<y<1时，Compression(y)=sin{y*(π/2)} （11）

当y取其它值时，Compression(y)=sign(y) （12）

式中，Compression(y)表示对数值y进行能级加强处理后所获得的信号，sign(y)表示取数值y的符号；y利用公式y=η*(x/σ)计算得到，该式中的x表示伪随机扫描信号在各采样时间位置上的样点值，σ表示伪随机扫描信号的均方根，η是常数；

所述的η常数为0.1-0.5。

10）对Num组能级加强信号进行高切滤波，得到Num组高切滤波信号；

所述的高切滤波公式利用（3）至（5）高切滤波函数。

11）对Num组高切滤波信号进行相关运算，处理过程与步骤4）相同；

12）使用相同时间点位置上各Num组高切滤波信号的之间的互相关最大值除以每组信号的零时刻自相关最小值得到Num组高切滤波信号的最大相对互相关值Rmax_i ^b；

此步骤得到的最大相对互相关值大小能够说明各组高切滤波信号之间的互相关能量强弱，进而为最终生成弱互相关能量的可控震源伪随机扫描信号提供定量评价依据。通过将步骤12）得到的最大相对互相关值与步骤5）所得到的最大相对互相关值进行比较，能够为步骤13）的执行提供依据。

13）对Num组高切滤波信号进行迭代修改，生成最终的Num组可控震源伪随机扫描信号；

所述的迭代修改是将步骤12）计算得到的相对互相关值Rmax_i ^b与步骤5）计算得到的相对互相关值Rmax_i ^a进行比较，当Rmax_i ^b<Rmax_i ^a时，并且正、负听时间范围内的最大相对互相关值低于门槛值，则将步骤10）中得到的Num组高切滤波信号重新输入至步骤3），进入下一次循环迭代的信号修改处理，若低于门槛值，则跳出循环，生成最终的Num组可控震源伪随机扫描信号；当Rmax_i ^b≥Rmax_i ^a时，再确定步骤7）中的β常量是否按照扫描间隔Δβ取过介于β_min至β_max之间的所有值，若否，则利用Δβ对β进行修改，并将步骤10）中得到的Num组高切滤波信号重新输入至步骤3），进行下一次循环迭代的信号修改处理，若是，则说明信号间的互相关值已经无法进行进一步的压制，跳出循环，即可生成最终的Num组可控震源伪随机扫描信号。

所述的门槛值为0.01至0.1。

所述的门槛值最佳为0.025。

14）将求得的扫描信号以Pelton和Sercel可控震源电控系统文本方式到导入相应的可控震源电控系统中，为可控震源地震勘探提供一种共振风险小、低频信息丰富、扫描信号之间的互相关能量弱且振幅谱与传统线性扫描信号相近的可控震源伪随机扫描信号，进而提高可控震源勘探质量。

附图说明

图1是二元线性反馈移位寄存器示意图；

图2是Gold序列族中用于生成可控震源伪随扫描信号的四组Gold序列码的部分截断；

(a)第1组Gold序列码；

(b)第2组Gold序列码；

(c)第3组Gold序列码；

(d)第4组Gold序列码；

图3为用于对初始Gold序列进行振幅谱整形的目标振幅谱；

图4为通过使用匹配滤波器进行过振幅谱整形后的四组振幅谱整形信号的部分截断；

(a)进行振幅谱整形后的第1组振幅谱整形信号；

(b)进行振幅谱整形后的第2组振幅谱整形信号；

(c)进行振幅谱整形后的第3组振幅谱整形信号；

(d)进行振幅谱整形后的第4组振幅谱整形信号；

图5为振幅谱整形信号的自相关与互相关；

(a)图4(a)所示振幅谱整形信号的正时刻自相关；

(b)图4(a)所示振幅谱整形信号的负时刻自相关；

(c)图4(a)与图4(b)所示振幅谱整形信号的正时刻互相关；

(d)图4(a)与图4(b)所示振幅谱整形信号的负时刻互相关；

图6为对Gold序列码进行振幅谱整形处理后所生成的振幅谱整形信号之间的最大相对互相关值；

(a)各振幅谱整形信号在正时间轴上的最大相对互相关值；

(b)各振幅谱整形信号在负时间轴上的最大相对互相关值；

图7为用于对各组振幅谱整形信号的互相关和自相关进行加窗处理的时窗函数；

(a)正时刻的部分时窗；

(b)负时刻的部分时窗；

图8为采用时窗函数进行加窗处理后的振幅谱整形信号自相关与互相关；

(a)图5(a)所示自相关序列的加窗处理结果；

(b)图5(b)所示自相关序列的加窗处理结果；

(c)图5(c)所示互相关序列的加窗处理结果；

(d)图5(d)所示互相关序列的加窗处理结果；

图9为进行迭代修改过的高切滤波信号之间的最大相对互相关值；

(a)各迭代修改过的高切滤波信号在正时间轴上的最大相对互相关值；

(b)各迭代修改过的高切滤波信号在负时间轴上的最大相对互相关值；

图10为四组最终生成的可控震源伪随机扫描信号的部分序列；

(a)进行迭代修改后生成的第1组可控震源伪随机扫描信号；

(b)进行迭代修改后生成的第2组可控震源伪随机扫描信号；

(c)进行迭代修改后生成的第3组可控震源伪随机扫描信号；

(d)进行迭代修改后生成的第4组可控震源伪随机扫描信号；

图11为最终生成的其中一组可控震源伪随机扫描信号的振幅分贝谱；

图12为可控震源电控系统：Sercel系统（左）和Pelton系统（右）分别对应的可控震源信号文本格式；

图13为对本发明方法进行归纳,生成可控震源伪随机扫描信号流程图。

具体实施方式

以下结合附图进行详细说明。

本发明以能够生成m序列的本原多项式为出发点，通过选取m序列优选对，进而将生成的Gold序列族中的Gold序列码作为最终所要生成的可控震源伪随机扫描信号的初始原型。并通过采用循环迭代的方式，其中包括：振幅谱整形、互相关压制以及能级加强等手段来实现对初始Gold序列码的修改。从而得到若干组共振风险小、低频信息丰富、扫描信号之间的互相关能量弱且振幅谱与传统线性扫描信号相近的可控震源伪随机扫描信号。

本发明具体实施方案及步骤详述如下：

1）根据所要生成的Num组可控震源伪随机扫描信号长度，在二元线性反馈移位寄存器上生成Gold序列码，任意选出Num组Gold序列码供生成后续的可控震源伪随机扫描信号；

所述的Num组可控震源伪随机扫描信号各组信号长度相等。

所述的Gold序列码生成是：

首先获取阶数为n的全部本原多项式，再将本原多项式的各项系数设置为二元线性反馈移位寄存器开关，利用此二元线性反馈移位寄存器即可生成与各组本原多项式对应的m序列，图1为二元线性反馈移位寄存器的示意图，利用此移存器可以根据每组本原多项式分别产生一组m序列；此后，从众多m序列中，确定一组互相关值最小的m序列优选对（对于线性反馈移位寄存器所激发出来的各组m序列，当本原多项式的阶数大于或等于12时，一定存在着一组m序列优选对，它们的互相关值小于其它任意一组m序列优选对互相关值）；之后再对得到的m序列优选对进行模2和运算，可生成一组Gold序列；每当改变两组m序列相对位置时，即生成一组新Gold序列。

当本原多项式阶数为n时，生成伪随机扫描信号的采样点长度为2ⁿ-1，此时采用相对位移方式可以得到（2ⁿ-1）个Gold序列。再加上两组原有的m序列，共可以生成（2ⁿ+1）个Gold序列码。根据所需要生成伪随机扫描信号的数量，即可从所生成的众多Gold序列码中，任意挑选出Num组，生成后续的可控震源伪随机扫描信号。（如图2所示，图中Num=4，扫描信号采用间隔Δt=2ms，生成Gold序列码的阶数为n=15）。通过后续的处理步骤做进一步修改，最终将生成Num组可控震源伪随机扫描信号。

2）根据所要生成的Num组可控震源伪随机扫描信号频率域的采样间隔Δf，以及最大频率f_max，利用下式得到目标振幅谱：

Target(f)=g₁(f)g₂(f) （1）

此处，

$g_1\left(f\right)={\left(\frac{f}{F_1}\right)}^2[{(1-{\left(\frac{f}{F_1}\right)}^2)}^2+{\left(\sqrt{2}\left(\frac{f}{F_1}\right)\right)}^2{]}^{-0.5}---\left(2\right)$

并且

当F₂<f≤F₃时， $g_2\left(f\right)={\left(\frac{1+\cos \left(\left(\frac{f-F_2}{F_3-F_2}\right)\mathrm{\&pi;}\right)}{2}\right)}^2---\left(3\right)$

当f≤F₂时，g₂(f)=1 （4）

当f>F₃时，g₂(f)=0 （5）

公式（1）至（5）中，f表示频率采样位置，其间隔为Δf，单位为Hz，g₁(f)为低切滤波器，g₂(f)为高切滤波器，F₁表示低切频率值的1/2，F₂表示高切频率值，F₃表示高频振幅为零值的频率值，且F₃<f_max；

图3为采用公式（1）至（5）所生成的目标振幅谱，图中f_max=250Hz、F₁=2Hz、F₂=100Hz、F₃=150Hz。

3）应用目标振幅谱与匹配滤波器对Num组Gold序列码进行振幅谱整形，生成Num组振幅谱整形信号；

步骤3）振幅谱整形处理是一个循环迭代的过程，当第一次执行步骤3）时，处理的信号为Gold序列码，当第二次及以后执行步骤3）时，此步处理的信号为步骤10）生成的Num组高切滤波信号。

匹配滤波器作为地震信号处理中一种将原始数据信号匹配为目标信号的常用工具，可以将步骤1）中所有生成的Num组Gold序列码的振幅谱匹配成与目标振幅谱一致，从而得到进行过振幅谱整形后的Num组振幅谱整形信号（如图4所示，图中Num=4，信号采样间隔Δt=2ms）。

4）对得到的Num组振幅谱整形信号进行相关运算；

步骤4）中的相关运算包括各组扫描信号间的互相关运算和每组信号的自相关运算，运算在频率域中实现。

步骤4）中的相关运算是首先对两组数据进行快速傅里叶变换，得到其频域的复数值，之后在频域中，使用第一组数据的复数样点值乘以另一组数据复数样点值的共轭；此后，再对得到的频域相关值进行快速傅里叶反变换，得到时域的相关值（如图5所示，取10秒长度）。需要说明的是，在频率域进行相关后所得到的数据结果，其前半个数据长度表示正时刻相关值，其后半个数据长度则表示负时刻相关值。

5）使用相同时间点位置上各Num组振幅谱整形信号之间的互相关最大值除以每组信号的零时刻自相关最小值得到Num组振幅谱整形信号的最大相对互相关值Rmax_i ^a；

在步骤5）中，本发明计算得到的最大相对互相关值用Rmax_i ^a表示（如图6所示，其中图6(a)为正时间范围的数值显示，时间范围是0秒至10秒，图6(b)为负时间范围的数值显示，时间范围是-10秒至0秒）。对于Rmax_i ^a的每个时间点位置上的数值，用相同时间点位置上各组扫描信号之间的互相关最大值除以每组信号的零时刻自相关最小值得到。

此步骤得到的最大相对互相关值大小能够说明振幅谱整形信号之间的互相关能量强弱，进而为最终生成弱互相关能量的可控震源伪随机扫描信号提供定量评价依据。

6）在时间域中，对Num组振幅谱整形信号的互相关和自相关进行加窗处理，得到加窗后的相关序列；

加窗处理的作用是使得从负听时刻到正听时刻范围内的数据保持不变，而在此时间范围外的数据逐渐向0值过渡。

步骤6）的加窗处理使用的时窗函数如下：

当RL<t≤RL+T时，Window(t)=1/2*(1-cos(π*（RL+T-t)/T)) （6）

当-(RL+T)≤t<-RL时，Window(t)=1/2*(1-cos(π*(RL+T+t)/T)) （7）

当t取其它值时，Window(t)=1 （8）

公式中，Window(t)为加窗处理结果，T为过渡时间长度，RL为听时间长度，t为时窗函数的时间值。

将公式（6）至(8)所示的时窗函数（如图7所示，图中RL=8秒，T=2秒，图中的前半轴表示正时刻，后半轴表示负时刻）分别与各组扫描信号的互相关和自相关进行相乘，完成加窗处理，得到加窗后的相关序列（如图8所示，取10秒长度）。

7）利用公式（9）和（10）对Num组振幅谱整形信号进行重构，得到互相关压制后的Num组重构信号；

对于重构获得的任一组扫描信号可用如下公式获得其近似解：

$F\left\{\stackrel{\&OverBar;}{S_i}\right\}\&ap;F\left\{S_i\right\}-\mathrm{\&beta;}*\left({\mathrm{\&delta;}}_i\right)---\left(9\right)$

公式（9）中，

${\mathrm{\&delta;}}_i\&ap;F\left\{R(S_i,\stackrel{\&OverBar;}{S_1})\right\}/F\left\{R(\stackrel{\&OverBar;}{S_1},\stackrel{\&OverBar;}{S_1})\right\}*F\left\{\stackrel{\&OverBar;}{S_1}\right\}+$

$F\left\{R(S_i,\stackrel{\&OverBar;}{S_2})\right\}/F\left\{R(\stackrel{\&OverBar;}{S_2},\stackrel{\&OverBar;}{S_2})\right\}*F\left\{\stackrel{\&OverBar;}{S_2}\right\}+$

$.........\left(10\right)$

$+F\left\{R(S_i,\stackrel{\&OverBar;}{S_{\mathrm{Num}-1}})\right\}/F\left\{R(\stackrel{\&OverBar;}{S_{\mathrm{Num}-1}},\stackrel{\&OverBar;}{S_{\mathrm{Num}-1}})\right\}*F\left\{\stackrel{\&OverBar;}{S_{\mathrm{Num}-1}}\right\}$

$+F\left\{R(S_i,\stackrel{\&OverBar;}{S_{\mathrm{Num}}})\right\}/F\left\{R(\stackrel{\&OverBar;}{S_{\mathrm{Num}}},\stackrel{\&OverBar;}{S_{\mathrm{Num}}})\right\}*F\left\{\stackrel{\&OverBar;}{S_{\mathrm{Num}}}\right\}$

式中，F{}表示对括号{}内的变量进行频率域运算；函数R(A,B)表示对信号A与信号B进行相关运算，且对相关运算结果进行加窗处理；S_i(i=1,2,……,Num)表示重构前的振幅谱整形信号；(i=1,2,……,Num)表示进行互相关压制后的重构信号；β表示一组常量，根据其扫描间隔Δβ，对最小值β_min与最大值β_max之间的取值范围进行扫描搜索；

本发明在首次迭代中，首先根据一初始值β_ini，并按照后续步骤13）来确定β取值是否有效，从而确定是否需要对β进行修改并进入下一次迭代。

所述的扫描间隔Δβ取值为0.1或0.2。

所述的β取值为0.1～2.0之间。

根据公式（9）与（10）即可得到Num组互相关压制后的重构扫描信号。

8）应用目标振幅谱对Num组重构信号进行振幅谱整形，过程与步骤3）相同；

此步骤的目的是为了消除由步骤10）所造成的振幅谱改变。经步骤8）处理后，由步骤7）所得的Num组重构信号的振幅谱得到恢复。

9）采用以下计算对Num组重构整形信号进行能级加强，得到Num组能级加强信号；

步骤9)的目的是为了尽可能提高可控震源的出力，以使得平板与地面能够紧密地接触，从而在使用伪随机扫描信号进行地震勘探时，能够提高震源激发能量，得到信噪比稳定的地震资料。实际上，步骤9）中的信号重构过程会降低信号的均方根能量，而步骤9）将会对信号的均方根能量进行补偿，从而提高信号的激发能量。

当-1<y<1时，Compression(y)=sin{y*(π/2)} （11）

当y取其它值时，Compression(y)=sign(y) （12）

式中，Compression(y)表示对数值y进行能级加强处理后所获得的信号，sign(y)表示取数值y的符号；y利用公式y=η*(x/σ)计算得到，该式中的x表示伪随机扫描信号在各采样时间位置上的样点值，σ表示伪随机扫描信号的均方根，η是常数，其取值范围为0.1-0.5；根据公式（11）与（12）即可得到Num组能级加强信号。

10）对Num组能级加强信号进行高切滤波，得到Num组高切滤波信号；

由于步骤9）是一组非线性过程，因此在经过此步骤处理后，会产生扫描信号目标谱范围以外的高频噪音。所述的高切滤波公式利用（3）至（5）高切滤波函数。应用此函数对上一步生成的信号进行滤波处理后，即可得到Num组高切滤波信号。

11）对Num组高切滤波信号进行相关运算，处理过程与步骤4）相同；

12）使用相同时间点位置上各Num组高切滤波信号的之间的互相关最大值除以每组信号的零时刻自相关最小值得到Num组高切滤波信号的最大相对互相关值Rmax_i ^b；

在步骤12）中，本发明将计算得到的最大相对互相关值用Rmax_i ^b表示（如图9所示，其中图9(a)为正时间范围的数值显示，时间范围是0秒至10秒，图9(b)为负时间范围的数值显示，时间范围是-10秒至0秒）。

此步骤得到的最大相对互相关值的大小能够说明各组高切滤波信号之间的互相关能量强弱，进而为最终生成弱互相关能量的可控震源伪随机扫描信号提供定量评价依据。通过将步骤12）得到的最大相对互相关值与步骤5）所得到的最大相对互相关值进行比较，能够为步骤13）的执行提供依据。

13）对Num组高切滤波信号进行迭代修改，生成最终的Num组可控震源伪随机扫描信号；

所述的迭代修改是将步骤12）计算得到的相对互相关值Rmax_i ^b与步骤5）计算得到的相对互相关值Rmax_i ^a进行比较，当Rmax_i ^b<Rmax_i ^a时，并且正、负听时间范围内的最大相对互相关值低于门槛值，则将步骤10）中得到的Num组高切滤波信号重新输入至步骤3），进入下一次循环迭代的信号修改处理，若已低于门槛值，则跳出循环，生成最终的Num组可控震源伪随机扫描信号；当Rmax_i ^b≥Rmax_i ^a时，再确定步骤7）中的β常量是否按照扫描间隔Δβ取过介于β_min至β_max之间的所有值，若否，则利用Δβ对β进行修改，并将步骤10）中得到的Num组高切滤波信号重新输入至步骤3），进行下一次循环迭代的信号修改处理，若是，则说明信号间的互相关值已经无法进行进一步的压制，跳出循环，即可生成最终的Num组可控震源伪随机扫描信号（如图10所示，图中Num=4，扫描信号采用间隔Δt=2ms）。

所述的门槛值为0.01至0.1。

所述的门槛值最佳为0.025。

14）将求得的扫描信号以Pelton和Sercel可控震源电控系统文本方式到导入相应的可控震源电控系统中，为可控震源地震勘探提供一种共振风险小、低频信息丰富、扫描信号之间的互相关能量弱且振幅谱与传统线性扫描信号相近的可控震源伪随机扫描信号，进而提高可控震源勘探质量。

与Sercel系统对应的可控震源信号文本格式如图12（左）所示，与Pelton系统对应的可控震源信号文本格式如图12（右）所示。通过将本文导入至相应的可控震源电控系统中，利用电控系统对与之相对应的可控震源进行控制，即可激发出相应的可控震源扫描信号，从而进行可控震源伪随机采集，获取得地震记录。

本发明的整个操作流程如图13所示。

为验证本发明的可行性，采用下面的已知实例对本发明进行测试，详述如下：

实例

实例中设置生成可控震源伪随机扫描信号的阶数为15，生成Gold序列码的样点采样间隔为2ms，并要求同时生成4组可控震源伪随机扫描信号。图2为利用m序列优选对生成的4组Gold序列码，本实例即以此4组Gold序列码作为可控震源伪随机扫描信号的初始原型，用于生成最终的可控震源伪随机扫描信号。本方法在对4组Gold序列码进行修改时，构造了如图3所示的振幅谱，作为对4组Gold序列码进行振幅谱整形的目标振幅谱。从图3中可见，振幅值在4Hz至100Hz的范围内是水平的，且在2Hz位置上的振幅值为4Hz位置上振幅值的一半，这使得生成的伪随机扫描信号能够具有较为丰富的低频信息。图4为通过使用匹配滤波器进行振幅谱整形后的4组振幅谱整形信号。由于4组振幅谱整形信号是由具有弱互相关特性的Gold序列码得到，因此对于此处得到的4组振幅谱整形信号而言，其自相关旁瓣能量和互相关能量相对于扫描信号的零时刻自相关而言仍然较弱。图5所展示的是在正、负10秒范围内，第一组振幅谱整形信号的自相关，以及第一组振幅谱整形信号与第二组振幅谱整形信号的互相关。从图5(a)可以看出，振幅谱整形信号的自相关旁瓣很弱，零时刻峰值很强，将图5(b)所展示的互相关值与图5(a)所展示的自相关峰值对比可以看出，振幅谱整形信号的互相关值相对于零相刻自相关值而言，能量同样很弱，其中，零时刻自相关的峰值已接近0.4，而互相关值仅在0.01附近，二者能量差别较大。图6则更加直接地展示了各组振幅谱整形信号间的最大相对互相关值（即4组信号的互相关最大值除以信号的零时刻自相关最小值所得到的比值）。其物理意义为，每一采样点时刻位置上，4组振幅谱整形信号的互相关最大值与零相刻自相关的最小值之比。由图6可以看出，4组振幅谱整形信号的互相关与零时刻自相关的比值仍然较小。实际上，4组振幅谱整形信号间的互相关值还可以做进一步地压制，以使得在采用伪随机扫描信号进行施工时，邻炮干扰能够被进一步削弱。由于只有在正、负听时间范围内的互相关值会对采集得到的相关地震记录产生干扰，故只需针对处于正、负听时间范围内的信号互相关值进行压制。因此，在进行真正的互相关压制之前，首先需要使用如图7所示的时窗函数，对图5所示的相关序列进行加窗处理，得到如图8所示的加窗后的相关序列。之后，再通过信号的重构过程来实现信号间的互相关压制。此外，为更好的适应地震资料采集的需求，还需对压制互相关后的信号进行振幅谱整形、能级加强、高切滤波等后续处理。信号修改是一个循环迭代的过程，图10为对Gold序列码进行迭代修改后，最终生成的4组可控震源伪随机扫描信号。图9为4组可控震源伪随机扫描信号间的最大相对互相关值(即4组信号的互相关最大值除以信号的零时刻自相关最小值所得到的比值）。通过将图9与图6进行对比可以看出，在正、负听时间范围内（-8秒至8秒），比值由压制前的0.035附近下降到了压制后的0.01附近，互相关能量得到了大幅度的压制，从而有利于在进行野外采集时对邻炮干扰的削弱。图11为其中一组可控震源伪随机扫描信号的振幅分贝谱，其余3组可控震源伪随机扫描信号的振幅谱也与此组扫描信号接近，从图11可见，生成的可控震源伪随机扫描信号与目标振幅谱的差异不大，与传统的线性扫描信号振幅谱类似。此外，由于可控震源伪随机扫描信号的各频率能量成分并不是随时间呈单调变化趋势，因此大大降低了施工过程中引起周围建筑共振的风险。同时，由图9所展示的4组可控震源伪随机扫描信号间的最大相对互相关值可以看出，各组信号之间的互相关能量非常弱。图12为其中一组可控震源伪随机扫描信号所对应的可控震源Sercel电控系统与Pelton电控系统的文本格式。

Claims (12)

1.一种可控震源伪随机扫描信号生成方法，特征是通过以下技术步骤实现：

1）根据所要生成的Num组可控震源伪随机扫描信号长度，在二元线性反馈移位寄存器上生成Gold序列码，任意选出Num组Gold序列码供生成后续的可控震源伪随机扫描信号；

2）根据所要生成的Num组可控震源伪随机扫描信号频率域的采样间隔Δf，以及最大频率f_max，利用下式得到目标振幅谱：

Target(f)=g₁(f)g₂(f) （1）

此处，

$g_1\left(f\right)={\left(\frac{f}{F_1}\right)}^2[{(1-{\left(\frac{f}{F_1}\right)}^2)}^2+{\left(\sqrt{2}\left(\frac{f}{F_1}\right)\right)}^2{]}^{-0.5}---\left(2\right)$

并且

当F₂<f≤F₃时， $g_2\left(f\right)={\left(\frac{1+\cos \left(\left(\frac{f-F_2}{F_3-F_2}\right)\mathrm{\&pi;}\right)}{2}\right)}^2---\left(3\right)$

当f≤F₂时，g₂(f)=1 （4）

当f>F₃时，g₂(f)=0 （5）

公式（1）至（5）中，f表示频率采样位置，其间隔为Δf，单位为Hz，g₁(f)为低切滤波器，g₂(f)为高切滤波器，F₁表示低切频率值的1/2，F₂表示高切频率值，F₃表示高频振幅为零值的频率值，且F₃<f_max；

3）应用目标振幅谱与匹配滤波器对Num组Gold序列码进行振幅谱整形，生成Num组振幅谱整形信号；

4）对得到的Num组振幅谱整形信号进行相关运算；

5）使用相同时间点位置上各Num组振幅谱整形信号之间的互相关最大值除以每组信号的零时刻自相关最小值得到Num组振幅谱整形信号的最大相对互相关值Rmax_i ^a；

6）在时间域中，对Num组振幅谱整形信号的互相关和自相关进行加窗处理，得到加窗后的相关序列；

7）利用公式（9）和（10）对Num组振幅谱整形信号进行重构，得到互相关压制后的Num组重构信号；

对于重构获得的任一组扫描信号可用如下公式获得其近似解：

$F\left\{\stackrel{\&OverBar;}{S_i}\right\}\&ap;F\left\{S_i\right\}-\mathrm{\&beta;}*\left({\mathrm{\&delta;}}_i\right)---\left(9\right)$

公式（9）中，

${\mathrm{\&delta;}}_i\&ap;F\left\{R(S_i,\stackrel{\&OverBar;}{S_1})\right\}/F\left\{R(\stackrel{\&OverBar;}{S_1},\stackrel{\&OverBar;}{S_1})\right\}*F\left\{\stackrel{\&OverBar;}{S_1}\right\}+$

$F\left\{R(S_i,\stackrel{\&OverBar;}{S_2})\right\}/F\left\{R(\stackrel{\&OverBar;}{S_2},\stackrel{\&OverBar;}{S_2})\right\}*F\left\{\stackrel{\&OverBar;}{S_2}\right\}+$

$..............\left(10\right)$

$+F\left\{R(S_i,\stackrel{\&OverBar;}{S_{\mathrm{Num}-1}})\right\}/F\left\{R(\stackrel{\&OverBar;}{S_{\mathrm{Num}-1}},\stackrel{\&OverBar;}{S_{\mathrm{Num}-1}})\right\}*F\left\{\stackrel{\&OverBar;}{S_{\mathrm{Num}-1}}\right\}$

$+F\left\{R(S_i,\stackrel{\&OverBar;}{S_{\mathrm{Num}}})\right\}/F\left\{R(\stackrel{\&OverBar;}{S_{\mathrm{Num}}},\stackrel{\&OverBar;}{S_{\mathrm{Num}}})\right\}*F\left\{\stackrel{\&OverBar;}{S_{\mathrm{Num}}}\right\}$

式中，F{}表示对括号{}内的变量进行频率域运算；函数R(A,B)表示对信号A与信号B进行相关运算，且对相关运算结果进行加窗处理；S_i(i=1,2,……,Num)表示重构前的振幅谱整形信号；表示重构所获得的信号；β表示一组常量，根据其扫描间隔Δβ，对最小值β_min与最大值β_max之间的取值范围进行扫描搜索；

8）应用目标振幅谱对Num组重构信号进行振幅谱整形，过程与步骤3）相同；

9）采用以下计算对Num组重构整形信号进行能级加强，得到Num组能级加强信号；

当-1<y<1时，Compression(y)=sin{y*(π/2)} （11）

当y取其它值时，Compression(y)=sign(y) （12）

式中，Compression(y)表示对数值y进行能级加强处理后所获得的信号，sign(y)表示取数值y的符号，y利用公式y=η*(x/σ)计算得到，该式中的x表示伪随机扫描信号在各采样时间位置上的样点值，σ表示伪随机扫描信号的均方根，η是常数；

所述的η常数为0.1-0.5。

10）对Num组能级加强信号进行高切滤波，得到Num组高切滤波信号；

11）对Num组高切滤波信号进行相关运算，处理过程与步骤4）相同；

12）使用相同时间点位置上各Num组高切滤波信号的之间的互相关最大值除以每组信号的零时刻自相关最小值得到Num组高切滤波信号的最大相对互相关值Rmax_i ^b；

13）对Num组高切滤波信号进行迭代修改，生成最终的Num组可控震源伪随机扫描信号；

14）将步骤13）求得的扫描信号导入可控震源电控系统，作为可控震源地震勘探伪随机扫描信号进行勘探。

2.根据权利要求1所述的方法，特点是步骤1）所述的Num组可控震源伪随机扫描信号各组信号长度相等。

3.根据权利要求1所述的方法，特点是步骤1）所述的Gold序列码生成是：

首先获取阶数为n的全部本原多项式，再将本原多项式的各项系数设置为二元线性反馈移位寄存器开关，利用此二元线性反馈移位寄存器即可生成与各组本原多项式对应的m序列；此后，从众多m序列中，确定一组互相关值最小的m序列优选对，并对此m序列优选对进行模2和运算，可生成一组Gold序列；每当改变两组m序列相对位置时，即生成一组新Gold序列。

4.根据权利要求1所述的方法，特点是步骤3）振幅谱整形处理是一个循环迭代的过程，当第一次执行步骤3）时，处理的信号为Gold序列码，当第二次及以后执行步骤3）时，此步处理的信号为步骤10）生成的Num组高切滤波信号。

5.根据权利要求1所述的方法，特点是步骤4）中的相关运算包括各组扫描信号间的互相关运算和每组信号的自相关运算，运算在频率域中实现。

6.根据权利要求1所述的方法，特点是步骤4）中的相关运算是首先对两组数据进行快速傅里叶变换，得到其频域的复数值，之后在频域中，使用第一组数据的复数样点值乘以另一组数据复数样点值的共轭；此后，再对得到的频域相关值进行快速傅里叶反变换，得到时域的相关值。

7.根据权利要求1所述的方法，特点是步骤6）的加窗处理使用的时窗函数如下：

当RL<t≤RL+T时，Window(t)=1/2*(1-cos(π*（RL+T-t)/T)) （6）

当-(RL+T)≤t<-RL时，Window(t)=1/2*(1-cos(π*(RL+T+t)/T)) （7）

当t取其它值时，Window(t)=1 （8）

公式中，T为过渡时间长度，RL为听时间长度，t为时窗函数的时间值；

将公式（6）至(8)所示的时窗函数分别与各组扫描信号的互相关和自相关进行相乘，完成加窗处理。

8.根据权利要求1所述的方法，特点是步骤7）所述的扫描间隔Δβ为0.1～0.8之间。

9.根据权利要求1所述的方法，特点是步骤10）所述的高切滤波公式利用（3）至（5）高切滤波函数。

10.根据权利要求1所述的方法，特点是步骤13）所述的迭代修改是将步骤12）计算得到的相对互相关值Rmax_i ^b与步骤5）计算得到的相对互相关值Rmax_i ^a进行比较，当Rmax_i ^b<Rmax_i ^a时，并且正、负听时间范围内的最大相对互相关值低于门槛值，则将步骤10）中得到的Num组高切滤波信号重新输入至步骤3），进入下一次循环迭代的信号修改处理，若已低于门槛值，则跳出循环，生成最终的Num组可控震源伪随机扫描信号；当Rmax_i ^b≥Rmax_i ^a时，再确定步骤7）中的β常量是否按照扫描间隔Δβ取过介于β_min至β_max之间的所有值，若否，则利用Δβ对β进行修改，并将步骤10）中得到的Num组高切滤波信号重新输入至步骤3），进行下一次循环迭代的信号修改处理，若是，则说明信号间的互相关值已经无法进行进一步的压制，跳出循环，即可生成最终的Num组可控震源伪随机扫描信号。

11.根据权利要求10所述的方法，特点是所述的门槛值为0.01至0.1。

12.根据权利要求10所述的方法，特点是所述的门槛值最佳为0.025。