CN103792580B - 一种拖缆勘探导航前绘炮点的获取方法 - Google Patents

Abstract

本发明涉及海上石油地震勘探技术领域，特别涉及一种海上拖缆地震勘探前绘炮点的获取方法。包括：将测线起点和测线终点的大地坐标分别转换成对应的平面投影坐标；计算测线平面方位角；计算前绘炮点的平面投影坐标；将得到的前绘炮点的平面投影坐标转换成前绘炮点的大地坐标，判断是否结束。以及另一种方法，包括：计算测线大地线方位角；计算前绘炮点的大地坐标；将得到的前绘炮点的大地坐标转换成前绘炮点的平面投影坐标；判断是否结束。本发明的有益效果是简洁、高效，算法精度高，避免了截断误差的影响。

Description

一种拖缆勘探导航前绘炮点的获取方法

技术领域

本发明涉及海上石油地震勘探技术领域，特别涉及一种海上拖缆地震勘探前绘炮点的获取方法。

背景技术

海洋地震勘探中，拖缆地震勘探是主要的勘探方法，地震勘探船拖曳一条或多条电缆，随着枪阵激发，拖曳电缆采集地震波反射数据，船载定位设备采集导航定位数据，定位精度直接关系到地震作业施工效率以及采集精度。前绘炮点是海洋地震勘探中测线上根据工区配置所提前给出的该测线上的所有理论响炮点位置，作业前事先给出测线起始点和相邻炮点之间的距离，根据作业模式的不同，计算出测线上理论响炮点的精确位置，根据配置，以枪阵中心点或CMP点或船参考点到达前绘炮点响炮，预计响炮点与前绘炮点距离越接近，认为精度越高，反映到作业质量上，则作业质量越高。但前绘炮点的理论位置难以获取，对响炮精度难以掌握。

发明内容

（一）要解决的技术问题

本发明的一个要解决的技术问题是提供了一种效率较高、精度较高的拖缆勘探导航前绘炮点坐标的获取方法。

本发明的另一个要解决的技术问题是提供了另一种效率较高、精度较高的拖缆勘探导航前绘炮点坐标的获取方法。

本发明的还一个要解决的技术问题是提供了一种可选择的拖缆勘探导航前绘炮点坐标的获取方法，在使用前，用户其可以选择使用前面两种方法之一

（二）技术方案

对于拖缆勘探导航前绘炮点坐标的获取方法的技术主题，本发明是通过以下技术方案实现的：

对于一种方法的技术方案，一种拖缆勘探导航前绘炮点的获取方法，包括如下步骤：

步骤F1：在勘探区域内定义一条测线，所述测线的一端为测线起点，所述测线的另一端为测线终点，相应得到所述测线起点的大地坐标和所述测线终点的大地坐标；并定义即将进行作业的起始炮号、炮号递增值和炮间距；

步骤F2：将测线起点的大地坐标转换成测线起点的平面投影坐标，将测线终点的大地坐标转换成测线终点的平面投影坐标；

步骤F3：根据得到的测线起点的平面坐标和测线终点的平面投影坐标，计算测线平面方位角；

步骤F4：根据得到的测线平面方位角、所述起始炮号、所述炮号递增值、所述炮间距和测线起点的平面投影坐标，计算前绘炮点的平面投影坐标；

步骤F5：判断前绘炮点到测线终点的距离是否大于所述炮间距，如果前绘炮点到测线终点的距离大于等于所述炮间距，则继续推算下一个前绘炮点的平面投影坐标；如果前绘炮点到测线终点的距离小于所述炮间距，则当前炮点即为所述测线上最后一个前绘炮点，计算结束。

进一步，所述步骤F2具体包括：

根据公式

