-
海底控制点三维坐标的确定通常采用GNSS(global navigation satellite system)结合声学观测量交汇的方法[1-6], 该方法首先由GNSS确定测量船的位置, 然后测量船在海面上按计划航线航行并开展水下声学测距, 利用距离交会的原理确定海底控制点的位置[7-8]。其中, 精确计算各观测历元测量船底换能器到海底应答器的斜距值是实现高精度海底控制点坐标测定的关键。
海底控制点定位中水下声学测距通常使用球面波, 所以声线的初始入射角未知, 而初始入射角是声线跟踪法计算斜距值的必要参数[9-10], 必须首先精确获得。传统的方法是由维·亚·科吾基娅和亚·伊·索罗金提出的, 它借助于其他手段获取的参考深度[11]、实测声速剖面和实测传播时间迭代计算初始掠射角, 然后根据最终迭代值所对应的水平位移和垂直位移计算船底换能器到海底应答器的斜距值; 文献[12]的方法是以换能器到应答器的大地线为约束条件, 不断调整初始入射角, 当声线跟踪法计算的水平位移与大地线的差值小于某一限差时, 所对应的初始入射角为所求值; 文献[13]的方法是根据应答器的概略坐标搜索本征声线, 从而确定初始入射角[14-16], 具体实现中以应答器和换能器间的垂直距离为约束条件进行声线跟踪, 以声线的到达位置和应答器概略坐标的距离为判别条件, 利用弦截法求得精确的初始入射角[17]。在计算过程中, 文献[12-13]都要用到应答器的概略坐标, 而由于概略坐标存在一定的偏差, 因此必然会对初始入射角的计算产生相应的影响。
为改善迭代算法的稳定性和准确性, 本文提出了利用牛顿法迭代计算初始入射角的方法, 首先介绍了传统方法的原理和迭代过程; 然后根据在深度和声速剖面一定的条件下传播时间随初始入射角的变化关系建立了迭代方程, 并给出了迭代计算的过程; 最后通过深海和浅海两组实测声速剖面开展了仿真实验, 对两种方法的优劣进行了对比和分析。
-
声线修正的基本原理是声线跟踪法, 在初始入射角和传播时间已知的条件下, 声线脚印位置相对于声源的水平位移和垂直位移可通过如下方法计算:首先, 忽略声速水平梯度的变化, 对水体进行垂直分层, 水层i内声线的水平位移xi和传播时间ti的计算如下[18]:
$$ {x_i} = \frac{{{{[1 - {{(p{c_i})}^2}]}^{1/2}} - [1 - {{\left( {p{c_{i + 1}}{)^2}} \right]}^{1/2}}}}{{p{g_i}}} $$ (1) $$ {t_i} = \frac{{{\rm{arcsin}}\left( {p{c_{i + 1}}} \right) - {\rm{arcsin}}\left( {p{c_i}} \right)}}{{pg_i^2{\rm{\Delta }}{z_i}}} \cdot {\rm{ln}}\left[ {\frac{{{c_{i + 1}}}}{{{c_i}}}} \right] $$ (2) 式中, p是snell常数; ci是水层i的表面声速; Δzi是水层i的厚度; gi是水层i内的声速梯度。
根据式(1)和式(2)以实测传播时间为判别条件, 采用深度追加水层的方法完成所有水层的计算, 具体步骤如下[19]:
1) 从换能器表面开始追加水层, 根据式(1)和式(2)计算声线在各层中的传播时间ti和水平位移xi;
2) 若实测的传播时间为T0, 累加各层的传播时间为$T = \mathop \sum \limits_{i = 1}^N {t_i}$, 通过比较T0和T的值判断是否要追加水层。若T < T0, 则需要追加水层, 重复步骤1);若T > T0, 则不需再追加水层, 多追加部分的厚度${\rm{\Delta }}z'$和水平位移${\rm{\Delta }}x'$计算如下:
$$ \left\{ {\begin{array}{*{20}{c}} {{\rm{\Delta }}z' = \left( {T - {T_0}} \right) \cdot {c_b} \cdot {\rm{cos}}{\theta _b}}\\ {{\rm{\Delta }}x' = \left( {T - {T_0}} \right) \cdot {c_b} \cdot {\rm{sin}}{\theta _b}} \end{array}} \right. $$ (3) 式中, cb和θb分别是最后一层的声速和入射角。获得多追加部分的厚度${\rm{\Delta }}z'$和水平位移${\rm{\Delta }}x'$后, 声线经过整个水体的水平位移x和垂直位移z分别为:
$$ \left\{ {\begin{array}{*{20}{c}} {z = \mathop \sum \limits_{i = 1}^N {z_i} - {\rm{\Delta }}z'}\\ {x = \mathop \sum \limits_{i = 1}^N {x_i} - {\rm{\Delta }}x'} \end{array}} \right. $$ (4) 通过式(4)可计算整个水层的水平位移和垂直位移, 进而根据勾股定理得到声源到脚印位置的斜距值。然而海底控制点定位的声波采用球面波, 初始入射角未知, 而初始入射角是声线跟踪法的必要条件, 因此需要迭代计算精确的初始入射角值。
-
传统方法迭代计算初始入射角示意图如图 1所示, O是换能器的位置, P是应答器的位置, θs是初始入射角的真值, zT是应答器与换能器的深度差。若初始入射角存在偏差δθs, 根据声线跟踪法的原理计算得到的应答器位置为P′, x和z分别为P′相对于换能器O的水平位移和垂直位移, 传统方法的思路是先利用已知观测量实测的传播时间和参考深度计算初始入射角的改正数, 再利用获得的改正数修正初始入射角的初值。
图 1 传统方法计算初始入射角的示意图
Figure 1. Schematic Diagram of Calculating Starting Incidence Angle by the Traditional Method
传统方法计算初始入射角改正数的方式是假设声线为直线, 即图 1中的线段OP′和OP为直线, 同时认为线段PP′相对于线段OP′和OP非常小, 因此近似认为$\angle OP{\rm{'}}P$为直角, 同时得到角α和β相等, 若OP′的斜距值为ρ, 线段PP′的斜距值近似认为是δθ·ρ, 根据图中的几何关系可得:
$$ \frac{x}{\rho } = \frac{{{z_T} - z}}{{{\rm{ \mathsf{ δ} }}\theta \cdot \rho }} $$ (5) 由式(5)可得初始入射角的改正数近似为:
$$ {\rm{ \mathsf{ δ} }}\theta \approx \frac{{{z_T} - z}}{x} $$ (6) 式中, 换能器和应答器之间的深度差zT无法精确获得, 通常用参考深度代替。
至此整理得到传统方法迭代计算初始入射角的步骤如下:
1) 由换能器的坐标(xT, yT)和应答器的概略坐标(x0, y0)计算初始入射角的初值θ0;
$$ {\theta _0} = {\rm{arctan}}\left( {\frac{{{x_T} - {x_0}}}{{{z_T} - {z_0}}}} \right) $$ (7) 2) 利用初始入射角的初值θ0和实测的传播时间T0, 根据声线跟踪法计算声线脚印位置相对于换能器的水平位移x和垂直位移z;
3) 根据式(6)计算初始入射角的改正数δθ, 修正初始入射角θ0, 得到新的初始入射角θ=θ0+δθ;
4) 返回步骤2), 将初始入射角的初值更新为θ, 重复步骤2)~4), 直到δθ小于某一限差, 终止计算, 最后一次更新的θ为所求值。
从以上步骤可以看出, 传统方法存在两方面的不足:一是在大入射角下, 水平位移x变大, 式(6)的δθ将减小, 这样容易导致迭代过程提前结束循环, 影响计算结果的准确性, 若是采用进一步减小限差的办法进行改善, 则又会大大影响计算的效率; 二是计算初始入射角的改正数δθ受初值θ0的影响, 定义迭代过程第一次得到的改正数δθ0为入射角改正数初值, 则当式(8)成立时, 表明迭代过程发散, 无法得到初始入射角值。
$$ \left| {{\theta _0} + {\rm{ \mathsf{ δ} }}{\theta _0} - {\theta _s}} \right| > \left| {{\rm{ \mathsf{ δ} }}{\theta _s}} \right| $$ (8) 式中, θs是真实的初始入射角; δθs是θ0和θs角度之差。式(8)表明传统方法存在迭代过程不收敛的情况。
-
由声线跟踪的基本原理可知, 在声速剖面和水深确定的条件下, 传播时间和初始掠射角具有唯一确定的函数关系, 若已知实测传播时间为T0, 可构建如下方程:
$$ T\left( \theta \right) = f\left( \theta \right) - {T_0} $$ (9) 式中, θ是初始入射角; T(θ)为对应的传播时间; f()代表声线跟踪算法, 所求的初始入射角等价于非线性函数T(θ)等于0时对应的θ值, 设θk是方程T(θ)=0时的近似根, 将函数T(θ)在θk处用泰勒级数展开, 并忽略二阶以上的项得到:
$$ T\left( \theta \right) \approx T\left( {{\theta _k}} \right) + T'\left( {{\theta _k}} \right) \cdot \left( {\theta - {\theta _k}} \right) $$ (10) 顾及到式(10), 可将方程式T(θ)=0近似表示为:
$$ T\left( {{\theta _k}} \right) + T'\left( {{\theta _k}} \right) \cdot \left( {\theta - {\theta _k}} \right) = 0 $$ (11) 式(11)为线性方程, 将其解记为θk+1, 得到初始入射角的迭代方程为:
$$ {\theta _{k + 1}} = {\theta _k} + \frac{{T\left( \theta \right)}}{{T'\left( \theta \right)}} $$ (12) 这就是利用牛顿法迭代计算初始入射角的原理, 其几何解释如图 2所示。
图 2 本文方法计算初始入射角的示意图
Figure 2. Schematic Diagram of Calculating Starting Incidence Angle by Our Proposed Method
θs是函数T(θ)等于0时横坐标的值, 为所求值, 若θ0是初值, 计算函数T(θ)在θ0处切线的斜率Tˊ(θ0), 该切线与横坐标的交点为θ1。依此类推, 计算f(θ)在θ1处的切线和该切线与横坐标的交点, 直到满足下式为止:
$$ \left| {{\theta _k} - {\theta _{k + 1}}} \right| < \varepsilon $$ (13) 式中, θk和θk+1分别是第k次和第k+1次计算初始入射角的结果; ε是某一设定的限差。因为函数T(θ)是单调递增的凹函数, 所以保证了切线的斜率始终大于零, 当初始入射角初值θ0大于所求值θs时, 从图 2中可以看出, 第k+1次更新的初始入射角值θk+1的数值始终在θk与θs之间, 当初始入射角初值θ0小于所求值θs时, 第一次迭代得到的初始入射角必然大于θs, 之后的计算过程将与初始入射角初值θ0大于θs的情况相同。从上述的分析可知, 初始入射角初值的选取对迭代计算的过程没有影响, 而且迭代的过程不会出现发散的情况, 因此牛顿法非常适合对非线性函数T(θ)的求解。
牛顿法迭代计算初始入射角的步骤如下:
1) 根据换能器的坐标(xT, yT)和应答器的概略坐标(x0, y0)计算初始入射角的初值θ0;
2) 计算函数T(θ)在θ0处切线的斜率Tˊ(θ0), 根据迭代方程(12)更新初始入射角的数值;
3) 判断式(13)的条件是否满足, 如果不满足, 返回步骤2);如果满足, 输出最后一次更新的初始入射角为所求值。
-
为对比分析两种迭代算法的有效性, 本文设计了两组对照实验, 两组对照实验的声速剖面如图 3和图 4所示。
图 3是某深海海域的实测声速剖面, 水深约为2 800 m, 图 4是某浅海海域的实测声速剖面, 水深约为56 m。仿真实验中设置初始入射角的范围为1°~89°, 每间隔1°选取一个初始入射角, 根据声线跟踪法分别计算在上述入射角和两组实测声速剖面条件下, 声线经过整个水层的水平位移和传播时间, 并把声线脚印位置当做海底控制点的真实位置, 把该传播时间当做实际的观测量。在不考虑参考深度偏差的情况下, 使用两种迭代方法分别计算不同初始入射角下的斜距值, 两种方法迭代所需的初值均设置为真实初始入射角加0.1°, 限差均设置为0.000 001°, 首先对比两种迭代算法在深海条件下程序的运行时间随初始入射角的变化, 结果如图 5所示。
图 5 两种方法计算初始入射角对应的运行时间
Figure 5. Running Time Corresponding to Each Starting Incidence Angle Calculated by Two Methods
从图 5可以直观地看出, 两种方法在初始入射角为1°和2°时计算效率均远远低于其他入射角。传统方法的计算效率有随初始入射角增大而降低的趋势, 而本文方法的计算效率趋势性变化不明显。当初始入射角小于78°时, 传统方法的计算效率高于本文的方法, 当初始入射角大于78°时, 本文方法的计算效率高于传统方法。在开始计算两种方法在不同实验条件下的测距误差之前, 首先仿真得出上述两条声速剖面下的T(θ)函数, 结果如图 6所示。
图 6 两条声速剖面下初始入射角和传播时间的关系
Figure 6. Relationship Between the Starting Incidence Angle and the Propagation Time Based on the sound velocity Profiles of Deep Sea and Shallow Sea
从图 6可以看出, 无论是深海还是浅海的声速剖面, 函数T(θ)均为单调递增的凹函数, 根据§1.3分析可知, 采用牛顿法求取初始入射角时迭代过程将不会出现发散或陷入死循环的状态。
图 7和图 8分别是两种方法在两条声速剖面下计算得到的测距误差与初始入射角的关系, 从计算结果中可以看出, 无论深海或是浅海的声速剖面, 传统方法计算得到的测距误差均随初始入射角的增大而发散, 尤其是当初始入射角超过60°时, 测距误差发散加剧, 主要原因是前文分析得到传统方法的第一点不足, 当初始入射角增大时, 水平位移增加, 改正数减小, 迭代过程容易提前结束循环, 导致测距误差随初始入射角的增加而增加, 同时, 在浅海声速剖面实验条件下, 当初始入射角为86°、87°、88°、89°时, 入射角改正数初值分别为-0.21°、-0.28°、-0.44°、-1°, 满足式(8)迭代过程发散的条件, 因此采用传统方法将无法得到所求值。相比而言, 本文方法计算的测距误差不随初始入射角的增大而发散, 在任何条件下的计算结果均很稳定。同时在浅海声速剖面大初始入射角的条件下, 传统迭代算法失效, 而本文方法依然可以获得准确的初始入射角值。
-
针对传统方法迭代计算初始入射角存在测距误差随初始入射角增加而发散和大入射角条件下迭代不收敛的问题, 本文提出了基于牛顿法的初始入射角迭代算法, 通过理论分析和实验验证得出了如下结论:
1) 在初始入射角初值和限差相同的情况下, 在绝大部分初始入射角范围内, 传统方法的计算效率高于本文的方法, 在少部分大入射角条件下, 本文方法的效率更高;
2) 传统方法计算得到的测距误差随初始入射角的增大而增加, 并且在大入射角条件下容易出现迭代过程不收敛的问题;
3) 本文方法计算得到的测距误差小于传统方法, 且不随初始入射角增加发生趋势性的变化, 也不受初始入射角初值的影响, 不会出现迭代过程不收敛的现象, 在任何角度下都可以获得准确的初始入射角值。
Calculating the Starting Incidence Angle by Iterative Method for Positioning Seafloor Control Points
-
摘要: 针对传统方法迭代计算初始入射角存在测距误差随入射角增加而发散和大入射角条件下迭代不收敛的问题, 提出了一种利用牛顿法迭代计算初始入射角的方法。首先介绍了传统算法的原理和步骤; 其次根据传播时间和初始入射角的函数关系建立了迭代方程; 最后通过深海和浅海两组仿真实验对两种方法的优劣进行了对比分析。结果表明, 当初始入射角小于78°时, 传统方法的计算效率高于所提方法; 除此之外, 当初始入射角大于60°时, 传统方法计算的测距误差开始发散, 而且在浅海仿真实验中, 当初始入射角大于85°时, 传统方法还出现了迭代不收敛的问题, 而所提方法在任何角度下均可收敛于真值。Abstract: Focuses on the problem of ranging error diverges with the increase of the starting incidence angle and non-convergence of iteration process at the large incidence angle by using the traditional method to calculate the starting incidence angle, this paper proposes a new method of calculating the starting incidenceangle by Newton method. Firstly, the principle and steps of traditional method are introduced. Secondly, the iteration equation is established according to the relationship between the propagation time and the starting incidence angle. Finally, the two methods are compared by two simulation experiments in deepsea and shallow sea. The results show that the computational efficiency of the traditional method is higher than our proposed method when the starting incidence angle is less than 78°. In addition, it is found out that when the starting incidence angle is greater than 60°, the ranging error calculated by the traditional method begins to diverge. In the shallow water simulation experiment, the traditional method also has the problem of non-convergence in the iteration process when the starting incidence angle is greater than 85°. In contrast, our proposed method can converge to the real value at any starting incidence angle.
-
Key words:
- seafloor control point /
- Newton method /
- ray-tracing method /
- starting incidence angle /
- ranging error
-
-
[1] Spiess F N, Chadwell C D, Hildebrand J A, et al. Precise GPS/Acoustic Positioning of Seafloor Reference Points for Tectonic Studies[J]. Physics of the Earth & Planetary Interiors, 1998, 108(2):101-112 http://www.sciencedirect.com/science/article/pii/S0031920198000892 [2] Gagnon K, Chadwell C D, Norabuena E. Measuring the Onset of Locking in the Peru-Chile Trench with GPS and Acoustic Measurements[J]. Nature, 2005, 434(7 030):205 http://europepmc.org/abstract/MED/15758997 [3] Bürgmann R, Chadwell D. Seafloor Geodesy[J]. Annual Review of Earth & Planetary Sciences, 2014, 42(1):509-534 http://www.irgrid.ac.cn/handle/1471x/956884 [4] Fujita M, Ishikawa T, Mochizuki M, et al. GPS/Acoustic Seafloor Geodetic Observation: Method of Data Analysis and Its Application[J]. Earth Planets & Space, 2006, 58(3):265-275 doi: 10.1186/BF03351923 [5] Chadwell C D, Spiess F N, Hildebrand J A, et al. Sea Floor Strain Measurement Using GPS and Acoustics[C]//Gravity, Geoid and Marine Geodesy. Berlin, Heidelberg: Springer, 1997: 682-689 [6] Watanabe S I, Ishikawa T, Yokota Y. Non-volcanic Crustal Movements of the Northernmost Philippine Sea Plate Detected by the GPS-Acoustic Seafloor Positioning[J]. Earth Planets & Space, 2015, 67(1):184-193 doi: 10.1186/s40623-015-0352-6 [7] Zhao Jianhu, Zou Yajing, Zhang Hongmei, et al. A New Method for Absolute Datum Transfer in Seafloor Control Network Measurement[J]. Journal of Marine Science & Technology, 2016, 21(2):216-226 doi: 10.1007/s00773-015-0344-z [8] 孙文舟, 殷晓冬, 马卓希, 等.浅海海底控制点水平坐标的简化计算方法[J].海洋测绘, 2018, 38(4): 31-33 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=hych201804006 Sun Wenzhou, Yin Xiaodong, Ma Zhuoxi, et al. Simplified Method for Horizontal Coordinates of Seafloor Control Points in Shallow Water[J]. Hydrographic Surveying and Charting, 2018, 38(4): 31-33 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=hych201804006 [9] Chadwell C D, Sweeney A D. Acoustic Ray-Trace Equations for Seafloor Geodesy[J]. Marine Geodesy, 2010, 33:164-186 doi: 10.1080/01490419.2010.492283 [10] Pedersen D, Melvin A. Acoustic Intensity Anomalies Introduced by Constant Velocity Gradients[J]. The Journal of the Acoustical Society of America, 1961, 33(4):465-474 doi: 10.1121/1.1908693 [11] Ballu V, Bouin M N, Calmant S, et al. Absolute Seafloor Vertical Positioning Using Combined Pressure Gauge and Kinematic GPS Data[J]. Journal of Geodesy, 2010, 84(1):65 doi: 10.1007/s00190-009-0345-y [12] Chadwell C D, Sweeney A D. Acoustic Ray-Trace Equations for Seafloor Geodesy[J]. Marine Geodesy, 2010, 33(2-3): 164-186 doi: 10.1080/01490419.2010.492283 [13] Sakic P, Valérie B, Crawford W C, et al. Acoustic Ray Tracing Comparisons in the Context of Geodetic Precise Off-shore Positioning Experiments[J]. Marine Geodesy, 2018, 1:1-16 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=10.1080/01490419.2018.1438322 [14] Mercer J A, Felton W J, Booker J R. Three-dimensional Eigenrays Through Ocean Mesoscale Structure[J]. Journal of the Acoustical Society of America, 1985, 78(1):157-163 doi: 10.1121/1.392552 [15] 王恕铨, Farmer D.求三维本征声线的一种新方法[J].声学学报, 1992, 17(2): 155-157 http://www.cqvip.com/QK/90012X/19922/862977.html Wang Shuquan, Farmer D.A New Method of Three-Dimensional Eigen-Ray-Tracing[J]. Acta Acustica, 1992, 17(2): 155-157 http://www.cqvip.com/QK/90012X/19922/862977.html [16] 张维, 杨士莪, 汤云峰, 等.不平整海底环境下的浅海本征声线求解方法[J].哈尔滨工程大学学报, 2011, 32(12):1 544-1 548 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=hebgcdxxb201112004 Zhang Wei, Yang Shi'e, Tang Yunfeng, et al. A Method of Seeking Eigen-Rays in Shallow Water with an Irregular Seabed[J]. Journal of Harbin Engineering University, 2011, 32(12):1 544-1 548 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=hebgcdxxb201112004 [17] 王百合, 冯西安, 黄建国, 等.一种分层海洋中求取本征声线的新方法[J].微处理机, 2006, 27(1):63-65 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=wclj200601021 Wang Baihe, Feng Xi'an, Huang Jiaguo, et al. A New Method for Eigenrays Searching in Stratified Ocean[J]. Microprocessors, 2006, 27(1):63-65 http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=wclj200601021 [18] 陆秀平, 边少锋, 黄谟涛, 等.常梯度声线跟踪中平均声速的改进算法[J].武汉大学学报·信息科学版, 2012, 37(5): 590-593 http://ch.whu.edu.cn/article/id/209 Lu Xiuping, Bian Shaofeng, Huang Motao, et al. An Improved Method for Calculating Average Sound Speed in Constant Gradient Sound Ray Tracing Technology[J]. Geomatics and Information Science of Wuhan University, 2012, 37(5): 590-593 http://ch.whu.edu.cn/article/id/209 [19] 赵建虎, 刘经南.多波束测深及图像数据处理[M].武汉:武汉大学出版社, 2008 Zhao Jianhu, Liu Jingnan. Multi Beam Sonnding and lmage Data Processing[M].Wuhan:Wuhan University Press, 2008 -