本科毕业论文-—近震定位及其应用.doc

桂林理工大学本科毕业设计·论文

目 次

摘要………………………………………………………………………………………Ⅰ
Abstract………………………………………………………………………………… Ⅱ
1 前言 ………………………………………………………………………………… 1
2 近震的理论基础………………………………………………………………………2
2.1 地震学的基本名词和概念……………………………………………………2
2.2 地震的分类……………………………………………………………………2
2.3 近震的主要震相………………………………………………………………3
3 正演计算（模型试算）…………………………………………………………… 5
3.1 基本原理………………………………………………………………………5
3.1.1 坐标变换………………………………………………………………5
3.1.2 正演计算………………………………………………………………6
3.2 正演模型的建立………………………………………………………………7
4 反演计算………………………………………………………………………………8
4.1 基本原理………………………………………………………………………8
4.1.1 当速度V未知的初定方法……………………………………………8
4.1.2 当速度V已知的初定方法……………………………………………10
4.2 模型反演结果…………………………………………………………………11
4.3 误差分析………………………………………………………………………11
5 实例剖析………………………………………………………………………………13
6 结束语…………………………………………………………………………………17

参考文献………………………………………………………………………………… 18
附录 近震定位程序………………………………………………………………………19

26

1 前 言
四川汶川地震和青海玉树地震发生时的地动山摇、房倒屋塌、生死分离的那种凄凉、悲惨的场面，又一次的让我们感受到自然灾害的威力，更是让灾区的人们心惊肉跳，恐慌难安，甚至谈“地震”色变，这就是地震给人们的最直观的印象。一次强烈地震的发生，常伴随着地面变形和地层错动，其破坏力是相当大的。主要表现为：大型建筑物破坏，普通民房破坏，山崩地裂，人畜伤亡。如今地震发生越来越频繁，2010年也被“尊称”为国际地震年。任何事情都具有双面性，地震对我们的生命安全造成了很大的威胁的同时也为专家学者们的研究地球提供了丰富的资料。目前人类对地球内部结构的认识主要是来自地震学研究，因此地震学也就成为地球科学中的一个重要学科。
由于天然地震具有很大的能量，它所产生的地震波可以穿透很大的深度，传播很远的距离，而且当遇到地球深部的波阻抗差异界面时产生反射波，其中，天然地震的中的近震资料为研究地球深部结构提供了一个重要的信息资源和途径。
利用地震记录进行定位始源于欧洲和日本。最初使用方位角方法，随后是几何作图法和地球投影法。《国际地震汇编》（ISS）最先采用最小二乘法进行计算修订震中。1961年，博尔特和威尔莫合作，改进了计算方案，并首先在ISS使用，随后，国际地震中央局与美国海岸和大地测量局先后使用。在我国，李善邦先生最初使用方位角和最小二乘进行地震的观测和定位。1953年，我国开始采用大量观测数据修订震中。目前我国的地震的定位方法兼作图法和计算法。其中近震的定位方法主要有石川法、和达法、高桥法、外心方位角法、假定发震时刻定位、等时量板法等十几种方法。计算机的近震定位方法主要分为速度未知的初定为方法和速度已知的初定为方法。
目前定位精度已达到较高的水平。如果有合理的台网分布和适用的走时表，对于Δ<100km的区域性事件，定位精度为1km左右；对于100-1000km的近震，可达2至5km；远震为5—10km。一般浅源地震（h<100km）的震源深度误差为深度值的10%左右。震源愈深，相对误差愈小；h>300km时，误差小于其深度的5%。发震时刻的误差为1/10—1/2s.

2 近震的理论基础
地震学主要研究地震的发生、地震波的传播及地球内部构造，地震波能穿透到地球内部，并能把地球内部的信息带回地面，具体说来，它主要是根据天然地震或人工地震的资料，运用物理学、数学及地质学的知识，来研究地震发生的状况及地震波传播的规律，以求进一步达到研究地震和预报地震的目的。同时地震是地球表层的振动，是地壳构造运动的一种形式，所以还利用地震波的传播特征来研究地壳和地球内部的构造。

