《地球重力学》-重力模型的确定与应用课件.ppt

1、1地球重力学2重力模型的确定与应用3地球重力模型地球重力模型及其作用地球重力场在地球物理学、海洋学和空间技术中占有特别重要的地位。它直接反映地球内部的密度分布。从地幔产生的长波信号，到大陆岩石圈和海底地壳的局部特征等，都反映在地球重力场中。用一组重力位系数来表示相应的地球重力场，称为地球重力场模型。它的作用可简单归纳为以下几点：卫星大地测量定位的精度取决于卫星定轨的精度，而全球重力场模型是精密定轨的基础。通过地球重力场模型及对地球外部重力场的分析，可为地球物理学和地质学提供地球内部结构和状态的信息。地球重力场模型可精确确定地球的扁率。各国的区域性坐标系与全球坐标系的精确转换，需要区域性大地水准

2、面资料，而大地水准面属地球重力场的一个等位面。大地测量观测是在地球重力场内进行的，数据的处理和归算要知道地球重力场。人造卫星、洲际导弹轨道的摄动与地球外部重力场密切相关。重力勘探是重力学原理在勘探地下资源方面的应用，根据局部重力场变化规律可以反推矿藏位置和范围。4地球重力模型的发展发展历程地球重力模型51.重力场展开与快速计算 6重力场的球谐函数展开7函数的勒让德展开( )ix正交展开方法设为一组归一化的正交多项式系（基函数），则且系数。这里的表示为函数与的内积。一般地取用有限展开，即一般会很小。有时用高阶展开式确定，也有时用小范围内的积分确定，也有的用一个内插值或一个经验值表示，在精度要求不

3、太高的情况下此项往往舍去。0( )( )kkkf xAx( ),( )kkAf xx0( )( )Nkkkf xAxff勒让德函数展开对一维问题。若，取为勒让德多项式，即取有11x ( )( )nnxP x0( )( )Nk kkf xA P x1121( )( )2kkkAf x P x dx8函数的勒让德展开勒让德函数展开对二维球面函数。若，；可用展开公式对二维勒让德函数的展开，一般称为（面）球函数展开。00200( , )cossin(cos )NnnmnmnmnmfAmBmP cos1( , )(cos )sin4nmnmnmmAfPdmB 9重力场的球谐函数展开基本公式 我们知道，地