$\left\{\begin{array}{c}x=\mathrm{FN}+0.9996\left[\begin{array}{c}S+\frac{l^{2}N}{2}\sin B+\frac{l^4}{24}\mathrm{NB}{\cos }^{3}B(5-t^2+{9\mathrm{\η}}^2+{4\mathrm{\η}}^2)+\\ \frac{l^6}{720}N\sin B{\cos }^{5}B(61-{58t}^2+t^2)\end{array}\right]\\ y=\mathrm{FE}+0.9996[\mathrm{lN}\cos B+\frac{l^{3}N}{6}{\cos }^{3}B(1-t^2+{\mathrm{\η}}^2)+\frac{l^{5}N}{120}{\cos }^5(1-{18t}^2+{14\mathrm{\η}}^2-{58\mathrm{\η}}^2{t}^2)]\end{array}\right.$

将测线起点的大地坐标转换成测线起点的平面投影坐标，将测线终点的大地坐标转换成测线终点的平面投影坐标；

其中，（x、y）对应待求点的平面投影坐标，B、L分别为待求点的大地坐标的大地纬度和大地经度；S是赤道到待求点的子午线弧长；l=L-L0，是待求点经度L相对中央子午线L0的经差；N是待求点的卯酉圈曲率半径；S和N可根据WGS-84椭球的参数求得，0.9996是UTM投影比例因子，t=tanB，η=e'cosB，e'是参考椭球的第二偏心率，东方向偏移FE=500000米；南北偏移FN北半球=0，FN南半球=10000000米；

其中在计算测线起点的平面投影坐标时，所述待求点为所述测线起点；在计算测线终点的平面投影坐标时，所述待求点为所述测线终点。

进一步，所述步骤F3具体包括如下步骤：

步骤F301：根据公式计算测线长度；

步骤F302：根据公式计算测线平面方位角；

其中A是测线平面方位角，（x₀，y₀）为所述测线起点的平面投影坐标，（x_e，y_e）为所述测线终点的平面投影坐标。

进一步，所述步骤F4具体包括：

根据公式 $\begin{array}{c}{x}_n=(n-n_0)D\cos A+x_0\\ y_n=(n-n_0)D\sin A+y_0\end{array}$ 计算前绘炮点的平面投影坐标；

其中，(x₀,y₀)是所述测线起点的平面投影坐标，(x_n,y_n)是所求第n炮的前绘炮点的平面投影坐标，n是炮号，n₀是所述起始炮号，A是测线平面方位角，D是所述炮间距。

进一步，所述步骤F4还包括：将得到的前绘炮点的平面投影坐标转换成前绘炮点的大地坐标。

进一步，所述将得到的前绘炮点的平面投影坐标转换成前绘炮点的大地坐标具体包括：根据公式

$\left\{\begin{array}{c}B=B_f-\frac{N_f\mathrm{tsn}{B}_f}{R_f}\left[\begin{array}{c}\frac{D^2}{2}-(5+{3T}_f+{10C}_f-{4C}_f^2-9e{\′}^2)\frac{D^4}{24}\sin B\cos B+\\ (91+{91T}_f+{298C}_f+{45T}_f^2-{252e\′}^2-{3C}_f^2)\frac{D^6}{720}\end{array}\right]\\ L=L0+\frac{1}{\cos B_f}\left[\begin{array}{c}D-(1+{2T}_f+V_f)\frac{D^3}{6}+\frac{l^{3}N}{6}{\cos }^{3}N(1-t^2+{\mathrm{\η}}^2)\\ +(5-{2C}_f-{3C}_f^2+{8e\′}^2+{24T}_f^2)\frac{D^5}{120}\end{array}\right]\end{array}\right.$

将得到的前绘炮点的平面投影坐标转换为前绘炮点的大地坐标；

其中，B、L为待求的前绘炮点的大地坐标，L0为中央子午线经度，由前绘炮点坐标(x_n,y_n)，椭球第一偏心率e，第二偏心率e'，该点在参考椭球上的卯酉圈曲率半径N，椭球长半轴a，椭球短半轴b等参考椭球参数算得。B_f，C_f，T_f，D，η等参数均是计算过程中的辅助参数。

对于另一种方法的技术方案，一种拖缆勘探导航前绘炮点的获取方法，包括如下步骤：

步骤S1：在勘探区域内定义一条测线，所述测线的一端为测线起点，所述测线的另一端为测线终点，相应得到所述测线起点的大地坐标和所述测线终点的大地坐标；并定义即将进行作业的起始炮号、炮号递增值和炮间距；

步骤S2：根据所述测线起点的大地坐标和所述测线终点的大地坐标，计算测线大地线方位角；

步骤S3：根据测线大地线方位角、测线起点的大地坐标、所述起始炮号、所述炮号递增值和所述炮间距，计算前绘炮点的大地坐标；

步骤S4：判断前绘炮点到测线终点的距离是否大于所述炮间距，如果前绘炮点到测线终点的距离大于等于所述炮间距，则继续推算下一个前绘炮点的大地坐标；如果前绘炮点到测线终点的距离小于所述炮间距，则当前炮点即为所述测线上最后一个前绘炮点，计算结束。

进一步，所述步骤S2具体包括如下步骤：

步骤S201：将所述测线起点的大地坐标和所述测线终点的大地坐标分别投影到以椭球中心为球心的辅助球的球面上，分别得到测线起点的辅助球坐标和测线终点的辅助球坐标；

步骤S202：在所述辅助球的球面上解算，得到测线起点的辅助球坐标至测线终点的辅助球坐标在所述辅助球上的辅助弧长和辅助方位角；

步骤S203：将得到的辅助弧长和辅助方位角归算到椭球面上，得到测线起点到测线终点的测线大地线长度和测线大地线方位角。

进一步，所述步骤S3具体包括如下步骤：

步骤S301：将测线起点的大地坐标、测线大地线方位角和所推前绘炮点与所述测线起点的大地线长度投影到以椭球中心为球心的辅助球的球面上；

步骤S302：在所述辅助球的球面上解算，得到前绘炮点的辅助球坐标；

步骤S303：将前绘炮点的辅助球坐标归算到椭球面上，得到前绘炮点的大地坐标。

进一步，所述步骤S3还包括：将得到的前绘炮点的大地坐标转换成前绘炮点的平面投影坐标。

进一步，所述将得到的前绘炮点的大地坐标转换成前绘炮点的平面投影坐标具体包括：根据公式

$\left\{\begin{array}{c}x=\mathrm{FN}+0.9996\left[\begin{array}{c}S+\frac{l^{2}N}{2}\sin B+\frac{l^4}{24}\mathrm{NB}{\cos }^{3}B(5-t^2+{9\mathrm{\η}}^2+{4\mathrm{\η}}^2)+\\ \frac{l^6}{720}N\sin B{\cos }^{5}B(61-{58t}^2+t^2)\end{array}\right]\\ y=\mathrm{FE}+0.9996[\mathrm{lN}\cos B+\frac{l^{3}N}{6}{\cos }^{3}B(1-t^2+{\mathrm{\η}}^2)+\frac{l^{5}N}{120}{\cos }^5(1-{18t}^2+{14\mathrm{\η}}^2-{58\mathrm{\η}}^2{t}^2)]\end{array}\right.$

将得到的前绘炮点的大地坐标转换成前绘炮点的平面投影坐标；

其中，x、y分别为待求点的平面投影坐标，B、L分别为待求点的大地纬度和大地经度，S是赤道到待求点的子午线弧长，l=L-L0，是待求点经度L相对中央子午线L0的经差，N是该点的卯酉圈曲率半径，S、N根据WGS-84椭球的参数求得，0.9996是UTM投影比例因子，t=tanB，η=e'cosB，e'是参考椭球的第二偏心率，东方向偏移FE=500000米；南北偏移FN北半球=0，FN南半球=10000000米；

这里的待求点为前绘炮点。

对于还一种方法的技术方案，一种拖缆勘探导航前绘炮点的获取方法，包括如下步骤：

判断用户选择的计算模式；如果用户选择的计算模式为平面投影模式，则执行步骤F；如果用户选择的计算模式为大地线模式，则执行步骤S；

其中所述步骤F具体包括：

步骤F1：在勘探区域内定义一条测线，所述测线的一端为测线起点，所述测线的另一端为测线终点，相应得到所述测线起点的大地坐标和所述测线终点的大地坐标；并定义即将进行作业的起始炮号、炮号递增值和炮间距；

步骤F2：将测线起点的大地坐标转换成测线起点的平面投影坐标，将测线终点的大地坐标转换成测线终点的平面投影坐标；

步骤F3：根据得到的测线起点的平面坐标和测线终点的平面投影坐标，计算测线平面方位角；

步骤F4：根据得到的测线平面方位角、所述起始炮号、所述炮号递增值、所述炮间距和测线起点的平面投影坐标，计算前绘炮点的平面投影坐标；

步骤F5：判断前绘炮点到测线终点的距离是否大于所述炮间距，如果前绘炮点到测线终点的距离大于等于所述炮间距，则则继续推算下一个前绘炮点的平面投影坐标；如果前绘炮点到测线终点的距离小于所述炮间距，则当前炮点即为所述测线上最后一个前绘炮点，计算结束。

所述步骤S具体包括：

步骤S1：在勘探区域内定义一条测线，所述测线的一端为测线起点，所述测线的另一端为测线终点，相应得到所述测线起点的大地坐标和所述测线终点的大地坐标；并定义即将进行作业的起始炮号、炮号递增值和炮间距；

步骤S2：根据所述测线起点的大地坐标和所述测线终点的大地坐标，计算测线大地线方位角；

步骤S3：根据测线大地线方位角、测线起点的大地坐标、所述起始炮号、所述炮号递增值和所述炮间距，计算前绘炮点的大地坐标；

步骤S4：判断前绘炮点到测线终点的距离是否大于所述炮间距，如果前绘炮点到测线终点的距离大于等于所述炮间距，则继续推算下一个前绘炮点的大地坐标；如果前绘炮点到测线终点的距离小于所述炮间距，则当前炮点即为所述测线上最后一个前绘炮点，计算结束。

（三）有益效果

与现有技术和产品相比，本发明有如下优点：

通过简洁、高效的算法分别实现了大地线模式和平面投影模式下的前绘炮点的坐标解算，具有较高精度；

在对前绘炮点的坐标计算过程中，所有点位的推算都由测线起始点开始计算，避免了截断误差的影响，且可以沿着测线向两端分别推算；

支持平面投影模式和大地线模式下前绘炮点坐标的计算。

附图说明

图1是本发明一种具体实施方式的流程示意图；

图2是工区、测线及前绘炮点的示意图。

具体实施方式

下面结合附图对本发明做进一步的说明：

分三个步骤：数据输入，选择放炮模式，坐标解算。该方法首先输入测线的起始点经纬度和炮间距，为了保证解算坐标的精度，需将秒数位保留至小数点后两位。然后选择计算模式，分别是大地线模式和平面网格模式两种。最后通过坐标解算计算前绘炮点的精确位置，涉及到通用横轴墨卡多正反投影变换，大地坐标正反算，及各种格式经纬度之间的相互转换等，实现了海洋地震勘探中前绘炮点的大地坐标及平面投影坐标的解算，人工干预少，解算精度高。

本发明具体由以下步骤组成，参考图1，详细说明如下：

步骤一：数据输入

主要输入起始点经纬度，起始炮号，炮号递增值及炮间距四个参数，其物理意义参见图2说明如下：测线一般是客户要求或提前给出的，作业人员按照客户给出的测线起始点坐标来进行勘探；起点终点是指定的，起始炮号也是指定的；测线的起点就是起始炮点，起始炮号已经指定，比如，起始炮号为1000，递增值为1，炮间距为25米，则下一炮为沿测线方向往前25米的位置，炮号为1001。所有炮点都是在测线上。

可以从起点向前往终点方向推算前绘炮点坐标，也可从起点向相反方向推算测线反向延长线的前绘炮点坐标。

步骤二：放炮模式选择

选择不同的放炮模式，会产生不同的前绘炮点坐标。在海洋地震勘探作业中最常用的两种放炮模式分别是大地线模式和平面投影模式。两种模式所采用的地球参考系均为1984世界大地坐标系（WGS-84），其原点位于地球质心，其椭球的几何参数为：椭球长半轴a6378137.0m，扁率f1/298.257223563,椭球的其他几个参数如椭球短半轴，偏心率等均可由这两个参数求得。

步骤三：坐标解算

坐标解算是前绘炮点坐标计算的关键部分，分为如下四大任务：大地坐标与平面坐标之间的相互转换，测线方位角计算，前绘炮点坐标计算及判断是否计算结束。

（1）大地坐标与平面投影坐标之间的相互转换：地震勘探船在勘探过程中，提供的测线起始点位和工区范围均为大地坐标，而最终前绘炮点坐标要求为大地坐标及其平面投影坐标。因此需要通过坐标转换与计算，求出测线起始点及两种模式下前绘炮点的平面投影坐标及大地坐标。在海洋地震勘探作业中最常用的投影方式是通用墨卡托投影（UTM投影），其坐标转换的公式如下：

由大地坐标（B,L）转换为平面投影坐标（x,y）的UTM投影正算公式（精度为±0.001米）：

$\left\{\begin{array}{c}x=\mathrm{FN}+0.9996\left[\begin{array}{c}S+\frac{l^{2}N}{2}\sin B+\frac{l^4}{24}\mathrm{NB}{\cos }^{3}B(5-t^2+{9\mathrm{\η}}^2+{4\mathrm{\η}}^2)+\\ \frac{l^6}{720}N\sin B{\cos }^{5}B(61-{58t}^2+t^2)\end{array}\right]\\ y=\mathrm{FE}+0.9996[\mathrm{lN}\cos B+\frac{l^{3}N}{6}{\cos }^{3}B(1-t^2+{\mathrm{\η}}^2)+\frac{l^{5}N}{120}{\cos }^5(1-{18t}^2+{14\mathrm{\η}}^2-{58\mathrm{\η}}^2{t}^2)]\end{array}\right.$

其中，B、L为该点的大地纬度和大地经度，S是赤道到该点的子午线弧长，l=L-L0，是该点经度L相对中央子午线L0的经差，N是该点的卯酉圈曲率半径，S、N根据WGS-84椭球的参数求得，0.9996是UTM投影比例因子，t=tanB，η=e'cosB，e'是参考椭球的第二偏心率，东方向偏移FE=500000米；南北偏移FN北半球=0，FN南半球=10000000米；

由平面投影坐标（x,y）转换为大地坐标（B,L）的墨卡多（UTM）投影反算公式（结果以“度”为单位，精度为±0.0001″）：

$\left\{\begin{array}{c}B=B_f-\frac{N_f\mathrm{tsn}{B}_f}{R_f}\left[\begin{array}{c}\frac{D^2}{2}-(5+{3T}_f+{10C}_f-{4C}_f^2-9e{\′}^2)\frac{D^4}{24}\sin B\cos B+\\ (91+{91T}_f+{298C}_f+{45T}_f^2-{252e\′}^2-{3C}_f^2)\frac{D^6}{720}\end{array}\right]\\ L=L0+\frac{1}{\cos B_f}\left[\begin{array}{c}D-(1+{2T}_f+V_f)\frac{D^3}{6}+\frac{l^{3}N}{6}{\cos }^{3}N(1-t^2+{\mathrm{\η}}^2)\\ +(5-{2C}_f-{3C}_f^2+{8e\′}^2+{24T}_f^2)\frac{D^5}{120}\end{array}\right]\end{array}\right.$

其中：

L0为中央子午线经度，B_f，C_f，T_f，D，η和N_f等参数均是计算过程中的辅助参数，由平面投影坐标（x,y），椭球第一偏心率e，第二偏心率e'，该点在参考椭球上的卯酉圈曲率半径N，椭球长半轴a，椭球短半轴b等参考椭球参数算得。

更具体一些来说：

$N_f=\frac{(a^2/b)}{\sqrt{1+{e\′}^2*{\cos }^2{B}_f}}=\frac{a}{\sqrt{1-{e\′}^2*{\sin }^2{B}_f}}$

$B_f=\frac{a(1-e_2)}{{(1-e^2*{\sin }^2{B}_f)}^{3/2}}$

$e_1=\frac{1-b/a}{1+b/a}$

M_f=(X_N一FW)/k₀

T_f=tan²B_fC_f=e′²cos²B_f

$D=\frac{Y_E-\mathrm{FE}}{k_0{N}_f}$

其中e是椭球第一偏心率，e'是椭球第二偏心率，N是该点的卯酉圈曲率半径，a是椭球长半轴，b是椭球短半轴。

具体地，其它各量比如X_n、k₀、Y_e、N_f，R_f等，是本领域技术人员常用的变量，例如可参见吕志平、张建军、乔书波编著的大地测量学基础[M]解放军出版社，2005。

（2）测线方位角计算

前绘炮点坐标的推算需要起点坐标，炮间距，炮号及测线距离方位角四个参数，其中前三个参数都是已知的，因此测线方位角的计算是前绘炮点坐标解算的关键问题之一。

例如，测线起点P_s，已知其大地坐标为(B₀,L₀)，对应起始前绘炮号10000，炮号沿测线方向递增，炮间距为n米，测线终点P_e，已知其大地坐标为(B_e,L_e)。

在平面投影模式下：首先通过坐标转换将测线起始点P_s和P_e的大地坐标通过UTM投影正算转为平面投影坐标(x₀,y₀),(x_e,y_e)，然后在平面直角坐标系下求测线距离和方位角：

距离： $S=\sqrt{{(y_e-y_0)}^2+{(x_e-x_0)}^2}$

方位角： $A=\arcsin \left(\frac{y_e-y_0}{S}\right)$

其中，方位角需要进行象限判断后加常数才是真正的方位角。

在大地线模式下的测线计算需要进行大地主题反算，即已知P_s和P_e两点的大地坐标(B₀,L₀)和(B_e,L_e)，求解P_s至P_e点的大地线的长度S及A。本专利采用是贝塞尔大地主题反算法,贝塞尔大地问题解算的基本原理是首先建立起以椭球中心为球心，以任意长为半径的辅助球，并按以下几个步骤解算：

①将椭球面元素（P_s和P_e的大地坐标）投影到球面上；

②在辅助球球面上解算球面三角形，求得P_s至P_e点的在辅助球面上的弧长S'及方位角A。

③将球面元素（P_s至P_e点的在辅助球面上的弧长S'及方位角A）归算到椭球面上，求得P_s至P_e点的大地线的长度S及方位角A。

其中在方位角A在归算过程中保持不变。

下面，给出根据上述原理计算方位角的一个具体例子：

在平面格网模式下的测线方位角和距离的计算较为简单，在大地线模式下的测线计算需要进行大地主题反算，即已知P₁和P₂两点的大地坐标(B₁,L₁)和(B₂,L₂)，求解P₁至P₂点的大地线的长度及方位角，目前有多重进行大地主题计算的方法，本专利采用的是贝塞尔大地主题正反算,贝塞尔大地问题解算的基本原理是首先建立起以椭球中心为球心，以任意长为半径的辅助球，并按三个步骤解算：

①按一定条件将椭球面元素投影到球面上；

②在球面上解算大地问题；

③将求得的球面元素按投影关系换算为相应的椭球面元素。

其中，贝塞尔大地问题结算公式的三个投影条件是：

①是投影后球面上点的球面纬度等于椭球面上对应点的归化纬度；

②椭球面上两点间的大地线投影到辅助球面上为大圆弧；

③大地方位角投影后数值不变。

大地主题反算的整个流程为：

（a）椭球面元素投影到球面上

①由B计算μ：

${\tan u}_1=\sqrt{1-e^2}\tan B_1$

${\tan u}_2=\sqrt{1-e^2}\tan B_2$

②由经差l求队赢球面上的经差O：

l=λ+sinm[α'σ+β'sinσcos(2M+σ)]

其

$\mathrm{\α}\′=(\frac{1}{2}+\frac{e^2}{8}-\frac{1}{16}{k\′}^2)e^2-\frac{e^2}{16}(1+e^2){k\′}^2+\frac{3}{128}{e}^2{k\′}^4$

$\mathrm{\β}\′=\mathrm{\ρ}\′\′[\frac{e^2}{16}(1+e^2){k\′}^2-\frac{e^2}{32}{k\′}^4]$

k'²e²cos²m

σ、M、m，需迭代计算。

（b）解算球面三角形

①求σ：

cosσ=sinu₁sinu₂+cosu₂cosu₂cosλ

②求A₁,A₂：

${\tan A}_1=\frac{\sin \mathrm{\λ}}{{\cos u}_1{\tan u}_2-{\sin u}_1\cos \mathrm{\λ}}$

${\tan A}_2=\frac{\sin \mathrm{\λ}}{{\cos u}_2\cos \mathrm{\λ}-{\tan u}_1\cos u_2}$

（c）将球面元素归算到椭球面上

①A₁,A₂无需换算；

②将σ转换为S

$S=\frac{1}{\mathrm{\α}}[\mathrm{\ρ}-\mathrm{\β}\sin \mathrm{\σ}\cos (2M+\mathrm{\σ})+\mathrm{\γ}\sin 2\mathrm{\σ}\cos (4M+2\mathrm{\σ})]$

（d）判断象限

（3）前绘炮点坐标计算：

根据起点坐标、测线方位角、炮间距及炮号，即可推算出各炮点的坐标，为了避免截断误差对前绘炮点坐标计算的影响，每个前绘炮点坐标的计算均由起点坐标算起。

在平面投影模式下，前绘炮点的计算公式如下：

x_n(n-n₀)DcosA+x₀

y_n(n-n₀)DsinA+y₀

其中(x₀,y₀)是起点的平面投影坐标，(x_n,y_n)是所求第n炮的前绘炮点坐标，n是炮号，n₀是起始炮号，A是测线方位角，D是炮间距。

在大地线模式下，前绘炮点的计算均是在大地线上进行推算，需要采用大地主题正算的方法来计算，所谓大地主题正算，即已知P_s点的大地坐标(B₀,L₀)以及P_s至P_n点的长度S及方位角A，求P_n点的大地坐标(B_n,L_n)，其中S(n-n₀)D，D是炮间距。则大地主题正算的主要流程为：

①椭球面元素投影到球面上，将大地线的长度S及方位角A转化为球面元素（P_s至P_e点的在辅助球面上的弧长S'及方位角A）：

②解算球面三角形，求得P_n点对应的球面位置；

③将球面元素（P_n点对应的球面位置）归算到椭球面上，得到P_n点的大地坐标(B_n,L_n)；

其中在方位角A在归算过程中保持不变。推算出前绘炮点的大地坐标之后，还需对其进行UTM投影变换，以得到前绘炮点对应的平面投影坐标。

下面以一个例子说明上述大地主题正算的主要流程：

（a）椭球面元素投影到球面上：

1)由B₁求u₁：

${\tan u}_1=\sqrt{1-e^2}{\tan B}_1$

2)计算辅助量m和M：

$\tan M=\frac{{\tan u}_1}{{\cos A}_1}$

sinm=cosu₁sinA₁

3)将S转化为σ：

V₀DS，i=1,2,...

V_i=αS+βsinσ_i-₁cos(2M+σ_i-₁)+γsin2σ_i-₁cos(4M+2σ_i-₁)

$\mathrm{\α}=\frac{\mathrm{\ρ}\sqrt{1+{e\′}^2}}{\mathrm{\α}}(1-\frac{k^2}{4}+\frac{7}{64}{k}^4-\frac{15}{256}{k}^6)$

$\mathrm{\β}=\mathrm{\ρ}(\frac{k^2}{4}-\frac{1}{8}{k}^4+\frac{37}{512}{k}^6)$

$\mathrm{\γ}=\mathrm{\ρ}(\frac{1}{128}{k}^4-\frac{1}{128}{k}^6)$

k²e'²cos²m

其中V的单位为度，当时，停止迭代。

（b）解算球面三角形：

求A₂，u₂，λ：

${\mathrm{tsnA}}_2={\tan A\′}_2=\frac{\tan m}{\cos (M+\mathrm{\σ})}$

tanu₂-cosA₂tan(M+σ)

tanλ₁=sinmtanM=sinu₁tanA₁

tanμ₂=sinmtan(M+σ)=sinu₂tanA₂

λ=λ₂-λ₁

（c）将球面元素归算到椭球面上：

1)由u₂求B₂：

${\tan B}_2=\sqrt{1-e^2}\tan u_2$

2)求L₂：