图2.1 地震波的分类
在弹性介质内传播的纵波和横波，由于存在整个弹性空间，因而这些波统称为体波。相对于体波而言，在弹性界面附近还存在着另一类波动，从能量上而说他们只分布在弹性分界面附近，称为面波。面波又分为瑞雷波和勒夫波(如右图2.1所示)。为便于对地震的了解研究，现将其基本知识介绍如下。
2.1 地震学的基本名词和概念
1．震源：地球内部发生地震的地方称为震源（或称震源区）。理论上将震源看成一个点,而实际上是一个区。
2．震源深度：将震源看作一个点，此点到地面的垂直距离称为震源深度，一般用字母h表示。
3. 震中：震源在地面上的投影点称为震中（或称震中区）。同时，地面上受破坏最严重的地区叫极震区，理论上震中区和极震区是相同的，实际上由于地表局部对质条件的影响，极震区不一定是震中区。与震中相对的地球直径的另一端称为对震中；或称震中对点。
4．震中距离：在地面上，从震中到任一点沿大圆弧测量的距离称为震中距离。一般用字母⊿表示。它可用线距离表示，也可用地心所张的角来表示。
5．发震时刻：发生地震的时刻，一般用字母Ο或T0表示。以北京时间标出。它比格林威治时间早8小时。
2.2 地震的分类
为研究方便，按震动的性质，可分为天然地震、人工地震及脉动三类。对于天然地震，有下述分类：
（一）按震源深度分：
1．浅源地震：震源深度小于60公里的天然地震称为浅震；也称正常深度地震。
2．中源地震：震源深度在60公里至300公里之间的地震称为中源地震。
3．深源地震：震源深度大于300公里的地震称为深震。已记录到的最深地震
的震源深度约700公里。有时也将中源地震和深源地震统称为深震。
（二）按震中距分：
1. 地方震：震中距小于100公里的地震。
2. 近震：震中距小于1000公里的地震。
3. 远震：震中距大于1000公里的地震。
2.3 近震的主要震相
由于近震的震中距很小，其地震波的传播主要局限于地壳及莫霍面附近，主要震相有Pg/Sg（）、P*/S*、Pn/Sn、P*P/S*S、P*S/S*P、PmP/SmS及PmS/SmP等，如图2.2所示。

图2.2 近震震相及其射线路径示意图
（αc、βc分别为康拉德界面及莫霍面的临界角）

1. Pg/Sg震相
由于地球浅部存在强烈的横向及纵向上的速度变化，尤其是有沉积该层地区的垂向梯度的影响，地震波射线在地壳浅部10公里的范围内发生廻折而到达地表，并在震中附近数公里的范围，这种波束与直达波，被称之为廻折或潜波（Diving Wave），对于纵、横波来说，分别命名为Pg和Sg。
2. P*/S*震相
P*/S*分别为地震波在地壳内部的康拉德界面产生的折射纵、横波，它们的速度分别为6.3~7.1km/s及3.6~3.9km/s。应当指出，由于康拉德界面并不是全球普遍存在的界面，在一些地区观测不到。
3. Pn/Sn震相
Pn/Sn分别为地震波在地壳底部的莫霍界面产生的折射纵、横波，它们的速度分别为7.9~8.2km/s及4.4~4.8km/s。一般来说，Pn震相在莫霍面的临界反射以远地区总是已初至波的形式出现，与其它折射波一样，波至呈现为线性分布。
4. P*P/S*S及P*S/S*P震相
P*P/S*S分别为地震波在地壳内部的康拉德界面产生的反射纵、横波，P*S/S*P分别为地震纵/横波入射地壳内部的康拉德界面后发生转换而产生的反射横/纵波。

