地球物理学报  2015, Vol. 58 Issue (3): 1042-1058   PDF    
最小曲率位场分离方法研究
纪晓琳, 王万银, 邱之云    
长安大学重磁方法技术研究所, 长安大学地质工程与测绘学院, 长安大学西部矿产资源与地质工程教育部重点实验室, 西安 710054
摘要:位场分离是位场数据处理和解释中的重点和难点之一.本文给出了单步长非原位和原位两种最小曲率位场分离差分迭代格式,并利用Fourier频谱分析理论研究了这两种迭代格式的收敛性.通过研究表明,单步长非原位迭代格式不收敛,只有单步长原位迭代格式收敛,但单步长原位迭代格式受迭代方向选择的影响,随着迭代次数的增大其影响逐渐消失.根据单步长非原位迭代格式的频谱特点,提出了叠加步长非原位和原位迭代格式,同样利用Fourier频谱分析理论研究了叠加步长非原位和原位迭代格式的收敛性.通过研究认为,一维叠加步长非原位迭代格式收敛,但二维叠加步长非原位迭代格式不收敛;不论是一维或二维,其原位迭代格式均收敛.进一步的理论研究表明,非原位迭代格式的频率响应是一个实偶函数,而原位迭代格式的频率响应是一个复函数;单步长迭代格式的频率响应具有一定的周期性,而叠加步长迭代格式的频率响应无周期性特征;叠加步长迭代格式比单步长迭代格式的收敛性好.
关键词单步长迭代     叠加步长迭代     原位迭代格式     非原位迭代格式     收敛性     频率响应    
The research to the minimum curvature technique for potential field data separation
JI Xiao-Lin, WANG Wan-Yin, QIU Zhi-Yun    
Gravity & Magnetic Institute of Chang'an University, College of Geology Engineering and Geomatics, Chang'an University, Key Laboratory of Western China's Mineral Resources and Geological Engineering, Ministry of Education, Chang'an University, Xi'an 710054, China
Abstract: Potential field data separation is one of the difficult and important problems in the potential field data processing and interpretation. Results of the potential field data separation have direct affection to the potential field data explanation. The minimum curvature technique has been widely used in the extension, blank padding and data gridding of the potential field data. We studied the frequency response and the convergence when using the minimum curvature technique for the separation of the potential field data, and take out a method for separation of the potential field data.
The minimum curvature difference iteration scheme has two types: 1-D and 2-D. We used the minimum curvature basic difference iteration scheme to construct the non in situ and in situ iteration of the single step. Studied these two iteration schemes using the Fourier spectrum analysis theory, get their frequency response and studied their convergence of the iteration schemes through their frequency response feature. Popularize the Fourier spectrum analysis theory to the 2-D minimum curvature difference iteration schemes, studied the frequency response feature and the convergence of the single step non in situ and in situ iterations. Took out the superposition step non in situ and in situ iterations schemes according to the frequency response feature of the single step non in situ iteration, and studied their frequency response feature and the convergence using the Fourier spectrum analysis theory.
Research of the Fourier spectrum analysis theory shows that: 1-D single step non in situ iteration cannot be used for separation of the potential field data because it has high-frequency amplification and its iteration is not convergent. But the 1-D superposition step non in situ iteration can be used for separation of the potential field data since it has good low-pass performance and its iteration is convergent. Neither single step nor superposition step, frequency response of their non in situ iteration schemes are a real even function. Frequency response of the 1-D from left to right or 1-D from right to left single step and superposition step in situ iteration are a low-pass complex function. Frequency response of the 1-D alternately iteration single step and superposition step in situ iteration are a low-pass real even function. The three types of the in situ iteration schemes are all convergent. Further theoretical study shows that when the iteration number is bigger than a bounded number, frequency response of the three types single step and superposition step in situ iteration schemes are tend to be the same. 2-D single step and superposition step non in situ iteration schemes cannot be used for separation of the potential field data because they have high-frequency amplification and their iteration are not convergent, but the convergence of the superposition step non in situ iteration is better than the single step non in situ iteration's. Frequency response of the 2-D non in situ iteration schemes are a real even function. 2-D single step and superposition step in situ iteration can be used for separation of the potential field data since they have good low-pass performance and their iteration is convergent. Frequency response of the 2-D in situ iteration scheme is a low-pass complex function. So through the comprehensive comparison, we considered that the convergence of the superposition step iteration scheme is better than the single step iteration's.
We studied the frequency response of the minimum curvature difference iteration schemes for potential field data separation using the Fourier spectrum analysis theory. We get the frequency response feature and the convergence of the single step and superposition step, non in situ and in situ iteration schemes. It was proved succeed and it can be popularized to the convergence research of the other difference iteration schemes.
Key words: Single step iteration scheme     Superposition step iteration scheme     In situ iteration scheme     Non in situ iteration scheme     Convergence     Frequency response    
1 引言

位场分离是位场数据处理和解释中的难点之一,也是重点之一.位场分离的效果直接影响位场解释的效果,目前位场分离方法有空间域和频率域两种类型,这里研究的最小曲率位场分离方法属于空间域中的一种.