L₂=L₁+l

l=O-sinm[α'σ+β'sinσcos(2M+σ)+γ'sin2σcos(4M+2σ)]

$\mathrm{\α}\′=(\frac{1}{2}+\frac{e^2}{8}-\frac{1}{16}{k\′}^2)e^2-\frac{e^2}{16}(1+e^2){k\′}^2+\frac{3}{128}{e}^2{k\′}^4$

$\mathrm{\β}\′=\mathrm{\ρ}\′\′[\frac{e^2}{16}(1+e^2){k\′}^2-\frac{e^2}{32}{k\′}^4]$

$\mathrm{\γ}\′=\mathrm{\ρ}\′\′\frac{e^2}{256}{k\′}^4$

k'²e²cos²m

（4）判断推算是否结束：

计算得到前绘炮点坐标之后，还需进行一次判断，以确定是否计算结束，判断的标准是当前绘炮点到测线终点的距离是否大于炮间距，若大于，则继续推算下一点坐标，若小于，则计算结束。

举个例子，推出当前炮点坐标后，比如1050炮的坐标，假设炮号递增，则下一炮炮号为1051，测线起始点不变，只是1051炮的炮点坐标相对测线起点的偏移量变了，同时相对终点的距离也变了。若1050炮离起点距离为1000米，炮间距25米，则下一炮离测线起点距离就为1025米，若是平面投影模式，则X偏移和Y偏移都随即要增加，根据这个偏移量即可推出1051点的XY坐标，大地线模式亦是如此，只不过对应的是经纬度的偏移，推出1051炮的坐标后，再计算1051炮到终点的距离，若大于炮间距，则继续推1052炮，若小于炮间距，则推算结束。

本发明的优点在于：

本发明通过简洁、高效的算法分别实现了大地线模式和平面投影模式下的前绘炮点的坐标解算，具有较高精度。

在对前绘炮点的坐标计算过程中，所有点位的推算都由测线起始点开始计算，避免了截断误差的影响，且可以沿着测线向两端分别推算。

另外，本领域技术人员还可理解的是，大地线模式在计算过程中一直使用的是经纬度，而平面投影模式一直计算使用的都是平面坐标，那么可以使每种模式下的前绘炮点都要有大地坐标系和平面投影系这两种坐标系下的坐标，也就是将每个模式下的坐标可以互相转换。

Claims (9)

1.一种拖缆勘探导航前绘炮点的获取方法，其特征在于，包括如下步骤：

步骤F1：在勘探区域内定义一条测线，所述测线的一端为测线起点，所述测线的另一端为测线终点，相应得到所述测线起点的大地坐标和所述测线终点的大地坐标；并定义即将进行作业的起始炮号、炮号递增值和炮间距；

步骤F2：将测线起点的大地坐标转换成测线起点的平面投影坐标，将测线终点的大地坐标转换成测线终点的平面投影坐标；

步骤F3：根据得到的测线起点的平面坐标和测线终点的平面投影坐标，计算测线平面方位角；

步骤F4：根据得到的测线平面方位角、所述起始炮号、所述炮号递增值、所述炮间距和测线起点的平面投影坐标，计算前绘炮点的平面投影坐标；

步骤F5：判断前绘炮点到测线终点的距离是否大于所述炮间距，如果前绘炮点到测线终点的距离大于等于所述炮间距，则继续推算下一个前绘炮点的平面投影坐标；如果前绘炮点到测线终点的距离小于所述炮间距，则当前炮点即为所述测线上最后一个前绘炮点，计算结束。

2.根据权利要求1所述的拖缆勘探导航前绘炮点的获取方法，其特征在于，所述步骤F3具体包括如下步骤：

步骤F301：根据公式计算测线长度；

步骤F302：根据公式计算测线平面方位角；

其中A是测线平面方位角，(x₀，y₀)为所述测线起点的平面投影坐标，(x_e，y_e)为所述测线终点的平面投影坐标。

3.根据权利要求1所述的拖缆勘探导航前绘炮点的获取方法，其特征在于，所述步骤F4具体包括：

根据公式计算前绘炮点的平面投影坐标；

其中，(x₀,y₀)是所述测线起点的平面投影坐标，(x_n,y_n)是所求第n炮的前绘炮点的平面投影坐标，n是炮号，n₀是所述起始炮号，A是测线平面方位角，D是所述炮间距。

4.根据权利要求1所述的拖缆勘探导航前绘炮点的获取方法，其特征在于，所述步骤F4还包括：将得到的前绘炮点的平面投影坐标转换成前绘炮点的大地坐标。

5.一种拖缆勘探导航前绘炮点的获取方法，其特征在于，包括如下步骤：

步骤S1：在勘探区域内定义一条测线，所述测线的一端为测线起点，所述测线的另一端为测线终点，相应得到所述测线起点的大地坐标和所述测线终点的大地坐标；并定义即将进行作业的起始炮号、炮号递增值和炮间距；

步骤S2：根据所述测线起点的大地坐标和所述测线终点的大地坐标，计算测线大地线方位角；

步骤S3：根据测线大地线方位角、测线起点的大地坐标、所述起始炮号、所述炮号递增值和所述炮间距，计算前绘炮点的大地坐标；

步骤S4：判断前绘炮点到测线终点的距离是否大于所述炮间距，如果前绘炮点到测线终点的距离大于等于所述炮间距，则继续推算下一个前绘炮点的大地坐标；如果前绘炮点到测线终点的距离小于所述炮间距，则当前炮点即为所述测线上最后一个前绘炮点，计算结束。

6.根据权利要求5所述的拖缆勘探导航前绘炮点的获取方法，其特征在于，所述步骤S2具体包括如下步骤：

步骤S201：将所述测线起点的大地坐标和所述测线终点的大地坐标分别投影到以椭球中心为球心的辅助球的球面上，分别得到测线起点的辅助球坐标和测线终点的辅助球坐标；

步骤S202：在所述辅助球的球面上解算，得到测线起点的辅助球坐标至测线终点的辅助球坐标在所述辅助球上的辅助弧长和辅助方位角；

步骤S203：将得到的辅助弧长和辅助方位角归算到椭球面上，得到测线起点到测线终点的测线大地线长度和测线大地线方位角。

7.根据权利要求6所述的拖缆勘探导航前绘炮点的获取方法，其特征在于，所述步骤S3具体包括如下步骤：

步骤S301：将测线起点的大地坐标、测线大地线方位角和所述前绘炮点与所述测线起点的大地线长度投影到以椭球中心为球心的辅助球的球面上；

步骤S302：在所述辅助球的球面上解算，得到前绘炮点的辅助球坐标；

步骤S303：将前绘炮点的辅助球坐标归算到椭球面上，得到前绘炮点的大地坐标。

8.根据权利要求5所述的拖缆勘探导航前绘炮点的获取方法，其特征在于，所述步骤S3还包括：将得到的前绘炮点的大地坐标转换成前绘炮点的平面投影坐标。

9.一种拖缆勘探导航前绘炮点的获取方法，其特征在于，包括如下步骤：

判断用户选择的计算模式；如果用户选择的计算模式为平面投影模式，则执行步骤F；如果用户选择的计算模式为大地线模式，则执行步骤S；

其中所述步骤F具体包括：

步骤F1：在勘探区域内定义一条测线，所述测线的一端为测线起点，所述测线的另一端为测线终点，相应得到所述测线起点的大地坐标和所述测线终点的大地坐标；并定义即将进行作业的起始炮号、炮号递增值和炮间距；

步骤F2：将测线起点的大地坐标转换成测线起点的平面投影坐标，将测线终点的大地坐标转换成测线终点的平面投影坐标；

步骤F3：根据得到的测线起点的平面坐标和测线终点的平面投影坐标，计算测线平面方位角；

步骤F4：根据得到的测线平面方位角、所述起始炮号、所述炮号递增值、所述炮间距和测线起点的平面投影坐标，计算前绘炮点的平面投影坐标；

步骤F5：判断前绘炮点到测线终点的距离是否大于所述炮间距，如果前绘炮点到测线终点的距离大于等于所述炮间距，则则继续推算下一个前绘炮点的平面投影坐标；如果前绘炮点到测线终点的距离小于所述炮间距，则当前炮点即为所述测线上最后一个前绘炮点，计算结束；

所述步骤S具体包括：

步骤S1：在勘探区域内定义一条测线，所述测线的一端为测线起点，所述测线的另一端为测线终点，相应得到所述测线起点的大地坐标和所述测线终点的大地坐标；并定义即将进行作业的起始炮号、炮号递增值和炮间距；

步骤S2：根据所述测线起点的大地坐标和所述测线终点的大地坐标，计算测线大地线方位角；

步骤S3：根据测线大地线方位角、测线起点的大地坐标、所述起始炮号、所述炮号递增值和所述炮间距，计算前绘炮点的大地坐标；

步骤S4：判断前绘炮点到测线终点的距离是否大于所述炮间距，如果前绘炮点到测线终点的距离大于等于所述炮间距，则继续推算下一个前绘炮点的大地坐标；如果前绘炮点到测线终点的距离小于所述炮间距，则当前炮点即为所述测线上最后一个前绘炮点，计算结束。