CC++如何实现字符串的佩里格拉姆属性检测？

本文共计2683个文字，预计阅读时间需要11分钟。

C++实现Perigram属性，描述地震波瞬时特征物理量：瞬时振幅、瞬时相位及瞬时频率（三瞬时参数），不仅可直接用来研究震源特性、构造等，还可反演。

C/C++实现Perigram属性

通常描述信号瞬时特征的物理量有：瞬时振幅、瞬时相位、及瞬时频率（“三瞬参数”），地震波的瞬时参数不仅可以直接用来研究岩性、构造等，而且也能够反演介质的品质因数等参数。在研究非平稳信号时，瞬时参数尤为重要。

假设原始信号为\(x(t)\)，通过Hilbert变换，将实信号转变为复信号\(S(t)=x(t)+iy(t)\)，并提取瞬时振幅、瞬时相位、瞬时频率三个参数。本文档将介绍使用C/C++代码，由三瞬属性得到Perigram以及Perigram×Cosine of Phase这两种属性的过程。

1 地震属性、意义及计算方法

(1)瞬时振幅 Instaneous Amplitude

(/振幅包络Amplitude Envelope/反射强度Reflection Strength)

瞬时振幅或反射强度是地震信号在某一时刻的总能量的平方根，它是地震轨迹的包络线。利用Hilbert变换提取复信号的虚部\(y(t)\)，瞬时振幅根据下式计算：

\[A(t)=\sqrt{x(t)^{2}+y(t)^{2}} \]

瞬时振幅代表了声阻抗（或反射系数）的差别，因此能够识别亮点、气藏聚集位置、层序边界、不整合面、大的岩性和沉积相变化、孔隙度或其它岩性参数的空间变化、断层的侧向变化等。其值总为正。

(2)瞬时相位 Instaneous Phase

瞬时相位是地震剖面上同相轴连续性的量度，无论能量的强弱，它的相位都能显示出来。其计算式为：

\[\theta(t)=\tan ^{-1}[\frac{y(t)} {x(t)}] \]

当波在各向异性的均匀介质中传播时，其相位是连续的；当波在有异常存在的介质中传播时，其相位将在异常位置发生显著变化，在剖面图中明显不连续。因此，利用瞬时相位能够较好地对地下分层和地下异常进行辨别，能够用来识别断层、尖灭、层序边界等。 瞬时相位数据呈现不连续的锯齿状，其值域范围为\(-180°\sim180°\)。

(3)瞬时频率 Instaneous Frequency

瞬时频率表示瞬时相位随时间的变化率：

\[\omega(t)=\frac{d \theta(t)}{d t} \]

瞬时频率能够在横向上有效判别地震资料的相关性，薄层的岩性和厚度的变化能够引起瞬时频率的变化。一般高瞬时频率反映了薄的泥岩沉积，低的瞬时频率反映富砂沉积，因此可以应用瞬时频率反映碎屑岩的沉积环境 。

(4) 反射强度交流分量 Perigram

Perigram是去除了直流分量的反射强度，反射强度低频（直流）分量定义为给定时间窗口\(\Delta T\)内反射强度的平均值：

\[\bar{A}(t)=\frac{1}{\Delta T} \int_{(t-\Delta T / 2)}^{(t+\Delta T / 2)} A(\tau) d \tau \]

当取一个极大的窗口时，\(A(t)\)接近直流分量，然后将其从反射强度中减去低频分量，得到Perigram：

\[g(t)=A(t)-\bar{A}(t) \]

Perigram 与反射强度的用途基本相同。其能够加强反射强度峰值异常，并且更适合分析和处理，因为其有正负。

(5)反射强度交流分量与相位余弦之积

公式表示为：

\[G(t)=g(t)\cos \theta(t) =[A(t)-\bar{A}(t)]\cos \theta(t) =x(t)[1-\frac{\bar{A}(t)}{{A}(t)}],\quad if\quad g(t)>0\\ G(t)=0,\quad if\quad g(t)\le 0\\ \]