最小曲率方法已经在位场数据处理方面得到了广泛应用,如离散数据的网格化(Briggs,1974Swain,1976Smith和Wessel,1990Carlos和Silva,1995Xiong,1999Michael等,2005王万银和邱之云,2011)、网格数据的扩边和补空处理(王万银等,2009张志厚等,2013).但在位场分离方面的研究很少.1991年,Mickus首次将Briggs(1974)提出的最小曲率迭代格式应用于美国东部弗吉尼亚地区布格重力异常的分离,成功的识别出了Richmond盆地的范围.通过地震反射剖面、磁测剖面、可控源音频大地电磁以及测井资料等进行约束,构建了2.5D重力模型.并反演计算了Richmond盆地中部基底的起伏变化,根据该基底起伏变化,选择网格间距为1km进行重力异常分离,并与6阶次多项式趋势分析方法的位场分离结果进行了对比.通过对比认为,最小曲率方法比6阶次多项式趋势分析方法的位场分离效果好.本文将给出单步长和叠加步长的非原位迭代和原位迭代最小曲率位场分离差分迭代格式,并利用Fourier频谱分析理论研究其收敛性. 2 1-D最小曲率位场分离方法及收敛性

1-D最小曲率方法的基本差分迭代格式为(王万银等,2009)

这里s为节点序号,M为节点数,u(s)为第s个节点上的位场值.若网格间距为Δx,则上式可进一步表示为

根据该基本差分迭代格式可以构造出非原位迭代和原位迭代两种类型的位场分离差分格式.下面分别给出这两种类型的差分迭代格式,并利用Fourier频谱分析理论研究其收敛性.

2.1 非原位迭代格式及其收敛性

若用lx=l×Δx表示迭代步长(以后称迭代步长为l),用u(k)(l)(sΔx)表示sΔx点迭代步长为l的第k次迭代结果.则非原位迭代格式就是利用第k-1次的迭代值u(k-1)(l)(sΔx)来计算第k次的迭代值u(k)(l)(sΔx).非原位迭代又分单步长和叠加步长两种,下面分别给出这两种非原位迭代格式,并利用Fourier频谱分析理论研究其收敛性.

2.1.1 单步长非原位迭代格式及其收敛性

根据基本差分迭代格式得到的单步长非原位迭代格式为


其中u(0)(l)(sΔx)=u(sΔx)为初始位场值,K为最大迭代次数.u(K)(l)(sΔx)为迭代结果,称为区域场;u(sΔx)-u(K)(l)(sΔx)为迭代剩余,称为剩余场.

对(2)式施行Fourier变换,并利用Fourier变换的位移性质

得到(2)式的Fourier变换结果为

其中H(1)1(f,l)表示迭代步长为l时的单步长非原 位迭代格式完成1次迭代的频率响应,可进一步写为

根据采样定理,u(sΔx)的截止频率为fc=1/(2Δx),则 2πf×lx = 2πf×l×Δx ≤lπ.现在来研究这个频率响应H1(1)(f,l)的特点.

(1)奇偶性及虚实性:因为H1(1)(-f,l)=H1(1)(f,l),而且虚部为0,所以H1(1)(f,l)是一个实偶函数.由于是偶函数,故只讨论2πf×lx≥0的函数特点.

(2)周期性:因为H1(1)(f,l)是由余弦函数构成,变量2πf×lx的周期为2πn.又由2πf×lx≤lπ得出,n=l/2.当l=2n时有n个周期,当l=2n+1时有n+0.5个周期.

(3)零值点位置:令H1(1)(f,l)=1- 1-cos(2πf×lx)2=0,推出cos2πf×lx=1- 3/2,得H1(1)(f,l)的零值点位置为2πf×lx=[2nπ±arccos(1-√ 3/2)],(n=0,1,2,…). 又由于2πf×lx≤lπ和反余弦函数的主值区间为0-π得出,在2πf×lx≥0时共有l个零值点,第l个零值点位置为2πf×lx=[2π[l/2]-(-1)l×arccos(1-√ 3/2)],l=1,2,…,l,这里[l/2]表示l/2取整数.由此得出第一个零值点位置为2πf×lx=arccos(1-√ 3/2),即随着迭代步长l的增大,第一个零值点频率f0= 1 2π×lx ×arccos(1-√ 3/2) 逐渐缩小,并与迭代步长l成反比.

(4)最大值:由H1(1)(f,l)的表达式得到其最大值为

可以看出,H1(1)(f,l)的最大值超过1,故H1(1)(f,l)是一放大滤波器,因此不能用来进行位场分离.

(5)极值点位置:由极值的必要条件(2πf×lx)]=0推出极值点的位置为2πf×lx=mπ.又由2πf×lx≤lπ得出,m=0,1,2,…,l.故在2πf×lx≥0时共有(l+1)个极值点.在该极值点上,其二阶导数为

由于2πf×lx≤lπ,当2πf×lx=(2j+1)π(j=0,1,…,[(l-1)/2]) 时为极小值点,共有[(l-1)/2]+1个极小值点,且极小值为-5/3;当2πf×lx=2jπ(j=0,1,…,[l/2])时为驻点,实际上就是极大值点,共有[l/2]+1个极大值点,且极大值为1,在该极大值点附近H1(1)(f,l)的变化率较小,因此才会表现出其二阶导数为0的情况.由于H1(1)(f,l)存在小于0的情况,故在滤波时存在某些频率的反相.由此看出,随着迭代步长的增大,极值点个数增加,不同频率的放大和压制作用不同.

若进行K次迭代,则其频率响应H1(k)(f,l)为

由于H1(k)(-f,l)=H1(k)(f,l),且虚部为0,所以H1(k)(f,l)也是一个实偶函数.H1(k)(f,l)的周期性与H1(1)(f,l)的相同.H1(k)(f,l)的零值点位置和个数与H1(1)(f,l)的相同.H1(k)(f,l)的极值点位置和个数与H1(1)(f,l)有所不同.当K为奇数时,H1(k)(f,l)的极值点位置和个数与H1(1)(f,l)的相同,极大值都为1;但极小值不同,H1(k)(f,l) 的极小值是H1(1)(f,l)极小值的K倍,为(-5/3)K<0.当K为偶数时,H1(k)(f,l)的值大于等于0,故其极小值点位置和个数与H1(1)(f,l)的零值点位置和个数相同,极小值为0;但H1(1)(f,l)的极小值和极大值位置都是H1(k)(f,l)的极大值位置,极大值为(5/3)K>0或1,个数为(l+1)个.