3 正演计算（模型试算）
测定震源位置及发震时刻最常用的方法是将非线性问题化为线性问题，然后利用最小二乘法原理迭代，全主元消元法等求解。通过反复修定，以得到震源参数的最佳值。此方法对近震经纬度的确定比较准确的。
近震的定位要确定四个参数，即震源的空间坐标（，，h）及发震时刻t。为此，我们需要知道台网各个台站的坐标，以及至少五个台的地震波到时。由台网提供的到时数据主要是指P 波的到时，因为各种续至波的相位通常是难以识别的。如果地震发生在台网内，那么只利用P波到时定位并不困难。如果地震发生在台网以外，这种答案是不确定的。为了进行参数的求解，一般要对地壳平均速度做出假设。
为了对此问题进行研究，有必要建立近震定位的正反演模型。对于正演模型，其
基本思路是通过设定震源初始参数（经度、纬度、高程、发震时刻及均匀介质中的波速）计算出地震波从震源到各台站的走时，进而求出各台站接收到地震波的时刻；这些数据就作为反演的数据来源。此模型的建立主要涉及两个关键问题：坐标的变换、走时计算。
3.1 基本原理
3.1.1 坐标变换
对于近震及地方震的计算，采用平面坐标比球面坐标方便。为此需要将经纬度换算成局部的直角坐标系。这里需要指出，为了处理数据方便，一般是将国内基准台网及区域台网的全部台站的经纬度及高程均输入计算机内。但是，在作具体计算时，只用到其中一部分台站，也就是只需要把这部分台站的经纬度（球面坐标）化为平面坐标。
因而，直角坐标系的原点是浮动的，每次都选来尽可能的靠近震中区，以减少变换过程中带来的误差。在计算程序中，坐标原点确定在接收到该次地震信号各台站经纬度的平均值0及0上，x轴向东，y轴向北，如图3.1所示。

图3.1 坐标转换示意图
于是，地面上任一点（其经纬度为、）相对于坐标原点O（0,0）的直角坐标、的换算公式为：
（3.1）
其中参数:

而参数和中的常数又为：
==0.0066935165 (地球偏心率)
=6378.160 km (地球长半径)
= 6356.778km (地球短半径)
经换算为平面坐标系后，就可以用计算机确定震源参数。但这样得到的震中位置仍是局部的、坐标。因此，必须化为经纬度，输出。将直角坐标值，化为经纬度，的换算公式为：
(3.2)
3.1.2 正演计算
如果已知震源及台站坐标；设震源坐标,发生时刻，均匀介质中波速；有n个台站，其中任一台站坐标，相应直达波到时如图3.2所示。

图3.2 正演计算示意图
由物理知识可知： ，然而实际中，台站及震中坐标是以经、纬度（地理坐标）表示的，在计算前首先要把地理坐标（经、纬度）转换成大地坐标（直角坐标）。
3.2 正演模型建立
假定某地震事件发生后被8个台站检测到(图3.3) ，这8个台站的经纬度分别是：
（96.000000， 28.000000）、（95.000000， 29.000000）、（97.000000， 29.000000）（94.000000， 30.000000）、（98.000000， 30.000000）、（95.000000， 31.000000）
（97.000000， 31.000000）、（96.000000， 32.000000）
通过给定两个模型：
Mod1的震中位置为：（96.000000， 30.000000）
Mod2的震中位置为：（95.500000， 29.500000），见图3.3

图3.3 台站分布图
正演的结果：
Mod1：t= 44.386244、29.397192、29.397192、38.647902、38.647902、29.523923
29.523923、44.386244
Mod2：t= 34.676012、14.846849、31.128296、31.261550、49.893035、34.715948、
44.367994、56.320175
对于mod1而言，震中在其中心，它与（96.000000，28.000000）台站相距20，有换算关系10 111.19，且速度v=5.0可知此正演结果与实际情况相符合。Modd2也是如此。
4 反演计算
对于反演模型，主要目的就是利用地震波在各接收台站的到时，结合台站经纬度，通过数学推导而确定震源参数。实际计算表明，只要输入数据正确，初定的震源位置及发震时刻足够精确，就可供地震速报之用。其具体方法如下。在此主要利用前面正演模型所得数据，作为反演之用，进而推导出震源参数。
4.1 基本原理
4.1.1 波速V未知的初定为方法
此方法通常用在台数n>4，震中距小于170km，已知直达纵波P到时，而且台站
方位和震中距分布合理的情况下。
设接收到某次地震直达波到时的台站有n个（n≥4）。其中任一台站的坐标为（xi,yi）相应的直达波到时为Ti。待求的震中位置为x,y,震源深度z,发震时刻为T（图4.1）；按照均匀地壳模型，可列出下列方程：

