位场球谐分析的基本理论.docVIP

2007.5 63第二章 位场球谐分析的基本理论球谐分析是卫星重力和磁场数据分析解释主要工具。地球外部的重力场可以表示成球谐级数形式，同时，通过不同阶次的球谐展开，可以对地球重力场进行分析，以达到显示地球重力场特征，进而研究地球重力异常各种成因的目的。本章主要介绍论文所涉及的球谐分析的一些基本理论。2.1 位场拉普拉斯方程的解可以证明，地球外部引力场是调和的，其满足拉普拉斯方程。在不同的坐标系中，拉普拉斯方程有不同的形式，为了便于讨论，本节分别对直角坐标系和球坐标系中拉普拉斯方程的求解进行讨论。2.1.1 直角坐标系中拉普拉斯方程的解假设V表示地球引力位，其为直角坐标系中空间点(x, y, z)的调和函数，则有2V=0，即 (2-1)用分离变量法解方程，令代入方程（2-1）则有式中X, Y, Z分别为X, Y, Z对x, y, z的二次导数。解方程，得V(x, y, z)的一般表达式其中为待定常数，根据边界条件来确定。2.1.2 球坐标系中拉普拉斯方程的解在球坐标系中，引力位V可表示为空间点(r, q, l )的函数，即V(r, q, l )，其中r为点的坐标径向距离，q 为余纬度，l 为经度，如图2-1所示。引力位V的拉普拉斯方程可表示成 (2-2)图2-1 球坐标系同理，可以用分离变量法解方程(2-2)，即令代入方程(2-2)，则有 (2-2a)上式中只第三项与有关，则第三项是一常数。令 (2-3)其中m是整数，可解得代入方程(2-2)，则有 (2-4)从上式可以看出，左边第一项只与有关，第二项只与有关，则两项都是常数。令 (2-5)这里n为整数。解方程(2-5)，可得当r 时，V有限，可知C=0. 将(2-4)和(2-5)代入方程(2-2a)，则有 (2-6)方程(2-6)解的形式为其中为连带勒让德函数，其中m是整数，是常参数。当m, n取不同的整数值时，方程(2-3)和(2-5)的解是特解，考虑到位场的叠加性质，将所有这些特解累加起来，的一般表达式可写成 (2-7)其中和是常参数。可以证明，地球外部引力位的球谐表达式可以写成 (2-8)其中和是通常所称的球谐系数。2.2 勒让德函数与连带勒让德函数2.2.1 勒让德与连带勒让德函数的一般形式及其递推公式勒让德函数球谐函数的核心组成部分。根据上一节方程(2-6)，把解代入，可得当m0时，把x = cosq 代入，可得到勒让德方程，即 （2-9）解方程(2-9)，可得勒让德函数表达式，即显然，勒让德函数的级数表达形式为 (2-10)其中为中最大的整数值。如图2-2所示，勒让德函数具有以下性质：(1) 当是偶数时是偶函数，是奇数时是奇函数;(2) 在1，+1的区域内有个零点。图2-2 (a) 函数P0(x)P6(x)在区间-1,1上的图形；(b) 函数P6(x)在球面上的图形，其中x=cosq ，灰色和白色分别代表正负区域。当m 0时，把x = cosq 代入，方程(2-9)变成连带勒让德方程，即解方程可得连带勒让德函数微分表达式为连带勒让德函数级数展开形式可写成 (2-11)利用(2-10)及(2-11)式可以计算勒让德函数。然而，计算高阶勒让德多项式时不方便的，需要导出其递推公式。由 (2-12) (2-13) (2-14)由此可知，把m0代入(2-13)即可得到(2-12)，因此(2-13)是一个通式。由此可得，与勒让德多项式导数有关的四个基本递推公式为 (2-15) (2-16) (2-17) (2-18)其中表示对多项式的求导。与连带勒让德函数有关的四个基本递推公式为 (2-19) (2-20) (2-21) (2-22)当时，则令。2.2.2 球函数的规格化勒让德函数与连带勒让德函数在-1，1区间都是非规格化的正交系。可以证明，勒让德函数和连带勒让德函数随着增大函数值也增大，对于连带勒让德函数来所，特别是当时，更是如此，这对实际应用是不方便的。为了避免这种现象，总是将它规格化，即使得这里是规格化后的连带勒让德函数。规格化前后有下列关系为式中为克罗内克符号，当k = 0时值为2，当k 0时值为1.假如有一个函数，它在整个单位球面上（）函数值时已知，如果展成N阶球谐函数，则有 (2-23)其中 上式即为规格化前后的球谐系数之间的互换关系。2.3 球谐变换 从前面的讨论可知，若已知球谐系数，就可以确定球谐函数。但如果已知球谐系数，也可以通过积分获得球谐系数，这种互换的过程通常称为球谐变换。2.3.1 球谐系数的求取根据方程(2-23)，球函数展开成球谐级数，其中区间为0, 2p ，q 区间为0, p ，其表达式为 (2-23a)其中为规格化连带勒让德函数，为规格化球谐系数。假设为两个与有关的参数，根据求取傅立叶级数系数的原理，有 (2-24) (2-25)再利用勒让德转换，得 (2-26) (2-27)积分式(2-24), (2-25)和(2-26), (2-27)就是利用球谐函数确定球谐系数与的过程，此过程称为球谐正变换。2.3.2 球谐函数的求取用球谐系数与计算球谐函数的过程，是与上述用球谐函数计算球谐系数的过程相反，被称位球谐逆变换，若将在球面作N阶球谐展开，其步骤为：(1) 求函数和 (2-28) (2-29)(2) 根据已获得的和，利用傅立叶级数展开式可求取球谐函数，即 (2-30)显然，球谐系数计算中，勒让德转换是十分重要的。实现勒让德转换方法有许多种，常见的有：标准法、直接法、快速多极法等。2.3.3 勒让德转换实现假设一个离散化的球函数，其截断阶次分别为，其余纬度必须满足. 首先求取满足的，可采用牛顿二分法计算。若是等间距的，求球谐系数第一步即采用公式(2-24)和(2-25)，第二步即采用公式(2-26)和(2-27)。离散化的求解球谐系数公式(2-26)和(2-27)可写为： (2-31) (2-32)其中是加权，区间为-1，1。勒让德转换可在高斯网格上进行转换，也可以在任意网格上进行转换。首先来看高斯网格上求解过程。求解权函数的步骤如下：由 (2-33) (2-34)用矩阵表示，即 (2-35) (2-36)其中，是维加权矩阵，是维矩阵，即勒让德综合表达式为。与连续勒让德函数一样，离散化勒让德函数在高斯积分下也是正交的，即则可反推得，即有 (2-37) (2-38) (2-39)从(2-39)式中选择满足的，矩阵是一个对角矩阵，其余各项为零，即有按照勒让德相关的推导公式可以推出其它计算的表达式。下表为截断阶次N=8,m=1时系数之差，其中coeff，经过公式计算出，然后通过公式计算新系数newcoeff。表21 截断阶次N=8,m=1时系数之差coeffnewcoeffCoeff- newcoeff111.554E-1522-2.665E-1533-8.882E-1544-7.105E-1555-4.441E-1566-2.665E-1577-6.217E-15等角网格上的勒让德转换满足纬度数据是截断阶数的2倍，即令：。有 (2-40) (2-41)其中权函数表2-2为截断阶次N=8,m=1时系数之差，其中coeff，经过公式计算出，然后通过公式计算新系数newcoeff。表22 截断阶次N=8,m=1时系数之差coeffnewcoeffCoeff- newcoeff11-2.068E-0822-2.558E-0833-6.67E-0844-4.25E-0855-1.271E-0766-3.07E-0877-2.739E-0787.9634350.0365655从表2-1和2-2可以看出，相同截断阶次N=8，而m=1时，所求系数与原来系数之差所处的数量级不在同一水平上，高斯网格法所求精度为，而等角网格法所求精度为，特别是最后一个系数相差较大。因此高斯网格法优于等角网格法。第三章 卫星重力异常处理方法研究由于重力异常所反映的密度不均匀体来自于上自地表，下至上地幔的多种因素叠加的结果，此外，测量数据还存在着各种干扰和误差；而我们所研究的对象往往是相对具体的、明确的。所以，在开展解释工作之前，需要对重力异常数据进行适当的处理。重力异常资料处理主要包括三个方面的内容：压制误差与干扰、异常区分、异常转换，其目的都是为突出研究对象的信息（信号）、压制其它方面的影响。对于局部（或一定范围上区域）重力资料而言，可以将重力观测面近似为平面，所有资料处理问题都可作为直角坐标系中问题进行理论推导和计算。然而，对于更大区域内的，乃至全球问题，重力观测面已不能近似为平面，而是一个球面问题，需要用球坐标进行表达。这时，建立在直角坐标系的数据分析与处理方法已不能适用。3.1 场源外部重力异常的表达3.1.1 直角坐标系中的重力异常重力异常的实质就是目标地质体与其周围岩石密度差异所形成的“剩余质量”对测点处所产生的附加引力在重力方向上的分力（或投影），因此，异常重力可以通过对场源积分获得，如图3-1所示。异常引力位积分表达式可写成 （3-1）式中dm 为剩余质量积分元，Dt (x,h,z ) 为剩余密度或密度差。根据定义，重力异常为异常重力位沿重力方向的导数，在直角坐标系中，重力方向可视为Z轴方向，即重力异常积分表达式为 （3-2）当地质体为均质时，其密度为常数，即(3-2)式可写成 （3-3）图3-1 直角坐标系中场源积分示意图根据(3-2)或(3-3)式，若已知地质体或剩余密度的的空间分布，便可计算出空间任意一点的重力异常。通常情况下，由于地质体或剩余密度是未知的，无法利用(3-2)或(3-3)式计算空间重力异常，但通过地面上测量到重力异常值进行外推。按照斯托克斯定理，地面以上空间任意一点的引力值可以通过对包围场源的曲面积分而获得。根据位场边值问题狄义赫利问题的解，若已知边界面上重力值异常，便可以给出场源外部空间任意一点的重力异常表达式。假设地面是水平的，且与XOY平面重合，如图3-2所示，不难得出平面以上空间任意一点重力异常，即 （3-4）式中S为无限平面，且有z0，即P点在XOY平面以上。理论上，确定空间重力异常，需要对整个无限平面上的异常场值进行积分。图3-2 直角坐标系中平面积分示意图3.1.2 球坐标系中的重力异常在球坐标系中，同样可以给出(3-3)或(3-3)式类似的异常引力位积分表达式，即 （3-5）当异常体为均质时，有 （3-6）显然，当异常体可近似为一个点源时，同样可以给出场源外部重力异常的积分表达式。但当场源不能近似为点源时，在给出具体异常体形态之前，各处的引力方向无法确定，以致不能直接给出重力异常的表达式。如图3-3所示，根据斯托克斯定理，若已知地球面表的（异常）引力位，同样可以给出地球球面（异常）引力位的积分表达式，即 （3-7）上式称为泊松积分公式。图3-3 球坐标系中球面积分示意图就地球重力场而言，重力位W与正常重力位U之差被称为扰动位。物理上，扰动位是指由对地球内部物质密度偏离“正常”密度（或平均密度）而成形的剩余质量以及地球表面相对地球椭球面起伏的物质共同产生的引力位。对于地球外部任一点P (r, q, l )，扰动位T 可用球谐函数表达21，即 （3-8）式中R表示地球椭球长半轴， 分别为规格化的球谐系数，为规格化的连带勒让德函数。对于卫星重力数据而言，R可近似为地球平均半径。根据Molodensky 边值问题的边界条件，具有球体近似的重力异常球谐表达式为 （3-9）可见重力异常球谐表达式与（3-8）式的扰动位表达式具有相类似的形式。利用球函数的正交性和规格化的定义，可根据球面S上重力值异常Dg(q, l )，计算球谐系数21，即 （3-10）需要说明，(3-9)式没有考虑垂向偏差，将大地水准面近似为地球椭球面。从(3-10)式可知，球谐系数和可通过已知的地球表面上的异常Dg(q, l )积分来确定。若S 为已知空间某一球面（如所有卫星轨道确定的球面等），则球面半径变成r ，(3-10)式中球面重力异常应改为Dg(r,q, l )。3.2 傅立叶波数域谱与球谐函数谱的特征3.2.1 傅立叶变换及其性质傅立叶变换(FourierTransformation)被广泛运用于各种信号处理与分析。对位场资料而言，直角坐标系中的许多数据分析与处理的解析表达，都可以映射到傅立叶频率域中。也可用频谱(波谱)表示，在频率域(波数域)中表达一般比空间域简单。若把位场异常的空间变化视为周期无限大的周期函数，于是可把位场异常分解成为各种频率的谐波，这些谐波幅度不是随时间变化，而是随空间变化。各种频率的谐波又具有各不相同的振幅和初位相，因此我们可以把位场异常看成是这种谐波所组成。这些谐波的振幅和初位相是频率的函数，此种关系分别称为振幅谱和相位谱，它们又统称为频谱。我们把位场异常转换成频谱来进行解释的方法称为频谱分析法。由于这里信号不是以时间尺度而是空间尺度进行描述，所以，在后面叙述中将“频率”改为“波数”。用一组空间谐波来表示位场异常，在数学上称为傅里叶展开。假设在一条长为2L的剖面上测得异常为f (x)，f (x)是以2L为周期的周期函数。最简单的空间波为正弦波(或余弦波，无本质差别)可写为 f (x) = A sin(wx+ )，其中A为振幅；为初位相；w为角波数(与角频率相当)，w =2 f， f为波数(与频率相当)，f =1/2L，2L为波长（与周期相当）。任一复杂的异常f (x)可以由不同频率的简单正弦波叠加而成，故f (x)可表示为 （3-11）经展开后得傅立叶级数，即 （3-12）其中ak, bk称为傅立叶级数系数，它们与振幅Ak及k之间的关系为如果在一条长为2L的剖面上测得异常为f (x)，则可确定傅立叶系数，即 （3-13）由此可以计算出Ak及k (k=0, 1, 2, 3,)。 数列Ak和k分别称为f (x)的振幅谱和相位谱。应用欧拉公式，(3-12)式可写成复数形式： （3-14）式中 （3-15）上面我们假定了异常f (x)是以2L为周期的函数，这样做往往与实际情况不符合。一般f (x)在剖面2L以外是为零的，即这时x的变化范围就可以是(-, +)，即为非周期函数。由于异常f (x)绝对可积，且在任何一个区间内都是有界的，只有有限个不连续点和有限个极值点，故可表示成傅里叶积分。由傅氏级数的复数形式(3-14)式出发，考虑到2L；经变换后，到傅里叶积分形式 （3-16）其中F (w)称为f(x)的谱函数，即 （3-17）由此称f (x) F (w )为的傅立叶变换，记作F f (x)，称F (w) f (x) 为傅立叶逆变换记作F -1 F (w )。考虑到w =2 f，(3-16)和(3-17)式可表示成 （3-18）由于F (w)是复函数。有实部和虚部，故其模A (f ) 称为振幅谱，其幅角f (f ) 称为相位谱，即 （3-19）由(3-18)式说明，异常f (x) 可分解为无穷多个频率f 谐波叠加，这些谐波的振幅和初相位由(3-18)表示。以下是几个傅立叶变换的性质：(1)叠加性若异常为N个异常叠加而成，即即其谱也可叠加，则有 （3-20）(2)微分定理若异常函数f (x)的谱函数为 F (f )，对异常函数求导，其谱的形式为 （3-21）(3)褶积对于两个任意函数f1 (x)和f2 (x)，其褶积函数f (x)定义为 （3-22）褶积函数f (x)的谱为 （3-23）以上讨论的是一维空间异常场的问题。实际上，我们更多地需要研究二维甚至三维空间的问题。对于二维空间异常f (x, y)，可以证明其傅立叶变换为 （3-23）式中u, v 分别为x, y方向上的波数，F (u, v)为f (x, y)的二维谱。对于三维空间异常f (x, y, z)，则傅立叶变换一般表达式可写成28 （3-24）通常，我们总是假设观测面为水平地面，异常数据是在XOY平面上，三维空间异常是利用二维空间异常f (x, y)进行外推而。所以，来由于，三维空间异常f (x, y, z)可以通过二维空间异常f (x, y)换算得到。相应的傅立叶变换也可以通过这个途径获得，其过程如下：由于场源外部位场异常f (x, y, z)满足拉普拉斯方程，则有 （3-25）对于f (x, y, 0)同样满足上式。由叠加性质(3-20)式和微分定理(3-21)式可知 （3-26）令 w = (u2+v2)1/2，并考虑(3-25)式，则(3-26)式可写为因此，位场异常f (x, y, z)的谱可表示为F (u, v, 0) e 2 w z (z0)。空间异常f (x, y, z)的三维傅立叶变换由(3-25)式变成 （3-27）显然，位场异常f (x, y, z)对x, y, z的导数及空间z = -h平面上的谱分别为 （3-28）对(3-28)式右边作逆变换便可得到异常对各方向上的导数和空间延拓结果，这就是利用傅立叶变换波谱分析法对异常进行分析处理。与(3-14)和(3-15)式类似，二维傅立叶级数形式为 （3-29）其中 -Lx x Lx ，-Ly y Ly 为函数f (x, y)定义域，其复数形式为 （3-29a）3.2.2 球函数谱的概念全球的重力场或重力异常场通常是通过球谐函数表达的，如(2-8)和(3-9)式。对于球函数(2-23a)式，球谐系数可根据球面场值计算出来，如(2-24), (2-25)和(2-26), (2-27)。这与二维位场异常傅立叶级数展开形式(3-29a)式所给出的形式类似，因此，球函数谱的概念也可以类似地给出29。对于二维傅立叶谱F (u, v)而言，当v不变时，F (u, v)相当于函数f (x, y)在x方向上谱，而当u不变时，F (u, v)相当于函数f (x, y)在y方向上谱28，所对应的振幅谱Ax (v), Ay (u)以及径向振幅谱A分别为 （3-30）用级数形式给出，即有 （3-31）其中分别为在x, y方向上和径向上的振幅谱。根据球函数理论29，函数f (, )可展开成 （3-32）其中为球谐函数，Cm,n为复球谐系数，即 （3-33）式中am,n, bm,n为球谐系数。现在讨论一下球函数的几何意义。从前面的讨论可知，n阶球函数可以写成 （3-34）图3-4 连带勒让德函数 (k=0,1,2,7)曲线簇下面我们分三种情况来讨论。1) 先讨论Pn (cos)。根据Pn (x)的性质可知，它在-1, +1区域内具有n个零点，即有n个不等的实根，因此它将球面划分成n+1个环带，但各环带的余纬度角间距不相等，其间距取决于Pn (x)的零值间距。在这些环带上Pn (x)的数值正负交错。对于偶阶的Pn (x)来说，因为它是偶函数，所以它的数值对称于赤道。对于奇阶的Pn (x)来说，因为它是奇函数，所以它的数值表现为南北半球不对称。由于Pn (x)具有这样的性质，因此称它为带球函数，相应的系数an称为带系数，实际上，。2) 再来讨论下面连带球函数对于n = m的情况。 或 由于连带球函数是由勒让德连带函数与的正弦(或余弦)函数的乘积所组成，当n = m时，勒让德连带函数只在0和1处为零，如图3-4；由此，连带球函数值是否为零，取决于sin(n) 或cos(n)函数何处取零值。显然，sin(n) 或cos(n)在的区间0, 2内有2n个零值，使连带球函数为零；对应着的这些经线将球面划分为2n个扇区，各扇区经度角间距相等。在扇区内sin(n) 或cos(n)是正负交错的；这种类型的连带球函数称为扇球函数，相应的系数称为扇系数。图3-5 连带球函数的几何意义示意图3) 现在讨论第三种情况，即对于n m (0 m n) 的连带球函数情况。根据勒让德多项式的性质，方程有n - m个实根，使函数为零，如图3-4；这相当于在nm个卯酉圈上使连带求函数为零，如此将球面划分成nm+1个正负相间的环带，此外，根据前面的讨论，又有2n个子午圈将球面划分为2m个正负相间的扇区，如此经纬网将球面划分为若干个球面四边形，形如“田”字状，但在两极地带为球面三角形；各个球面四边形（或球面三角形）之间正负交错，这种连带球函数称为田球函数(或格球函数)，相应的系数称为田系数。以上三种情况的几何意义见图3-5，图中阴影处为负值区。综上所述，第n阶球函数共有2n+1项，其中有一个带函数，两个扇函数以及2(n-1)个田函数。如果把函数f (, ) 从0阶展开到n阶，则其项数共有1+3+5+(2n+1) = (n+1)2项。从上述球函数的几何意义可以更清楚的看出，所谓将函数展开成球谐级数，就是用具有上述几何意义的各种简单的周期性函数的迭加来近似地表示它。各种球函数前的系数分别称为带、扇和田球函数系数，它们的大小分别确定了各种球函数在带、扇以及田之间正负交错量的数值大小。由此可见，各种球谐系数与(3-31)式给出傅立叶级数系数类似，都是对应阶次上的幅值。若令(3-33)式给出球谐复系数cm, n，的模为振幅30，则有 （3-35）序列即为函数的球谐振幅谱。3.2.3 球函数谱与波长从二维傅立叶波数域的谱和球谐函数谱的讨论中不难发现，直角坐标系中的二维傅立叶波数域的谱反映了平面上不同波长的振幅或能量，而球谐振幅谱则反映了球面上不同角距的振幅或能量。 （3-36）或 （3-37）这里k (k = 0, 1, 2, 3, )，(k/Lx)和(k/Ly)代表两个方向的波数，k/Lxy可视为复合波数。波数越大，则波长越小。对于球谐函数来说，不同阶次的球谐函数，分割的网格不同，阶次愈高，分割得愈密。从网格的纬向分布来看，由卯酉线分割成带的数目由n阶的零值点数决定，其中不包括两极处的零值点。当m = 0时，有n个零值点，球面被分成n+1个带，而当m = 1, 2 , , n 时，的零值点数为nm，带的数目为n-m+1，如图3-4所示。显然，最大带数为n+1。由于这样分割下各带带宽并不相等，但可以想象，南北半球带宽的分布是对称的。对于不同地m，带宽极距间隔近似为/(n-m+1)，其变化范围在/(n+1) 之间。这里，最小间隔与最大间隔分别为/(n+1) 和。从网格的纬向分布来看，球函数中的谐波函数cos (m) 或sin (m) 的零值点数决定着子午线分割扇区的数目。当m = 0时，没有零值点，即不分割；当m = 1, 2 , n 时，零值点数为1 2n，其分割的扇区数目也对应为1 2n，可见，扇区最大数目为2n。但每个扇区的经度角距间是相等的，即为 /m，其变化范围为 /n 。显然，最小间距与最大间距分别为 /n和2。就n阶球谐系数而言，在卯酉线方向与子午线方向的分割角距有不同的几何尺度，所表达的含义也有区别。具体讨论如下：1) 随着卯酉分割线偏离赤道越远，相邻子午分割线之间的球面距离越小，其构成的分割网格球面面积也就越小，在两极地区则演变为球面三角形且面积最小，各处网格面积显然不相等；2) 当n确定后，随着m的增大，卯酉线分割密度逐渐减小，即网格距逐渐增大，而子午线分割密度逐渐增大，即网格距逐渐减小。这表明随着m的增大，球谐系数既包含了子午线方向波长由小到大的成分，又包含了卯酉线方向波长由大到小的成分，体现了两种不同的效应。对比傅立叶波数域谱与波数(或波长)的关系可知：1) 傅立叶(离散)谱是随波数(或频率)而变化，波数的倒数即为波长，也可以说是随波长而变化；而球谐谱是按阶次数(或球面角度分割数量)来定义的，随阶次不同而变化；2) n阶傅立叶(离散)谱(模)仅代表着n倍基波数（或1/ n基波长）的振幅成分；而n阶球谐谱既包含了n阶所对应空间尺度(弧长)的振幅成分，又包含了n阶中0 n次所对应空间尺度(弧长)的振幅成分；3) 阶数n和次数m是球谐谱的两个变量，在半径为R的球面上，阶次数(几何地)对应着分割角间隔或球面距离(弧长)，因此，可以将球谐系数视为弧长的函数。研究n阶m次球谐函数的波长，则需要研究三个方面的问题：其一，如何构建波长序列，即给出一系列由基波长至最大波长的集合；其二，需要从卯酉线和子午线两个方向来考虑，即建立描述球面二维空间尺度，以便研究球面各向异性问题；其三，如何利用波长序列所张的空间，刻画不同尺度球面的函数特征。首先来看子午线圈上的半圆周分割问题。前面的讨论已经提到了午线方向的分割问题，即最小分割极距角为 (n+1)，因为n阶勒让德函数的零值卯酉线可把圆周分成(n+1)段，此时m = 0；而最大分割弧段是半圆周，即R，此时m = n；当0 m n时，随着m增大，弧段长度由最小变到最大，形成一个与m, n有关的序列。由于零值卯酉线并不是以等余纬度间隔划分球面，即划分的带宽不同，所对应的弧段长度也不同，但可以近似地将零值卯酉线看作以等余纬度间隔划分球面，即等分的弧段长度为R (n-m+1)，这对球谐谱的分析不产生实质性的影响。如此弧段长度仅随n, m不同而异。相同的n-m之差会有相同的弧段长，但由此可形成一个弧段长序列，由大到小依次为R, R/2, R/3, R/4, R/(n+1)。(3-38)式给出了午线方向弧段长度的计算方法，即 （3-38）其中为球谐系数在子午线方向弧段长。现在来讨论卯酉圈上的圆周分割问题。从前面的讨论可知，对于确定的n和m，午线方向弧段长度是可以唯一确定的，但卯酉圈上弧段长度不能被唯一地确定。球面上的卯酉圈被谐波函数零值子午线分割成弧段，在半径为R的球面上，当m = 0时，有|e-im| = 1，球面上没有被分割，理论上，球面函数为球带函数，弧段长为2R。当0 m n时，球面被球函数零值子午线分割成2m个相等的扇区，赤道线上被分割的弧段长度为R/m，为相同m的最大弧段。由于午线圈已被分成了n-m+1个正负相间带，则造成在不同的余纬度上切割的卯酉圈弧段长不同，逾向两极接近，相邻的零值子午线之间距离愈小，则卯酉圈上被分割的弧段也就愈小，减小程度取决于弧段所处的余纬度正弦函数值sin。可见，最小弧段应该位于两极最近的弧段，即(R/m)sin/(n+1)。对于确定的n和m，每个扇区(无论正负)都有n-m个长度不同的弧段。为了便于分析，可以考虑将扇区内各弧段进行平均，以各弧段平均长度作为相应的几何尺度，即但由于子午方向最小弧段仅为R/(n+1)，而按照上述平均弧度长计算，就150阶球谐而言，最小平均弧度长可小于200m，这显然与实际的分辨率不符。因此，对于分割网格边长小于子午方向最小弧段长的信号成分，可以作为弧段长等于或大于R/m来考虑。这样，卯酉圈上被分割的弧段长仅与m有关。随着m增大，弧段长度由最大变到最小，形成一个与m, n有关的序列。(3-39)式给出了弧段长度计算方法，即 （3-39）其中为球谐系数在卯酉线方向弧段长。如图3-6所示，和分别为n阶m次田球函数网格的尺度，它们构成了球面上是两个空间尺度序列，且彼此正交。图3-6 球田函数网格边弧长示意图平面直角坐标中，离散数据分布在以Dx和Dy为间距的网格上。在讨论离散傅立叶波数域谱时，最大波数为数据间隔的倒数(1/Dx或1/Dy)，最小波数为数据区域边长倒数1/(M-1)Dx或1/(N-1)Dy，M, N分别为两个方向上的数据点数，它们分别对应于两个方向上平面数据的最短波长成分和最长波长成分，数据中任何波长的信息，都与谱中相应波数的成分互映。对于球面而言，n阶m次球函数间球面划分成间距为R/(n+1)和(R/n)sin/(n+1)的网格，由球函数几何意义可知，球谐系数为对应网格中的幅值，因此，弧段序列和可以用于刻画球面函数各种波长的信息，作为分析球面函数的几何尺度。3.3 球谐谱滤波3.3.1 傅立叶波数域滤波的特点把重力异常数据视为输入信号，设计某种数据分析模型或算法，通过对输入信号进行处理得出结果（输出信号），以达到某种目的。这种过程也被称为数据滤波。对于位场数据，在空间域中，数据滤波可以通过褶积来实现，如最小二乘平滑、滑动平均等。数据滤波也可以通过傅立叶变换，在波数域中用异常数据的谱与滤波器的谱函数相乘来完成，如3.2.1所述。这种方法更方便、更快捷。(a) 原始重力异常(b) 异常频谱图(c) 分离后的深源场(d) 分离后的浅源场图3-7 傅立叶波数域中匹配滤波区分浅源场与深源场应用实例.在波数域中，由于某些波长的异常信号反映了一定的地质目标，可以根据处理的目的选择或设计滤波器，如利用汉宁窗过滤掉高频（短波长）的成分、利用匹配滤波器将浅源场（短波长成分）与深源场（长波长成分）信号分开、设计提取一定波长范围信号的带通滤波器，等等。图3-7展示了在傅立叶波数域中利用数字滤波技术处理某地重力异常的应用实例。除了滤波之外，在傅立叶波数域中还可以进行异常转换，详见3.4、3.5节。3.3.2 球谐谱波长滤波及实现针对m, n的球谐谱滤波技术已经被应用于许多方面31-36，但依据阶次的滤波，难以给出几何的意义。而依据波长（或波长倒数）滤波，则有利于对处理结果的分析以及滤波器的设计。傅立叶波数域滤波是通过滤波器函数对异常场的谱进行处理，在离散傅立叶变换中则是相当于对傅立叶级数的系数进行处理。与离散傅立叶滤波相似，球谐谱的滤波也是对球谐系数进行处理，其处理方式与傅立叶波数域中类似。球谐函数关于n的振幅谱K(n)可以通过球谐系数计算。例如，根据(3-9)式给出的重力异常球谐式，有 （3-40）其中am,n和bm,n都是规格化后的球谐系数。图3-8球谐谱关于阶次n滤波原理示意图图3-8展示了球谐谱关于阶次滤波过程中设置截断高阶、去除低阶和保留带阶的方法，图中虚线为对滤波因子处理后的滤波器。若将(3-35)式中的Km,n视为m, n的函数，即为二维振幅谱，用波长表示即可写成 （3-41）这里 Ds, Ds 分别记作波长，为卯酉方向及子午方向上的弧长的2倍，弧长可按(3-38)和(3-39)式计算。而径向波长则为。(a) 球谐谱关于波数滤波的原理示意图(b) 球谐谱关于波长滤波的原理示意图图3-9球谐谱关于波数或波长滤波的原理示意图在图3-9(a)和(b) 描绘了球谐谱关于波数和波长滤波过程中设置“高通”、“低通”和“带通”滤波器的方法；按照通常的习惯，图(a)中将横坐标用1/Ds r表示。球谐谱滤波因子h(1/Ds r)或h(Ds r)的选择，可根据需要来设计。“高通”、“低通”和“带通”滤波器的设计如图3-9中所示，即在希望保留的波长或波数范围内滤波因子h(1/Ds r)或h(Ds r) 取1，而在希望截去的信号所在的波长或波数范围内h(1/Ds r)或h(Ds r)取0，为了减小“吉布斯效应”，可以在滤波因子的跳变处进行“镶边”处理，如图3-8和3-9中的虚线。3.4 球谐重力异常导数换算在重力场异常资料分析和处理中，异常导数换算具有特殊的意义。例如，在研究局部地区卫星重力异常特征问题时，许多异常特征往往不清晰，需要进行一些处理，对异常进行求导，可以突出更多的异常细节，有助于对产生异常原因的理解和异常的解释。按照求坐标的三个方向，球面重力异常导数换算分为三种情况，即径向导数、纬向导数和经向导数。具体讨论如下。3.4.1 球谐重力异常球面径向导数换算重力异常对径向求导相当于对向径r的求导，它反映了重力异常在半径方向上的梯度变化特征。由地球重力异常的球谐表达式(3-9)，可以推导出重力异常对半径方向的一阶导数表达式，即， (3-42)由此可知，重力异常对半径方向的一阶导数与相比较重力异常，只是将重力异常球谐系数乘上一个因子 (-n/ r)，其它参数没有改变。由此，可以对k阶导数进行类推，即重力异常对半径的k阶导数表达式为 (3-43)不难得知，重力异常对半径方向的k阶导数，将重力异常球谐系数乘以因子。3.4.2 球谐重力异常球面纬向导数换算重力异常对纬向求导即为对余纬度方向弧长sq 求导，它可以反映了重力异常在纬度方向上的梯度变化特征。由地球重力异常的球谐表达式(3-9)，可以推导出重力异常对纬向一阶导数表达式，即 (3-44)其中连带勒让德函数对余纬度的一阶导数的计算可以用(2-11)及(2-22)式给出的方法。重力异常对余纬度方向的一阶导数与重力异常相比较，就是将球谐系数乘上一个因子，连带勒让德函数换成连带勒让德对余纬度的导数。同理，重力异常对余纬度方向的k阶导数表达式为 (3-45)其中连带勒让德函数对余纬度的k阶导数的计算方法可以由(2-22)式递推求得，其二阶导数的计算具体计算方法为 (3-46)3.4.3 球谐重力异常球面经向导数换算重力异常对纬向求导类似，重力异常对经向求导即为对经度方向弧长sl 求导，它可以反映了重力异常在经度方向上的梯度变化特征。由地球重力异常的球谐表达式(3-9)，可以推导出重力异常对经度方向的一阶导数表达式，即 (3-47)显然，与重力异常球谐式(3-9)相比，对经度方向的一阶导数球谐系数由变成。如果把球谐系数写成复数形式，则有其中i表示单位虚数。重力异常对经度方向的一阶导数则可表示为 (3-48)由此可知，重力异常对经度方向的一阶导数与重力异常相比较，就是球谐系数复数形式乘上。重力异常对经度方向的k阶导数表达式为 (3-49)3.5 球谐重力异常空间延拓位场空间延拓是位场资料处理另一种重要的方法，它可以通过换算，展示场源以外位场的空间分布特征。位场空间延拓通常分为向上延拓和向下延拓。向上延拓相当于远离观测面，随着延拓高度的增加，场的局部特征逐渐“淡化”，而区域性的特征逐步得以突出；向下延拓则相反，随着延拓高度的降低，区域性的特征逐渐“淡化”，场的局部特征逐渐突出。理论上，