通过上面的频率响应特征分析认为,由于在某些频率范围内,H1(1)(f,l)和H1(k)(f,l)的幅值均大于1,即存在振幅放大作用,故(2)式的迭代是不收敛的.随着迭代次数的增加,高频放大作用增强.

图 1为步长分别为1,2和3时完成1次单步长迭代的频率响应,其中横坐标为fx×Δx,是一个无量纲之量,其最大值为0.5(图 2图 8中的横坐标与此相同).由图 1看出,H1(1)(f,l)具有一定的周期性,当l=1时,有半个周期;当l=2时,有一个周期;当l=3时,有一个半周期.H1(1)(f,l)共有l个零值点,当l=1时,有一个零值点;当l=2时,有两个零值点,第一个零值点是l=1时第一个零值点的1/2;当l=3时,有3个零值点,第一个零值点是l=1时第一个零值点的1/3.H1(1)(f,l)共有(l+1)个极值点,当l=1时,有一个极大值点和一个极小值点;当l=2时,有两个极大值点和一个极小值点;当l=3时,有两个极大值点和两个极小值点;极大值为1,极小值为-5/3,并且在极大值附近,H1(1)(f,l)的水平变化率较小.

图 1 迭代次数为1时不同迭代步长非原位迭代频率响应(实部)
蓝色线条,绿色线条和红色线条分别表示步长为1,2和3.
Fig. 1 The frequency response(real part)of the non in situ iteration scheme of different single steps when the number of iteration is 1
Blue line,green line and the red line are represent for the single step are 1,2 and 3 respectively.

图 2为步长为4时完成1次、2次和3次单步长迭代的频率响应.由图 2看出,偶次迭代的频率响应特征相同,奇次迭代的频率响应特征相同.偶次迭代频率响应值大于等于0;奇次迭代的频率响应值有正有负.不同迭代次数频率响应的零值点位置相同,即零值点位置随迭代步长变化,不随迭代次数变化.奇次迭代的极小值点位置和极大值点位置与1次迭代的相同;偶次迭代的极小值点位置与1次迭代的零值点位置相同,偶次迭代的极大值点位置与1次迭代的极小值和极大值点位置相同.随着迭代次数的增加其频率放大作用迅速增强.

图 2 迭代步长为4时不同迭代次数非原位迭代频率响应(实部)
蓝色线条,绿色线条和红色线条分别表示迭代次数为1,2和3.
Fig. 2 The frequency response(real part)of the non in situ iteration scheme of the different iteration numbers when the single step is 4
Blue line,green line and the red line are represent for the number of iteration are 1,2.
2.1.2 叠加步长非原位迭代格式及其收敛性

图 1来看,不同迭代步长的频率响应具有一定的周期性,而且正、负相伴.如果把不同迭代步长的迭代结果求和再取平均值是不是会削弱其对某些频率的放大作用?为此,利用不同迭代步长的迭代结果相加再求其平均值来构成新的非原位迭代格式,可以表示为

该式称为叠加步长为L(最大迭代步长)的非原位迭代格式.若对其施行Fourier变换可得到1次迭代的频率响应为

那么K次迭代的频率响应H(K)2(f,L)为

由(7)和(8)式可以看出,叠加步长非原位迭代格式的频率响应也是一个实偶函数.由于该函数较复杂,用解析法较难讨论其函数特点,故用图形来加以说明.图 3为叠加步长分别为2,3和4时完成1次非原位迭代的频率响应,图 4为叠加步长为4时完成1次、2次和3次非原位迭代的频率响应.由图 3看出,随着叠加步长L的增大其第一个零值点的频率逐渐趋向于0,且无明显周期性特征,零值点的个数与叠加步长L成正比关系;并随着叠加步长的增大高频压制能力增强.由图 4看出,随着迭代次数的增加其零值点的位置不变,但高频压制能力明显增强.偶次迭代的频率响应大于等于0,奇次迭代的频率响应有正、有负.随着迭代次数的增加,大于第一个零值点的频率响应逐渐趋向于0,只有小于第一个零值点以内的频率响应得以保留.因此,迭代次数决定压制高频的强度,叠加步长决定了滤波的截止频率.由此可见,叠加步长迭代格式是收敛的,是可以用来进行各种数据处理的.

图 3 迭代次数为1时不同叠加步长非原位迭代频率响应(实部)
蓝色线条,绿色线条和红色线条分别表示叠加步长为2,3和4.
Fig. 3 The frequency response(real part)of the non in situ iteration scheme of the different iteration steps when the superposition number is 1
Blue line,green line and the red line are represent for the number of iteration are 1,2 and 3 respectively.

图 4 叠加步长为4时不同迭代次数非原位迭代频率响应(实部)
蓝色线条,绿色线条和红色线条分别表示迭代次数为1,2和3.
Fig. 4 The frequency response(real part)of the non in situ iteration scheme of the different iteration number when the superposition step is 4
Blue line,green line and the red line are represent for the number of iteration are 1,2 and 3 respectively.

通过以上对非原位迭代格式的研究表明,单步长非原位迭代格式具有高频放大作用,迭代不收敛,不能用来进行位场分离;而叠加步长非原位迭代格式具有很好的低通性能,迭代收敛,能够用来进行位场分离.不论是单步长或叠加步长,其迭代格式的频率响应均为实偶函数.