图4.1 震源计算机定位原理图
（4.1）
i=1,2,…，n.
式中为地震波在地壳中的平均传播速度，因区域不同而异。
式（4.1）为非线性方程组，为便于直接求解，首先应将其线性化，变为线性方程组。为此将（4.1）式展开得
(4.2)
i=1,2,… ,n.
对于第一个台，（i=1）可以写出
(4.3)
将(4.3)式中i分别代以2，3，…，n与(4.3)式相减得
（4.4）
共有n-1个线性方程式。(4.4)式以消去了深度参数z，只有x、y、T、V四个未知数。由于n-1一般都大于4，故(4.4)式所表示的线性方程个数多于未知数个数，即构成所谓的超定方程组

…… …… ……
（4.5）
为了表达简洁起见，将此超定方程组用矩阵形式表达
(4.6)
其中，对于n个台站的数据,

其中U=V2， W=V2T。
当地震台站数目n小于5时，上述方程组为不定方程，有多解；大于5时，方程组为超定方程，通过变换
（A’为A的转置）
使得超定方程组变为正定方程组，从而计算出震源参数，这就是最小二乘法原理。
只要解出此方程组(4.6)，相应震源参数也可求得。即：
震中位置: （x，y）
介质的波速: V＝
发震时刻: T=W/U.
又由：
任取一台站(i)数据代入此式便得震源深度：

4.1.2 当波速V已知的初定为方法
因为所研究区域位于青藏高原，其地壳的厚度较大，所以近震的深度基本在上地
壳的范围之内，其速度大约为v=5.8km/s.当速度为已知值后，在上述方法中的线性方程组就会减少一个未知数，所以（3.5）变化之后 ，
（4.7）
其中

W=V2T。
只要解出此方程组(4.7)，相应震源参数也可求得。即：
震中位置: （x，y）
发震时刻: T=W/ V2.
同理，再由求出震源的深度。
这两种方法的主要流程为：

图4.2 震源计算机定位流程图

4.2 模型反演结果
为了更好的检验程序，在正演t的结果上均加了6s，下面为各个反演方法所得到的结果：
1.当介质的速度V未知时，两个模型反演出的地震参数分别为：
mod1:震中(96.000000,29.999999) t0=6.000025 v=5.000002 h=9.998836
mod2:震中(95.490894,29.495549) t0=5.797093 v=4.994111 h=20.333211
2.当介质的速度V作为已知条件（V=5.0）时，两个模型反演出的地震参数分别为：
mod1:震中(96.000000,29.999999) t0=6.000000 h=9.999972
mod2:震中(95.490701,29.495359) t0=5.874814 h=18.794602
4.3 误差分析
从数据的反演来看，深度的误差很大，这主要是因为震源深度对地震波到时的读取误差反应极其敏感，是地震学中最难准确测定的参数之一，各种方法对于震源深度的估计都具相当程度的不确定性，影响着人们对震源过程的认识。具体来说，震源深度的精度误差受震中距、到时误差和速度模型（地壳模型）三个因素制约，而且这些因素对震源深度的影响是非线性的。
1.当地震波传播速度一定时，震源深度的误差随着震中距或台站位置的增大和走时残差的增大而增大。
2.走时残差一定时，震源深度误差随着震中距的增大和地震波速 度的增大而增大。
3.当速度已知，走时残差一定时，越浅的地震，定位误差可能越大。这一结论与许多深源地震的定位结果相一致。
尽管震源深度是地震学中很难精确测定的参数之一，各种方法测定震源深度的结果不同，但可以将不同的测定震源深度的结果进行对比分析，可能会改善对震源深度的测定精度。地震发生后，立即在震源区布设流动观测地震台站（网）是修正主震震源深度精度的有效方法。