4、球的引力位可以表示为地球的扰动位可以表示为大地水准面高与重力异常同样可写为00),(),(nnmnmnmnmnmnSDRCrarGMV20*),(),(nnmnmnmnmnmnSDRCrarGMT20*),(),(nnmnmnmnmnmSDRCRN20*),(),() 1(nnmnmnmnmnmSDRCng说明：系数由地球重力模型给定，求和从2开始是满足“质量相等、质心重合”，尽管地球模型的不同，但求和教师有限的，面球谐函数的计算稍后讨论。10勒让德函数的计算方法 11勒让德函数的计算由大地水准面高的计算公式中面球谐函数为20*),(),(nnmnmnmnmnmSDRCRN( , )cos(c

5、os )nmnmRm P ( , )sin(cos )nmnmSm P 并且其中，。不难导得不同的递推公式： ,0,0,()!(cos )21(cos ),(cos )2(21)(cos )()!nnnknkn kPnPPnPn k21,1( )3 1Pxx0,0( )1Px 列递推公式：,1,2,(21)(21)21(1)(1)( )( )( )()()23()()nknknknnnn kn kPxxPxPxn k n knn k n k 行递推公式：,1,2222(2)(1)( )( )( )(1)()(1)()1n kn kn kxkn kn kPxPxPxn kn kn kn kx 1

6、2勒让德函数的计算正向横推法 其基本思路(见图)是利用固定列的勒让德函数递推公式计算较高阶的勒让德函数。求出部份0,1，2阶项，先求出纵向的0,1列（红线），其公式为2,1,1,( )( )( )nknknkn kn kPxaxPxbPx1,(23)(25)(2)(2)nknnan kn k ,(25)(1)(1)(21)(2)(2)n knn kn kbnn kn k 其中：自第3列开始，依次用下式从左向右递推（称为正向横推方法）：,2,1,222(1)()( )( )( )(2)(1)(2)(1)1n kn kn kxknknkPxPxPxnknknknkx 说明：后式中若n=k时首项收敛

7、末项为0，若k固定n时首项为0末项系数为1，n固定k不可能（kn）。13勒让德函数的计算正向纵推法 正向纵推方法是先求得对角线上元素，再向下求出其它的勒让德函数(见图)。先推导出对角线上的函数的公式：导得利用对角线上的函数值推导出紧挨对角线下方(同列)一函数的值，即 再用列递推公式求得每列余下的函数值： 122,1,121( )(1)( )2n nnnnPxxPxn,11,1( )21( )n nnnPxxnPx,1,2,(21)(21)21(1)(1)( )( )( )()()23()()n knknknnnn kn kPxxPxPxn k n knn k n k 说明：后式中若k固定n时首

8、项系数为2x末项系数为-1，若n且k时首项系数为0末项系数为1。（若2x1方收敛）14勒让德函数的计算反向横推 该方法是正向纵推方法的改进(见右图)。先推导出对角线上的函数的公式：导得利用对角线上的函数值推导出紧挨对角线下方(同列)一函数的值，即 再用列递推公式求得每列余下的函数值： 122,1,121( )(1)( )2n nnnnPxxPxn,11,1( )21( )n nnnPxxnPx说明：后式中若k固定n时首项系数为0末项系数为-1，若n且k时首项系数为0末项系数为1。,1,2222(2)(1)( )( )( )(1)()(1)()1n kn kn kxkn kn kPxPxPxn

9、kn kn kn kx 15展开级数的快速计算 16级数的快速计算常规计算方法首先，在实际计算过程当中，大地水准面高的计算公式中的地球位系数的个数不可能是无穷多个，它总是一个有限值，即写为在上式中我们不难发现，计算某上点的大地水准面高需要多次计算三角函数，可以证明需计算三角函数计算量十分庞大。为避免这种重复，交换求和的顺序，即max*20(cossin)(cos )NnnmnmnmnmNRCmDmP maxmax*0(cossin)(cos )NNnmnmnmmn mNRCmDmP再整理为上式中不但三角函数的计算量较小，而且它对勒让德函数的计算量也没有增加。)(cossin)(coscosma

10、xmaxmax0*nmNmnnmnmNmNmnnmPDmPCmRN)sin),(cos),(max0mmBmmARNNm17级数的快速计算基本公式与分析1979年,利兹（Rizos）描述了二维均匀间距的地理网格上用位系数计算大地水准面高或重力异常的计算机技术，常称为利兹法。他编制的程序,是为一个由经纬度所限定的区域中,计算已确定了经纬度间隔的网格点值而设计的。这个计算区,可以是局部区域,也可以是全球。大地水准面的高的计算公式为：级数的快速计算方法之一：格网分裂法*2020cossin(cos )( )( )nnnmnmnmnmnmnmnmNRCmDmPRDP 由格网可以看出,同纬度一行点的值相

11、同,同经度一列点的值也相同。(cos )nmP( )nmD18级数的快速计算级数的快速计算方法之一：格网分裂法经度的分裂将经度按等间距写为，是格网最左列的经度原点。由于：0q 0*00coscoscossinsinnmnmnmCmCmmqCmmq 00sinsincoscossinnmnmnmSmSmmqSmmq *000000coscossinsincossinnnmnmnmnmnmnmNRmqmCmSmqmSmCP max120( , )cos( , )sinNmNRQ mmqQ mmq 其中Q仅与纬度相关。方法评述方法的优点是先由径度原点计算每一纬度点上的系数Q，再代入公式计算大地水准面

12、高。这样对纬度而言，没有重复计算。利用递推性质计算三角函数。具体公式由读者自选导出。计算速度结果表明，利兹方法比常规方法快近100倍。可见利兹方法比常规方法具有更高的效率。但必须指出，利兹方法只适用于计算等间距地理网格点的重力量，若计算单点的重力量，则这两种算法的计算速度将不会有太大的差别。 19级数的快速计算基本原理将公式写成为如下形式级数的快速计算方法之二：矩阵求和法化简：将勒让德函数的（纵向）递推性质写为矩阵形式有maxmax*0cossinNNnmnmnmmn mNRm Cm SPmax0cossinNmmmNRAmBm*Nmnm nmn mACP*Nn nn mAC PTAC P12

13、0nn nn nPa Pb P0DPP1100()TTTAC D PPDC*Nmnm nmn mBSP*Nn nn mBS PTBS P矩阵形式20方法评述该方法比较稳定，减少了有效数字的损失。用该方法计算局部区域或全球网格点的重力量时，比常规方法快20倍。该方法的另一个优点是还可以计算重力场的其它量。 级数的快速计算级数的快速计算方法之二：矩阵求和法化简（续）：记，此时由有 0TAP Y1()TYDCTD YC1mmm mAPY同理得：11112nn nn nnYa Yb YC其中2mmm mBPY其中22212nn nn nnYa Yb YS1(21)(23)sin(1)(1)nnnan

14、mnm125 (1)(1)21 (2)(2)nnnmnmbnnmnm矩阵化为分量21级数的快速计算级数的快速计算方法之三：格网对称法原理与方法为了避免重复计算，节约内存，提高计算速度，我们利用缔合勒让德函数相对于赤道的奇偶对称性，结合地理格网的一些特点，导出了一种高效率的计算对称于赤道区域的重力量的方法，称之为“格网对称法”。下面，介绍该算法。由利兹法中交换求和顺序后得到的大地水准面高的计算公式max120( , )cos( , )sinNmNRQ mmQ mm公式的进一步简化关于三角函数的计算，按三角函数递推关系式对称性运用关于对称性的利用。考虑到勒让德函数相对于赤道的奇偶对称性，在南纬为负

15、的情况下，我们有 cos( +2) )2cos cos(1) ) cosmmm sin(2) )2cos sin(1) ) sinmmm(sin()( 1)(sin )n mnmnmPP 22级数的快速计算级数的快速计算方法之四：快速富氏变换法方法的原理将大地水准面高的计算公式写为方法评述利用快速富氏变换法的原理构造的计算大地水准面高的计算方法的计算速度，在计算全球的大地水准面高时，其计算速度比格网分裂法（Rizos）快近1.5倍，比矩阵求和（Clenshaw）法快7倍多。 1*20( ,)cossin(sin)knijnmjnmjnminmNRCmDmP 11*0(,)(sin)cos()s

16、in()(,)( 1)(sin)kkijnminmnmn mijmn mnmiNPRCmjDmjNP *1(sin)iknnmnmiin mnmnACPDB *1( 1)(sin)iknn mnmnmiin mnmnACPDB10(,)cossin(,)iikijnniiijmnnNABRmjmjNAB 232.地球重力模型的解算 24地球重力模型的解算可以证明以下公式称为展平因子且。不难证明 平均值及展平因子对重力异常的级数式在方形区域上积分（取平均）得重力异常的平均值为其中平均值所在区域的面积，也可能取为一个与之相似园。即*20(1)( , )( , )nnm nmnm nmnmgnC R

17、D S *20(1)( , )( , )nnmnmnmnmnmgnCRdDSd 02 21 120212100sinsin2 (1 cos)()(coscos)dd dd d 2000( , )2(cos )(,)nmnnmRdPR 2000( , )2(cos )(,)nmnnmSdPS 002(cos )sinnnqPd 10100(cos)(cos)(21)(1 cos)nnnPPqn*000020(1)(,)(,)nnnmnmnm nmnmgnqCRD S 园形区域与方形区域的比值25地球重力模型的解算积分的化简将球面重力异常数据的网格形式划分为若干小块，以C为例:地球模型的基本计算公

18、式由重力异常的球谐公式得系数计算式 *20(1)( , )( , )nnm nmnm nmnmgnC RD S 14(1)nmnmnmnmRCgdnSD11111( , )sin4(1)4(1)jiijjinmnmnmijiijijnCgRdgRd dnnq 11sinsin1(cos )sin4(1)iijjnmnmijiijnmmCgPdnmq 1sinsin14(1)nmijjnmijiijnmmCIPgnmq26地球重力模型的解算其中：拉让德函数积分的计算下述积分的计算方法很多。 运用的递推性质可以导得11(cos )sin( )iinmiitinmnmtIPPdPt dt 122,2

19、,2( , )1( , )1(1, )1jnmjtnmnmtnA n mtIPIPA n mPnA nmn12,2,( ) (1)11jnnjtnnn ntnB n B ntIPIPPnn00jIPt 221013()2jjIPtt121113(1cos)2jjttIPttt27地球重力模型的解算若干说明地球重力模型的计算除了工作量大外，还有许多问题需要考虑。如零阶项由于计算公式是从2阶项开始的，必须对重力异常及其相关资料进行必要的约束，以满足“质量相等、质心重合”要求。椭球改正由于重力异常观测值进行重力数据处理后是归算到大地水准面上，而地球模型的计算中重力异常资料是在积分曲面（球面）上的，二

20、者存在差异；因此要考虑该差异的影响，该影响称之为椭球改正。大气改正在正常重力与位系数计算公式中使用的参数一般是采用国际组织提供的参数，而IUGG提供的GRS80等系统中均包含有大气质量，为求得地球表面上可用的地球重力模型，则该部份质量必须去掉。多种资料的联合处理由于地球表面的复杂性，它表面上的测量方法与手段各异。为保证原始资料的精确度，一般不要求事先将它们转换成重力异常。例如，在陆地上重力资料较为丰富，这其中随着地区不同而资料的精度也不一致，在广阔的沙漠、原始森林、崇山峻岭等地区重力的测量也是极为困难，因此该地区的资料可能是垂线偏差或者是高程资料。而海洋区的重力测量数据寥寥无几。随着卫星技术的

21、发展，现今海洋上均以卫星测高资料已有较大的发展。另外，卫星重力（梯度）、GPS、传统的水准资料、VLBI（甚长基线）、激光测星（月）等等，它们都将为地球重力模型的确定能提供必要的信息。不能忘记的是，先前的任何地球重力模型也是好的基础背资料。将这许多资料进行联合处理必将得到好的结果，但处理方法的选择必是一项十分重要的工作。283.重力学中积分的计算方法 29在地球重力学中，积分的应用十分广泛，现简要介绍几种常见的积分应用。三维体积分的应用三重体积分主要用于计算地球的引力位及其相关量。质体的引力位如计算均质地球的重力场分布，计算层层质量分布的地球引力场分布；也可用于计算旋转时的地球重力场分布；重力

22、归算如层间改正（即无限平面层的引力），地形改正（局部地形起伏的引力），重力均衡（计算柱体对柱体顶部或柱体底部的引力）；重力正反演在重力正反演的计算中，正演的计算在一些规则的物质形状下是可以直接进行积分的，而对一般形状下大多数积分也只能采用数值积分的办法（见下节）；在反演的计算中，能直接求解的是少之又少，如均质球等。面积分的应用这时的面积分包括二维平面积分和球面积分甚至一般的曲面积分。二维平面积分在地球引力位的计算中常常被使用，主要原因有三：一是当质量层厚度远小于平面层大小时，质量层常常被当成平面处理；二是研究地球的断层面或物质不连续面时；三是由于位理论中的等效原理的原因。层面引力位正演计算中，

23、用于计算厚度较小的物质层产生的引力位，还可用于研究断层（密度分界面）产生的引力位；单层位质面引力位常用的有单层引力位、双层引力位、多极子引力位等。多用于理论证明与位的高阶梯度的计算；大地水准面计算计算大地水准面高的司托克斯公式是一个球面积分，类似地还有一些相关量也可以用司托克斯公式的“变换”得到；地球重力模型（位系数）的计算理论上它也是一个球面积分。 积分在重力学中的应用30积分在重力学中的应用重力学中积分的计算方法重力学中的一些积分计算都比较困难，其主要原因是积分都是高维积分。本节介绍积分在重力学中的计算方法，如直接积分法（投影法与截面法）、变量替换法（柱坐标、球坐标）方法等。 坐标选择法其

24、基本思想是选择坐标系使得积分简单。如取坐标原点位于质心或位于计算点，也可取计算点位于某一坐标轴上。例如均质（密度为常数）球体的引力位的计算。在球坐标下，选球心为坐标原点，且使计算点在轴上，则有 如果计算点不在轴上，则积分的计算就比较麻烦223000111110441sin()3RGGGMVGrdr d dr drRrrlrrr 0sindl10112,(sin2,(rrrdlrrr当时)当时)31积分在重力学中的应用等效原理方法等效原理是物理学中一些重要的物理结果，也可能是一些经过数学论证成立的结论，在地球重力学的计算中常常被使用；如二度体问题的正演计算等。现以一个简单的例子来说明。 均质球体

25、外部的重力引力分布设均匀球体的密度为，半径为。利用等效原理：质球体的产生的重力可等效为一等质量的质点产生的的重力。 建立以球心为原点、计算点在轴上的坐标系。由于球体的质量故它产生的引力位为与一质点产生的引力位相等，即这里G为万有引力常数，r=z为计算点至球心的距离。不难看出，等效原理使得可直接写出引力位公式，计算点置于坐标轴上使得计算更加简单。 343vMd d dR GMVr32积分在重力学中的应用直接积分法直接积分法是对积分中一部分容易积分的部份进行先行积分（也称为逐步积分法）。积分计算的技巧很多，如对称性、无关性的利用等等。不能一一叙述，仅举两个例子。圆柱的引力求密度为、半径为R、高为h

26、的圆柱过圆心垂线上圆柱外的引力。取坐标原点于园柱顶端中心。设待求点距圆心距离为a，由对称性可知： 。而采用柱坐标得 若半径R时为层间改正，即得0 xyFF222223220 0 02()()hRzGzaFrdrdGhRaRh arza 2zFG h 33积分在重力学中的应用位理论方法在重力学中有许多积分难以计算，特别是当积分区域较为复杂时尤为困难。例如，若地球的密度为常数，问地球外空间一点处的引力为多大？当考虑地球的复杂形状时，引力位的积分是难以精确计算的。若存在一种地球内部的密度分布使得（地）球表面的引力位为常数，问（地）球外空间任一点处的引力位为多大？ 积分（方程）反演方法的讨论对均值问题

27、，反演（反求质体内部的场源密度）有时是可能简单实现的。这里均值是指密度为常数的情况。均值问题是重力正反演中的计算最简单的问题。在大多数情况下，积分是可进行的（但并不是说反演结论是绝对正确的，问题的非唯一性可能还存在）。 化积分为微分有时还利用位的微分性质，将积分化为微分进行讨论。常见的讨论方法有，微分方程的差分方法、微分方程的有限元-边界元方法等。 34积分在重力学中的应用积分的积分变换法 质体的引力位由引力位的计算公式（直角坐标形式）为求V的逆变换 222( , , )( , , )()()()d d dV x y zGxyz 求R的逆变换后可得问题的解 作积分变换 222( , , )(

28、, , )*GV x y zx y zxyz( , , ),( , , )( , , )0,( , , )x y zx y zx y zx y z222( , , )GR x y zxyz( , , )( , , ) ( , , )V u v wu v w R u v w44( , , )VGGm x y z 2 2224 () ( , , )4uvw VuvwG m 222( , , )( , , )()G uv wV uv wuvw222( , , )()GR u v wuvw1222( , , )( , , )()Gu v wV x y zFuvw35积分在重力学中的应用积分的积分变换法

29、 大地水准面高大地水准面高的计算公式是一个曲（球）面积分。先考虑平面形式的计算公式。为方便使用快速积分变换，计算大地水准面高的司托克斯公式这是二维平面形式的计算公式。可得问题的解 司托克斯公式的平面近似公式为：( )4RNg Sd2( )(,)PPRSS x x y yl2R ddxdy1( ,)( ,)( ,)4PPPPPPN x yg x yS x yR 11( ,)( ,)( ,)4PPPPPPNx yFF g x yF S x yR 36积分在重力学中的应用积分的数值积分方法网格分区方法对于区域大的积分，特别是一些远处信息对计算点值影响较小的积分，往往采用分区积分的办法。具体办法是：远

30、区采用近似积分的办法，如采用经验公式的办法、舍去的办法、模型的办法等等；近区采用实际资料的格网数值方法。 最后一式可称为半数值方法（即部分网格法部分积分法）。球面积分的处理方法地球形状研究中常常用到球面积分20 0( , )( ,)( , , )cosQkgKdd 001 1( , )( ,)( , , )cosJIijijiijjiQkgK 110001 1( , )( ,)( , , )cosjijiJIijjiQkgKd d 37积分在重力学中的应用积分的数值积分方法三维积分的处理方法在重力正反演的计算中，若质体庞大并且物质密度变化不剧烈时，其远区域的影响可用常密度处理，近区域采用网格积

31、分方法，即精细分区由于计算点附近的重力信息（质量影响）对计算点的贡献远远大于远处的重力信息的作用，因此计算点附近应该具有更加细致精确的重力信息。在上述计算中，一般不统一采用同一规格大小的观测资料，而是选取由远而近越来越精细的网格划分，即把近区再划分为多级分区；如 区、区、区等等。注意：要注意球坐标下球（园）形区域与直角坐标下网格区域的链接问题，尽量避免区域的重复或遗漏。01( , )iiiiQGdGdGlll 001120151 138积分在重力学中的应用积分的数值积分方法奇异积分的处理方法一般而言，积分有确切的计算值。但在数值计算时，函数可能变为奇异而无法求出，这种情况往往是不真实的。 球面

32、积分法对奇异点的邻域，采用球面极坐标，若则先作局部积分即可消去奇异现象。三维积分中奇异积分的处理方法与此相类似。 局部拟合法对奇异积分的问题还可以应用观测资料的局部拟合方法处理。其基本思想是将计算点附近几个（四个或九个）网格（块）上的积分上的观测值拟合成一个多项式代入积分中再进行积分即可。1( , , )( , )KK 394.线性重力问题的求解方法 40重力反问题 重力学中的反（演）问题在地球重力学中，反（演）问题很多。如已知地表观测值求场源（质量）分布、推求断层，利用地表重力异常寻求地球内部的剩余密度等等。 一般地，场源函数与地表物理场之间存在着积分关系（也可以表示为其它关系）为其中，为场

33、源的作用区域，称G为积分核函数（又称格林函数），它刻划了“场”g与“源”m间的“关系”。 ( )( , ) ( )Qg PG P Q m Q d 若场函数g为（地表观测的）已知函数，积分核函数G根据场与场源的物理实际给出，求解场源的物理问题称为反演问题。反（演）问题就是求解上述积分方程问题。 有时为叙述方便，将它写成为内积形式或算子形式。若存在“逆”算子使得m=G-1g，则称m为算子方程的（反演问题意义下的、广义的）解。一般情况下，积分算子与微分算子互为逆算子，如牛顿引力位的积分与引力位的调和性质。 ( )( , ),( )g PG P Q m QgGm41适定性与最优解反演问题的适定性对一般

34、的反演问题而言，其解是不适定的。例如：存在着“非0”的场源m，使得Gm=0，称m为“零解”，亦称“0场源”，即无法观测的场源。此时，若m1是g的场源，则m1+m也是的场源，即场源不唯一。 反演问题不适定的另一个特点是解的不稳定性，若观测值增加充分小的，积分的解是的函数，若解随减小而增大则称为不稳定（请参阅相应的参考书）。 反演解的不唯一性极大地阻碍了反演问题的求解，也阻碍了地球重力学的发展。对此，如同数学中求解解不唯一的问题一样，增加约束条件使得问题具有唯一解。增加的条件必须符合物理实际，又要使得问题的求解成为可能（解唯一、方便解算）。由于所增加条件的不同，产生了多种反演问题的解决办法。在一般

35、情况下可表述为： 在条件Gm=g下，满足 的解称为“最优解”。min ( )f mmin ( )f m反演问题解的存在性问题更是经常遇到，如当观测数据多于待球求量时问题往往无解。如：线性代数中方程个数多于未知量的个数时方程无解。此时，人们往往寻求一个查对而言“比较好的解”。在一般情况下，约束条件可表示为“差异最小”，即这里，称m*为问题的最佳解，称为范数，意义为“距离” 。*mingGm42问题的线性化方法几点说明从推导过程可以看出： 解是不唯一的。即有限个观测值不可能确定无限个系数（精确模型）。对近似模型，不同的观测值确定不同的近似模型。 对于任一组有限（个）观测值，总可以找到一个解使得满足

36、。 随着观测值个数的不同或观测数据间的相关性不同，方程的求解方法很多，一般寻求一个“最佳最优解”，即最好的解。线性化方法设为一组归一化的正交多项式系（基函数），则一维问题变为 系数待定。上式代入g=Gm中得系数gk为已知且未知量的系数为线性的。 1( )( )kkkm xmx1( )( )kkkg xm gx上式中无限多个系数需要无限多个观测值去求定，这在实用中是不可能的。在实际应用中，总是用有限展开方法作的近似，即设则以上推导仍成立。再用足够多的观测值可求得所有的系数。 1( )( )Nkkkm xmx43适定问题的最优解适定问题的唯一解若使用N个无关的观测值去求定N个系数m，则方程写为：由

37、观测值的无关性可得矩阵G可逆，即问题有唯一解。求解该线性方程的方法很多，如消去法、迭代方法、共轭梯度法等等。 gGm44超定问题的最优解超定问题的最佳解若使用多于N个（如NP个）无关的观测值去求定N个系数，则该问题（在线性代数中）无解。该问题在数学上称为超定问题，可以寻求一个与问题比较接近的“解”，该“解”称为超定解，也称为最佳解。 何为“接近”？一般的定义是“距离”最小。设为超定解，即式中为范数（距离），常用的是L2范数： minegGm221Niiff在该范数的意义下，可以求得Gm=g的解：事实上，由与即可求得上式。该方法一般称最小二乘方法（欧氏几何意义下的距离最小），其解也称为最小二乘解

38、。 1()TTmG GG g() ()TegGmgGm()0TeGgGmm 45欠定问题的最优解欠定问题的最优解欠定问题是当线性方程的个数少于未知数的个数时，方程变为有很多解。在很多解中寻求一个符合某种特性的解的方法是最优化方法，求得的解称为最优解。“符合某种特性”在一般情况下往往选取“最小”（距0最小）性质，即要求 mindm在L2范数的意义下，问题变为：。这是一个二次规划问题。可用最优化方法求解。 min ( )st.Tf mm mGmg引入拉格朗日乘子，将上式写为 min ( , )()TTf mm mg Gm使用乘子法求得：2m=TG与Gm=g，将前式中的m代入后式可得 ，其中系数矩阵

39、为可逆方阵，故可求得，再代回前式求得解，即 2TGGg12()TGGg1()TTmGGGg46混定问题的最优解混定问题超定问题条件过多，欠定问题条件过少；它们的求解都比适定问题困难。无论是哪一类问题都存在求逆矩阵的问题，当矩阵行列式的值较小（接近于零）时，方程的求解是十分困难的（求解不稳定）。究其原因是相关性较强，即尽管约束（条件）十分多，但约束之间的独立性较差，形成了形式上为超定实际上为欠定的局面，俗称为混定问题。一般解法是求2min( )() ()TTmgGmgGmm m即得。选择阻尼（或加权）因子使得方程有解 21()TTmG GIG g2()TTG GI mG g2当因子充分大时，问题

40、以欠定部份为主，反之以超定部份为主。增加阻尼因子的方法也称为正则化方法。 247线性问题的最优解线性问题超定问题条件过多，欠定问题条件过少；它们的求解都比适定问题困难。无论是哪一类问题都存在求逆矩阵的问题，当矩阵行列式的值较小（接近于零）时，方程的求解是十分困难的（求解不稳定）。现将它们的解汇总如下：超定问题的解 混定问题的解 21()TTmG GIG g适定问题的解 欠定问题的解 1()TTmG GG g1()TTmGGGg1mGg485.重力场的逼近方法 49司托克斯逼近方法 50司托克斯逼近方法大地水准面的计算方法介绍大地水准面是物理大地测量学乃至其它许许多多科学中的一个非常重要的基础性

41、概念。计算大地水准面高常用的有两种方法，即（1）司托克斯积分公式方法，或称为重力异常方法（2）球谐级数方法，或称为位系数方法这两种方法各有特色：司托克斯公式方法公式简单、计算方便，它的弱点是需要全球的重力异常资料；这一点是难以保证的，因为海洋上的重力测量十分困难。而球谐级数方法中一旦有了位系数、之后其计算也十分易于实现；但该方法的不足之处是，位系数总是有限的，不可能无穷下去；并且当较大时，位系数的计算也十分困难；况且球谐函数的计算中收敛性、稳定性也值得深入研究。重力场的逼近方法就是寻求一种顾及以上两种方法中各自优点的方法。以下就分别讨论这种方法。 dgSRN420,nnmnmnmnmnmSBR

42、ARN51司托克斯逼近方法司托克斯公式的逼近方法区域的划分 在司托克斯公式中，积分区域是对全球进行积分，由于全球资料的不足而设法改为对局部区域（地区）的积分，而对计算点影响较小的远区域则尽量使用其它的资料。为此，我们将地球按计算点影响划分在远区域和近区域两部分（如图），那么，积分可以写为 0144RRNSgdSgd远区函数的展开公式 0000SS 021(cos )2n nnnSQ P00( )(cos )sin( )(cos )sinnnnQSPdSPd 52司托克斯逼近方法远区的公式大地水准面高的远区影响为 020021sin(cos )442nnnRRnNSgd dQ Pgd dPgng

43、nn)(cos412NNgQRNnnn100)(2 20)(240nnngQRdgSRN大地水准面高的公式大地水准面高等于远区域与近区影响的和大地水准面高可以由近区域的重力异常积分（第一积分项）和远区域的位系数求和（第二求和项）相加求得。从而避免了远区域重力异常资料的不足的影响。 53司托克斯逼近方法实用计算公式 位系数总是有限的，常常取为有限值。可以证明：舍去的部份（称为截断误差）是一个小量。关于积分的计算下节讨论。 NgQRdgSRNNnnnmax020)(2400( )(cos )sin( )(cos )sinnnnQSPdSPd *0(1),nnnmnmnmnmmgnC RD S 54

44、一般的逼近方法 55一般的逼近方法公式的改化 引入一个已知函数，将司托克斯公式改写为区域的划分 将积分区域按为界分为远区域和近区域，则上式可改写为)(S dgSRdgSSRN)(4)(4 dgSSRdgSRdgSSRN10)()(4)(4)(4210NNNN远区的作用 56一般的逼近方法远区的作用远区的处理I在中取容易得到可得远区域的影响为远区的处理II 由展开式dgSRN)(410)(cos212)(nnnPQnS212nnngQRN0sin)(cos)(dPSQnn00000(0)21( )() (cos )( )( )()2nnnnSQPSS 0sin)()(sin)()(00dPSSd

45、PSQnnn212nnngQRN202)(2nnngQRN20021)()(2nnnngQQRNNN57一般的逼近方法函数的讨论在一般情况下，函数取为分段函数的形式此时的大地水准面高为这里的为待定函数，常常选为)(S*00( )(0)( )0()SS*( )S)()(cos)(212)(00*mmnnnSPQnSNNNN100)()(4*0dgSSRNnmnmknkknngXQkQQRN200*00*1)(212)()(210*00)(212)(2mnnnkkmkngXQkQRN58一般的逼近方法几个特殊的逼近方法司托克斯逼近方法如前面的介绍，司托克斯逼近方法的公式为并且函数为连续函数的逼近方

46、法在司托克斯逼近方法中，函数为间断函数，这与司托克斯函数的性质严重不符；为改善函数在处连续状况，取。那么 0000SS20)(2nnngQRN)(S)()(0*SS20200021)(2)()(2nnInnnnngQRgXSQRNN)()(cos)(11)(1000nnnInPPSnQQ2021)(2nnIngQRNN59一般的逼近方法几个特殊的逼近方法函数为光滑函数的逼近方法取，；则该方法也称为带边界约束条件的逼近方法。)()(0*SS)()(00*SS)(sincoscos)()()(00001*SSSS10000000sin)()(sincos)(nnnIInXSXSSQQ2021)(2

47、nnIIngQRNN60一般的逼近方法几个特殊的逼近方法有限展开的逼近方法该方法是取司托克斯级数的有限项构成函数，即并且mnnmPnnSS20*)(112)()(mknknIIInXkkQQ211220221)(2222nnIIInmnngQRgnRNN0)()(400dgSSRNmnmnnmkngXkkQnRN20201112)(1221020112)(2mnnnmkngXkkQRN61一般的逼近方法几个特殊的逼近方法最优化的逼近方法选取函数使得它为特别需要强调的是，这样选取的函数是在的意义下的最优逼近。并且mnnnmPQnSS00*)()(212)()(0minsin)()(2*dSS102021)(2)(2mnnIVnmnnngQRgQRNN计算公式为0)()(40dgSSRNmnmnngQRN21)(210200)(212)(2mnnnmknngXQkQRN62谢谢！