2.2 原位迭代格式及其收敛性

原位迭代格式就是利用第k-1次的迭代值u(k-1)(l)(sΔx)和第k次的迭代值u(k)(l)(sΔx)来计算第k次的迭代值u(k)(l)(sΔx).原位迭代又分单步长和叠加步长两种,下面分别给出这两种原位迭代格式,并利用Fourier频谱分析理论研究其收敛性.

2.2.1 单步长原位迭代格式及其收敛性

对于1-D数据来讲,原位迭代具有方向性,共有三种迭代格式.一种是从左到右(从小序号到大序号)进行迭代,另外一种是从右到左(从大序号到小序号)进行迭代,还有一种是从左到右和从右到左交替进行迭代.假设迭代从左到右进行(从小序号到大序号),则其迭代格式(称为第一种迭代格式)为

但如果迭代是从右到左进行(从大序号到小序号),则对应的迭代格式(称为第二种迭代格式)为

若先从左到右,后从右到左交替进行迭代,则对应的迭代格式(称为第三种迭代格式)为

(9)、(10)和(11)式均为单步长原位迭代格式.为了分析其收敛性,先对(9)式施行Fourier变换,并根据Fourier变换的位移性质得到

其中H3(1)(f,l)表示步长为l时完成(9)式1次原位迭代的频率响应,可写为

A=[4cos(2πf×lx)-cos(2πf×2lx)],
B=[4sin(2πf×lx)-sin(2πf×2lx)],

H3(1)(f,l)是一个复函数,由于A是偶函数,B是奇函数,所以H3(1)(f,l)的实部是偶函数,虚部是奇函数,振幅是偶函数.因此,只讨论2πf×lx≥0的情况.

根据采样定理,u(sΔx)的截止频率为fc=1/(2Δx), 则2πf×lx≤lπ.若l=1,则当2πf×lx=0时,A=3,B=0,H3(1)(f,l)=1;当时,;当2πf×lx=π时,A=-5,B=0,H3(1)(f,l)=-5/11<1.可以看出,H3(1)(f,l) 实部的最大值为1,虚部的最大值为16/41,振幅的最大值为1,为一低通滤波器,即(9)式是收敛的.若进行K次迭代,则其频率响应 H3(K)(f,l)为

该频率响应H3(K)(f,l)也是一个复函数,其实部为偶函数,虚部为奇函数,振幅是偶函数.与H3(1)(f,l)特点基本相同,为一低通滤波器,并且随着迭代次数的增加,其低通性能更加明显.图 5为迭代步长分别为1,2,3时完成(9)式1次原位迭代的频率响应,图 6为迭代步长为4时完成(9)式1次、2次和3次原位迭代的频率响应.由图 5看出,该频率响应具有一定的周期性特征,随着迭代步长l的增大,周期个数增大;频率响应实部的零值点个数与迭代步长l相同,而且第一个零值点的频率随着迭代步长l的增大逐渐趋向于0.也就是说,迭代步长控制着不同频率振幅的压制强度.由图 6可以看出,随着迭代次数的增加,其频率响应也具有一定的周期性特征,其基本特征受迭代步长l的控制,但随着迭代次数的增加,其高频压制能力增强.由此可见,随着迭代步长的增大和迭代次数的增加,单步长原位迭代的低通性能增强.
图 5 迭代次数为1时不同单步长原位迭代(从左到右)频率响应
(a)实部;(b)虚部;(c)振幅;蓝色线条,绿色线条和红色线条分别表示步长为1,2和3.
Fig. 5 The frequency response of in situ iteration scheme ofthe different single steps(from left to right)when the number of iteration is 1
(a)Real part;(b)Imaginary part;(c)Amplitude. Blue line,green line and the red line are represent for the single step are 1,2 and 3 respectively.

图 6 单步长为4时不同迭代次数原位迭代(从左到右)频率响应
(a)实部;(b)虚部;(c)振幅;蓝色线条,绿色线条和红色线条分别表示迭代次数为1,2和3.
Fig. 6 The frequency response of in situ iteration schemeof the different iteration numbers(from left to right) when the single step is 4
(a)Real part;(b)Imaginary part;(c)Amplitude. Blue line,green line and the red line are represent for the number of iteration are 1,2 and 3 respectively.

为了研究不同迭代方向对迭代结果的影响,对(10)式施行Fourier变换得到其频率响应为

可以看出,H4(1)(f,l)也是一个复函数,其实部为偶函数,虚部为奇函数,振幅为偶函数,并且H4(1)(f,l)与H3(1)(f,l)为一对共轭函数.即H4(1)(f,l)和H3(1)(f,l)除相位移不同外,其他特点均相同.其实部小于等于1,虚部小于1,为一低通滤波器,即(10)式也是收敛的.若进行K次迭代,则其频率响应H4(K)(f,l)为

该频率响应也是一个复函数,其实部为偶函数,虚部为奇函数,与H4(1)(f,l)的特点相同,也为低通 滤波器,并且随着迭代次数的增加,其低通性能更加明显.

对(11)式施行Fourier变换,得到两次迭代的频率响应为

该频率响应为实偶函数,其实部为偶函数,而且小于等于1,为一低通滤波器,即(11)式是收敛的.若进行K次迭代,则其频率响应H5(K)(f,l)为

该频率响应也是一个实偶函数,与H5(2)(f,l)的特点相同,也为一低通滤波器,并且随着迭代次数的增加,其低通性能更加明显.

当迭代达到一定次数后,H4(K)(f,l)和H3(K)(f,l)的虚部趋向于0,此时H5(K)(f,l)、H4(K)(f,l)和H3(K)(f,l)的实部相等,即有下式成立为