5 实例剖析
所研究区域位于青藏高原南迦巴瓦地段，该处台站的分布如图5.1。

5.1 区域台站及示例震源分布图

1．台网记录的地震事件里该近震记录
根据美国国家高级地震系统（ANSS）提供的全球地震事件目录，可知该事件的位置（96.6050,29.5730），发震时刻为2003年第230天的09:03,该地震事件总有38道记录，震中距范围在0.50-30。经过上述方法反演后得到的位置为：（95.6102,29.4248），反演的结果与实际经纬度分别相差（0.0052,0.1382）。下图分别是利用geotool软件在选定地震事件前后，通过屏幕抓图得到的图像。其中，红色的p标注的时刻是geotool软件算出来的该地震记录道的起跳时刻，是理论到时。黑色p标注处的时刻是经人工判断后标记的该地震道记录的起跳时刻，是实际到时。
图5.2中，a为定位前的地震记录图像，b为定位后的图像。由图可看出反演后的结果。反演后理论到时和实际的起跳时刻仍然存在几秒的误差，因为模型毕竟是理想状态下的情况，而实际的资料并不是那么完美，在实际操作中很难得到极其精确的台站到时，一方面存在人为因素，如时间拾取的精度；另一方面也存在机器设备的物理因素。而且由于地表的覆盖层厚度，松散度等等的不同等等都会引起波至到时的误差，所以利用地震资料到时时，还要经过判断和人工调整。在该地震事件中，各个记录道起跳刻的最小误差为2s，最大的误差为5s。
桂林理工大学本科毕业设计·论文

a： 近震定位前 b： 近震定位后
图 5.2 实例一近震定位前（a）和定位后（b）
2．台站网记录的地震事件里没有近震记录时
当台站网记录的地震事件里没有近震记录时，利用上述反演方法进行地震的定位，可确定该地震事件发生的经度和纬度为：
（96.7124,25.7474）。将得出其经度、纬度、发震时刻、深度等信息，按照标准格式写入origin事件，利用geotool软件对该事件就行处理后，结果如下图5.3所示，其中最小的误差为0. 1s，最大的误差为1.3s。

a:近震定位前 b：近震定位后
图5.3 实例二近震定位前（a）和定位后（b）
为了进一步验证此种定位的正确与否，再加一例辅证。同样采用上述反演方法定位后，可确定该地震事件发生的位置为（96.7124,25.7474）。利用geotool软件的p波波至标定后的结果如下图5.4，其最小的误差为0. 09s，最大的误差为3s。

a:近震定位前 b：近震定位后
图5.4 实例二近震定位前（a）和定位后（b）

桂林理工大学本科毕业设计·论文

6 结束语
通过本次毕业设计，使我了解了近震震源参数定位的方法和原理，从正演和反演的理论知识落脚到实际资料的处理。解决此问题的基本思想就是将非线性问题化为线性问题，然后利用最小二乘法原理迭代，全主元消元法等求解，通过反复修定，以得到震源参数的最佳值。
此次设计的目的在于借助震源的计算机初步定位这个问题，进一步了解一般科学计算研究所采用的正、反演计算方法及实例应用。
本次毕业设计是在王教授悉心指导下，经过不断的学习和修改完成的。王教授渊博的学识，丰富的实践经验，高瞻远瞩、敏锐的科学眼光，将是我永远学习的楷模；他乐观的心态、谦逊的为人、朴实的生活态度，令我深深敬佩；他严谨细致、一丝不苟的作风将是我以后学习、工作中的榜样，他循循善诱的教导和不拘一格的思路给予我无尽的启迪。谨致以衷心的感谢和崇高的敬意。
本次毕业论文的撰写也得到众多老师和同学的关心和帮助。在此一并表示衷心的感谢。祝愿他们身体健康，工作顺利，事业上取得更大成功。　　
由于本人水平有限，本篇论文及程序有待进一步改进，请各位老师批评指正。

