㈠ 《崩坏学园2》修改器和修改方法是什么
打开葫芦侠,进入游戏,来到地狱芽衣这,按葫芦侠选择数值搜索,输入9.0000.
搜索9.0000之后有些人可能会搜出很多数据,然后点继续搜索9.0000。
然后修改6555.0000.,如果不想升级那么快可以改小点。(修改65=1W经验)
操作流程:
1.十倍攻击修改:
打开烧饼修改器,在打开游戏,选择2级的枪械大湿或者散弹女王的徽章。(PS:徽章在1关第五张图会掉。没有的同学多刷几次。)然后进图,常规搜索1023530972改成1120390349锁定就可以了。以后都可以十倍攻击。除非你大退游戏。
2.十倍经验修改:
在地图外搜1066192077改1092616192。然后进生存难度的图。(PS:只是改了生存图的经验。普通图无效。)打一次升一级多。 注意低调使用。如果怕封号可以改五倍经验。搜索1066192077改1084887584。
㈡ X射线谱的数据处理
由射线探测器得到的谱脉冲信号经放大器和多道分析器,进入微机系统。计算机通过软件对X射线能谱的常规数据进行处理,从而对样品的含量做定性和定量分析。谱数据处理包括:原始谱数据的光滑;自动寻峰及确定峰位的能量:待测元素的定性分析:峰边界道的确定;峰面积计算等内容。
(1)原始谱数据的光滑预处理
由于测量过程的统计性原理误差,使得X射线谱形式带有统计涨落的特点。所获取的X射线谱统计误差更为明显。在寻峰时,计数的统计涨落可能被误认为是一个谱峰,因此需要对谱先进行光滑预处理。本系统选用二阶多项式五点光滑法。具体计算公式为:式中,ni-2、ni-1、ni+1、ni+2为待光滑的第i道光滑前及光滑后的计数。在进行谱线光滑时,可以重复数次,以达到需要的光滑效果。
(2)自动寻峰?
从获取的X射线谱中找到峰位并换算成相应的能量是X射线谱定性分析的基础。本软件系统选用比较法作为计算机自动寻峰的方法。具体方法如下如果第i道计数满足以下不等式:则认为峰位在i-1,i,i+l道中,再从这三道中选出计数最大的道址即为峰位。k为找峰阈值,一般取值为1—1.5。
(3)系统的能量刻度?
能量刻度是确定本测量系统的道址m与X射线能量E之关系。即:经过自动寻峰得到特征X射线的峰位m,代入(3)式即可求解出该峰位对应的X射线能量。经与元素特征X射线能量库的逐一检索、查证,即可确定待测元素的种类。
(4)峰边界道的确定?
确计算特征峰的净峰面积是定量分析的依据。为此,必须根据实际情况确定特征峰的边界道址。由于采用的半导体探测器具有很好的能量分辨率,对于较高含量的元素,一般都能得到清晰的谱峰。因此,我们采用谱峰的全宽度做为确定本底计数的依据。
(5)净峰面积计算
本系统选用沃森法计算净峰面积。在峰边界L、H内取点F作为计算本底的宽度(F又称为面积因子),β1,β2为左、。
㈢ 什么是奇异谱分析方法
奇异谱分析
主成分分析( PCA, Principal Component Analysis) , 也称为经验正交函数( EOF, Emp irical Orthogonal Function) ,可以由多维的时间序列中获取时间序列的主要成分, 是常用的多元统计分析方法之一, 主要将多个彼此相关的指标变换为少数几个彼此独立的综合指标即主成分, 并要求主成分能反映原始数据的几乎全部信息, 其中, 常用于对一维的时间序列进行分析的方法称为奇异谱分析( SSA, Sin2 gular spectrum analysis) 。 奇异谱方法( SSA) 是一种特别适合于研究周期振荡行为的分析方法, 它是从时间序列的动力重构出发, 并与经验正交函数相联系的一种统计技术, 是EOF分解的一特殊应用。分解的空间结构与时间尺度密切相关, 可以较好地从含噪声的有限尺度时间序列中提取信息, 目前已应用于多种时间序列的分析中。SSA的具体操作过程是, 将一个样本量为n的时间序列按给定嵌套空间维数(即窗口长度) 构造一资料矩阵。当这一个资料矩阵计算出明显成对的特征值, 且相应的EOF几乎是周期性或正交时, 通常就对应着信号中的振荡行为, 可见SSA 在数学上相应于EOF 在延滞坐标上的表达。
对给定的X1 , X2 , �6�8, Xn的时间序列, 给定嵌套维数M, M <N /2, 建立时滞矩阵: S 为对称阵且主对角线为同一常数, 称为Toep litz矩阵, 其特征值为: (4) 式即为序列{ xi } 的奇异谱, 称对X的奇异值运算为奇异谱分析。
最大的特征值对应的特征向量看为第一阶模式, 第二大的特征值对应的特征向量看为第二阶模式, 依次类推。第一阶模式代表了信号的最大变化趋势, 第二阶模式代表了与第一阶模式无关的剩余信号量的最大变化趋势, 依此类推。在实际分析过程中, 通常只选取前面的低阶模式进行分析。计算出S的特征值λk 和相应的特征向量, 序列的SSA展开为: 式中, i = 1, 2, �6�8N 2M + 1; j = 1, 2, �6�8M。Ekj为时间EOF, 记为T2EOF, aik为时间主分量, 记为T2PC: SSA最重要的应用是重建成分RC (Reconstruction ) 。由第k个的T2EOF和T2PC重建xi 的成分记为xki , 公式为 同时也可对各分量进行重建, 用于对原信号的分析。SSA可以提取具有显着振荡行为的信号分量, 并可选择若干有意义的分量进行序列重建。其中低频信号的重建分量, 显示了原始序列的主要演变特征。
㈣ 重、磁数据常规处理方法
所谓常规数据处理是指在重、磁数据处理中经常使用到的位场转换或滤波处理,如上延、求导数、化极等。下面对本次研究中使用的常规数据处理方法做一简述。
空间域异常处理、转换的基本公式均可以写成如下的褶积形式
东北地球物理场与地壳演化
式中:fa、fb分别表示原始异常处理和转换前、后的异常;为权函数。不同的处理、转换间只是它们的权函数不同而已。
由于上述两式与电子电工学中描述滤波器滤波特性的卷积的形式完全相同,因此异常的处理和转换又称为异常的滤波。由卷积定理有
东北地球物理场与地壳演化
式中: 分别为原始异常和处理、转换后异常的波谱; 为权函数φ的波谱(也称为处理与转换因子、波数响应、滤波算子)。
将(4-5)、(4-6)、(4-7)式比较可知,空间域的处理和转换是褶积运算,而波数域是乘积运算。而且,波数谱的连乘可以完成连续的多种变换。因此数域的转换方法要简单得多。随着电子计算机的广泛应用,特别是快速傅里叶变换算法的问世,使区域重磁资料数据处理中的波数域方法成为主要的方法。
在本次研究工作中,数据处理的计算工作是在波数域中进行的。十分明显,为实现波数域的异常处理、转换,必须已知或设计出处理、转换因子 。
根据计算,重、磁异常波谱公式是一些独立因子的乘积,其通式为
东北地球物理场与地壳演化
式中:A为参数因子,只与地质体的剩余密度或剩余质量有关;B为参数因子,只与地质体的磁化强度或磁矩有关;H(ϖ,h)为深度因子,仅与地质体的深度有关;S(ϖ,a,b)为水平尺寸因子,仅与地质体的水平尺寸有关;L(ls,ms,ns,mt,nt)为方向因子,只与磁化强度和地磁场的方向有关;D(ϖ,ξ,η)为位移因子,它是由于坐标原点任选而增加的因子。
(一)解析延拓
根据异常波谱特征的计算,可知无限延深直角棱柱体异常波谱深度因子H=e(-hϖ)。若原始异常体的深度为h1,解析延拓后异常体的深度为h2,则有
东北地球物理场与地壳演化
式中:Δh=h2-h1。
于是 因子延拓因子应为
东北地球物理场与地壳演化
式中:ϖx、ϖy分别为与x轴、y轴的对应的圆波数, 为径向圆波数。
上延时Δh>0,下延时Δh<0。上延可压制高波数成分(即突出低波数成分),属低通滤波;下延可放大高波数成分(即突出高波数成分),属高通滤波,但对低波数成分无压制作用。
(二)导数计算
重、磁异常的导数计算广泛应用于异常的处理和解释。原因在于:①异常的导数在不同形状的地质体上有不同的特征,有助于对异常的解释和分类;②异常导数可以突出反映浅部地质因素,而压制区域深部地质因素的影响,在一定程度上可以划分不同深度和大小异常源产生的叠加异常;③在利用欧拉反褶积对重、磁异常进行构造反演计算时,要用到异常水平一阶导数和垂向一阶导数。
1.垂向导数
由于重、磁异常函数f(x,y,z)的n阶垂向导数可以用下列公式表示
∂nf(x,y,z)/∂n=-∂nf(x,y,h)/∂hn,因此由深度因子H=e(-hϖ)可以求出异常的n阶垂向导数因子。
东北地球物理场与地壳演化
显然一、二阶垂向导数因子分别为:
东北地球物理场与地壳演化
重、磁异常垂直导数可放大高波数成分(即突出高波数成分),但对低波数成分有压制作用。
2.水平导数
由傅里叶变换(FT)的微分性质可知,沿x和y方向的异常水平导数因子分别为
东北地球物理场与地壳演化
式中: 当n取1,即为一阶水平导数的转换因子。如果异常f(x,y,h)对任意水平方向l的导数为: (其中,α为l与x轴的夹角)。依FT的微分性质可得到异常的方向导数因子。
由上式可知,异常的水平导数可突出某一方向上异常特征(或构造线),如α=45°时,能突出135°方向的构造线。
(三)化地磁极
球体总强度磁异常T的谱Δ 为:
东北地球物理场与地壳演化
式中:qs=j(lscosθ+mssinθ)+ns,qt=j(ltcosθ+mtsinθ)+nt;而ls,ms,ns为磁化强度M的方向余弦;lt,mt,nt为地磁场T0的方向余弦;θ为径向圆波数的极角。 为引力位的谱;G为万有引力常数;ρ为密度;ϖ为圆频率;μ0为真空导有磁率。
令qt=nt=1时,即垂直磁异常的谱为:
东北地球物理场与地壳演化
再令qs=ns=1时,即化极后垂直磁异常的谱为:
东北地球物理场与地壳演化
比较以上各式,便可得到化极转换因子为:
东北地球物理场与地壳演化
从(4-18)式看出,化极因子与ϖ无关,因此磁异常化极无滤波作用。另外,从(4-18)式可知,磁异常化极需已知磁化强度方向。然而剩磁与感磁方向不一致时,磁化强度的方向是难以确定的,尤其在大面积磁测资料处理时,区内磁体很多,更无法了解它们的磁化强度方向,因此往往假设磁化强度的方向与地磁场方向一致。另外,实践中还认为,在研究区范围内的地磁场方向是相同的。这种假设在测区不大的情况下对结果的影响较小。图4-2显示出化磁时磁化倾角对结果的影响。根据这一计算结果,若测区的南北纬度差在10°之内,用统一的地磁场方向余弦来作化磁极运算,其结果受影响不大。因此,对一般成矿预测为目的的区域性资料研究问题并不突出。但当作深部地质研究或大面积区域地质研究时,就要注意这种情况,有些学者已开始了测点地磁场方向余弦各异时的化磁研究工作。
图4-2 不同磁化方向化磁极后的曲线对比
(四)谱分析方法
谱分析方法作为重、磁异常数据处理、转换的重要方法,有着广泛的应用。利用径向平均对数能谱分析可以估算重、磁场源的平均深度,为进一步的处理和解释提供基础信息。
下面对径向平均对数能谱分析和平均深度估算的原理进行简介:
经计算可知,球体重力异常的波谱为:
东北地球物理场与地壳演化
则球体重力异常功率谱为:
东北地球物理场与地壳演化
对数径向功率频谱为:
东北地球物理场与地壳演化
式中A=2πGρv=2πGm。上式表明,球体重力异常对数径向功率谱与径向圆波数中呈线性关系,见图4-3a,故可利用lnE(ϖ)的拟合直线斜率求解出球体中心深度:
东北地球物理场与地壳演化
式中:f为波数,而ϖ=2πf。
另外,由球体磁异常功率谱也可计算深度。把引力位谱( =2πGme-hϖ/ϖ)代入(4-15)式,经整理可求得球体垂直磁化(qs=1),垂直磁异常(qt=1)的功率谱为:
东北地球物理场与地壳演化
式中B=2πμ0VM/4π=2πμ0m/4π;m为磁矩;V为球体体积;M为磁化强度。
东北地球物理场与地壳演化
上式说明,球体垂直磁化垂直磁异常的对数径向功率谱与圆波数呈非线性关系(图4-3b)。但是,在高波数段近于线性关系,可用(4-21)式计算深度。
图4-3 对数径向功率谱
(五)重、磁对应分析
基于泊松定理发展起来的重磁异常对应分析方法,是重磁数据综合解释的重要方法,能对重磁异常的相关性进行定量研究,有效地将重磁信息进行综合,对重磁资料定性地赋予地质意义,并突出地质目标的反映,为重磁资料的地质解释提供有用的信息,特别是在强磁性火山岩解释中具有重要的作用。
重磁异常对应分析方法的基本理论如下:
由同一场源引起的重力异常和磁异常间的关系可以简单地用泊松方程描述。当垂直磁化时,泊松方程可表示为:
东北地球物理场与地壳演化
式中:Δz⊥为垂直磁化的垂直磁异常;M为场源磁化强度;G为万有引力常数;Δρ为场源剩余密度;Δg为重力异常; 为重力异常的垂向一阶导数。
上式表明垂直磁化的垂直磁异常与重力场的垂向一阶导数满足线性关系,而且拟合直线的截距为零。
由于原始资料不可避免地存在某些干扰因素,通常进行重磁异常的线性回归分析时,选用如下稍加推广的泊松方程:
东北地球物理场与地壳演化
式中:b为斜率;A为截距。
将Δz⊥与 作线性回归分析则可得到斜率b与截距A的估计值。两个离散序列的相关导数可以由下式求得:
东北地球物理场与地壳演化
式中:Cxy(k) 为两个离散序列x(t)={x1,x2,…,xn}和y(t)={y1,y2,…,yn}的相关函数,k为延迟时间。当x(t)=y(t)时,称为自相关函数Cxx(k)或Cyy(k)。
计算处理时,给定适当大小的分析窗口,将窗口内各点垂直磁化磁异常和重力异常的垂向一阶导数进行最小二乘线性回归,求得中心点的相关系数R、斜率b和截距A。
相关系数R反映了在给定窗口内重磁异常的线性相关程度,即宏观地反映了重磁异常的“同源性程度”。相关系数绝对值接近于1的窗口区间重磁异常的“同源性”较好,它们或者同源、或者都离场源较远、或者同处异常的拐点等。其中R接近+1时,重磁异常正相关;R接近-1时,重磁异常负相关。当R绝对值较小时,重磁异常相关性差,重磁异常可能不同源,或存在邻近异常干扰,或是存在方向不同于地磁场的强剩磁磁性体等。
斜率b反映了所有场源泊松比的加权平均值,称为广义泊松比。只有在重磁异常同源的前提下,回归所得的斜率b才有意义。仅由b不能直接确定M和Δσ,但若在解释中结合其他地质、地球物理信息,就能从中获得关于物性分布的有用信息,从而为进一步的定量解释提供依据。
截距A反映了实测资料中的长波长成分,它主要反映重磁异常数据的背景变化。在重磁异常完全同源的理想情况下,A=0。
由于重磁异常对应分析是对场源之间的相关系数进行定性和半定量研究的方法,它能分离和鉴别不同类型的异常,从而勾画出与异常场源相对一致的地质单元和构造分区,不相关区说明重磁异常不同源或存在邻近异常干扰。
(六)欧拉反褶积与构造反演
欧拉反褶积方法使用欧拉(Euler)齐次关系,对经方向谱分析过的数据快速估计重、磁场源的位置和深度,是一种既能够利用重磁网格数据,又对剖面数据有效地确定地质体位置(边界)和深度的定量反演方法(Reid等,1990)。这种方法并不需要已知地质信息(密度、磁化率等)的控制。使用该方法可以将位场及其梯度以及场源位置之间的关系用欧拉齐次方程表示,而场源的不同形状即地质构造的差异则表现为方程的齐次程度,就是所谓的地质构造指数,地质构造指数或齐次程度实质上表现了场随离开场源距离的衰减率。模型研究和应用实例表明,这个方法对确定断层、磁性接触带、岩脉、喷出岩体等构造位置或勾绘它们的轮廓有较高的精度。
位场的欧拉方程是由Thompson推导的。首先建立一个直角坐标系,取观测平面为z=0,z轴向下为正,x轴指北,y轴指东。考虑在此坐标系中的任一函数f(x,y,z),如果
东北地球物理场与地壳演化
则称函数是n阶齐次的。此外可证明,如果f(x,y,z)是n阶齐次的,则满足下列方程
东北地球物理场与地壳演化
此偏微分方程称为欧拉齐次方程,或称欧拉方程。
对于位于(x0,y0,z0)的点磁源,在观测平面上任一点(x,y,z)处的总磁场强度具有如下形式:
东北地球物理场与地壳演化
式中 N=1,2,3,……。G不依赖于(x,y,z)。对于(3-23)这样的函数,其欧拉方程可写成
东北地球物理场与地壳演化
方程(4-29)是n=-N阶齐次的。三个坐标方向的梯度值可以利用空间域或波数域的一般位场变换计算出来。如果梯度值通过观测获得,直接用于方程(4-29)则更可取。
方程(4-29)虽然是根据磁源异常推导的,但对重力异常也同样适用。该方程用于平面网格重、磁异常数据的反演计算。如果假定方程中横向梯度∂ΔT/∂y为零,则可得到适用于剖面数据计算的方程。这对于众多走向方向不变的二度情况很显然就是这样。
齐次度N被定义为“构造指数”,它是重、磁异常场源深度变化“陡缓”的量度。特定的地质构造具有特定的衰减率(即:构造指数)。例如:倾斜断层的磁场、水平薄岩脉的磁场按线性的规律变化,构造指数就为1。表4-1列出了构造指数对应的地质构造。
利用不同坐标点(x,y,z)上的场值ΔT及其三个方向上的梯度值以及方程(4-29)组成的线性方程组,最后可以解出未知变量x0、y0、z0,进而确定构造形迹及位置。
表4-1 欧拉构造指数表
但是,直接用方程(4-29)及其变换的二度形式解决构造问题,会使解的精度极不可靠和不稳定。主要原因有如下几个方面:
(1)很难知道磁场ΔT的绝对水平,区域场或邻近磁异常的影响几乎总是存在的。
(2)根据线性方程组与系数的关系,较低的构造指数才会有较好的深度估计值。但大多数磁异常是偶极性的,有较高的构造指数。同时又有许多线性构造的指数接近于零而使反褶积发散。
(3)实测异常是多种构造指数特征的复杂叠加,很难用一些简单模型来模拟,亦很难将具有线性特征的构造识别与分离出来。
为克服上述三个方面的问题,釆用下面的一些办法:
从观测数据中消除偏差是通过网格数据进行窗口计算解决的。对网格数据假定异常在方程(4-29)求值的窗口范围内有一常量偏差,观测值为
东北地球物理场与地壳演化
这里B是常数。从方程(4-30)中解出ΔT,代入(4-29),整理得
东北地球物理场与地壳演化
如果构造指数小于0.5,即构造指数接近零时,这样可能造成对深度值的过低估计。为此需要提供一个补偿值A,使得欧拉齐次方程在构造指数较低时写为
东北地球物理场与地壳演化
式中:A是与场幅值有关的一个参数。不同的构造形体有不同的A,A可以通过将已知的某一构造的解析式代入欧拉方程(4-29)而求出。
图4-4 重力异常欧拉反褶积计算结果示意图
图4-4中的曲线为重力异常等值线,圆圈为反演解的构造位置。圆圈直径的大小代表了不同的构造深度。
(七)重、磁人机交互剖面正反演
该项技术的优点是便于将重、磁异常的处理、转换方法得出的结果和其他地质、地球物理方法获取的先验信息输入到模型里,形成初始模型。并且根据计算结果和实际重、磁异常的差异,随时方便地修改模型,直观地监督和指导正反演过程。重、磁人机交互剖面正反演流程见图4-5。
图4-5 人机交互正反演流程图
1.重力人机交互正反演技术
重力人机交互正反演技术(Gamble,1979)主要是依据A截面为多边形的二度体重力异常计算方法来实现的。通过对初始模型计算出的重力效应与测线上的布格重力异常进行对比,不断修正模型,直至达到计算出的重力效应与测线上的布格重力异常之差满足预定精度。重力人机交互正反演流程见图4-5。
图4-6 二度体A截面
A截面为多边形的二度体重力异常计算方法:
假设二度体的剩余密度为σ,以计算点作为坐标原点,x轴与二度体走向垂直,z轴铅垂向下(图4-6)。若n边形第k个顶点的坐标为(ξk,ζk),其中k=1,2,...,n。则(ξk,ζk)与(ξk+1,ζk+1)两个顶点连线上ξ与ζ有如下关系:
东北地球物理场与地壳演化
引用解Δg正问题的基本公式,首先对ξ求积分,得
东北地球物理场与地壳演化
式中:s为多边形的A截面积;l为A截面的周长。
将式(4-33)代入(4-34)得
东北地球物理场与地壳演化
对上式积分可得到如下结果
东北地球物理场与地壳演化
或写成下面的形式
东北地球物理场与地壳演化
(4-36)式或(4-37)可以编成计算机程序,用以计算A截面为任意多边形的二度体的重力异常,进而可以进行重力人机交互正反演。在具体编程计算时应注意以下几个问题:
(1)因多边形的边数为n,故ξn+1=ξ1,ζn+1=ζ1;
(2)(4-36)和(4-37)两式是假设计算点位于原点时导出的,因此,当任意计算点P(x,y)的重力异常时,式中的ξk和ζk应以ξk-x和ζk-z来代替;
(3)在(4-36)和(4-37)式中,反正切函数的取值范围应在-π到π之间,即当ξk+1>ξk时,反正切函数在0到π之间取值;反之,则在-π到0之间取值。
2.磁法人机交互正反演技术
磁法人机交互正反演技术主要是依据A截面为多边形的二度体磁力异常计算方法来实现的。其基本思想同重力人机交互正反演技术相一致。
由于V2=Δg,所以根据公式(4-37)可以求出引力位的二阶导数
东北地球物理场与地壳演化
将(4-38)和(4-39)式代入下面二度体的磁异常公式,就可以利用该式进行磁力人机交互正反演计算。
东北地球物理场与地壳演化
式中:μ0为真空的磁导率;MS为有效磁化强度;is为有效磁化倾角;I0为地磁场倾角;A'为x轴与磁北的夹角。在具体编程序上机计算时应注意的问题等方面与重力人机交互方法相同(图4-5)。
㈤ 检验检测新方法确认的内容包括哪些
方法确认具体来讲包括了标准方法的证实和非标方法的确认两个方面。
我们先来看看方法确认和方法证实的目的是什么:
非标准方法确认目的:
该非标准方法能否合理、合法使用;
标准方法证实目的:
实验室是否有能力按标准方法开展检测、校准活动。
接下来我们看看标准方法证实和非标方法的确认应该如何做:
标准方法证实:
从人、机、料、法、环、测几个方面去证实实验室有能力满足标准方法的要求,有能力开展检测、校准活动。
证实的内容要从七方面去做:
(1) 对执行新标准所需的人力资源的评价,即检测、校准人员是否具备所需的技能及能力;必要时应进行人员培训,经考核后上岗;
(2) 对现有设备适用性的评价,诸如是否具有所需的标准、参考物质,必要时应予补充。
(3) 对设施和环境条件的评价,必要时进行验证。
(4) 对物品制备,包括前处理、存放、辅助试剂等各环节是否满足标准要求的评价。
(5) 对作业指导书、原始记录、报告格式及其内容是否适应标准要求的评价。
(6) 对新旧标准进行比较,尤其是差异分析与比对的评价。
(7) 按标准要求进行完整模拟检测,出具完整结果报告。
标准方法证实应有相关的文件规定及其证实的记录,标准方法变更后应重新证实。
非标准方法确认:
一个非标方法的确认,在文件中要包括以下内容:
a) 方法适当的标识;
b) 方法所适用的范围;
c) 检测或校准样品是什么类型,以及对样品的描述;
d) 被测参数的范围;
e) 方法对仪器、设备的要求,包括仪器设备关键技术性能的要求;
f) 需要用到的的标准物质;
g) 方法对环境条件的要求,对环境稳定周期的要求;
h) 操作步骤,包括:
—样品的标志、处置、运输、存储和准备;
—检测、校准工作开始前需要进行的检查;
—检查设备工作是否正常,需要时,使用前之前对设备进行校准和调整;
—结果的记录方法;
—安全注意事项;
i) 结果接受(或拒绝)的准则、要求;
j) 需记录的分析数据;
k) 不确定度评定。
非标方法的技术确认,需要从五个方面确认:
(1)使用参考标准或标准物质进行比较;
(2)与其他方法所得的结果进行比较;
(3)实验室间比对;
(4)对影响结果的因素作系统性评审;
(5)根据对方法的理论原理和实践经验的科学理解,对所得结果不确定度进行的评定。
技术确认要尽可能全面,并需有确认记录。
㈥ 基于不同数据源的土地利用变化遥感动态监测方法
李翔宇 樊彦国
(中国石油大学地球资源与信息学院,山东东营,257061)
摘要:本文从所拥有的遥感数据源的可能情况出发,分别介绍了各种情况下利用遥感进行土地利用变化动态监测的方法,分析了其优势和劣势。
关键词:遥感;土地利用变化;动态监测;方法
1 引言
我国是一个人多地少的国家,土地是我们赖以生存的资源。建立土地动态监测系统以快速准确地提供各类土地资源面积及其分布、土地资源动态变化状况及土地资源生态环境信息是十分必要的,这样可以保证我国在科学翔实的资料基础上对土地资源进行科学的规划及合理的利用,实现土地资源的可持续健康发展。可是传统的统计或实地调查方式,耗时耗力,劳民伤财,并且难以适应土地利用的快速变化,而遥感可以提供及时准确且覆盖面广的地面影像资料,并且周期短、信息量大,通过后期的分析、处理、比较,可以使人们迅速准确地掌握土地利用变化的详细信息,即实现土地利用的动态监测。现在,遥感技术已成为进行土地利用变化动态监测的重要手段。
基于遥感影像的土地利用变化监测方法大致可分为两类:光谱直接比较法和分类结果比较法。多数变化提取算法属于前一种,主要包括影像差值法、比值法、主成分分析法和变化矢量分析法等,这些算法直接通过两时相数据的光谱差异确定变化发生的区域,但不能得出变化图斑的类型;后一种方法通过对各自时相的数据进行土地利用分类,通过对两个分类结果的比较提取变化信息,但其精度受两时相数据分类精度的制约。实际操作中可以根据所持有数据源的不同而采用相应的方法。
2 基于单一传感器的土地利用变化监测方法
2.1 基于单一传感器多时相遥感影像
当遥感数据源为单一传感器但可以获得多时相遥感影像时,可以考虑以下几种方法。
2.1.1 单变量图像差值法[1]
单变量图像差值法比较简单,是使用最广泛的一种探测方法。它是将两个时相的遥感图像按波段进行逐像元相减,从而生成一幅新的代表二时相间光谱变化的差值图像。辐射值的显着变化代表了土地覆盖变化,在差值图像中接近于零的像元就被看做是未变化的,而那些大于或小于零的像元表示其覆盖状况发生了某种变化,从而设定适当的阈值就可以把变化信息提取出来。
2.1.2 图像比值法[1,2]
比值处理被认为是辨识变化区域相对较快的手段。它是对于两个时相多谱段数据中同名像元的光谱灰度值施以除法运算。显然,经过辐射配准后,在图像中未发生变化的像元其比值应近似为1,而对于变化像元而言,比值将明显高于或低于1。比值法可以部分地消除阴影影响,突出某些地物间的反差,具有一定的图像增强作用。
2.1.3 图像回归法[1]
图像回归法是首先假定时相Ⅰ的像元值是另一时相Ⅱ像元值的一个线性函数,通过最小二乘法来进行回归,然后再用回归方程计算出的预测值来减去时相Ⅰ的原始像元值,从而获得两时相的回归残差图像。
2.1.4 植被指数差值法[2]
植被指数差值法是用近红外与红光波段间的比值(植被指数)代替原始波段作为输入数据进行差值运算来生成变化图像。由于植物普遍对红光强烈吸收和对近红外光强烈反射,因此红光和近红外波段之间的比值有利于提高光谱差异。
2.1.5 主成分分析法[3]
(1)差异主成分法 两时相的影像经纠正、配准之后,先对影像作相差取绝对值处理,从而得到一个差值影像。差值影像作主成分变换之后的第一分量应该集中了该影像的主要信息,即原两时相影像的主要差异信息。这个分量可以被认为是变化信息而被提取出来,从而生成变化模板,作为指导下一步变化类型确认和边界确定的参考信息。
(2)多波段主成分变换 由遥感理论可得知,地物属性发生变化,必将导致其在影像某几个波段上的值发生变化,所以只要找出两时相影像中对应波段上值的差别并确定这些差别的范围,便可发现土地利用变化信息。在具体试验中将两时相的影像各波段进行组合,成一个两倍于原影像波段数的新影像,对该影像作主成分变换。由于变换结果前几个分量上集中了两个影像的主要信息,而后几个分量则反映出了两影像的差别信息,因此可以抽取后几个分量进行波段组合来产生出变化信息。一般说来,在上述多波段主成分变换之后,采用0、1、2分量进行波段组合能较好地反映出新旧时相影像的变化部分。
(3)主成分差异法 本方法和差异主成分方法所不同之处在于影像作主成分变换与差值处理的顺序不一样。要求先对两时相的影像作主成分变换,然后对变换结果作差值,取差值的绝对值为处理结果。在实际的试验中,两时相影像作主成分变换后相差的第一分量已经涵盖了几乎所有的变化信息。因此,可以认为这一分量属于影像的变化信息。
2.1.6 变化向量分析法[1]
由于多时相遥感数据中任一像元矢量都可用多维测量空间中的一个点来表示(空间的维数等于原始波段数),通过对不同时相下的同名像元矢量进行相减所得到的变化矢量就可以用于描述该像元第一时相 t1 到第二时相 t2 期间在多维空间中所发生的位置变化。其中变化矢量的模代表了变化的强度,而方向则指示了发生变化的类型。设时相 t1、t2 图像的像元灰度矢量分别为 G=(g1,g2,…,gk)T 和H=(h1,h2,…,hk)T,则变化矢量为:ΔG=G -H。ΔG 包含了两幅图像中所有变化信息。变化强度由变化矢量的模||ΔG||决定,||ΔG||越大,表明图像的差异越大,变化发生的可能性越大。因此,提取变化和非变化像元,可根据变化强度||ΔG||的大小设定阈值来实现,即像元||ΔG||超过某一阈值时,即可判定为土地利用类型发生变化的像元;而变化的类型,可由ΔG的指向确定。
这种方法利用多频段信息,在提取变化位置的同时可以得到变化类型信息,是一种较理想的算法。当然,要用好变化向量分析法还取决于分析过程中变化/未变化阈值是否取值合理以及相关分类方法是否适当。
2.1.7 分类后比较法
分类后比较法是对两期遥感影像进行监督或非监督分类,然后比较在各图像系列同一位置上的分类结果,进而确定土地利用类型变化的位置和所属类型。该方法可直接获得变化类型信息,但如何选择合适的分类方法提高分类精度是准确获得变化信息类型的关键。
2.1.1至2.1.6均属于光谱直接比较法,此方法对变化比较敏感,可以避免分类过程所导致的误差,但需要进行严格的辐射标准化,排除大气状况、太阳高度角、土壤湿度、物候等“噪声”因素对图像光谱的影响,由于目前对各种干扰(尤其是物候)导致的辐射差异的校正方法仍不成熟,因此,只能通过选择同一传感器、同一季相的数据来尽可能减小“噪声”。同时光谱直接比较法只注重变化像元的提取,而不能提供变化中土地类型的转化信息(如地类属性)。与之相对照,分类后比较法对辐射纠正要求相对较低,适用于不同传感器、不同季相的数据的比较,同时该方法不仅可以提供变化信息,而且还能够给出各时期的土地利用类型信息。但这种方法的最终精度受到影像分类精度的限制,而且它对影像的全部范围都要进行分类计算而不管它们是否已经发生变化,这样无疑大大增加了变化信息检测的计算量。
在目前的土地利用遥感监测研究中,结合光谱直接比较法和分类后比较法的混合动态监测方法逐渐受到重视,并有了一些成功的案例研究。Jenson 通过对湿地变化的动态监测研究表明:先利用光谱直接比较探测变化区,再进行图像分类确定变化类型的混合法是一种非常有效的变化检测方法[4];Macleod和Congalton的研究也表明以差值法为基础的混合动态监测法优于传统分类后比较法[5]。这样可以集两者之所长,取得更好的监测效果。
2.2 基于单一传感器单时相遥感影像
无论是光谱直接比较法还是分类后比较法都是基于多个时相的遥感影像来进行土地利用变化监测。而当前期遥感影像无法或者难以获得的情况下,依靠后期的单时相遥感影像与前期的土地利用现状图也可以进行动态监测,这就是采用将土地利用现状图叠加在遥感图像上的方法来监测土地利用变化情况[6]。具体说来,是利用土地利用现状图中不变的明显地物标志(如线状地物交叉点)作为控制点对遥感图像进行配准,然后将土地现状图叠加再校正后的遥感图像上,检查各图斑是否吻合,若图斑的角点有偏移,则发生变化。可通过遥感图像辨识当前的土地利用类型,而土地利用现状图含有先期的土地利用类型信息,所以可以比较容易地辨识土地利用类型的变更情况,并可测算出变化图斑的面积。若其中有不能确定的图斑,可以辅以外业调查,以提高监测精度。
3 基于多源遥感的土地利用变化信息监测方法
不同传感器都具有各自的优势,获得的图像各有所长,如美国陆地卫星(Landsat)TM图像光谱信息丰富;法国SPOT卫星图像具有全色通道而空间分辨率高;SAR图像不受光照条件的影响而且几乎不受大气和云层的干涉,可用于探测地物的复介电常数和表面的粗糙度等等。利用不同传感器的多源遥感影像进行融合,可以使其优势互补,在此基础上的土地利用变化动态监测已成为国际遥感界研究的主题之一。以TM影像和SPOT影像为例,目前应用多光谱TM和全色SPOT数据融合的方法主要有LAB变换、HIS变换、线性复合与乘积运算、比值运算、BROVEY 变换、高通滤波变换(HPH)和主成分分析(PCA)等方法[7],经上述算法融合后的图像可以有效地同时保留SPOT高分辨率图像的精细纹理和TM多光谱图像的丰富色彩信息,从而有利于提高图像的空间分辨率和光谱分辨率,为发生变化的地类图斑的提取提供良好的数据源基础。
3.1 光谱特征变异法[8]
针对基于多源遥感的土地利用变化监测,变化信息的提取方法除了2.1所述方法之外还可以选择光谱特征变异法。
同一地物反映在SPOT影像上的信息是与其反映在TM影像上的光谱信息一一对应的。因此作TM和SPOT影像融合时,才能如实地显示出地物的正确光谱属性。但如果两者信息表现为不一致时,那么融合后影像的光谱就表现得与正常地物有所差别,此时就称地物发生了光谱特征变异(例如同一位置,前期在遥感影像上呈现为绿色的麦地,后期新修道路在影像上呈现较亮的灰度,那么叠加之后会呈现一条绿色的道路,与正常地物相异),这部分影像在整个的影像范围内是不正常和不协调的,这些地物可以通过影像判读的方法勾绘出来,这种变化信息提取的方法具有物理意义明显、简洁的特点。但是经过试验发现,发生光谱特征变异的地物在几何尺寸上要足够的大才能被人工目视发现。此外,该方法的效率还受到被监测区地物光谱特性的限制。
3.2 变化信息提取方法的选择
根据土地利用动态监测项目所获取的数据源,可将遥感数据组合分为下述几种类型,针对不同的类型要采取相应的方法以获取较好的效果。
3.2.1 具有两时相的 TM 和 SPOT 数据
这种情况是最好的。在该条件下,先对两时相的数据以某一纠正后的TM或SPOT影像(首先处理TM还是SPOT视数据的具体情况而定,原则是利于TM和SPOT数据的配准融合处理)为参考分别作纠正和配准处理,为保留并结合原始数据中纹理信息和光谱信息要融合相对应的TM和SPOT影像,在两时相融合影像的基础上采用主成分差异的方法来提取变化信息。另外还可以用新时相的 SPOT 影像与旧时相的 TM 影像进行融合生成光谱特征变异影像来指导发现变化的区域。
3.2.2 具有两时相的 TM 和一个时相的 SPOT 数据
在此数据源的基础上,首先仍对某一时相的TM或SPOT数据作纠正处理,然后将其他时相的TM和SPOT数据都统一以这个纠正后的TM (SPOT)为参考影像作影像到影像的纠正和配准。之后,选择光谱特征变异的方法来寻找大部分的变化信息,借助于两时相的TM影像确认变化;此外,利用主成分分析的办法对两时相的TM数据进行处理,得到变化信息模板,将模板叠置在判读影像上补充单一方法进行变化提取的遗漏。
3.2.3 具有两时相的 SPOT 和一个时相的 TM 数据
通常,前面的数据预处理纠正配准部分同3.2.2相同,然后对其中交错时相的TM和SPOT数据进行融合得到光谱特征变异影像,借助于两时相的SPOT数据发现影像中纹理信息的变化,从而辅助提取影像中的变化信息部分。除此之外,两时相的SPOT影像数据理论上说,可以直接作比较得到变化的部分,但是由于成像条件的不同,这样直接比较的方法会导致产生很多伪变化信息,干扰了真正变化部分的提取。因此,首先要对原始SPOT影像进行去噪及辐射校正等预处理,然后才能用来提取变化的信息。
3.2.4 具有单时相的 SPOT 影像和另一时相 TM 影像的数据
首先要对SPOT和TM数据进行纠正处理,然后利用纠正后的SPOT和另一时相TM影像融合得到光谱特征变异影像,并以此作为判读变化信息的主要参考数据。此外,单时相的SPOT数据可以作为新增波段加入到原始的 TM 数据中去进行主成分分析来提取变化的信息,辅助发现漏判的变化图斑。
利用遥感进行土地利用动态监测的方法非常多,这些方法各有自己的优势和劣势,实际工作中,要针对所拥有的数据源的情况,综合各方面要求来选择合适的方法,也可以综合几种方法取长补短以达到更好的监测效果。至于如何更有效地识别土地变化的类型以及如何提高分类的精度仍有很大的研究空间。
参考文献
[1]卢珏.土地利用动态监测变化信息提取算法评估[J].湖北农学院学报,2002,22 (5):394~396
[2]张银辉,赵庚星.试论土地利用遥感动态监测技术方法[J].国土资源管理,2001,18 (3):15~18
[3]杨贵军,武文波,陈步尚,夏春林.土地利用动态遥感监测中变化信息的提取方法[J].东北测绘,2003,26 (1):18~21
[4]Jensen J R,Cowen D J,Narumalani S,et al.An evaluation of coast watch change detection protocol in South Carolina [J].Photogram metric Engineering and Remote Sensing,1993,59 (6):1039~1046
[5]Macleod R D,Congalton R G.A quantitative comparison of change-detection algorithms for monitoring eelgrass from remotely sensed data [J].Photogram metric Engineering and Remote Sensing,1998,64 (3):207~216
[6]吴连喜,严泰来,张玮,薛天民,程昌秀.土地利用现状图与遥感图像叠加进行土地利用变更监测[J].农业工程学报,2001,17 (6):156~160
[7]张炳智,张继闲,张丽.土地利用动态遥感监测中变化信息提取方法的研究[J].测绘科学,2000,25 (3):46~50
㈦ 常用的数据分析方法有哪些
常见的数据分析方法有哪些?
1.趋势分析
当有大量数据时,我们希望更快,更方便地从数据中查找数据信息,这时我们需要使用图形功能。所谓的图形功能就是用EXCEl或其他绘图工具来绘制图形。
趋势分析通常用于长期跟踪核心指标,例如点击率,GMV和活跃用户数。通常,只制作一个简单的数据趋势图,但并不是分析数据趋势图。它必须像上面一样。数据具有那些趋势变化,无论是周期性的,是否存在拐点以及分析背后的原因,还是内部的或外部的。趋势分析的最佳输出是比率,有环比,同比和固定基数比。例如,2017年4月的GDP比3月增加了多少,这是环比关系,该环比关系反映了近期趋势的变化,但具有季节性影响。为了消除季节性因素的影响,引入了同比数据,例如:2017年4月的GDP与2016年4月相比增长了多少,这是同比数据。更好地理解固定基准比率,即固定某个基准点,例如,以2017年1月的数据为基准点,固定基准比率是2017年5月数据与该数据2017年1月之间的比较。
2.对比分析
水平对比度:水平对比度是与自己进行比较。最常见的数据指标是需要与目标值进行比较,以了解我们是否已完成目标;与上个月相比,要了解我们环比的增长情况。
纵向对比:简单来说,就是与其他对比。我们必须与竞争对手进行比较以了解我们在市场上的份额和地位。
许多人可能会说比较分析听起来很简单。让我举一个例子。有一个电子商务公司的登录页面。昨天的PV是5000。您如何看待此类数据?您不会有任何感觉。如果此签到页面的平均PV为10,000,则意味着昨天有一个主要问题。如果签到页面的平均PV为2000,则昨天有一个跳跃。数据只能通过比较才有意义。
3.象限分析
根据不同的数据,每个比较对象分为4个象限。如果将IQ和EQ划分,则可以将其划分为两个维度和四个象限,每个人都有自己的象限。一般来说,智商保证一个人的下限,情商提高一个人的上限。
说一个象限分析方法的例子,在实际工作中使用过:通常,p2p产品的注册用户由第三方渠道主导。如果您可以根据流量来源的质量和数量划分四个象限,然后选择一个固定的时间点,比较每个渠道的流量成本效果,则该质量可以用作保留的总金额的维度为标准。对于高质量和高数量的通道,继续增加引入高质量和低数量的通道,低质量和低数量的通过,低质量和高数量的尝试策略和要求,例如象限分析可以让我们比较和分析时间以获得非常直观和快速的结果。
4.交叉分析
比较分析包括水平和垂直比较。如果要同时比较水平和垂直方向,则可以使用交叉分析方法。交叉分析方法是从多个维度交叉显示数据,并从多个角度执行组合分析。
分析应用程序数据时,通常分为iOS和Android。
交叉分析的主要功能是从多个维度细分数据并找到最相关的维度,以探究数据更改的原因。
㈧ 常用的数据处理方法
前面所述的各种放射性测量方法,包括航空γ能谱测量,地面γ能谱测量和氡及其子体的各种测量方法,都已用在石油放射性勘查工作之中。数据处理工作量大的是航空γ能谱测量。
(一)数据的光滑
为了减少测量数据的统计涨落影响及地面偶然因素的影响,对原始测量数据进行光滑处理。消除随机影响。
放射性测量数据光滑,最常用的光滑方法是多项式拟合移动法。在要光滑测量曲线上任取一点,并在该点两边各取m个点,共有2m+1点;用一个以该点为中心的q阶多项式对这一曲线段作最小二乘拟合,则该多项式在中心点的值,即为平滑后该点的值。用此法逐点处理,即得光滑后的曲线,光滑计算公式(公式推导略)为
核辐射场与放射性勘查
式中:yi+j、为第i点光滑前后的值;为系数;为规范化常数。
五点光滑的二次多项式的具体光滑公式为
核辐射场与放射性勘查
如果一次光滑不够理想,可以重复进行1~2次,但不宜过多重复使用。
光滑方法,还有傅里叶变换法,以及多点平均值法,多点加权平均值法等。
使用那种方法选定之后,一般都通过编程存入计算机,进行自动化处理。
图7-2-1是美国东得克萨斯州一个油田上的航空γ放射性异常中的两条剖面图(A-B和B-C)。经过光滑处理后,低值连续,清晰明显,与油田对应的位置较好。说明四个油藏都在铀(w(U))和钾(w(K))的低值位置。
图7-2-1 美国东得克萨斯油田航空γ放射性异常剖面图
(二)趋势面分析方法
趋势分析主要反映测量变量在大范围(区域)连续变化的趋势。在原始数据中常含有许多随机误差和局部点异常,直观反映是测量曲线上下跳动或小范围突变。使用趋势分析处理是为了得到研究区域辐射场的总体分布趋势。
趋势面分析,实质上是利用多元回归分析,进行空间数据拟合。根据计算方法不同,又可分为图解法趋势面分析和数学计算法趋势面分析。图解法趋势面分析的基本思路是对观测数据采用二维方块取平均值法,或滑动平均值法计算趋势值。方块平均值法是对每一方块内的数据取平均值,作为该方块重心点的趋势值。滑动平均值法是设想一个方框,放在测区数据分布的平面图上,把落在方框内的测点数据取平均值,记在方框中心上,最后得到趋势面等值图。一般讲做一次是不够的,需要如此重复3~9次。一般都有专门程序可供使用(不作详述)。如图7-1-14(a)为原始数据等值图,中间有许多呈点状高值或低值分布,经过四次趋势面分析之后可以清楚地看出三个低值异常区。
计算法趋势面分析是选定一个数学函数,对观测数据进行拟合,给出一个曲线。拟合函数常用的有多项式函数,傅里叶级数,三角函数以及指数函数的多项式函数等。目前以二维多项式函数应用最多。
(三)岩性影响及其校正分析
不同岩石、不同土壤中放射性核素含量是有差别,有的相差还比较大,有的相差甚至超过10%~20%。这是油田放射性测量的主要影响因素。
一个测区可能出现不同土壤分布,把不同放射性水平的土壤上测量结果校正到同一水平(叫归一化方法)是非常重要的工作,主要有下面三种方法。
1.确定土壤核素含量的归一化方法
利用γ能谱测量资料,根据测区地质图或土壤分布图,分别统计总道的总计数率和铀、钍、钾含量的平均值。然后进行逐点校正,即逐点减去同类土壤的平均值,其剩余值即为异常值。
核辐射场与放射性勘查
式中:分别为第 i类土壤中测点 j的总计数和铀、钍、钾含量。分别为i类土壤的平均总计数和铀、钍、钾的平均值。分别为扣除各类土壤平均值后的剩余值,即为各测点不同土壤校正后的归一化的油田的放射性异常。根据需要可以用来绘制平面剖面图或等值线图,即为经过不同岩性(土壤)校正后的油田放射性异常图。
这个方法的缺点是计算工作量较大。
2.用钍归一化校正铀、钾含量
对自然界各种岩石中的钍、铀、钾含量的相关性研究(D.F.Saundr,1987),发现它们的含量具有很好的相关性(表7-2-2);而且随岩性不同含量确有相应的增加或减小,据此可以利用钍的含量计算铀和钾的含量。钍有很好的化学稳定性,钍在地表环境条件下基本不流失。因此,利用钍含量计算出来的铀、钾含量,应当是与油藏存在引起的铀、钾
表7-2-2 几种岩石的钍、铀、钾含量
异常无关的正常值。用每点实测的铀、钾,减去计算的正常值,那么每个测点的铀、钾剩余值(差值)应当是油气藏引起的异常值。这样就校正了岩性(土壤)变化的影响。
对于航空γ能谱测量的总道计数率,也同样可以用钍含量(或计数率)归一化校正总道计数率,效果也非常好。
具体方法如下。
1)对铀、钾的归一化校正。
2)根据航空γ能谱测量或地面γ能谱测量数据,按测线计算铀、钍、钾含量。根据岩石(土壤)中钍与铀,钍与钾的相关关系(表7-2-1),认为铀和钍存在线性关系,钾和钍存在对数线性关系,于是建立相应的拟合关系式。
核辐射场与放射性勘查
式中:A、B、A′、B′为回归系数(对每个测区得到一组常数);wi(Th)为测点i实测的钍含量;w点i(U)、w点i(K)为i点由钍含量计算的铀、钾含量。
计算每个测点的铀、钾剩余值:
核辐射场与放射性勘查
式中:wi(U)、wi(K)为测点i的实测值。剩余值Δwi(U)和Δwi(K)为油藏引起的异常值。
南阳-泌阳航空γ能谱测区,测得的钍、铀、钾含量,按钍含量分间隔,计算其平均值,列于表7-2-3。根据此表中数据,由(7-2-7)和(7-2-8)式得:
核辐射场与放射性勘查
表7-2-3 南阳-泌阳航空γ能谱计算的钍、铀、钾
3)对总道γ计数率的归一化校正。钍比较稳定,可以认为与油气藏形成的放射性异常无关。经研究得知,原岩的总道计数率(I点i)与钍含量的对数值存在近似的线性关系,即
核辐射场与放射性勘查
根据γ能谱实测数据求得实测i点的总道计数率(Ii)与I点i的差值:
核辐射场与放射性勘查
即为消除岩性影响的,由油气藏引起的γ总计数率异常值。
图7-2-2 钍归一化校正岩性影响的结果
图7-2-2为任丘双河油田,两条测线(1100线和11010线)。用钍归一化法,消除岩性影响的结果。油田边界高值和油田上方低值,除钾11010线外都比较明显清晰。与已知油田边界基本一致。
㈨ 叠前地震数据重建方法研究
霍志周
(中国石化石油勘探开发研究院,北京 100083)
摘 要 地震勘探的目的是为了获得地下构造的精确成像。由于人为因素和环境原因,地震数据在空间方向上往往是不规则采样或缺失采样的,因此经常需要在空间方向对缺失的地震数据进行重建。最小范数傅立叶重建方法是基于估算非规则采样地震数据傅立叶系数的方法,一旦准确求得这些系数,就可以通过傅立叶反变换将地震数据重建到任何合适的空间位置。该方法的主要优点是既可以处理规则采样数据有空道的情况,也可以处理非规则采样的数据;该方法的缺点是无法重建含空间假频以及含空隙过大的地震数据。针对含空间假频的地震数据重建问题,本文通过将最小范数傅立叶重建方法和多步自回归方法相结合,较好地克服了最小范数傅立叶重建方法的缺点。通过对不同的理论和实际地震数据算例的验证,表明了该重建方法的有效性和实用性。
关键词 地震数据重建 最小范数反演 傅立叶变换 多步自回归
Research on Pre-stack Seismic Data Reconstruction Method
HUO Zhizhou
(Exploration and Proction Research Institute,SINOPEC,Beijing 100083,China)
Abstract The objective of exploration seismology is to obtain an accurate image of the subsurface.Due to human-related reasons and environmental circumstances,more often than not the seismic data can be irregularly sampled or missing sampled in spatial direction.Therefore,it often needs to reconstruct missing seismic data along spatial direction.Fourier reconstruction with minimum norminversion is based on estimating the Fourier coefficients that describe the irregularly sampled seismic data,and once these coefficients have been obtained, seismic data can be reconstructed on any suitable spatial location via inverse Fourier transformation.The main advantages of Fourier reconstruction are flexible,as it can not only handle regularly sampled data with gaps,but also can handle irregularly sampled data.The disadvantage of this method is that the method can’t handle spatially aliased seismic data and seismic data with large gaps.In this article,for reconstruction question of spatially aliased seismic data,Fourier reconstruction with minimum norminversion and multi-step autoregressive method is combine.This method overcomes the shortcomings of the Fourier reconstruction method.Several different theoretical and practical seismic data would be reconstructed using multi-step autoregressive method,that prove the effectiveness and practicality of this method。
Key words seismic data reconstruction;minimum norm inversion;Fourier transforms;multistep autoregressive
众所周知,地震数据的采集严重影响地震数据最终的成像结果,而地震数据采集中很常见的一个问题就是地震数据沿着空间方向是非规则采样或是稀释采样的。地震数据在空间方向上稀疏采样的原因主要是出于经济因素的考虑,稀疏采样比较经济,但意味着采集到较少的数据,而且会导致地震数据中含有空间假频,尤其是在3D地震勘探中。引起地震数据在空间方向上非规则采样的原因主要有:地表障碍物的存在(建筑物、道路、桥梁等)或地形条件因素(禁采区和山区、森林、河网地区等)、仪器硬件(地震检波器、空气枪、电缆等)问题引起的采集坏道以及海洋地震数据采集时电缆的羽状漂流等。在地震数据处理过程中,非规则采样和稀疏采样不但会引起人为误差,而且会对基于多道技术的DMO、FK域滤波、速度分析、多次波衰减、谱估计和波动方程偏移成像等方法的处理结果带来严重的影响,因此通过对原有的地震数据进行重建,使其包含的地球物理信息更加真实地反映地下地质体的地球物理特征,使得后续地震数据处理能够更好地满足对复杂地质构造进行精细刻画的要求,为油气勘探提供更有效的指示和帮助等具有重要的现实意义[1,2]。
基于傅立叶变换的地震数据重建方法不需要地质或地球物理假设,只要求地震数据是空间有限带宽的,并且计算效率高。傅立叶重建方法利用最小二乘反演估算非规则采样数据的傅立叶系数,如何更好地估算傅立叶系数是该方法的核心。一旦傅立叶系数被正确估算出来,数据可以重建到任意采样网格上。Duijndam等[3]将傅立叶重建方法应用于非规则采样地震数据的规则化上,并成功解决了参数选择等一系列问题。Hindriks和Duijndam[4]将该方法扩展到3D地震数据重建中。Liu和Sachhi[5]提出了最小加权范数插值的傅立叶重建方法,该带限重建方法利用自适应谱加权范数的正则化项来约束反演方程的解,将数据的带宽和频谱的形状作为带限地震数据重建问题的先验信息,因此得到了比传统的带限数据傅立叶重建方法更好的解,但没有给出好的反假频方法。Zwartjes和Sachhi[6]提出了使用非二次型正则化项的稀疏约束傅立叶重建方法,以改善地震数据含较宽的空道时的重建效果,并较好地解决了含有空间假频的地震数据的重建问题。傅立叶重建方法不但可以重建规则采样的地震数据,而且可以重建非规则和随机采样的地震数据,但是不能很好地重建含有空间假频的地震数据。
本文对基于最小范数解的傅立叶地震数据重建方法的研究分析,通过最小二乘反演方法得到傅立叶域的系数来进行地震数据重建。为了改进最小范数傅立叶重建方法不能重建空道间距过大的地震数据和无法重建含有空间假频的地震数据的缺点,本文采用了最小范数傅立叶重建方法和多步自回归方法相结合的思想进行地震数据重建,该方法不但能重建空道间距大的地震数据,而且可以重建含有空间假频的地震数据。
1 最小范数傅立叶重建方法
傅立叶重建是从非规则采样数据上恢复信号的一种方法,它是基于采样定理的,也就是说一个带限的连续信号能够从规则采样数据中恢复。如果非规则采样信号的平均采样率超过Nyquist采样率,则非规则采样的信号也可以重建。在规则采样的情况下,离散傅立叶变换是正交变换。但是当采样是非规则时,傅立叶变换的基函数不再是正交的,这就意味着直接用离散傅立叶变换计算傅立叶系数将产生误差。利用最小二乘反演计算傅立叶系数就是一种补救措施[7]。
假设数据是在空间方向上是不规则采样的,每个采样点的位置分别为[x0,…,xn,…,xN-1]。使用真实的采样位置和采样间隔的中点法则,非规则采样数据的离散傅立叶变换可由以下离散求和的形式表达:
油气成藏理论与勘探开发技术(五)
上式为非均匀离散傅立叶变换。其中,空间采样间隔△xn定义为:
油气成藏理论与勘探开发技术(五)
在波数域规则采样意味着数据在空间域是周期性的,所以 X为非规则采样数据的长度。如果直接用NDFT(Non-uniform Discrete Fourier Transform)计算波数,则由于采样非规则而会引起极大的误差,因此实际计算时通常采用最小二乘反演来计算波数。
首先定义由规则采样波数计算任意空间位置采样数据的数学变换,把它当作正演模型。假设带限数据的波数域带宽为[-M△k,M△k],在波数域规则采样,△k为空间波数采样间隔,则由波数域重建任意空间位置xn的离散傅立叶反变换为
油气成藏理论与勘探开发技术(五)
记系数矩阵为 不规则采样数据为dn=P(xn,ω),待求的规则波数为
油气成藏理论与勘探开发技术(五)
则将公式(3)写成矩阵形式为油气成藏理论与勘探开发技术(五)
在实际的地震数据处理中,由于数据可能不完全是带限的,所以部分空间波数成分会超出定义的频带范围,这些超出的成分构成了上述正演模型的误差和噪音,因此在上式中需要噪声项:
油气成藏理论与勘探开发技术(五)
Duijndam等[3]通过最小二乘反演估计得到非规则采样数据d(xn,t)的空间波数 从非规则采样数据向量d中计算出未知的规则采样的傅立叶系数向量 可以归结为求解一个不适定线性反演问题,需要对其进行正则化,借助一些先验信息构建出合适的解。可以使用任何所需的参数估计技术,首先我们假设噪音n=N(0,Cn)和先验信息
油气成藏理论与勘探开发技术(五)
都是高斯分布的,噪音的协方差矩阵为Cn,其平均值为零。利用贝叶斯参数反演方法通过寻找后验概率密度函数油气成藏理论与勘探开发技术(五)
的最大值来进行反演,其中 是似然函数, 表示模型向量的先验分布。分别满足
油气成藏理论与勘探开发技术(五)
油气成藏理论与勘探开发技术(五)
求 的最大后验概率解转化为求下面目标函数的最小化解,建立目标函数
油气成藏理论与勘探开发技术(五)
最小化目标函数得:
油气成藏理论与勘探开发技术(五)
这里, 为计算要得到的规则采样波数,AH为矩阵A的共轭转置矩阵, 为先验模型的协方差矩阵。
下面我们对(9)式进行简化。首先对于地震数据,通常没有先验模型信息,因此 一般没有理由假设空间波数之间的相关性,所以 是对角阵,通常的形式为 是先验模型的方差。准确地表达噪音的协方差矩阵Cn是不现实的,因为关于噪音详细的信息是未知的。Duijndam等[3]给出的噪音协方差矩阵为Cn =c2W-1,c是常数;W为权系数组成的对角阵,即W=diag(△xn)。根据离散傅立叶变换理论,应选择△k≤2π/X,这里X=∑n△xn,为数据的长度,即X=xN-1-x0,则(9)式变为
油气成藏理论与勘探开发技术(五)
其中, 称为阻尼因子。λ可以通过L-curve或者广义交叉验证(GCV)方法确定,最佳的选取方法是[4]:
油气成藏理论与勘探开发技术(五)
式中:F为用户给定的常数,表示期望的数据信噪比值。但在实际地震数据重建过程中,λ一般取AHWA矩阵主对角元素的1%。
方程(10)的解称为最小范数解,也称为阻尼最小二乘解,该重建方法称为最小范数傅立叶重建方法(Fourierreconstruction with minimum norminversion,FRMN)[8]。通常非规则采样时,式(10)的系数矩阵AHWA为病态的Toeplitz矩阵。当不加权矩阵W时,AHA形成的Toeplitz矩阵病态程度受非规则采样数据之间的致密程度控制。非规则采样地震数据中地震道靠得越近,间距△x越小,则Toeplitz矩阵的条件数就越大,求解越困难;加上权系数矩阵W后,AHWA形成的Toeplitz矩阵病态程度受各数据之间的最大空隙△xa的大小控制,△xa=max(△xn)。系数矩阵AHWA的条件数与最大空隙△xa的关系如下[7]:
油气成藏理论与勘探开发技术(五)
由上式可见,最大空隙△xa越大,矩阵AHWA病态程度越大,求解方程时就越难以收敛。如果定义空间Nyquist采样间隔为
油气成藏理论与勘探开发技术(五)
则当△xa≥3△xNyq时,系数矩阵AHWA已经无法保证迭代收敛[3]。也就是说当非规则采样地震数据的空隙太大时,不能得到满意的重建效果。这是傅立叶重建方法的固有弊病。
方程(10)实际求解时一般在频率域逐频率求解。在求解方程时,由于低频部分只需要很小的波数带宽就能完整重建数据,因此求解方程(10)的规模小,求解相对容易;而高频部分则需要较大的波数带宽,因此求解式(10)中的未知数多,求解需要更多的计算时间,而且解也不稳定。因此,利用最小范数傅立叶方法重建的地震数据低频部分有较高的精度。
2 多步自回归方法
自回归模型(预测滤波器)在信号处理领域具有广泛的应用,它是一种模拟信号演化的技术[9]。自回归模型可以应用于信号预测和噪音消除[10]、地震道内插[11,12]以及参数频谱分析[13]等方面。t-x域的线性同相轴变换到f-x域是复正弦函数,该函数可以通过自回归算子来模拟。Spitz[11]和Porsani[12]提出了自回归的重建方法,成功地解决了规则采样含空间假频地震数据的插值问题,这些方法是利用低频信息来恢复数据的高频部分。但这种方法只适用原始地震数据是空间规则采样的情况,而且只能用于加密插值。
多步自回归方法(multistep autoregressive,MSAR)[14]是对Spitz单步预测方法的拓展,使其应用范围从只能进行道加密插值扩展到能对不规则缺道地震数据进行插值重建。假设地震数据包含有限个线性同相轴,由N个等间距的地震道组成,部分地震道是缺失的。首先将地震数据从时间域变换到频率域,在f-x域,地震数据可以用向量x(f)表示,xT(f)=[x1(f),x2(f),x3(f),…,xN(f)],其中只有M道数据是已知的。分别用n={n(1),n(2),n(3),…,n(M)}和m={m(1),m(2),m(3),…,m(N-M)}表示已知数据和未知数据(缺失道)的下标,目标是从xn(f)中恢复出xm(f)。
由L个近似线性的同相轴构成的地震数据在f-x域可表示为
油气成藏理论与勘探开发技术(五)
式中:△x和△f分别表示空间域和频率域采样间隔;pj表示第j个线性同相轴的斜率;Aj表示振幅。对于每个频率成分f,上式表明在f-x域每个线性同相轴都可以用复谐波函数来表示。考虑当△x′=α△x,△f′=△f/α时,得到:
油气成藏理论与勘探开发技术(五)
此外,通过自回归模型的形式,可将L个谐波函数的叠加表达为
油气成藏理论与勘探开发技术(五)
其中P(j,n△f)表示预测滤波因子。同样的,对于△x′和△f′,有
油气成藏理论与勘探开发技术(五)
比较表达式(15)、(16)和(17),可得:
油气成藏理论与勘探开发技术(五)
该式即为多步自回归方法的基础。它表明在频率轴上,对于预测滤波器的每个成分都是可预测的。这就意味着,如果已知某些频率的预测滤波器,可以预测得到其他频率的预测滤波器。也就是说,我们可以从傅立叶方法重建得到的无空间假频的低频成分的预测滤波器中提取高频成分的预测滤波器,进而重建得到缺失地震道的高频成分。
假设用最小范数傅立叶方法重建得到的低频数据的频率范围为f∈[fminr,fmaxr],在f-x域线性同相轴向前和向后预测的多步预测滤波器可以由下列方程组确定:
油气成藏理论与勘探开发技术(五)
式中:*表示复共轭;L表示预测滤波器的长度;Pj(f)表示预测滤波器。这些方程对应一种特殊类型的自回归模型,向前自回归方程(19)和向后自回归方程(20)是通过每次向前和向后跳α步来实现的。通过自回归方程(19)和(20)可以计算出在α步时的预测滤波器Pj(f)。参数α=1,2,…,αmax是步长因子,用于从频率f中提取频率αf的预测滤波器。由于步长因子是一个正整数,很显然低频部分为数据重建算法提供了重要的信息。步长上限αmax依赖于地震道数N和预测滤波器的长度L,该参数由下式给出
油气成藏理论与勘探开发技术(五)
这里[.]表示取整数部分。
当用多步自回归方法从已重建的低频数据x(f)中计算出高频数据x(f′)的预测滤波器时,同Spitz插值方法相似,可以通过已知的数据和预测滤波器重建出缺失的数据。向前和向后自回归重建方程为
油气成藏理论与勘探开发技术(五)
设地震数据中含有L个不同斜率的线性同相轴,地震数据的有效频带范围为[fmin,fmax],含空间假频的不规则道缺失的地震数据的重建实施步骤为:(1)首先将原始地震数据变换到f-x域,用最小范数傅立叶方法重建无空间假频的低频段[fminr,fmaxr]的地震数据,得到低频段地震数据,其中fminr=fminr。对于不含空间假频的有限带宽信号而言,FRMN重建得到的地震数据精度较高;(2)运用方程(19)和(20),从低频段[fminr,fmaxr]中提取高频成分的预测滤波器Pj(f′);(3)利用已知道数据和预测滤波器Pj(f′)重建缺失的地震数据;(4)最后将重建后的地震数据反变换回t-x域。遇到复杂地震数据时,同相轴可能不满足线性假设,可将地震数据划分成多个小时空窗,分窗口进行重建。综上所述,从无空间假频低频段[fminr,fmaxr]数据中提取缺失数据高频成分f′=αf的预测滤波器,然后利用已知数据和预测滤波器计算缺失数据的高频成分,最终完成多步自回归重建。
3 理论数据算例
为了验证多步自回归算法的有效性,本节中我们将该算法应用于理论数据,进行缺失道的重建以及加密插值。第一个理论数据如图1(a)所示,是由7个不同斜率的线性同相轴组成,其f-k谱含有严重的空间假频(如图1(c)所示)。共有81道,道间距为5m,时间采样间隔为2ms,采样点数为901。图1(b)是从原始数据中随机抽去了40%的地震道后得到的数据。图1(d)是图1(b)对应的f-k谱。从图1(d)中可以看出,由于地震道的缺失而导致f-k谱上产生严重的噪音。
图1 多步自回归法理论算例
图2 最小范数傅立叶重建方法与多步自回归法的理论联合应用(一)
图2(a)是利用FRMN方法重建出的低频数据,其f-k谱如图2(c)所示。重建出的低频数据被MSAR算法用于提取预测滤波器来重建数据的高频部分。对于数据低频端的预测滤波器是通过预测滤波器的外推来估计。通过FRMN + MSAR方法重建后的完整数据如图2(b)所示,其对应的f-k谱如图2(d)所示,与原始数据的f-k谱(图1(c))相对比,几乎完全一样,由采样缺失引起的噪音已被消除。与原始数据(图1(a))相对比,缺失的地震道被填充,线性同相轴的连续性也很好。
图3 最小范数傅立叶重建方法与多步自回归法的理论联合应用(二)
图4 图3中数据对应的f-k谱
图5 最小范数傅立叶重建方法与多步自回归方法的实际应用
为了进一步验证算法在复杂情况下的适用性,我们选取了Marmousi模型数据中的一个单炮数据(图3(a)),共有96道数据,道间距为25m,时间采样间隔为4ms,采样点数为750。随机抽去了其中的27道数据(图3(b)),用FRMN + MSAR方法对该数据进行重建,图3(c)显示的是用FRMN方法重建的低频段的数据,图3(d)显示的是用FRMN+MSAR方法重建的完整单炮数据。由于模型很复杂,所以原始单炮数据的f-k谱有空间假频的存在(图4(a))。图4(b)是图3(b)对应的f-k谱,可以看出含有严重的噪音。图4(c)和图4(d)分别是3(c)和图3(d)对应的f-k谱。重建后的数据f-k谱中的噪音消除了,缺失的道也得到了填充,而且同相轴也保持很好的连续性。
图6 图5中数据对应的f-k谱
4 实际数据算例
本节我们将对实际数据进行重建,以验证FRMN +MSAR方法的适用性。选取一个共偏移距地震剖面的部分数据(图5(a)),总共有201道,道间距为12.5m,时间采样间隔为2ms。随机抽去其中30%的地震道(图5(b))进行重建,图5(c)展示的是FRMN方法重建的低频段的数据,图5(d)展示的是FRMN+MSAR重建的完整数据。图6(a)、图6(b)、图6(c)和图6(d)分别是图5(a)、图5(d)、图5(c)和图5(d)对应的f-k谱。可以看出,重建前后数据f-k谱的变化很小。重建后数据的缺失道得到了恢复,且同相轴连续,重建的结果接近于原始数据。
5 结论
本文在最小范数傅立叶重建方法的基础上,结合多步自回归方法进行含空间假频地震数据的重建。多步自回归方法是对Spitz方法的拓展,也是基于近似线性同相轴的假设。因此在处理复杂地震数据的时候一般难以满足这个假设,这时可采用小时空窗的方法来进行计算,在小时空窗中可以认为满足近似线性的假设。但是时空窗太小会使数据量不足,反而会导致重建的结果不好或可能无法重建。众所周知,为了能够求解大多数的地球物理问题,必须基于某些假设条件。一般在处理实际数据时,都是部分地违背这些假设的。事实上,对于中等程度弯曲的同相轴本方法同样能取得比较理想的重建结果,说明本文的重建方法具有很好的稳定性。实际上,对于含有大间距空道的地震数据,该方法同样取得了较好的重建结果。通过对一些理论数据和实际数据进行重建实验,验证了本文中重建方法的有效性和实用性。另外,地震数据的重建效果同原始数据的复杂程度以及谱的性质、缺失地震道的数量及位置和缺失道间距的大小等多方面原因有关,需要进一步研究这些因素对重建算法的影响。
参考文献
[1]Eiken O,Haugen G U,Schonewile M A,and Duijndam A J W.A proven method for acquiring highly repeatable towed streamer seismic data[J].Geophysics,2003,68(4):1303~1309.
[2]Wever A,Spetzler J.Criteria for source and receiver positioning in time-lapse seismic acquisition:74th Ann.Internat.Mtg.,SEG,Expanded Abstracts,2004:2319~2322.
[3]Duijndam A J W,Schonewille M A,and Hindriks C O H.Reconstruction of band-limited signals, irregulary sampled along one spatial direction[J].Geophysics,1999,64(2):524~538.
[4]Hindriks K,Duijndam A J W,Reconstruction of3 -D seismic signals irregularly sampled along two spatial coordinates[J].Geophysics,2000,65(1):253~263.
[5]Liu B,Sacchi M.Minimum weighted norminterpolation of seismic records[J].Geophysics,2004,69(6):1560~1568.
[6]Zwartjes P M ,Sacchi M D.Fourier reconstruction of nonuniformly sampled,aliased seismic data[J]. Geophysics,2007,72(1):21~32.
[7]Feichtinger H G,Grochenig K,Strohmer T.Efficient numerical methods in non -uniformsampling theory[J].Numerische Mathematik,1995,69:423~440.
[8]Zwartjes P M,Fourier reconstruction with sparse inversion:Ph.D.thesis,Delft University of Technology,2005.
[9]Takalo R,Hytti H,Ihalainen H.Tutorial on univariate autoregressive spectral analysis[J].Joural of Clinical Monitoring and Computing,2005,19(6):401~410.
[10]Canales L L.Random noise rection:54th Ann.Internat.Mtg .,SEG,Expanded Abstracts,1984: Session:S10.1.
[11]Spitz S.Seismic traces interpolation in f-x domain[J].Geophysics,1991,56(6):785 ~794.
[12]Porsani M J.Seismic trace interpolation using half-step prediction filters[J].Geophysics,1999,64(5):1461~1467.
[13]Marple S L.Digital spectral analysis with applications.Englewood Cliffs,New Jersey:Prentice-Hall Inc,1987.
[14]Naghizadeh M,Sacchi M D.Multistep autoregressive reconstruction of seismic records[J]. Geophysics,2007,72(6):111-118.