发明内容
为了克服现有技术中的上述和其他缺点,本发明的目的在于提供一种改进的多子波分解和重构方法。
为了实现上述目的,提供了一种地震道多子波分解和重构方法,包括如下步骤:(a)读入地震道数据,利用读入的地震道数据初始化剩余信号;(b)建立基础原子库D,利用建立的基础原子库D,初始化动态原子库AD;(c)计算动态原子库AD中的每个原子与剩余信号的互相关值;(d)针对动态原子库AD中的每个原子,依次与剩余信号匹配出相关值较大,且范围不重叠的多个原子,形成有效原子集,同时生成新的动态原子库AD;(e)根据原子的互相关率大小,从所述有效原子集中提取匹配原子,并利用提取出的匹配原子计算新的剩余信号;(f)利用步骤(e)的结果判断是否达到迭代终止条件,如果没有达到,则重复步骤(c)、(d)、(e)进行迭代,直至满足终止条件,如果达到了迭代终止条件,则进行下面的步骤;(g)保存所提取的原子信息;(h)利用所述原子信息,结合井信息和已知储层信息,进行地震道数据的重构。
其中,在步骤(e)中,通过从剩余信号中依次减去所提取出的匹配原子,来计算新的剩余信号。
其中,在步骤(a)中,初次读入地震道数据时,通过将地震道数据值赋给剩余信号来初始化剩余信号。
步骤(d)包括如下步骤:将每个原子的最大互相关值除以整个原子库中所有原子的互相关值中的最大互相关值并求绝对值,得到该原子的互相关率;判断每个原子的互相关率是否大于设定的第一阈值;利用互相关率大于设定的第一阈值的那些原子形成新的动态原子库AD。
步骤(d)包括如下步骤:针对动态原子库中的每个原子,将互相关值按照降序排序,从大到小选择不超过第一预定数量的时间范围不重叠的原子;将上述步骤中所获得所有原子按照互相关值降序的方式排列,从大到小选择不超过第二预定数量的时间范围不重叠的原子,构成有效原子集。
其中,在步骤(e)中,从所述有效原子集中提取互相关率大于第二阈值的匹配原子,并计算新的剩余信号。
其中,在步骤(f)中,通过计算剩余信号的能量与地震道数据的能量的比值来判断是否达到迭代终止条件。
其中,在步骤(f)中,通过判断提取的匹配原子的个数是否达到设定个数来决定是否终止迭代。
其中,在步骤(f)中利用提取出的匹配原子或新的剩余信号来判断是否达到迭代终止条件,如果没有达到迭代终止条件,则更新动态原子库和剩余信号后,返回步骤(c)。
其中,在步骤(g)中,保存所提取的原子的起始时间、相关系数、频率、相位相关的信息。
其中,在步骤(h)中,基于能量特征选择性重构、基于频率特征选择性重构,或者能量和频率特征综合选择性重构地震道信号。
其中,据频带范围、或者原子能量大小次序、或频率和能量的综合准则来重构地震道。
因此,根据实施例的多子波分解和重构方法可以更快速、准确的对地震道进行分解,且基于频率和能量特征综合筛选的重构方法,增加了能量和频率综合QC的手段,而重构剖面的实时动态QC,有利于用户做出重构参数合理选择的判断。
具体实施方式
在描述本发明的实施例之前,解释一下本发明中提到的一些术语的概念。这些术语在现有技术中已经为本领域技术人员所公知和使用。
原子:自然界中“原子”的概念可大体理解为组成物质的基本单元。在本发明中,原子的概念也大体类似——即组成地震信号的基本单元,可假定地震信号为不同频率、不同振幅、不同时间、不同相位的原子构成,原子也可称为子波,包括有雷克子波、Morlet子波、Gaussian子波等不同类型。
子波:即为原子的意思,根据用户选择的子波类型、主频序列、相位序列,就可构建一个原子库。假设用户选择的子波为雷克子波,它的公式为这样根据不同的主频f就可以得到相应的子波序列。如图7所示,示出了主频分别为10(左上)、20(右上)、30(左下)、40(右下)HZ的雷克子波序列。
以下,参照附图来详细说明本发明的实施例。参照图3和图4,将详细描述根据实施例的多子波分解和重构的获取方法。
首先,以本领域公知的方法进行野外勘探采集,以得到原始采集的数据,然后,经过地震资料处理,得到能进行多子波分解和重构的地震数据。地震数据可以是叠加前的数据或叠加后的数据。水平叠加是针对地震勘探过程中,多个激发点和多个接收点在地下的某个位置有反射,将不同接收点收到的来自地下同一反射点的不同激发点的信号,经动校正后叠加起来。
在根据本发明实施例的多子波分解和重构方法中,可以为整道的地震数据,也可以为部分采样点构成的地震数据。
下面将具体地描述根据示例实施例的多子波分解和重构方法。
图3示出了根据本发明示例性实施例的改进的多子波分解和重构方法。
在步骤S10,读入地震道信号,用读入的地震道数据将剩余信号初始化。每一道地震道信号都是以一个离散序列的数组形式存在的,假定为S(i)(i=1,2...1000),其中,i表示地震道数据的时间点或采样点。剩余信号是指每次匹配出原子后信号的剩余值,也是一个数组,假定为R(i)(i=1,2...1000)。将剩余信号初始化的意思就是:在j=1时,即第一次迭代时,将地震道信号的值一一赋给剩余信号,即R(j)=S(i),这个仅限于步骤S10一开始的时候,之后剩余信号就会发生变化了,每当通过迭代提取出一个原子,就要将对应的这个原子的振幅值序列在对应的时间从剩余信号中减掉。除了第一次迭代时,剩余信号=原始信号,其他时候的剩余信号都是减去提取的原子后的信号,是一个动态的即时更新的过程。
本实施例中的剩余信号可由流程图4中的S448步骤计算出。即,在步骤S448中,每次匹配出一个原子Aip,都将从当前剩余信号的对应的时间点上将Ajp的振幅值减去,从而更新剩余信号,大体示意图如图6所示。
在步骤S20,建立一个初始的基础原子库D,初始化动态原子库AD。在该步骤中,基础原子库的建立,可根据自定义的起始频率、终止频率和步长、相位等参数设定。假设选择的子波类型为雷克子波,起始频率、步长和终止频率分别为10HZ、1HZ、90HZ,则原子库里总共有81个雷克子波,频率分别为10HZ,11HZ,12HZ......90HZ。将该原子库作为初始的基础原子库D,同时也作为初始的动态原子库AD。其中,原子库的建立过程中的起始频率和终止频率的设定主要通过用户对地震数据的频谱分析或历史经验而来,比如说通过频谱分析,该区域频带范围主要为10-80HZ,那么就可将原子库的起始和终止频率设为10-80或相当范围即可,而步长的设置则需通过分解所花时间及分解精度综合考虑而来。此外,步骤S10和步骤S20的顺序可以互换,也可以同时执行。
在步骤S30,计算动态原子库AD里每个原子与剩余信号R(j)的互相关值序列。在该步骤中,互相关值序列长度由信号长度与原子长度决定,针对每一个原子,用户可以选取前m个互相关值A(k)m较大的原子作为最优原子,这样,整个动态原子库将有个最优原子,其中,k表示原子库中的原子的序号。
在步骤S40,一次性从个最优原子中筛选出n个互相关值较大且范围不重叠的原子,从而构建出有效原子集,并更新动态原子库AD。在该步骤中,n为用户设定的自然数,意为一次迭代提取的最多原子数目,可以根据用户所需分解精度来设置,同时受最终提取阈值以及是否重叠的限制。
原子重叠与否的标准就是看在时间上是否重叠,假定信号为S(i)(其中,i=1,2...1000),相关值最大的原子起始时间为第100个点,该原子的时窗长度为220,那么该原子时间域就是从100到319,如果相关值第二大的原子起始时间为300,那么这两个原子的范围就重叠了。
在步骤S50,根据有效原子集中的原子的互相关率的大小,从中提取匹配原子,更新剩余信号。
在步骤S60,利用剩余信号能量或已提取的匹配原子个数,判断是否达到迭代终止条件,如果没有达到,则重复步骤S30、S40、S50,直至满足终止条件。在该步骤中,终止条件有多个,地震道剩余能量小于设置的阈值的时候,或者提取的原子个数大于设定的个数时。地震道能量=各序号点振幅值的平方之和。假设S(i)={0.2,0.4,0.8,-0.2......},那么地震道能量就等于Sum=(0.2)2+(0.4)2+(0.8)2+(-0.2)2+......。剩余能量就是剩余信号的能量,每次分解出某个或多个原子后,从地震信号中减去该原子后,获得新的剩余信号,剩余能量等于更新后的剩余信号的振幅平方之和。
在步骤S70,保存所提取的原子信息到分解信息存储文件(起始时间、相关系数、频率、相位等)。
步骤S80,读取分解信息存储文件,结合井信息,有针对性的进行多子波筛选重构。在该步骤中,可根据频带范围、或者原子能量大小次序、或频率和能量的综合准则来重构地震道,且可实时动态的进行重构剖面QC(质量控制),有利于用户做出重构参数合理选择的判断。除了根据用户地质目标、井信息来有选择性的重构外,还可根据剖面的形态、波组特征等信息来判断重构的参数是否合理。
下面,结合附图4详细描述实施附图3中的方法的示例性流程图。
假设所需处理的地震道总时长为5s,以2ms采样率进行采样,则共有2501个采样点,地震道数据可表述为S(i)(i=1、2、3、......、2501)。
在步骤S410中,读入地震道数据,初始化剩余信号。该步骤与步骤S10相同。
在步骤S412中,生成基础原子库,初始化动态原子库AD。该步骤与步骤S20相同。
在步骤S414中,依次将动态原子库AD(j)(其中,j表示本轮迭代中第j次迭代)中的所有原子与剩余信号R(j)做互相关计算。具体地,依次选取第k个原子(在本示例中,1≤k≤81),将第k个原子与剩余信号R(j)求互相关,得到一个互相关序列。即,针对原子库中的每个原子,获得一个互相关序列,因此,在本示例中共有81个互相关序列。
在步骤S416中,针对每个原子,采用快速降序排序法,选出相关值最大、且相关范围不重叠的m个原子(m为用户自定义的自然数,例如,1≤m≤信号长度/子波长度),舍弃其他原子(S418)。这样整个动态原子库将匹配出个最优原子(步骤S420),并基于选取的个原子进行下面的步骤,以构建原子集以及更新动态原子库AD。但是,并不是每个原子一定能选出m个有效原子,也可能少于m个。
在上述步骤S414中,是将动态原子库AD中的所有原子均与同一个剩余信号求互相关。在步骤S416中,针对每个原子计算出的互相关值序列,选取前m个较大的互相关值。假如针对每个原子选取前5个较大的互相关值,则81个原子的总的互相关值为405个,即,在步骤S420中获得个最优原子。下面,描述构造原子集、生成动态原子库的详细步骤。
在步骤S420中,获得个最优原子。然后,在步骤S422中,对个最优原子进行快速降序排序,判断哪个原子的互相关值绝对值最大。提取出互相关值绝对值最大的那一个原子作为主原子P(步骤S424),然后对剩余的原子执行步骤S426。
在步骤S426,从S422中得到的个原子中(即,主原子P除外),依次判断是否进入第二类原子S,主要是根据该原子与本轮迭代已提取出来的原子(即主原子P以及筛选出的第二类原子S)在时间上是否重叠。具体地说,从中选取互相关值最大的那个原子,判断该原子所对应的时间与主原子P的时间是否重叠,如果不重叠,则加入第二类原子S,否则丢弃。然后,再次选择互相关值第二大的原子,判断该原子的互相关值与主原子P以及已经放入第二类原子S中的原子在时间上是否重叠,依次类推,直到选取出总共n个原子(主原子P+第二类原子S)。最后,将主原子P、第二类原子S共同作为本次迭代提取的有效原子,构建出有效原子集。
这里A(k)m,k表示第k个原子,m表示为原子库里序号为k的那个子波与剩余信号匹配出m个最优原子。举例说明:假定主频为10HZ,序号为3的原子,在时间为100的地方匹配出一个最优原子,在时间为500的地方匹配出一个最优原子,那么k=3,m=2。
在对降序排序的互相关值进行选取时,先选择互相关值最大的原子,然后选择互相关值第二大的原子,判断是否与互相关值第一大的原子重叠,是则有效,否则舍弃。然后判断互相关值第三大的原子,判断与前面选出的有效的原子是否重叠,如果不重叠,则保留,否则舍弃。如果选择出的较大互相关值对应的原子的个数达到了n,则停止该操作。其中,n的选取一方面需要考虑分解精度的要求,另一方面根据信号长度而来,同时在程序计算中,最终提取的有效原子个数将小于或等于n,是因为将受到最终提取阈值的限制。n的选取不宜过大,否则分解精度降低;也不宜过小,否则分解速度降低。
在上述步骤中,除根据原子范围重不重叠的准则来判断是否提取原子外,还能根据地震道能量是否减小的原则来判断是否提取该原子。地震道能量=各序号点振幅值的平方之和。假设S(i)={0.2,0.4,0.8,-0.2......},那么地震道能量就等于Sum=(0.2)2+(0.4)2+(0.8)2+(-0.2)2+......。判断能量是否减少的步骤如下:每次匹配出原子后,就从S(i)里将该原子减掉,然后计算能量,看是否小于上次迭代出的Sum,若减小,则提取该原子。
同时,在步骤S432中,针对动态原子库中的原子,获得每个原子的最大互相关绝对值。在步骤S434中,将动态原子库AD里的每个原子所对应的最大互相关值除以主原子P的最大互相关值计算每个原子的互相关率。
互相关值即为原子与信号互相关计算后得到的一个序列,互相关值有正负之分,如果两个子波完全相同,只是相位反转,相关值即为-1,如图8所示。最大互相关值就是每个原子和剩余信号的最大相关值,即为该序列中绝对值最大的那个互相关值;互相关率是一个相对概念,即为单个原子对应的最大互相关值除以所有原子中的最大互相关值,再求绝对值。
假设原子库里有81个原子,那么每个原子都会有一个最大互相关绝对值,假设存在多个相等的最大值时,任选一个即可(此处比较范围为同一个原子在不同时间的相关值,步骤S432)。也就是说,在步骤S420中,比较范围是个原子,从而得到主原子。而互相关绝对值最大的原子比较范围是81个原子的最大相关值,即,在步骤S432和S434中,比较范围是81个原子,针对每个原子,从不同时间的多个互相关值中选出最大的那个,除以主原子P的最大互相关值,从而计算出互相关率。
在步骤S436中,判断每个互相关率是否大于用户设定的动态原子库阈值,如果大于,则将该原子放入新的动态原子库AD,否则舍弃该原子。从而,在步骤S438中获得新的动态原子库AD(j+1),用于在下一次迭代中替换步骤S414中使用的动态原子库AD(j)。
动态原子库阈值也是个基准值,在界面上由用户设定或默认而来,动态原子库AD通过以下公式判断:
(每个原子的最大互相关值绝对值/所有原子的最大互相关值绝对值)>=动态原子库阈值。
假设10HZ对应的原子互相关率(互相关率=10HZ原子对应的最大互相关值绝对值/所有原子的最大互相关值绝对值)为0.5,设定的动态原子库阈值小于或等于0.5(动态原子库阈值越小,则分解精度越高,计算速度越慢,一般默认为0.32),那么10HZ这个原子就可认为是有效的原子,加入到动态原子库里,依次类推,计算出原始原子库里其他主频的原子是否有效,有效则加入到动态原子库里。其目的就是动态地减小原子库的规模。类似于比赛选手的选拔,假如共有400人参与选拔,就将初选目标定为上一次考试成绩大于80分的,这样选择范围就大大缩小,提高效率。
在上述步骤中构建原子集并更新动态原子库AD之后,在步骤S442中,计算各个原子的最大互相关率,即,每个原子的最大互相关值除以主原子的最大互相关值。在步骤S444中,依次从得到的原子集里选择互相关率最大的匹配原子Ajp,其中,j表示为第j次动态原子库与剩余信号求互相关,p表示本轮迭代的第p个原子。所谓匹配原子就是被匹配出的有效原子,也可理解为被选择的原子。然后在步骤S446中判断该原子Ajp的互相关率是否大于最终提取阈值。如果大于,则进行步骤S448,提取该有效原子Ajp,并计算新的剩余信号R(j+1)p,剩余信号R(j+1)p的计算是将迭代前采用的剩余信号减去所提取的原子Ajp后得到的剩余信号。然后进行步骤S450的判断。如果步骤S446的判断结果为否,则进行步骤S458,用新的动态原子库AD(j+1)更新动态原子库AD(j)。
通过步骤S424和S430,已经得到了一个有效原子集(原子集里共有n个原子,用Ajp表示其中一个原子),因此,步骤S444中的“依次”的概念就是从这个n个原子的原子集合里依次选择一个互相关率最大的原子。假如第一次选择的是互相关率最大的原子A1(可以判断出,第一次被选择的原子是主原子P,计算出的互相关率为1),则下一次被选择的原子是除了被选择的原子P之外剩余的原子中互相关率最大的原子(可以判断出,第二次被选择的原子为首先进入第二类原子S中的那个原子),以此类推。
对每一个原子而言,所谓互相关率就是把该原子的最大互相关值除以整个原子库里的最大互相关值,然后求绝对值,定义这个比值为互相关率。在步骤S446中,原子库中互相关值最大的原子为主原子P,将所提取出的原子Ajp的最大互相关值除以主原子P的最大互相关值并求绝对值,从而得到该原子Ajp的互相关率。假如互相关率大于用户给定的某个值(称为最终提取阈值,该阈值主要控制信号的分解精度,值越大则精度越高,但是相对的计算效率越低),那么就有效,则提取该原子Ajp,否则无效。
在步骤S450,判断是否达到迭代终止条件。如果达到了迭代终止条件,则终止迭代,分解完成(S460),保存所提取的原子信息到分解信息存储文件(起始时间、相关系数、频率、相位等);如果未达到迭代终止条件,进行步骤S452。
步骤S450中,采用的迭代终止条件可以有两个,第一个是可以通过判断所提取的原子是否达到了预定数量,例如,达到了200个,则可以终止迭代。第二个是通过判断剩余信号的能量与地震道数据的能量比值是否小于预定的阈值,例如,小于0.05,则可以终止迭代。
如果没有达到迭代终止条件,则在步骤S452判断本轮主原子和第二类原子是否已经提取完毕,如果为否,则返回步骤S444,进行下一次提取。如果已经提取完毕,则进行步骤S454。
在步骤S454,判断该原子Ajp的互相关率是否大于重新修改的原子库阈值。如果大于重新修改原子库阈值,则进行到步骤S458,采用在步骤S438中生成的动态原子库AD(j+1)更新当前的动态原子库AD(j),即令AD(j)=AD(j+1),并将第一剩余信号R(j)更新为第二剩余信号R(j+1)后,返回步骤S414,进行下一轮迭代,重新产生新的动态原子库并构建新的有效原子集;如果步骤S454的判断结果是小于重新修改的原子库阈值,则进行到步骤S456。在步骤S456中将动态原子库更新为第一步生成的基础原子库D,即81个完整的原子,并将第一剩余信号更新为第二剩余信号R(j+1)p,然后回到步骤S414,进行下一轮迭代,重新产生新的动态原子库并构建新的有效原子集。重新修改的原子库阈值也是由界面而来,可为设定的默认值,也可由用户给定。
如果在步骤S446中,判断的结果是Ajp的互相关率不大于最终提取阈值,也进行到步骤S458。
在完成分解后,可以读取存储的分解信息存储文件,结合井信息,有针对性的进行多子波重构(S462),重构方式有三种:①基于频率特征选择性重构;②基于能量特征选择性重构;③能量和频率特征综合选择性重构。
下面,详细解释一下上述三种重构方式。
分解后,原子信息里存储有主频、能量等信息,这样用户重构时可以通过判断子波的频率进行重构,也可根据能量进行重构,也可根据频率和能量的交集进行重构。
如图9所示,假如用户只希望选择主频为10-30HZ的原子进行重构,就将8HZ、14HZ、24HZ、26HZ这四个子波叠加即可。
假如用户想根据能量最大的1-4重构,那么就将26HZ、14HZ、32HZ、50HZ的加起来即可。假如用户想将主频10-30HZ、能量1-4的进行重构,那么就是以上两种情况的交集,即将26HZ和14HZ两个子波加起来。
在本发明的技术方案中,设置了原子库阈值、最终提取阈值、重新修改原子库阈值等参数,这些阈值是根据用户对地震资料处理时精度与速度的平衡选取而来,以上三个阈值中,最终提取阈值和重新修改的原子库阈值越大,则精度越高,分解速度越慢,而原子库阈值则是越小分解精度越高,分解速度越慢。本发明将由以上三个阈值组成的基本框架流程称为原子库阈值技术。同时,每次原子库与剩余信号互相关值计算后,提取了多个原子,本发明将该技术称为多原子提取技术。
如上所述,根据实施例的基于原子库阈值技术和多原子提取技术进行地震道的分解,可以极大地减少子波分解迭代次数,提高了计算速度,并可以对地震道进行相对准确的多子波分解(效果对比见图5中的表格)。
根据本发明的实施例采用了多原子提取技术,每次迭代提取出多个原子。这样,单次迭代所提取出来的能量明显增大,总迭代次数就会显著下降,运算效率明显提高;在匹配追踪的相关性计算中,先采用大间距进行滑动求相关,再在相关值相对较大的样点序号附近进行精细匹配,也大大降低了相关次数。
根据本发明的实施例采用了原子库阈值技术,建立了动态原子库,有效地减小了原子字典的冗余度,缩减了单次迭代的运行时间。