参考文献

[1] Seth Stein,Michael Wysession. An Introduction to Seismology Earthquakes,and Earth Structure.Blackwell Publishing, 2003.
[2] 李善邦，中国地震. 地震出版社,1981.
[3] 单娜琳,程志平,刘云祯.工程地震勘探. 冶金工业出版社，2006.
[4] 姜枚,王有学,钱辉.造山的高原—青藏高原及其邻区的宽频地震探测与地壳上地幔结构.地质出版社，2009.
[5] 熊章强,方根显.浅层地震勘探.地震出版社，2002.
[6] 大港油田科技丛书编委会.地震勘探资料处理和解释技术.石油工业出版社，1999.
[7] 时振梁,张少泉,赵荣国等.地震工作手册.地震出版社，1990.
[8] 傅淑芳,刘宝城,李文艺.地震学教程（上、下册）. 地震出版社，1980.
[9] 彭国伦.Fortran 95 程序设计.中国电力出版社， 2002.

附录 程序
program test
implicit none
integer :: i, j, k, nstn, istn
real*8 :: x0, y0, v0, t0, xx, yy, coef(1000,5)
real*8 :: delt(1000), solut(4,5), x(1000), y(1000), t(1000)
real*8 :: ph0, la0, phi, lai ,v,h, MinLat, MaxLat, MinLon, MaxLon character :: filename*100

write(*,*) "input filename"
read(*,'(a)') filename
!read the picks for all stations
open(1,file=filename)
read(1,*) nstn
x0=0;y0=0
do istn=1,nstn
read(1,*) x(istn),y(istn), t(istn)
x(istn)=x(istn)*atan(1.00)/45
y(istn)=y(istn)*atan(1.00)/45
x0=x0+x(istn)
y0=y0+y(istn)
end do
la0=x0/nstn
ph0=y0/nstn
close(1)

goto 10
!calculating sythetic data
la0=99*atan(1.00)/45
ph0=27*atan(1.00)/45
open(2,file='mod13.dat')
write(2,'(I2)') nstn
v=5.0
h=10.0
do istn=1,nstn
lai=x(istn)
phi=y(istn)
call lp2xy(lai,phi,la0,ph0,xx,yy)
t(istn)=sqrt(xx**2.00+yy**2.00+h**2.00)/v
write(2,'(6f15.6)') x(istn)*45.00/atan(1.00), y(istn)*45.00/atan(1.00), t(istn)+6
end do
write(2,'(4f20.6)') la0*45.00/atan(1.00),ph0*45.00/atan(1.00),v,h
close(2)
！stop

10 continue
!(lai,phi)---->(x,y)
do istn=1,nstn
lai=x(istn)
phi=y(istn)
call lp2xy(lai,phi,la0,ph0,xx,yy)
x(istn)=xx
y(istn)=yy
end do
!calculating coefficents for the linear equations
do istn=2,nstn
coef(istn,1)=x(istn)-x(1)
coef(istn,2)=y(istn)-y(1)
coef(istn,3)=(t(istn)**2-t(1)**2)/2.00
coef(istn,4)=-t(istn)+t(1)
delt(istn) =(x(istn)**2-x(1)**2+y(istn)**2-y(1)**2)/2.00
end do

!least squar-root method to get the solution
! dt(1) de(1,1) de(1,2) de(1,3) de(1,4) de(1,5)
! dt(2) de(2,1) de(2,2) de(2,3) de(2,4) de(2,5)
! .
! .
! .
! dt(n) de(n,1) de(n,2) de(n,3) de(n,4) de(n,5)

