|
/ c: i" A% c; k% r3 F
《测绘学报》 , i9 l% H7 k* \, ~* D5 S
构建与学术的桥梁 拉近与权威的距离
& f) U5 @ Z5 h' z/ l 海底控制点定位的半参数平差模型法 + E" q! Z# E: l |
孙文舟1 , 殷晓冬1, 暴景阳2, 曾安敏3
3 b0 c: q$ y L. R) Q0 {- S+ ] 1. 海军大连舰艇学院军事海洋与测绘系, 辽宁 大连 116018; a- b4 e/ I$ Z5 J' |! h
2. 武汉大学测绘学院, 湖北 武汉 430079; 4 B2 ?) K' ]2 a: ~1 Q5 }4 l
3. 地理信息工程国家重点实验室, 陕西 西安 710054
. G$ v* h- j4 u @ G/ | 收稿日期：2018-04-27；修回日期：2018-09-10 / }8 ] ?# ]/ \ w/ e) c( J
基金项目：国家重点研发计划（2016YFB0501701）；国家863计划（2013AA122501）；国家自然科学基金（41474015；41874016）
7 E- O' r9 R) K7 @" t5 u 第一作者简介：孙文舟(1991-), 男, 博士生, 研究方向为海洋大地测量。E-mail:1519374228@qq.com
/ N+ W) q" j h" q9 F- y, b 摘要：针对圆走航确定海底控制点三维绝对坐标存在垂直解精确度差的问题，本文提出了基于半参数平差模型处理测距系统误差的方法，首先分析了声速变化对测距误差的影响，其次在圆走航观测方程的基础上建立了水下控制点三维坐标解算的半参数平差模型，分别给出了观测时长为测距误差长周期项整数倍与非整数倍时的解算方法。仿真试验表明，当内波存在且观测时长不足测距误差长周期项的整数倍时，该方法相比于差分法和最小二乘法可明显提高水下控制点垂直解的精确度。 / C& [5 V& `% i" G3 P! N
关键词：半参数平差模型 圆走航 控制点坐标定位 距离交会定位 测距系统误差
+ f# a9 G5 A# G4 O- c# C7 n Semi-parametric adjustment model methods for positioning of seafloor control point 2 G& H8 V5 j, |: @' @3 R. @ I% L/ n
SUN Wenzhou , YIN Xiaodong1, BAO Jingyang2, ZENG Anmin3 % C- G: J% t% n1 |5 K6 @& `
Abstract: This paper focuses on solving the problem of seafloor control point absolute positioning with low vertical accuracy based on the survey ship sailing circle. The method of dealing with the systematic error based on semi-parametric adjustment model was proposed. Firstly, the influence of sound speed change on ranging error is analyzed. Secondly, a semi-parametric adjustment model for determining three-dimensional coordinates of underwater control points was established. And respectively proposed solutions under two different conditions, the observation duration is an integral multiple or non-integer multiple of the long-period term of the ranging error. Simulation experiment results show that this method can obviously improve the accuracy of vertical solution of seafloor control point compared with difference technique and least square method when internal waves exist and observation duration is less than an integer multiple of the long-period term of the ranging error. , e P3 d6 Z; S: q2 P& c' `* @
Key words: semi-parametric adjustment model sailing circle positioning of seafloor control point intersection positioning model systematic ranging error " N! u# f" E' e! }7 e6 Q
高精度确定海底控制点的三维坐标是建立海底大地控制网的关键环节[1]。船载GNSS结合水下声学测距是实施水下声标(控制点)交会定位的有效手段[2-10]，影响这种方法定位精度最主要的因素是声速在时间和空间上复杂变化所引起的测距误差。采用圆走航的方式可以保证水面观测点相对于海底控制点对称分布[11]，对提高控制点水平坐标的精确度有重要帮助，然而垂直方向的解依然较差且不稳定。
. s+ \3 Q M* [; Q8 e 为了提高海底控制点垂直解的精确度，国内外学者进行了不同的尝试。文献[12]提出了利用观测历元之间差分的方法，这种方法可以有效地消除声速误差长周期项的影响，但是短周期项依然存在。文献[13]利用传播时间建立了声速误差模型，提出了在500 m水深情况下，不依赖声速剖面的水下静态目标三维坐标快速求解方法。文献[14]认为声速误差项为常数，通过南海的试验验证了该方法相比于传统方法的有效性。近年来，压力传感器的测深精度不断提高，又有学者提出了深度约束下的海底控制点坐标确定方法[15-17]，进一步提高了三维坐标解的稳健性和精度。
; G' j: E7 g. D+ v 上述方法都是对声速变化引起的测距误差中可以参数化的部分进行建模改进。根据文献[2, 18]对水下声学测距误差的研究，由内波引起的测距误差形态复杂，难以用少量的参数表达，而引入参数过多又往往会导致法方程的病态甚至秩亏。因此，本文尝试利用半参数平差模型进行控制点坐标的解算，以期提高垂直解的精确度。 ! v. t& N* @0 \4 `: h
1 海底控制点坐标半参数平差模型解算方法1.1 声速变化引起的测距误差分析
. I$ m3 u7 `4 D: K' x9 T v$ L" ? 影响声速变化的主要因素是太阳辐射和海水运动[19]。声速剖面的变化包括日变化、季节变化以及附加之上的随机扰动，声速剖面的季节变化主要与太阳直射点在南北回归线之间周期性的运动有关，而声速剖面日变化主要与太阳高度角的日变化有关，由外力扰动引起的海洋内波造成水体微团垂向运动，使得上下层海水混合引起温度的变化从而导致海水声速变化[20]。文献[2]对声学测距误差的研究表明，由声速变化引起的测距系统误差项包含两部分，长周期项与短周期项，其中长周期误差项与潮汐的影响有关，而短周期误差项则主要与海洋内波现象有关。之后文献[18]在夏威夷岛附近进行了相关试验，再一次验证了这个结论，并且认为声速主要受温度变化影响，盐度变化对其影响微乎其微[18]。 ; H, R* d5 T; z5 Q1 h% F& U
为进一步分析文献[2, 18]的试验结论，研究声速对测距误差的影响规律，本文进行了如下试验，首先利用文献[21]的标准方程构建了背景声速剖面[21-22]  ! `& @ n5 `9 x3 a6 ~9 M: ~
(1)
0 w6 N) E# s2 x 式中，C1是声道轴处的声速；z1是声道轴的深度；ε=7.4×10-3是扰动系数；B=1300为声道厚度尺度；z是水层深度。
- o B3 @7 Z2 |% Q/ h5 N z3 s 利用表层温度的变化规律以及热传导方程构建声速时空变化的部分[23] 
3 Q g; |" E. e( {8 y (2) x8 E6 J! ?: Y6 R8 w
式中，Δc0是声速变化的振幅；ω=2π/T是声速变化的角速率；T为周期；A是涡动热传导系数；ρ是海水密度；φ是初始相位角。
) p$ i Q3 {) C- a6 R2 R( c& w 式(2)可表示声速剖面在不同时间不同深度的声速日变化值，通过式(1)和式(2)可构造一个稳定变化的声速剖面。 , n1 s: ?5 S7 `2 m( U
在海底控制点坐标解算的过程中往往采用某一固定的声速剖面，用上述方法构建的声速剖面簇的平均值作为固定声速剖面，计算在水深为3000 m，声速采样间隔为5 m的条件下，声速变化所引起的测距误差。为了不引入入射角偏差而导致的误差，仿真试验采用垂直波束，设表层温度日变化的振幅为0.2℃，每小时采样一次，共采样24次，表层海水声速变化与测距误差的计算结果如图 1所示。 
5 I d' X$ E- k6 q3 d1 v 图 1 表面声速变化与测距误差曲线Fig. 1 Surface sound velocity variation and ranging error curve
4 l0 Z& ]. z& ^$ c6 p 温度的日变化只影响表层海水声速，图 1(a)是对表层海水等间隔取样得到的海水声速平均值在时间上的变化，图 1(b)是各采样时间所对应的测距误差的变化。从图中可以看出，测距误差的变化规律与声速的变化规律基本相似，周期相同但相位存在差异。
! N; K0 [' i3 V* I 文献[2, 18]的试验结果显示，测距误差长周期项的变化与潮汐的变化规律相近。这种现象是由于潮波的振动产生了日潮内波或半日潮内波，在潮内波的影响下海水产生垂向运动，从而导致了上层温度较高的海水与下层温度较低的海水进行热量交换，引起了海水温度的变化，温度变化是引起海水声速变化最主要的因素，其周期与潮内波的周期相同。测距误差短周期项的变化与海洋内波有关，与潮内波的作用机理相同，海洋内波造成了海水垂向运动，引起了海水温度的变化进而导致声速的变化。相比于声速长周期项变化，短周期项变化较为复杂，因此与其对应的测距误差短周期项随时间变化同样具有复杂的形式，难以进行参数化建模。 $ e- O" X5 Y+ w9 ^5 L. a n
1.2 圆走航距离交会定位模型
" i) p7 q) h, U# ~, U' _7 m- n 水下控制点定位的观测方程通常表示为 
: x0 l+ B1 \. a (3)
6 U1 M6 w/ B4 h3 I 式中，ρi是第i时刻换能器到海底应答器的距离观测值；f(Xi, Xo)是两者之间的直线距离，Xi是第i时刻换能器的位置，Xo是海底应答器的位置；δρdi是由应答器电路延迟引起的系统性误差，是可计算的改正量；δρvi是声速变化及声线弯曲引起的系统性误差，为主要误差；εi是第i时刻测量的随机噪声。
8 S1 `" `; I* y3 g 若δρdi已通过外部设备进行改正，则式(3)线性化展开得[14](4)
% D: z ^1 D! m( n 式中，ρi0是第i时刻距离观测量的近似值；li=(xi-x0)/ρi0, mi=(yi-y0)/ρi0, ni=(zi-z0)/ρi0；dx、dy、dz是三维坐标的改正数，式(4)的误差方程可写为 
6 `6 U- f! J" z6 ]. L% O (5)
, T1 Y( s; G% I% q& v$ p1 O 若进行了n次观测，可将式(5)写成矩阵形  # q" n2 `' `% t8 U3 i6 P; c
(6)
" X6 a* Y, o: ?- e 式中，X=[dxdydz]T;L=[ρ10-ρ1ρ20-ρ2… ρn-ρn0]T; δρv=[δρv1 δρv2 … δρvn]T是由声速变化引起的测距误差向量。
" O( B% @ V5 [: w 由上文分析可知，由声速变化引起的测距误差变化复杂，尤其是测距误差的短周期项，难以对其进行参数化的建模，若将每次观测的系统误差δρv都作为未知参数，则此时未知参数的个数为n+3，而观测方程的个数为n，因此式(6)无法得到未知参数的唯一解。本文采用半参数模型进行控制点三维坐标的解算。 4 a) j$ T* i% N" N% ~
1.3 控制点坐标解算的半参数模型
& h1 c7 q7 P+ X 根据式(4)，可写出如下矩阵形式的观测方程 (7) 6 y x! ?' A) Z# i8 N7 s6 m
式中，S=[δpv1 δpv2 … δpvn]T，是由声速变化引起的测距系统误差的n维未知向量。式(7)即为海底控制点坐标解算的半参数模型[24-25]。
- A6 y: ]# R: d' d* h 式(7)的误差方程是 : z! _& V0 Y& f L: |5 q/ M! A" ]
(8)
+ ` r$ D+ c/ ?0 [7 ^ 由于未知参数多于方程的个数，因此无法直接利用最小二乘原理进行求解，需要选择新的平差准则[22] 
% T3 ]2 @: }7 Y (9) " |6 ^% {/ F7 {6 ]7 H$ o5 o! E
式中，α是光滑因子，对两参数V和S起平滑作用；R是一个给定的正定矩阵，本文中选定R=GTG，其中  + J/ K. Z, {# z: W. n3 n
同时，若令所有历元声速变化引起的测距误差和为0，则又可设置一个限制条件  1 d1 h* y; F% M9 |. {
(10)
- V7 i! E8 r& D4 n8 q* I' t7 f 利用式(8)-式(10)可构造拉格朗日函数 5 d6 a: H; f$ ~6 D4 z* o: y7 a6 W
(11)
* X3 Y: i! R1 k7 p7 M 根据所构造的拉格朗日函数即可推导出式(8)-(9)的解，这里直接给出推导结果 
! e& ?+ i- @' @2 f# a8 g$ q6 d. d (12)
) p+ {% V* {/ {4 Y" g3 e 将计算得到的坐标改正数加上所设定坐标初始值即可获得海底控制点的三维绝对坐标。
c& ^+ i, y; Z; g4 f 需要说明的是，满足式(10)的条件是观测时长必须是测距误差长周期项周期的整数倍，否则很难保 
' M, n2 Q. w6 L+ E3 U! D/ _* p ，即测距误差之和为0这一约束条件。当观测时长不足长周期项周期的整数倍时，应当首先利用式(12)计算S向量，分析测距误差的变化规律，然后对声速误差长周期项进行参数化处理，可表示为 
0 S+ I6 Y9 O/ K" J" N: w4 K (13)
9 A' B0 g H8 R( L* s 式中，δpvl是测距误差的长周期项；δpv0是长周期项测距误差的振幅；T是测距误差的长周期项的周期；φ是初始相位。对长周期项参数化处理后，短周期项仍然认为不可参数化，同样根据式(12)计算坐标改正数，加上坐标初始值获得控制点三维绝对坐标。 1 P: K! i4 O: Y4 T7 P
2 仿真试验分析
% m% g1 ?% B& ]( j& m+ l. u 为了验证半参数平差模型解算海底控制点坐标的有效性，本文进行了如下仿真试验：水深设置为3000 m，应答器布放于海底，测量船以应答器为中心，半径为3000 m的圆航迹航行，速度为6节(约3 m/s)，每10 s实施一次水声测距。解算的坐标系采用笛卡儿坐标系，应答器为原点，x轴指向东，y轴指向北，z轴指向天顶。测距误差的仿真采用与文献[12]相同的方法，表示为  7 Z; k8 V- H3 i$ m4 p
(14) : k/ J) X2 Z' y, W7 ?" A
式中，误差共包含4项，第1项和第2项分别是声速变化引起的短周期项和长周期项，第3项是测区相关性误差，第4项为随机性误差。参照文献[2]在北太平洋的试验和文献[18]在夏威夷岛附近的试验，测距误差中的长周期项与潮汐变化规律近似，周期为12 h，振幅约20 cm，由内波引起的测距误差短周期项的周期从几十分钟到几个小时不等，振幅约为12 cm，因此令c1=12 cm、c2=20 cm、c3=2 cm，Tw=20 min，随机误差满足方差为5 cm的高斯分布，t0是初始时刻，t是任意时刻，x是t时刻测量船的三维坐标，x′是海底应答器的三维坐标。除此之外，设置潮汐的周期为12 h，振幅为5 m，海面波浪的周期为20 s，振幅为2 m，测量船水平方向的定位精度为5 cm，垂直方向定位精度为10 cm。
" l" V6 G. T# x) n8 C 试验中将本文的方法与其他两种方法进行对比，方法1是直接对观测值利用最小二乘进行解算，方法2是采用文献[12]所提出的差分算法进行解算，本文方法为方法3。
8 t- C3 k" J) e E; @& ?* ~' _ 第1次仿真试验的采样数为4320，采样时长刚好为测距误差长周期项的一个周期，试验中分别考虑存在内波和不存在内波两种情况，表 1给出了不同方法解算所得控制点坐标偏差的统计参数。 2 o% n9 w# C* `0 e: z" n, |! O% j
表 1 不同方法解算所得控制点坐标偏差的统计参数Tab. 1 Statistic parameters of underwater control points coordinate biases by different methods 3 R* C. i1 w- K5 n6 i
2 Y0 d( r& Z; m9 t; F cm方法|dx||dy||dz|meanSTDMSEmeanSTDMSEmeanSTDMSE无内波方法10.097 70.079 60.015 30.105 20.070 70.015 62.528 20.070 36.396 0方法20.076 30.046 10.007 70.065 50.050 40.006 61.261 60.644 21.965 3方法30.051 90.043 20.004 40.090 00.065 80.012 02.523 50.068 06.372 3有内波方法10.102 40.094 20.018 60.063 70.057 90.007 12.525 00.063 16.379 1方法20.065 90.066 70.008 40.071 40.036 80.006 314.745 21.770 0220.268 9方法30.085 00.057 20.010 20.079 40.055 00.009 12.534 90.062 56.429 4
: g; V% b& ]7 Q" c
7 m# A( V; V+ ~ Q" O- Y# G8 M 对比水平方向和垂直方向坐标偏差的精确度(MSE)可以发现，声速误差对控制点垂直解的精确度影响较大，而对水平解的精确度影响小，这与圆走航测距定位消除水平方向误差的设计相吻合，也被国内外研究所证实。对比两个试验采用差分算法计算的垂直坐标偏差的精确度，可验证文献[12]的结论，即差分可以消除测距误差中长周期项的影响，但是短周期项的影响依旧存在。当存在内波时，水平方向依旧可以达到较高精确度，但是垂直方向解的精确度会明显变差。将3种方法进行横向比较，当不存在内波时，差分算法是3种方法中垂直解精确度最高的，高于半参数模型和最小二乘法，而在实际的海洋环境中，更普遍的情况是存在海洋内波，此时，差分算法垂直解的精确度降低，而最小二乘法与半参数模型解的精确度较高。
# X4 E& c: z! H4 E 第2次试验的采样数为720，采样时长为测距误差长周期项周期的1/6，试验中同样分别考虑存在内波与不存在内波两种情况，不同方法解算控制点坐标偏差的统计参数见表 2。
( s8 T/ y0 V3 T- D+ N 表 2 不同方法解算所得控制点坐标偏差的统计参数Tab. 2 Statistic parameters of underwater control points coordinate biases by different methods - F, v8 }; Z* @/ i! v1 ]
, e; m s% o- y) x% Z9 u
cm内波状态方法|dx||dy||dz|meanSTDMSEmeanSTDMSEmeanSTDMSE无内波方法10.296 90.128 70.103 00.215 40.111 20.057 516.023 10.238 9256.791 7方法20.201 70.130 70.056 10.174 20.106 50.040 51.679 52.315 45.286 6方法30.289 10.208 60.122 70.432 60.219 70.230 62.225 91.957 98.404 6有内波方法10.466 70.268 60.282 70.209 40.214 80.085 415.916 50.181 1253.365 7方法20.449 50.160 10.225 10.239 80.123 00.071 117.395 59.809 2389.202 0方法30.460 80.151 40.232 90.228 00.165 90.076 84.25 812.867 625.532 6
* s* B, K5 u7 I" k: a( y9 X" h
! {* P) j2 g; Q- s' V) Q4 w 将试验2与试验1垂直坐标偏差的精确度进行对比可得，无论是否存在内波，当观测时长不是测距误差长周期项的整数倍时，最小二乘法垂直解的精确度都会降低。差分方法的结论与前一试验一致，整体精确度降低，这主要与观测数据的数量减少有关。半参数模型的垂直解则依旧可以得到较高的精确度。图 2-图 5中测距误差是指采用文献[12]的方法仿真得到的真实测距系统误差，而误差估值是指半参数模型在不同试验条件下对每次测距误差的估值，可以看出半参数模型可以较准确地描述测距误差的变化规律，对比图 2和图 4、图 3和图 5，可以发现在内波存在试验中，半参数模型对测距误差估值的结果明显发散，精度明显下降，表现为图 3中测距误差与误差估值的中误差分别为0.17和0.22 cm；图 5中测距误差与误差估值的中误差分别为0.10和0.18 cm。对比表 2中半参数模型在存在和不存在内波两种条件下垂直坐标偏差的精确度，可以得出内波存在时垂直解误差的精确度由6.43增加到25.53，证明内波会影响半参数模型解算的精确度，但依旧小于最小二乘法的253.37和差分算法的389.20。从表 2中还可以看出，两种条件下，虽然最小二乘法垂直解的精度高于半参数模型，但由于其准确度较低，所以其精确度远小于半参数模型。第2次试验对测距误差长周期项进行了参数化处理，结果显示，无内波时对振幅的估计为20.02 cm，而当内波存在时，对振幅的估计为15.92 cm，相对于所设定的值20 cm，估值的准确性明显降低，对应垂直方向坐标解的误差也增大，这个结果表明，当测距误差同时包含参数项和非参数项时，非参数项误差会影响对参数项估值的准确性。 
5 l: c( T; c( u 图 2 半参数模型对测距误差的估值(无内波)Fig. 2 Ranging error estimated by semi-parametric model (no internal waves) 
* f; g* W5 A; E0 V4 v9 \% X5 \ 图 3 半参数模型对测距误差的估值(有内波)Fig. 3 Ranging error estimated by semi-parametric model (with internal waves) 
- g# d; q# P1 i ^) H+ s3 Y0 x 图 4 半参数模型对测距误差的估值(无内波)Fig. 4 Ranging error estimated by semi-parametric model (no internal waves) 
1 s1 y' o# f* L/ ~ 图 5 半参数模型对测距误差的估值(有内波)Fig. 5 Ranging error estimated by semi-parametric model (with internal waves) 1 K0 p: N: p/ l7 E! ^. N
在控制点三维坐标解算的过程中，最小二乘法和差分算法要想达到最优解，都需要较为苛刻的条件，最小二乘法要求观测的时长为测距误差长周期项的整数倍，而差分算法要求不存在短周期项的测距误差，前者显然实际工程应用价值较低，因为测量船连续圆走航12 h甚至24 h难度大，而后者所需要不存在内波现象的海洋环境条件更加难以满足。因此半参数模型的优势则明显突出出来，它不但可以取得较为精确的坐标解算值，而且对客观环境的依赖度相对较低。 * O3 y& q! s$ ^
3 结论
8 x' K1 }. r B y& b% z+ u3 s 针对传统方法确定海底控制点三维坐标垂直解精确度差的问题，本文提出了利用半参数平差模型处理测距系统误差的方法，通过分析和试验得出了以下结论：
/ m- U% I& V1 I9 I+ W5 y$ r( A/ Z- r (1) 测距误差的周期性变化与声速周期性变化有关。海水通过涡动热传导和内波所引起对流热传导，将上层海水热量向下层海水传递，引起了温跃层声速的变化，从而产生了测距误差周期性的变化规律。
6 I. p. @- h. k" @7 d9 V (2) 由海洋内波等随机海洋现象引起的测距误差短周期项(非参数项)不会影响控制点坐标水平解的精确度，但会影响控制点坐标垂直解的精确度。 : c' g6 `; o G8 J
(3) 利用半参数平差模型可以有效处理由声速变化引起的测距系统误差，从而得出优于传统方法的海底控制点三维坐标解算结果。 $ ]& v' D: g, j& o7 b" ?2 j
【引文格式】孙文舟, 殷晓冬, 暴景阳, 等. 海底控制点定位的半参数平差模型法. 测绘学报，2019，48(1)：117-123. DOI: 10.11947/j.AGCS.2019.20180187
' |/ a K! u9 U" Y: T 
; h# J5 _' Y4 T 精
7 Z, `% L2 C) |7 s3 k 彩 P% B, z {" T, B
回
. \# g5 H& v; B4 L5 m$ ` 顾 ; D/ w! u8 X4 Z5 v& n9 V6 W
测绘学科2019年第1期网刊发布 8 w, M: T; u# U8 Y/ C
测绘学科18人上榜，第64批中国博士后基金面上项目出炉，
& H; c- M3 g( j3 |6 M, F 174人！四川、陕西、黑龙江、海南四大省测绘局招聘地学人才！ 2 z! P1 p; |$ @& F: q
重磅！中共自然资源部党组关于深化科技体制改革提升科技创新效能的实施意见 * W1 D' K$ \, C8 d' a
Nature排行榜：中国地学院所TOP 50！ ) q8 r$ N* F% o+ G% e4 V: W
李德仁院士：论遥感卫星商业化的道路—如何从事业走向商业 , `5 B+ E/ T( d/ @
2018年最新自然指数公布，哪些高校表现优异？ 7 ^9 g. z8 K( \! N" @
2018年中科院地学类期刊分区名单（完整版）！ * ?/ N& W- w" M! Y7 p8 H3 ~
重磅！关于2019年度国家自然科学基金项目申请与结题等有关事项的通告 - G, m. M) s6 \ z5 Q3 M9 i' r* ]
权威 | 专业 | 学术 | 前沿
) k7 V) B; x n. j* |8 o: f 微信投稿邮箱 | song_qi_fan@163.com
1 `+ y$ z; x$ F% }2 } 微信公众号中搜索「测绘学报」，关注我们，长按上图二维码，关注学术前沿动态。
7 R& r( J5 h3 \ 欢迎加入《测绘学报》作者QQ群： 297834524
# T2 U) w2 e) ~ N+ h9 F7 p/ c 进群请备注：姓名+单位+稿件编号 , f- D. [4 R! o; ^6 |
) N% `3 r2 ~ h' t; C+ J5 t
7 x- d& C% D& R5 l1 w" z1 l1 f1 s6 n9 @
. \ A" z$ O, ^
| 
