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);
}
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);
}
}
#include "convolution.h"
void hilbert(int n, float x[], float y[]);
#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);
//需要用到卷积函数，见https://www.cnblogs.com/GeophysicsWorker/p/16399209.html
//convolution.cpp
}
#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);
#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);
}
#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.https://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.https://doi.org/10.1111/j.1365-2478.1986.tb00479.x
6 PostStack™ Family Software Reference Manual.Landmark Software & Services.