当Perigram>0时，该属性为反射强度交流分量与相位余弦的乘积，否则为0。换言之，当Perigram值为正时，Perigram×相位余弦的值等于输入数据\(x(t)\)乘以一个小于1的比例；而当Perigram值为负时，值为0，对应低反射强度。该属性可用于突出强振幅与连续相位，同反射强度的用途相似。例如，当周围低能反射体被减弱至0时，与含气砂有关的部分将被明显点亮。

2 代码实现 (1) Perigram

void Compute_perigram(float* dataInput, int nt, float* Perigram){ float *Imaginary = NULL; //定义变换后得到的虚部 float *Amplitude_envelope = NULL; //振幅包络 float *DC_component = NULL; //直流分量 float sum = 0.0f; Imaginary = (float*)calloc(nt, sizeof(float)); Amplitude_envelope = (float*)calloc(nt, sizeof(float)); DC_component = (float*)calloc(nt, sizeof(float)); hilbert(nt, dataInput, Imaginary); //希尔伯特变换 for (int i = 0; i < nt; i++) { Amplitude_envelope[i] = sqrt(pow(dataInput[i], 2) + pow(Imaginary[i], 2)); //实部、虚部平方和开根号求瞬时振幅 sum += Amplitude_envelope[i]; } float Avg = sum/nt; for (int i = 0; i < nt; i++) { DC_component[i] = Avg; //直流分量为全局平均 Perigram[i]=Amplitude_envelope[i]- DC_component[i]; //Perigram } free(Imaginary); free(Amplitude_envelope); free(DC_component); } (2) Perigram*Cosine of Phase

void Compute_perigramXcosine(float* dataInput, int nt, float* PerigramXcos) { float *Imaginary = NULL; //定义变换后得到的虚部 float *Amplitude_envelope = NULL; //振幅包络 float *DC_component = NULL; //直流分量 float sum = 0.0f; Imaginary = (float*)calloc(nt, sizeof(float)); Amplitude_envelope = (float*)calloc(nt, sizeof(float)); DC_component = (float*)calloc(nt, sizeof(float)); hilbert(nt, dataInput, Imaginary); //希尔伯特变换 for (int i = 0; i < nt; i++) { Amplitude_envelope[i] = sqrt(pow(dataInput[i], 2) + pow(Imaginary[i], 2)); sum += Amplitude_envelope[i]; } float Avg = sum / nt; for (int i = 0; i < nt; i++) { DC_component[i] = Avg; if ((Amplitude_envelope[i] - DC_component[i]) > 0) { //如果g(t)大于0 PerigramXcos[i] = dataInput[i] * (1 - DC_component[i] / Amplitude_envelope[i]); } else { //g(t)小于0，那乘积余弦就赋为0 PerigramXcos[i] = 0; } } free(Imaginary); free(Amplitude_envelope); free(DC_component); } 主函数

int compute_type = 1; //选择计算属性类别，当前为Perigram //遍历每一道 for (int itrace = 0; itrace < ncdp; itrace++) { fread(&traceheader, sizetraceheader, 1, fpinput); fread(dataInput, nt * sizeof(float), 1, fpinput); if (compute_type==1) //Perigram { float* Perigram = NULL; Perigram = (float*)calloc(nt, sizeof(float)); Compute_perigram(dataInput, nt, Perigram); fwrite(&traceheader, sizetraceheader, 1, fpoutput); fwrite(Perigram, nt * sizeof(float), 1, fpoutput); free(Perigram); } else { //Perigram*cos of phase float* Perigram_cos = NULL; Perigram_cos = (float*)calloc(nt, sizeof(float)); Compute_perigramXcosine(dataInput, nt, Perigram_cos); fwrite(&traceheader, sizetraceheader, 1, fpoutput); fwrite(Perigram_cos, nt * sizeof(float), 1, fpoutput); free(Perigram_cos); } } 完整代码 I hilbert.h

#include "convolution.h" void hilbert(int n, float x[], float y[]); II hilbert.cpp

#include "hilbert.h" #define PI 3.14159265358 #define LHHALF 30 /* half-length of Hilbert transform filter*/ #define LH 2*LHHALF+1 /* filter length must be odd */ void hilbert(int n, float x[], float y[]) /***************************************************************************** Compute Hilbert transform y of x ****************************************************************************** Input: n length of x and y x array[n] to be Hilbert transformed Output: y array[n] containing Hilbert transform of x *****************************************************************************/ { static int madeh = 0; static float h[LH]; int i; float taper; /* if not made, make Hilbert transform filter; use Hamming window */ if (!madeh) { h[LHHALF] = 0.0; for (i = 1; i <= LHHALF; i++) { taper = 0.54 + 0.46*cos(PI*(float)i / (float)(LHHALF)); h[LHHALF + i] = taper * (-(float)(i % 2)*2.0 / (PI*(float)(i))); h[LHHALF - i] = -h[LHHALF + i]; } madeh = 1; } /* convolve Hilbert transform with input array */ convolve_cwp(LH, -LHHALF, h, n, 0, x, n, 0, y); //需要用到卷积函数，见www.cnblogs.com/GeophysicsWorker/p/16399209.html //convolution.cpp } III Perigram_Attributes.h

#include<stdio.h> #include<stdlib.h> #include<math.h> #include<string.h> #include"hilbert.h" void Compute_perigram(float* dataInput, int nt, float* Perigram); void Compute_perigramXcosine(float* dataInput, int nt, float* PerigramXcos); IV Perigram_Attributes.cpp

#include "Perigram_Attributes.h" void Compute_perigram(float* dataInput, int nt, float* Perigram){ float *Imaginary = NULL; //定义变换后得到的虚部 float *Amplitude_envelope = NULL; //振幅包络 float *DC_component = NULL; //直流分量 float sum = 0.0f; Imaginary = (float*)calloc(nt, sizeof(float)); Amplitude_envelope = (float*)calloc(nt, sizeof(float)); DC_component = (float*)calloc(nt, sizeof(float)); hilbert(nt, dataInput, Imaginary); //希尔伯特变换 for (int i = 0; i < nt; i++) { Amplitude_envelope[i] = sqrt(pow(dataInput[i], 2) + pow(Imaginary[i], 2)); //实部、虚部平方和开根号求瞬时振幅 sum += Amplitude_envelope[i]; } float Avg = sum/nt; for (int i = 0; i < nt; i++) { DC_component[i] = Avg; //直流分量为全局平均 Perigram[i]=Amplitude_envelope[i]- DC_component[i]; //Perigram } free(Imaginary); free(Amplitude_envelope); free(DC_component); } void Compute_perigramXcosine(float* dataInput, int nt, float* PerigramXcos) { float *Imaginary = NULL; //定义变换后得到的虚部 float *Amplitude_envelope = NULL; //振幅包络 float *DC_component = NULL; //直流分量 float sum = 0.0f; Imaginary = (float*)calloc(nt, sizeof(float)); Amplitude_envelope = (float*)calloc(nt, sizeof(float)); DC_component = (float*)calloc(nt, sizeof(float)); hilbert(nt, dataInput, Imaginary); //希尔伯特变换 for (int i = 0; i < nt; i++) { Amplitude_envelope[i] = sqrt(pow(dataInput[i], 2) + pow(Imaginary[i], 2)); sum += Amplitude_envelope[i]; } float Avg = sum / nt; for (int i = 0; i < nt; i++) { DC_component[i] = Avg; if ((Amplitude_envelope[i] - DC_component[i]) > 0) { //如果相减后仍大于0 PerigramXcos[i] = dataInput[i] * (1 - DC_component[i] / Amplitude_envelope[i]); } else { //相减后小于0，那乘积余弦也就赋为0 PerigramXcos[i] = 0; } } free(Imaginary); free(Amplitude_envelope); free(DC_component); } V main.cpp

#include<stdio.h> #include<stdlib.h> #include<math.h> #include<string.h> #include<iostream> #include<vector> #include<algorithm> #include"hilbert.h" #include"segy.h" #include"Perigram_Attributes.h" int main() { const char *filenameInput = "test.sgy"; const char *filenameOutput = "Output66.sgy"; bhed fileheader; segy traceheader; unsigned int nt = 0; unsigned int sizefileheader = sizeof(fileheader); unsigned int sizetraceheader = sizeof(traceheader); unsigned int ncdp = 0; long long size_file = 0; long long size_trace = 0; FILE *fpinput = NULL; FILE *fpoutput = NULL; float *dataInput = NULL; float *dataOutput = NULL; fpinput = fopen(filenameInput, "rb"); fpoutput = fopen(filenameOutput, "wb"); if (fpinput == NULL) { printf("can not open %s file\n", filenameInput); return false; } if (fpoutput == NULL) { printf("can not open %s file\n", filenameOutput); return false; } fread(&fileheader, sizefileheader, 1, fpinput); nt = fileheader.hns; _fseeki64(fpinput, 0, SEEK_END); size_file = _ftelli64(fpinput); size_trace = nt * sizeof(float) + sizetraceheader; ncdp = (size_file - (long long)sizefileheader) / size_trace; _fseeki64(fpinput, sizefileheader, SEEK_SET); fwrite(&fileheader, sizefileheader, 1, fpoutput); dataInput = (float*)calloc(nt, sizeof(float)); dataOutput = (float*)calloc(nt, sizeof(float)); /************以上为segy读写，不再赘述***********/ int compute_type = 2; //遍历每一道 for (int itrace = 0; itrace < ncdp; itrace++) { fread(&traceheader, sizetraceheader, 1, fpinput); fread(dataInput, nt * sizeof(float), 1, fpinput); if (compute_type==1) //Perigram { float* Perigram = NULL; Perigram = (float*)calloc(nt, sizeof(float)); Compute_perigram(dataInput, nt, Perigram); fwrite(&traceheader, sizetraceheader, 1, fpoutput); fwrite(Perigram, nt * sizeof(float), 1, fpoutput); free(Perigram); } else { //Perigram*cos float* Perigram_cos = NULL; Perigram_cos = (float*)calloc(nt, sizeof(float)); Compute_perigramXcosine(dataInput, nt, Perigram_cos); fwrite(&traceheader, sizetraceheader, 1, fpoutput); fwrite(Perigram_cos, nt * sizeof(float), 1, fpoutput); free(Perigram_cos); } } free(dataInput); free(dataOutput); fclose(fpinput); fclose(fpoutput); return 0; }

参考：

1 唐金炎. 基于地震属性参数的地震相识别方法研究[D].成都理工大学,2011.

2 高静怀,汪文秉,朱光明.小波变换与信号瞬时特征分析[J].地球物理学报,1997(06):821-832.

3 BASIC THEORY AND GEOLOGICAL MEANING OF SEISMIC ATTRIBUTES (linkedin.com)

4 B. Gelchinsky, E. Landa, and V. Shtivelman, (1985), "Algorithms of phase and group correlation," GEOPHYSICS 50: 596-608.doi.org/10.1190/1.1441935

5 SHTIVELMAN, V., LANDA, E. and GELCHINSKY, B. (1986), PHASE AND GROUP TIME SECTIONS AND POSSIBILITIES FOR THEIR USE IN SEISMIC INTERPRETATION OF COMPLEX MEDIA. Geophysical Prospecting, 34: 508-536.doi.org/10.1111/j.1365-2478.1986.tb00479.x

6 PostStack™ Family Software Reference Manual.Landmark Software & Services.