! de(1,1) de(2,1) ... de(n,1)
! de(1,2) de(2,2) ... de(n,2)
! de(1,3) de(2,3) ... de(n,3)
! de(1,4) de(2,4) ... de(n,4)
! de(1,5) de(2,5) ... de(n,5)

do i=1,4
do j=1,4
solut(i,j)=0
do k=1,nstn
solut(i,j)=solut(i,j)+coef(k,i)*coef(k,j)
end do
end do
end do

do i=1,4
solut(i,5)=0
do k=1,nstn
solut(i,5)=solut(i,5)+coef(k,i)*delt(k)
end do
end do

call gaussian(solut,4,5)
xx=solut(1,5)
yy=solut(2,5)
v0=sqrt(solut(3,5))
t0=solut(4,5)/solut(3,5)
h=0
！print *,xx,yy
do istn=1,nstn
h=h+sqrt(v0*v0*(t(7)-t0)*(t(7)-t0)-(x(7)-xx)**2-(y(7)-yy)**2)
end do
h=h/nstn

call xy2lp(xx,yy,la0,ph0,lai,phi)
print *,'------lai phi v0 t0 h-------'
print *, lai*45.00/atan(1.00),phi*45.00/atan(1.00),v0,t0, h
end program test

SUBROUTINE Gaussian(A,N1,M1)
implicit none
integer ::n1,m1,k1,k2,ik,i,j ,l, i1
REAL*8 :: A(n1,m1)
REAL*8 :: BMAX,T,EPS

EPS=0.
DO K1=1,N1
BMAX=0.
DO I=K1,N1
IF(BMAX-ABS(A(I,K1)).LT.0) THEN
BMAX=ABS(A(I,K1))
L=I
END IF
END DO
IF(BMAX.LT.EPS) STOP 4444
IF(L.NE.K1) THEN
DO J=K1,M1
T=A(L,J)
A(L,J)=A(K1,J)
A(K1,J)=T
END DO
END IF
T=1./A(K1,K1)
K2=K1+1
DO J=K2,M1
A(K1,J)=A(K1,J)*T
DO I=K2,N1
A(I,J)=A(I,J)-A(I,K1)*A(K1,J)
END DO
END DO
END DO
DO IK=2,N1
I=M1-IK
I1=I+1
DO J=I1,N1
A(I,M1)=A(I,M1)-A(I,J)*A(J,M1)
END DO
END DO
RETURN
END SUBROUTINE Gaussian

subroutine lp2xy(lai,phi,la0,ph0,xx,yy)
implicit none
real*8 :: ph0, la0, phi, lai, ee, a, r1,r2, r3, dl, dp, sinxcos, xx, yy
ee=0.0066935165
a=6378.160
sinxcos=sin(ph0)*cos(ph0)
R1=sqrt(1-ee*sin(ph0)**2.00)
r2=a/R1
r3=a*(1-ee)/R1**3
dl=lai-la0
dp=phi-ph0
xx=r2*dl*cos(ph0)
yy=r3*dp+0.50*r2*dl**2.00*sinxcos

return
end subroutine lp2xy

subroutine xy2lp(xx,yy,la0,ph0,lai,phi)
implicit none
real*8 :: ph0, la0, phi, lai, ee, a, r1,r2, r3, sinxcos, xx, yy
ee=0.0066935165
a=6378.160
sinxcos=sin(ph0)*cos(ph0)
R1=sqrt(1-ee*sin(ph0)**2.00)
r2=a/R1
r3=a*(1-ee)/R1**3
lai=xx/(r2*cos(ph0))+la0
phi=yy/r3-0.50*xx**2*tan(ph0)/(r2*r3)+ph0

return
end subroutine xy2lp

当速度给定时，在求解方程参数和未知数的量稍有变化。即
do istn=2,nstn
coef(istn,1)=x(istn)-x(1)
coef(istn,2)=y(istn)-y(1)
coef(istn,3)=(-t(istn)+t(1))*v**2
delt(istn)=(x(istn)**2-x(1)**2+y(istn)**