随着迭代次数的增加,H4(K)(f,l)→H3(K)(f,l)→H5(K)(f,l),即(9)、(10)和(11)式具有相同的收敛性,但当迭代次数较小时,三种迭代格式的结果会有一定的差别.

2.2.2 叠加步长原位迭代格式及其收敛性

仿照非原位迭代格式的做法,也可以利用不同迭代步长的迭代结果相加再求其平均值的迭代方法构建叠加步长原位迭代格式.对于(9)式所对应的叠加步长迭代格式为

对于(10)式所对应的叠加步长迭代格式为

对于(11)式所对应的叠加步长迭代格式为

对(21)式施行Fourier变换后得到1次迭代的频率响应为

那么K次迭代的频率响应H6(K)(f,L)为

可以看出,H6(1)(f,L)和H6(K)(f,L)都是复函数.图 7为叠加步长分别为2,3,4时完成(21)式1次原位迭代(从左到右)的频率响应,图 8为叠加步长为4时完成(21)式1次、2次和3次原位迭代(从左到右)的频率响应.由图 7图 8可以看出,随着叠加步长的增大和迭代次数的增加,其低通性能明显增强,而且无明显周期性特征.

图 7 迭代次数为1时不同叠加步长原位迭代(从左到右)频率响应
(a)实部;(b)虚部;(c)振幅;蓝色线条,绿色线条和红色线条分别表示叠加步长为2,3和4.
Fig. 7 The frequency response of in situ iteration schemeof the different superposition steps(from left to right) when the number of iteration is 1
(a)Real part;(b)Imaginary part;(c)Amplitude. Blue line,green line and the red line are represent for the superposition steps are 2,3 and 4 respectively.

图 8 叠加步长为4时不同迭代次数原位迭代(从左到右)频率响应
(a)实部;(b)虚部;(c)振幅;蓝色线条,绿色线条和红色线条分别表示迭代次数为1,2和3.
Fig. 8 The frequency response of in situ iteration scheme of the different iteration numbers(from left to right) when the superposition step is 4
(a)Real part;(b)Imaginary part;(c)Amplitude. Blue line,green line and the red line are represent for the number of iteration are 1,2 and 3 respectively.

仿照上述做法,对(22)式施行Fourier变换后得到1次迭代的频率响应为

那么K次迭代的频率响应H7(K)(f,L)为

可以看出,H7(1)(f,L)和H7(K)(f,L)都是复函数,并且H7(1)(f,L)和H6(1)(f,L)互为共轭.

同样,对(23)式施行Fourier变换后得到2次迭代的频率响应为

那么K次迭代的频率响应H6(K)8(f,L)为

可以看出,H8(2)(f,L)和H6(K)8(f,L)都是实偶函数.

当迭代达到一定次数后,H7(K)(f,L)和H6(K)(f,L) 的虚部趋向于0,此时H6(K)8(f,L)、H7(K)(f,L)和H6(K)(f,L)的实部相等.随着迭代次数的增加,(21)、(22)和(23)式具有相同的收敛性,但当迭代次数较小时,三种迭代格式的结果会有一定的差别.

图 9为某地一条东西向重力异常,最左边的点号为-64,最右边的点号为104,共有169个点.从-14号到12号这段区域(27个点)由于一些条件的限制没有采集到数据.图 9a给出了迭代步长为3,分别按(9)、(10)和(11)式完成4次原位迭代结果,图 9b给出了叠加步长为3,分别按(21)、(22)和(23)式完成4次原位迭代的结果.从图 9a可以看出,(9)式(从左到右)的迭代结果极值位置向原始位 场极值位置的右方偏移,(10)式(从右到左)的迭代结果极值位置向原始位场极值位置的左方偏移,而(11)式(先从左到右,再从右到左,交替进行)的迭代结果极值位置介于上述两个迭代结果之间,三者之间的偏移量相对较大.从图 9b可以看出,(21)式(从左到右)的迭代结果极值位置向原始位场极值位置的左方偏移,(22)式(从右到左)的迭代结果极值位置向原始位场极值位置的右方偏移,而(23)式(先从左到右,再从右到左,交替进行)的迭代结果极值位置基本位于上述两个之间,三者之间的偏移量相对较小.通过对比可以看出,叠加步长迭代结果相对较稳定,且偏移量较小.

图 9 迭代次数为4时原位迭代结果
(a)单步长为3时4次迭代结果;(b)叠加步长为3时4次迭代结果.
Fig. 9 The results of in situ iteration scheme when the number of iteration is 4
(a)Single step is 3;(b)Superposition step is 3.

通过以上对1-D原位迭代格式的频率响应分析研究认为,从左到右和从右到左的单步长及叠加步长原位迭代格式频率响应是一低通复函数,因此其迭代结果是收敛的,但具有一定的相位移;交替进行迭代的单步长和叠加步长原位迭代格式频率响应是一低通实偶函数,其迭代结果也是收敛的,而且无相位移.通过进一步比较认为,三种单步长和叠加步长原位迭代格式的频率响应均为低通滤波器,故其迭代收敛.当迭代到达一定次数后,三者的迭代结果趋于相同.叠加步长的原位迭代格式比单步长的原位迭代格式具有更好的收敛性或低通性. 3 2-D最小曲率位场分离方法及收敛性

2 -D最小曲率方法的基本差分迭代格式为(王万银等,2009)

