w0(n,:)=exp(-n.^(2*p)*k.^2/(2*N.^2))*(n.^p)/(N*sqrt(2*pi))是窗函数，k是t的离散取样意义。公式的话就是简单的在S变换里把f变成了f^p（f的p次方)
这里p是要赋值的，它是给定的一个数值。如果p=1那就是原S变换，如果p>1则窗函数随f变化呈指数增长，一般这个效果不值得应用。反而是0<p<1时，则窗函数不对称性和窗宽随f呈对数变化，f小时窗变化剧烈，f大时窗变化反而很小，会比较有用的。
冯强强1，2 温明明1，2 吴衡1，2
（1.广州海洋地质调查局 广州 510760；2.国土资源部海底矿产资源重点实验室 广州 510760）
第一作者简介：冯强强（1987—），男，硕士，助理工程师，现在从事海上地球物理探测方法的研究工作。
摘要 本文选取了两条不同的单道地震接收电缆（AAE接收电缆、GEO接收电缆）在野外采集到的原始数据，并利用快速傅立叶变换（FFT）以及S变换等方法，从数据本身、频率域、时频域方面进行了对比分析。结果表明接收段较短、且内部检波器组合更紧致的AAE电缆具有较高的浅层分辨率，更适合于浅水作业，而接收段较长，检波器较多，检波器间距稍大的GEO电缆可以采集到高信噪比的资料，更适合于深水作业。
关键词 单道地震 S 变换 时频域 检波器
1 引言
单道地震是解决海上工程地质问题以及海洋区域地质调查的一种重要手段［1，2］，相对于使用全波形反射的多道地震调查，单道地震对震源和检波器的偏移距离相对于水深和对反射倾角的要求很小，另外，由于其施工方法简单、配置灵活，处理简单、便于解释、正被广泛应用于全球的海洋地球物理调查中，如海域陆架区域地球物理勘查［3，4］，另外在新能源-天然气水合物的勘探［5，6］中也得到了应用。
单道地震接收电缆的好坏对于高质量的单道地震数据的采集有着重要的影响［7］，常见的接收电缆通常采用多个检波器组合的形式获取单道信号，因此检波器的灵敏度以及其组合方式等都会影响到采集数据的质量。另外，作业环境对于单道地震数据的采集也有着重要的影响［8］。
2 单道地震电缆设备以及数据比较分析
2.1 单道地震接收电缆性能比较
单道地震电缆内部的检波器通常是由压电材料制造，利用压电效应［9］将地震波引起的水压变化转变成电信号，然后通过电缆或者光缆等介质传输到船上的记录系统进行显示和存储。早期的海上地震勘探通常使用悬浮在固定或者缓慢漂流船只上的单个检波器来进行采集［10］。后来发展到特定排列方式组成的阵列拖缆拖在船尾，检波器被集成安置在充满油的柔韧塑料管中，保持与海水相适应的波阻抗，保证检波器与水的紧耦合。这一段包含检波器的塑料管就是本文所指的信号接收段，这种装置形式被设计成流线型来减小机器震动和水流引起的噪声。
本文所比较的两个电缆，一条是 Applied Acoustics Engineering公司生产的AAE20型单道接收电缆（以下简称为AAE电缆，图1左），其由20个检波器单元组成，间距0.15m，检波器段长度4.5m，其频率响应为145～7KHz（-3dB），总的灵敏度为-167 dB ref 1 v per。另外一条是Geo Marine Survey Systems公司生产的48检波器单元的Geo⁃Sense Mini⁃Streamers接收电缆（以下简称为GEO电缆，图1右），检波器间距0.4m，接收部分总长约20m，其检波器型号为AQ-2000，性能参数如下：灵敏度（@100Hz）为-201dB ref 1v per+/-1.5dB，频率响应为1Hz到10KHz（+/-2dB）。
图1 两种电缆实际形态（左边蓝色为AAE型电缆，右边红色为GEO型电缆）
Fig.1 The samples oftwo cables（AAE Cable on the left，GEO Cable on the right）
2.2 数据采集及比较分析
本次野外数据采集作业水域水深在250m左右，统一采用GEO-6 KJ电火花震源系统作为激发源，选取合适且大小相同的激发能量，同时为了避免船尾螺旋桨引起的尾流以及环境噪声［11］对数据产生影响，采用如下装置参数：接收电缆释放长度选择为45m，炮检距（接收电缆与震源缆之间的距离）约为10m。另外，采集参数如下：采样频率（Sample frequency）4000Hz，延时（Delay）300ms，记录长度（Record Length）500ms，因此每道数据有2000个采样点
本文从野外实测数据中选取了一组进行比较，该组测线分别由上述两种不同的电缆进行数据采集，间距约300m，地质情况类似，且作业时测线方向一致，海况及作业环境相差不大，因此具有较强的可比性。
图2和图3分别显示GEO电缆和AAE电缆的原始实测数据，从图中可以看出两条电缆都有良好的穿透能力以及地层分辨能力，所采数据也都能满足常规单道地震数据的采集规范。
由于作业区域水相对较浅，AAE电缆检波器接收段较短且接收段内检波器之间距离较小，因此从图3的地震剖面可以看出该电缆接收到的数据在浅部有较高的分辨率。相比之下，GEO电缆检波器多且间距较大，其地震剖面信噪比高，尤其深部地层反射信号更加清晰。
图2 GEO电缆采集到的数据
Fig.2 Data collected by GEO Cable
图3 AAE电缆采集到的数据
Fig.3 Data collected by AAE Cable
图4和图5为从上述剖面中各抽取的一道数据，由于前置放大器设置不同造成两条电缆信号回波强度不一致，并不影响数据的信噪比及分辨率。
图4 GEO电缆抽取一道数据
Fig.4 One channel extraction of GEO Cable
图5 AAE电缆抽取一道数据
Fig.5 One channel extraction of AAE Cable
另外，分别对所抽取数据进行快速傅立叶变换（FFT），观测其频率域特征，图6和图7分别为上述两道数据的振幅谱，可以看到两条电缆能量主要集中在0-1000Hz之间，但是单从频率域并不能对两条电缆的分辨率、穿透能力以及信噪比等参数进行直观分析。
图6 GEO抽取数据的振幅谱
Fig.6 The amplitude spectrogram of data extraction by GEO Cable
图7 AAE电缆抽取数据的振幅谱
Fig.7 The amplitude spectrogram of data extraction by AAE Cable
因此，本文引入了S变换对所抽取数据进行了分析。传统的傅立叶变换虽然能对信号的整体性质进行分析，但却难获得信号的一些瞬时信息，特别是对于非平稳信号，而这些局部性质却是分析这些非平稳信号的关键。时频分析技术可以同时从时间域和频率域分析信号的局部特征，因此在傅立叶分析的基础上，发展了一系列的时频域信号分析论，如短时窗傅立叶变换、小波变换等。
S变换是由Stockwell（1996）等提出的，是连续小波变换的一种扩展，它的原理是基于一个高度和宽度随频率变化的高斯窗。函数的S变换定义如下：
南海地质研究（2014）
其中t表示时间，f表示频率，τ是一个控制高斯窗函数 在时间轴上位置的一个参数，e-2πft是一个相位因子，起到相位校正作用，这是较小波变换的一个优点［13］。另外，其离散形式可以利用快速傅立叶变换FFT的结果来计算，因此计算效率比较高。信号h（t）的离散形式h［kT］，k=0，1，…，N-1（T为采样间隔）的离散傅立叶变换表示如下：
南海地质研究（2014）
其中n=0，1…，N-1，则h（kT）的离散S变换可表示为：
南海地质研究（2014）
其中：j、m、n=0，1…，N-1。
图8和图9分别为GEO电缆和AAE电缆所抽取数据的S变换显示，从中可以看出两条电缆都能够采集到有效的地层反射信号，海底以及近海底地层反射信号最强，深部地层的反射信号也明显被记录到。另外，前述的AAE电缆浅层分辨率高，GEO电缆信噪比高也在图上有直观的反映。
3 结论及分析
本文采集了两条不同的单道地震接收电缆的数据，通过观察其剖面、频率域特征、时频域等特征，对上述两条电缆的特性、适用范围等进行了分析比较。
图8 GEO电缆抽取数据的S变换域显示
Fig.8 The sample of S⁃transform domain of data extraction by GEO Cable
图9 AAE电缆抽取数据的S变换域显示
Fig.9 The sample of S⁃transform domain of data extraction by AAE Cable
AAE电缆其检波器接收部分是由20个检波器组合形成，可以有效的压制环境噪声［14］。此外，由于其检波器间距只有0.15m，因此整个接收段部分长度仅为4.5m，相对较短，这种装置形式对于浅水区有明显的优势，在满足质量要求的前提下具有更高的分辨率。而本文所使用的GEO电缆其接收段较长，且检波器组合达到48个，除了达到压制噪声的目的外还可以有效采集到来自深水区以及深部地层的反射波，避免丢失回声数据，另外其具有较宽的频率响应以及良好的灵敏度，使其能采集到高信噪比的资料。
检波器组合具有低频响应［7］，因此会对高分辨率数据采集有一定影响，但是检波器组合又是高信噪比数据采集必不可少的，因此在野外作业时要根据作业环境、接收条件、工作目的等选择合适的接收电缆。
致谢
感谢参与数据采集的全体成员，中国地震局地球物理研究所吴庆举研究员提供了本文的S变换程序，温明明高工、吴衡高工、牟泽霖、林桂柱在本文的资料分析中给予悉心指导！
参考文献
［1］李军峰，肖都，孔广胜，等.2004.单道海上反射地震在海上物探工程中的应用［J］.物探与化探，28（4）：365 368
［2］岳保静.2010.单道地震资料处理方法及应用［D］.硕士学位论文，青岛：中国科学院海洋研究所
［3］D Van Rooij，B De Mol，V Huvenne，et al.2003.Seismic evidence of current⁃controlled sedimentation in the Belgica mound province，upper Porcupine slope，southwest of Ireland［J］.Marine Geology，195（1-4）：31-53
［4］Nicolas Weber，Eric Chaumillon，Michel Tesson，et al.2004.Architecture and morphology of the outer segment of a mixed tide and wave⁃dominated⁃incised valley，revealed by HR seismic reflection profiling：the paleo⁃Charente River，France［J］.Marine Geology，207（1-4）：17-38
［5］张光学，黄永祥，陈邦彦.2003.海域天然气水合物地震学［M］.北京：海洋出版社
［6］吴志强，文丽，童思友，等.2007.海域天然气水合物的地震研究进展［J］.地球物理学进展，22（1）：218-227
［7］付清锋.2005.地震检波器新技术发展方向［J］.石油仪器，19（6）：1-4
［8］李丽青，梁蓓雯，徐华宁.2007.海上单道地震资料中多次波的衰减［J］.石油物探，46（5）：457-464
［9］吴毅强，邓忠华.1995.压电器件［M］.北京：电子工业出版社
［10］E J W Jones.1999.Marine Geophysics［M］.John Wiley＆Sons
［11］刘建勋.2007.提高海上单道反射地震记录信噪比和分辨率的方法技术［J］.物探化探计算技术，29（增刊）：116-122
［12］Stockwell R G，Manisinha，R P Lowe.1996.Localization of the complex spectrum：the S transform［J］.IEEE Trans.Signal Process，44：998-1001
［13］Stefano Parolai.2009.Denoising of Seismograms Using the S Transform［J］.Bulletin of the Seismological Society of America，99（1）：226-234
［14］陆基孟.1993.地震勘探原理（上、下册）［M］.北京：石油大学出版社
Comparison Between Two Different Single-Channel Seismic Streamers
Feng Qiang Qiang1，2，Wen Ming Ming1，2，Wu Heng1，2
（1.Guangzhou Marine Geological Survey，Guangzhou，510760；2.Key laboratory of Marine Mineral Reasources，MLR，Guangzhou，510760）
Abstract：In this paper we collect the raw data from two different singlechannel reflection seismic streamers（AAE streamers and GEO streamers）and analyze those data comparatively from frequency domain and time-frequency domain by using fast Fourier transform（FFT）and Stransform.The result shows that AAE streamers which have shorter receiving section and compact geophones can get higher shallow resolution and therefore are more suitable for shallow water.While GEO streamers are more suitable for deep water because of its longer receiving section，more geophones，and the ability to acquire high signal⁃to⁃noise ratio data with longer span between geophones.
Key word：Single⁃channel reflection seismic；S⁃transform；Time⁃frequency domain；Geophone
CS%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%产生Stripmap SAR的回波
clear all
clc
thetaT=0%T平台波束斜视角
thetaT=thetaT*pi/180%rad
thetaR=0%R平台波束斜视角
thetaR=thetaR*pi/180
c=3e8%光速
fc=1.5e9%载频
lambda=c/fc%波长
%%测绘带区域
X0=200%方位向[-X0,X0]
Rtc=3000
Rrc=3000
Rc=(Rtc+Rrc)/2
R0=150%距离向[Rc-R0,Rc+R0]
%%距离向(Range),r/t domain
Tr=1.33e-6%LFM信号脉宽1.33us (200m)
Br=150e6%LFM信号带宽 150MHz
Kr=Br/Tr%调频斜率
Nr=1024
r=Rc+linspace(-R0,R0,Nr)
t=2*r/c%t域序列
dt=R0*4/c/Nr%采样周期
f=linspace(-1/2/dt,1/2/dt,Nr)%f域序列
%%方位向(Azimuth,Cross-Range),x/u domain
v=100%SAR 平台速度
Lsar=300%合成孔径长度
Na=512
x=linspace(-X0,X0,Na)%u域序列
u=x/v
du=2*X0/v/Na
fu=linspace(-1/2/du,1/2/du,Na)%fu域序列
ftdc=v*sin(thetaT)
ftdr=-(v*cos(thetaT))^2/lambda/Rtc
frdc=v*sin(thetaR)
frdr=-(v*cos(thetaR))^2/lambda/Rrc
fdc=ftdc+frdc%Doppler调频中心频率
fdr=ftdr+frdr%Doppler调频斜率
%％目标位置
Ntar=3%目标个数
Ptar=[Rrc,0,1 %距离向坐标,方位向坐标,sigma
Rrc+50,-50,1
Rrc+50,50,1]
%%产生回波
s_ut=zeros(Nr,Na)
U=ones(Nr,1)*u%扩充为矩阵
T=t'*ones(1,Na)
for i=1:1:Ntar
rn=Ptar(i,1)xn=Ptar(i,2)sigma=Ptar(i,3)
rtn=rn+Rtc-Rrc
RT=sqrt(rtn^2+(rtn*tan(thetaT)+xn-v*U).^2)
RR=sqrt(rn^2+(rn*tan(thetaT)+xn-v*U).^2)
R=RT+RR
DT=T-R/c
phase=-pi*Kr*DT.^2-2*pi/lambda*R
s_ut=s_ut+sigma*exp(j*phase).*(abs(DT)<Tr/2).*(abs(v*U-xn)<Lsar/2)
end
%方位向fft
s_kt=fftshift(fft(fftshift(s_ut).')).'
%CS变换
kc=4*pi/lambda
kc=kc*ones(1,Na)
kx=fu/v
p_kx0=-sqrt(kc.^2-kx.^2)%相位项泰勒展开的系数函数
p_kx1=2*kc/c/p_kx0
p_kx2=-2.*kx.^2/c^2./p_kx0.^3
C_kx=-(c*p_kx1/2+1)
Ks_r=1-2*Kr*Rc.*p_kx2
Ks_kx_r=Kr/pi./Ks_r
r0=Rc
s2_ut=exp(j*pi*C_kx.*ones(Nr,1)*Ks_kx_r.*(t'*ones(1,Na)-2*r0*(1+C_kx)/c).^2)%设计的线性调频信号
S_cs=s_kt.*s2_ut
%距离向fft
S_kw=fftshift(fft(fftshift(S_cs)))
%距离向匹配滤波
w=2*pi*f
rmc_r=exp(j.*w*2*C_kx*r0/c).*exp(j.*w.^2/4/pi/Kr/(1+C_kx))
rmc_r=rmc_r'*ones(1,Na)
S_rmc=S_kw.*rmc_r
%距离向ifft
S_kt=fftshift(ifft(fftshift(S_rmc)))
d_kxr=4*pi/c^2*Kr*C_kx*(1+C_kx).*(Rc-r0).^2%CS变换带来的相位误差
S_kt=S_kt.*exp(-j*d_kxr)%消除相位误差
%方位向匹配滤波
FU=ones(Nr,1)*fu
H_kx=exp(j*pi/fdr*(FU-fdc).^2)%方位向压缩因子
I_ut=S_kt.*H_kx
I_ut=fftshift(ifft(fftshift(I_ut.'))).'
subplot(221)
G=20*log10(abs(s_ut)+1e-6)
gm=max(max(G))
gn=gm-40%显示动态范围40dB
G=255/(gm-gn)*(G-gn).*(G>gn)
imagesc(x,r-Rc,-G),colormap(gray)
grid on,axis tight,
xlabel('Azimuth')
ylabel('Range')
title('(a)原始信号')
subplot(222)
G=20*log10(abs(S_rmc)+1e-6)
gm=max(max(G))
gn=gm-40%显示动态范围40dB
G=255/(gm-gn)*(G-gn).*(G>gn)
imagesc(x,r-Rc,-G),colormap(gray)
grid on,axis tight,
xlabel('Azimuth')
ylabel('Range')
title('(b)距离向匹配滤波后频谱')
subplot(223)
G=20*log10(abs(S_kt)+1e-6)
gm=max(max(G))
gn=gm-40%显示动态范围40dB
G=255/(gm-gn)*(G-gn).*(G>gn)
imagesc(x,r-Rc,G),colormap(gray)
grid on,axis tight,
xlabel('Azimuth')
ylabel('Range')
title('(c)消除相位误差后频谱')
subplot(224)
G=20*log10(abs(I_ut)+1e-6)
gm=max(max(G))
gn=gm-60%显示动态范围40dB
G=255/(gm-gn)*(G-gn).*(G>gn)
imagesc(x,r-Rc,G),colormap(gray)
grid on,axis tight,
xlabel('Azimuth')
ylabel('Range')
title('(d)目标图象')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
RD
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%产生Stripmap SAR的回波
clear all
thetaT=0%T平台波束斜视角
thetaT=thetaT*pi/180%rad
thetaR=0%R平台波束斜视角
thetaR=thetaR*pi/180
c=3e8%光速
fc=1.5e9%载频
lambda=c/fc%波长
%%测绘带区域
X0=200%方位向[-X0,X0]
Rtc=3000
Rrc=3000
Rc=(Rtc+Rrc)/2
R0=150%距离向[Rc-R0,Rc+R0]
%%距离向(Range),r/t domain
Tr=1.5e-6%LFM信号脉宽 1.5us (200m)
Br=150e6%LFM信号带宽 150MHz
Kr=Br/Tr%调频斜率
Nr=512
r=Rc+linspace(-R0,R0,Nr)
t=2*r/c%t域序列
dt=R0*4/c/Nr%采样周期
f=linspace(-1/2/dt,1/2/dt,Nr)%f域序列
%%方位向(Azimuth,Cross-Range),x/u domain
v=100%SAR 平台速度
Lsar=300%合成孔径长度
Na=1024
x=linspace(-X0,X0,Na)%u域序列
u=x/v
du=2*X0/v/Na
fu=linspace(-1/2/du,1/2/du,Na)%fu域序列
ftdc=v*sin(thetaT)
ftdr=-(v*cos(thetaT))^2/lambda/Rtc
frdc=v*sin(thetaR)
frdr=-(v*cos(thetaR))^2/lambda/Rrc
fdc=ftdc+frdc%Doppler调频中心频率
fdr=ftdr+frdr%Doppler调频斜率
%％目标位置
Ntar=3%目标个数
Ptar=[Rrc,0,1 %距离向坐标,方位向坐标,sigma
Rrc+50,-50,1
Rrc+50,50,1]
%%产生回波
s_ut=zeros(Nr,Na)
U=ones(Nr,1)*u%扩充为矩阵
T=t'*ones(1,Na)
for i=1:1:Ntar
rn=Ptar(i,1)xn=Ptar(i,2)sigma=Ptar(i,3)
rtn=rn+Rtc-Rrc
RT=sqrt(rtn^2+(rtn*tan(thetaT)+xn-v*U).^2)
RR=sqrt(rn^2+(rn*tan(thetaT)+xn-v*U).^2)
R=RT+RR
DT=T-R/c
phase=pi*Kr*DT.^2-2*pi/lambda*R
s_ut=s_ut+sigma*exp(j*phase).*(abs(DT)<Tr/2).*(abs(v*U-xn)<Lsar/2)
end
%%距离压缩
p0_t=exp(j*pi*Kr*(t-2*Rc/c).^2).*(abs(t-2*Rc/c)<Tr/2)%距离向LFM信号
p0_f=fftshift(fft(fftshift(p0_t)))
s_uf=fftshift(fft(fftshift(s_ut)))%距离向FFT
src_uf=s_uf.*(conj(p0_f).'*ones(1,Na))%距离压缩
src_ut=fftshift(ifft(fftshift(src_uf)))%距离压缩后的信号
src_fut=fftshift(fft(fftshift(src_ut).')).'%距离多普勒域
%%二次距离压缩,距离迁移校正原理仿真
src_fuf=fftshift(fft(fftshift(src_uf).')).'%距离压缩后的二维频谱
F=f'*ones(1,Na)%扩充为矩阵
FU=ones(Nr,1)*fu
p0_2f=exp(j*pi/fc^2/fdr*(FU.*F).^2+j*pi*fdc^2/fc/fdr*F-j*pi/fc/fdr*FU.^2.*F)
s2rc_fuf=src_fuf.*p0_2f
s2rc_fut=fftshift(ifft(fftshift(s2rc_fuf)))%距离多普勒域
%%方位压缩
p0_2fu=exp(j*pi/fdr*(FU-fdc).^2)%方位向压缩因子
s2rcac_fut=s2rc_fut.*p0_2fu%方位压缩
s2rcac_fuf=fftshift(fft(fftshift(s2rcac_fut)))%距离方位压缩后的二维频谱
s2rcac_ut=fftshift(ifft(fftshift(s2rcac_fut).')).'%方位向IFFT
subplot(221)
G=20*log10(abs(s_ut)+1e-6)
gm=max(max(G))
gn=gm-40%显示动态范围40dB
G=255/(gm-gn)*(G-gn).*(G>gn)
imagesc(x,r-Rc,-G),colormap(gray)
grid on,axis tight,
xlabel('Azimuth')
ylabel('Range')
title('(a)原始信号')
subplot(222)
G=20*log10(abs(src_fut)+1e-6)
gm=max(max(G))
gn=gm-40%显示动态范围40dB
G=255/(gm-gn)*(G-gn).*(G>gn)
imagesc(fu,r-Rc,-G),colormap(gray)
grid on,axis tight,
xlabel('Azimuth')
ylabel('Range')
title('(b)距离多普勒域频谱')
subplot(223)
G=20*log10(abs(s2rc_fut)+1e-6)
gm=max(max(G))
gn=gm-40%显示动态范围40dB
G=255/(gm-gn)*(G-gn).*(G>gn)
imagesc(fu,r-Rc,-G),colormap(gray)
grid on,axis tight,
xlabel('Azimuth')
ylabel('Range')
title('(c)RMC后的RD域频谱')
subplot(224)
G=20*log10(abs(s2rcac_ut)+1e-6)
gm=max(max(G))
gn=gm-60%显示动态范围40dB
G=255/(gm-gn)*(G-gn).*(G>gn)
imagesc(x,r-Rc,G),colormap(gray)
grid on,axis tight,
xlabel('Azimuth')
ylabel('Range')
title('(d)目标图象')
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）