其中,,α2=2α2,α3=-4(1+α2),α4=-4(1+α22ΔxΔy 分别为x和y方向上的节点距;M和N分别为x方向和y方向上的节点个数;s和t分别为x和y方 向上的序号.下面给出2-D最小曲率非原位迭代和原 位迭代格式,并利用Fourier频谱分析理论研究其收敛性.

3.1 非原位迭代格式及其收敛性

若用u(k)(l)(sΔx,tΔy)和u(k-1)(l)(sΔx,tΔy)分别表 示步长为lx=l×Δx和ly=l×Δy时第k次和第k-1次迭代结果,则非原位迭代格式就是利用第k-1次的迭代值u(k-1)(l)(sΔx,tΔy)来计算第k次的迭代值u(k)(l)(sΔx,tΔy).非原位迭代又分单步长和叠加步长两种,下面分别给出这两种非原位迭代格式,并利用Fourier频谱分析理论研究其收敛性.

3.1.1 单步长非原位迭代格式及其收敛性

根据基本差分迭代格式得到单步长非原位迭代格式为

这里u(0)(l)(sΔx,tΔy)=u(sΔx,tΔy)表示初始位场值.

对(31)式施行Fourier变换,并利用2-D Fourier变换的位移性质

得到(31)式的Fourier变换结果为

这里H1(1)(fx,fy,l)表示步长为l完成(31)式1次单步长非原位迭代的频率响应,可写为

为了讨论方便,假设Δx=Δy,则α=1,α0=-1/20,α1=1,α2=2,α3=-8,α4=-8.上式可进一步写为

根据采样定理,u(sΔx,tΔy)的截止频率为fcx=1/(2Δx)和fcy=1/(2Δy),则 2πfx×lx ≤lπ,2πfy×ly ≤lπ.现在来研究这个频率响应H1(1)(fx,fy,l)的特点.

(1)奇偶性及虚实性:因为H1(1)(-fx,fy,l)=H1(1)(fx,fy,l),H1(1)(fx,-fy,l)=H1(1)(fx,fy,l),而且虚部为0,所以H1(1)(fx,fy,l)是一实偶函数.由于是偶函数,故只讨论2πfx×lx≥0和2πfy×ly≥0的情况.

(2)周期性:因为H1(1)(fx,fy,l)是由余弦函数构成.变量2πfx×lx的周期为2πn,由2πfx×lx≤lπ得出,n=l/2;当l=2n时有n个周期,当l=2n+1时有n+0.5个周期.同理,变量2πfy×ly的周期为2mπ,由2πfy×ly≤lπ得出,m=l/2;当l=2m时有m个周期,当l=2m+1时有m+0.5个周期.

(3)零值点位置:由于H1(1)(fx,fy,l)是一个二元函数,故其零值点位置较难讨论,可以根据频率响应图形来进行讨论.

(4)最大值:由H1(1)(fx,fy,l)的表达式得到其最大值为

可以看出,H1(1)(fx,fy,l)的最大值大于1,故H1(1)(fx,fy,l)是一放大滤波器.

(5)极值点位置:由极值的必要条件

推出极值点的位置为2πfx×lx=nπ和2πfy×ly=mπ.又由2πfx×lx≤lπ和2πfy×ly≤lπ得出,n=0,1,2,…,l和m=0,1,2,…,l.在2πfx×lx≥0方向上共有(l+1)个极值点,在2πfy×ly≥0方向上也共有(l+1)个极值点.仿照一维滤波响应的讨论可以得到:当2πfx×lx=(2j+1)π(j=0,1,…,[(l-1)/2])时为极小值点,共有[(l-1)/2]+1个极小值点,且极小值为-11/5;当2πfx×lx=2jπ(j=0,1,…,[l/2])时为驻点,实际上就是极大值点,共有[l/2]+1个极大值点,且极大值为1.当2πfy×ly=(2k+1)π(k=0,1,…,[(l-1)/2])时为极小值点,共有[(l-1)/2]+1个极小值点,且极小值为-11/5;当2πfy×ly=2kπ(k=0,1,…,[l/2])时为驻点,实际上就是极大值点,共有[l/2]+1个极大值点,且极大值为1.

若进行K次迭代,则其频率响应H1(k)(fx,fy,l)为

根据H1(1)(fx,fy,l)的性质容易得到H1(k)(fx,fy,l)也是一个实偶函数.H1(k)(fx,fy,l)的周期性与H1(1)(fx,fy,l)的相同.H1(k)(fx,fy,l)的零值点位置和个数与H1(1)(fx,fy,l)的相同.H1(k)(fx,fy,l)的极值点位置和个数与H(1)1(fx,fy,l)有所不同;当K为奇数时,H1(k)(fx,fy,l)的极值点位置和个数与H1(1)(fx,fy,l)的相同,极大值都为1,但极小值不同,H1(k)(fx,fy,l)的极小值是H1(1)(fx,fy,l)极小值的K倍,为(-11/5)K<0;当K为偶数时,H1(k)(fx,fy,l)的所有值均大于等于0,故其极小值点位置和个数与零值点位置和个数相同,但H1(1)(fx,fy,l)的极小值和极大值位置都是H1(k)(fx,fy,l)的极大值位置,极大值为(11/5)K>0或1,个数为(l+1)个.

图 10为迭代步长分别为1,2和3时完成1次单步长非原位迭代的频率响应,其中横坐标为fx×Δx、纵坐标为fy×Δy,二者均为无量纲之量,其最大值均为0.5(图 11图 13中的横坐标和纵坐标与此相同).由该图可以看出,单步长非原位迭代格式的频率响应具有周期性的特征,而且振幅谱大于1.因此单步长原位迭代格式(31)式不收敛,不能用于各种数据处理.

图 10 迭代次数为1时单步长非原位迭代频率响应(实部)
(a),(b)和(c)分别表示步长为1,2和3时的频率响应.
Fig. 10 The frequency response(real part)of the non in situ iteration scheme of the single step when the number of iteration is 1
(a),(b) and (c)are represent the frequency response for the single steps are 1,2 and 3 respectively.
3.1.2 叠加步长非原位迭代格式及其收敛性

仿照1-D最小曲率叠加步长非原位迭代格式,同样可以构造2-D最小曲率叠加步长非原位迭代格式,可以表示为

该式称为叠加步长为L(最大迭代步长)的非原位迭代格式.若对其施行Fourier变换可得到1次迭代的频率响应为

那么K次迭代的频率响应H6(K)2(fx,fy,L)为

由(36)和(37)式可以看出,叠加步长非原位迭代格式的频率响应也是一个实偶函数.图 11为叠加步长分别为2,3和4时完成1次和4次非原位迭代的频率响应.由图 11看出,随着叠加步长L的增大其第一个零值点的频率逐渐趋向于0,且无明显周期性特征.当叠加步长为2和3时,其频率响应振幅大于1.故2-D叠加步长非原位迭代格式不收敛,这 与1-D叠加步长非原位迭代格式的收敛性完全不同.

图 11 叠加步长非原位迭代频率响应(实部)
(a),(b)和(c)分别表示叠加步长为2,3和4时1次迭代的频率响应;(d),(e)和(f)分别表示叠加步长为2,3和4时4次迭代的频率响应.
Fig. 1 The frequency response(real part)of the non in situ iteration scheme of superposition step
(a),(b) and (c)are represent the frequency response for the superposition steps are 2,3 and 4 when the number of iteration is 1; (d),(e) and (f)are represent the frequency response for the superposition steps are 2,3 and 4 when the number of iteration is 4.

通过以上对非原位迭代格式的研究表明,不论是单步长非原位迭代格式或叠加步长非原位迭代格式均具有高频放大作用,迭代不收敛,不能用来进行各种数据处理及位场分离.只是叠加步长非原位迭代格式的收敛性相对于单步长非原位迭代格式要好.不论是单步长或叠加步长,其迭代格式的频率响应均为实偶函数.

3.2 原位迭代格式及其收敛性

原位迭代格式就是利用第k-1次的迭代值u(k-1)(l)(sΔx,tΔy)和第k次的迭代值u(k)(l)(sΔx,tΔy)来计算第k次的迭代值u(k)(l)(sΔx,tΔy).对于1-D数据来讲,原位迭代有三种迭代格式,分别是从左到右,从右到左和交替进行.但对于2-D数据来讲,原位迭代方式较多.通过对1-D原位迭代格式收敛性研究得知,不论是单步长或叠加步长原位迭代格式均收敛,虽受迭代方向选择的影响,但这种影响较小,当迭代一定次数后就不受迭代方向选择的影响.因此,在2-D数据原位迭代格式中,只研究从左到右、从下到上的迭代格式,称为第一种迭代格式.下面分别给出这种单步长和叠加步长非原位迭代格式,并利用Fourier频谱分析理论研究其收敛性.

3.2.1 单步长原位迭代格式及其收敛性

2 -D最小曲率单步长原位迭代格式为

对(38)式施行Fourier变换,并利用Fourier变换的位移性质得到(38)式的Fourier变换结果为

这里H3(1)(fx,fy,l)表示步长为l完成(38)式1次迭代的频率响应,可写为

由于该式比较复杂,难以用解析式来研究其收敛性,故利用图形来研究.不失一般性,假设Δx=Δy,则α=1,α0=-1/20,α1=1,α2=2,α3=-8,α4=-8.(40)式进一步写为

根据采样定理,u(sΔx,tΔy)的截止频率为fcx=1/(2Δx)和fcy=1/(2Δy),则 2πf×lx ≤lπ,2πf×ly ≤lπ.若l=1,当2πfx×lx=0和2πfy×ly=0时,H3(1)(fx,fy,l)=1+i0;当2πfx×lx=0和2πfy×ly= 时,H(1)3(fx,fy,l)= ,H(1)3(fx,fy,l)= <1;当2πfx×lx=0和2πfy×ly=π时,H(1)3(fx,fy,l)=;当2πfx×lx= 和2πfy×ly=0时,H(1)3(fx,fy,l)=,H(1)3(fx,fy,l)=<1;当2πfx×lx=π和2πfy×ly=π时,H(1)3(fx,fy,l)= .可以看出,该频率响应是一个低通响应,其振幅的最大值为1.

若进行K次迭代,则其频率响应H3(K)(fx,fy,l)为

图 12为迭代步长分别为1、2和3时1次单步长原位迭代的频率响应及迭代步长为3时4次原位迭代的频率响应.由该图可以看出,单步长原位迭代格式的频率响应依然表现出一定的周期性特征,随着迭代步长的增大,其低通性能增强,并随着步长的不同而有所变化;随着迭代次数的增加,其低通性能也增强.

图 12 单步长原位迭代频率响应
(a),(b)和(c)分别表示迭代步长为1时1次迭代的实部,虚部和振幅;(d),(e)和(f)分别表示迭代步长2时1次迭代的实部,虚部和振幅;(g),(h)和(i)分别表示迭代步长为3时1次迭代的实部,虚部和振幅;(j),(k)和(l)分别表示迭代步长为3时4次迭代的实部,虚部和振幅.
Fig. 12 The frequency response of in situ iteration scheme of the single step
(a),(b) and (c)are represent for the real part,the imaginary part and the amplitude when the single step is 1 and the iteration number is 1;(d),(e) and (f)are represent for the real part,the imaginary part and the amplitude when the single step is 2 and the iteration number is 1;(g),(h) and (i)are represent for the real part,the imaginary part and the amplitude when the single step is 3 and the iteration number is 1;(j),(k) and (l)are represent for the real part,the imaginary part and the amplitude when the single step is 3 and the iteration number is 4.

通过以上单步长原位迭代格式的频率响应特征分析可以看出,单步长原位迭代格式的频率响应为低通滤波器,故其迭代收敛,可以用于各种数据处理及位场分离.

3.2.2 叠加步长原位迭代格式及其收敛性

仿照叠加步长非原位迭代格式,可以构建叠加步长原位迭代格式.对于(35)式,叠加步长原位迭代格式表达式为

对其施行Fourier变换后得到1次迭代的频率响应为

那么K次迭代的频率响应H(K)4(fx,fy,L)为

图 13分别为叠加步长为2和3时1次迭代的频率响应图及叠加步长为3时4次迭代的频率响应图.由图 13可以看出,该频率响应的振幅小于等于 1,是一低通滤波器;叠加步长迭代的频率响应无周期性特征,而且随着叠加步长的增大其低通性能增强;随着迭代次数的增加,其低通性能增强更快.

通过以上对原位迭代格式的研究表明,不论是单步长原位迭代格式或叠加步长原位迭代格式均具有低通性能,迭代收敛,故可用来进行各种数据处理及位场分离.

图 13 叠加步长原位迭代频率响应
(a),(b)和(c)分别表示叠加步长为2时1次迭代的实部,虚部和振幅;(d),(e)和(f)分别表示叠加步长为3时1次迭代的实部,虚部和振幅;(g),(h)和(i)分别表示叠加步长为3时4次迭代的实部,虚部和振幅.
Fig. 13 The frequency response of in situ iteration scheme of the superposition step
(a),(b) and (c)are represent for the real part,the imaginary part and the amplitude when the superposition step is 2 and the iteration number is 1;(d),(e) and (f)are represent for the real part,the imaginary part and the amplitude when the superposition step is 3 and the iteration number is 1;(g),(h) and (i)are represent for the real part,the imaginary part and the amplitude when the superposition step is 3 and the iteration number is 4.
4 结论与建议

(1)本文利用Fourier频谱分析理论研究了最小曲率位场分离方法差分迭代格式的频率响应,得到了单步长迭代和叠加步长迭代、非原位迭代和原位迭代的频率响应特征及收敛性.实践证明这种研究方法是成功的,而且可以推广到其他差分格式的频率响应特征及收敛性研究中.

(2)本文给出了最小曲率位场分离方法的迭代格式,并通过Fourier频谱分析技术研究了其频率响应特征及收敛性,在实际应用中还需要进一步研究迭代次数和迭代步长的选择问题.通过研究认为,叠加步长迭代格式比单步长迭代格式的收敛性好、而且稳定性好.

(3)除单步长迭代和叠加步长迭代外,还可以选用不同单步长迭代组合进行迭代,即每次迭代采用不同的步长.对不同单步长迭代组合迭代也可以采用Fourier频谱分析技术研究其频率响应特征及收敛性.这里不再赘述.

致谢     感谢论文评审专家和潘作枢教授提出的修改意见,感谢本文编辑对文字的加工和修改.
参考文献
[1] Briggs I. 1974. Machine contouring using minimum curvature. Geophysics, 39(1):39-48.
[2] Carlos A M, Silva B C. 1995. Interpolation of potential-field data by equivalent layer and minimum curvature: A comparative analysis. Geophysics, 60(2):399-407.
[3] Swain C J. 1976. A fortran Ⅳ program for interpolating irregularly spaced data using the difference equations for minimum curvature. Computer & Geoscience, 1:231-240.
[4] Xiong L. 1999. Comparison of some gridding methods. The Leading Edge, 898-900.
[5] Michael D O, Richard S S, Marc A V. 2005. Gridding aeromagnetic data using longitudinal and transverse horizontal gradients with the minimum curvature operator. The Leading Edge, 142-145.
[6] Mickus K L, Aiken C L V, Kennedy W D. 1991. Regional-residual gravity anomaly separation using the minimum-curvature technique. Geophysics, 56(2):279-283.
[7] Smith W H F, Wessel P. 1990. Gridding with continuous curvature splines in tension. Geophysics, 55(3):293-305.
[8] Wang W Y, Qiu Z Y. 2011. The research to a stable minimum curvature gridding method in potential field data processing. Progress in Geophys. (in Chinese), 26(6):2003-2010,doi:10.3969/j.issn.1004-2903.
[9] Wang W Y, Qiu Z Y, Liu J L, et al. 2009. The research to the extending edge and interpolation based on the minimum curvature method in potential field data processing. Progress in Geophys. (in Chinese), 24(4):1327-1338,doi:10.3969/j.issn.1004-2903.
[10] Zhang Z H, Xu S Z, Yu H L, et al. 2011. Study of extending methods of iteration of downward continuation in potential field. Journal of Zhejiang University (Engineering Science), (in Chinese), 47(5):918-924,doi:10.3785/j.issn.1008-973X.
[11] 王万银,邱之云,刘金兰,黄翼坚,于长春,李焓. 2009. 位场数据处理中的最小曲率扩边和补空方法研究. 地球物理学进展,24(4):1327-1338,doi:10.3969/j.issn.1004-2903.
[12] 王万银,邱之云. 2011. 一种稳定的位场数据最小曲率网格化方法研究. 地球物理学进展,26(6):2003-2010,doi:10.3969/j.issn.1004-2903.
[13] 张志厚,徐世浙,余海龙,吴乐园. 2013. 位场向下延拓的迭代法的扩边方法. 浙江大学学报(工学版),47(5):918-924, ,doi:10.3785/j.issn.1008-973X.