Numerical Simulation and Analysis of Tsunami Impacts on the Coastal Area of South China
-
摘要: 应用海啸数值模拟模型（cornell multi-grid coupled Tsunami model，COMCOT）对华南沿海地区可能遭受的海啸进行模拟研究具有现实意义。基于COMCOT模型及马尼拉海沟潜在震源分段及其参数设置，模拟了不同震源引发多场海啸到达华南沿海地区的波高变化和到时特征。模拟分析发现，海啸到达中国福建省南部至广东省西部一带能够产生最大2.5 m波高的海啸波，海南岛东部沿岸也会产生最大1.5~2.0 m的海啸波，海啸传播到广西北部湾地区的最大波高在0.3 m以下；海啸波发生1.9 h后最先到达海南岛东部沿岸，2.6 h后到达香港、澳门及广东的汕尾一带地区，3.6 h后到达福建南部一带及琼州海峡地区，8.9 h以后到达北部湾地区。华南沿海不同地区海啸的首波波高和最大波波高对马尼拉海沟震源位置的敏感性存在差异性，其中广东沿海中段对马尼拉海沟震源位置最为敏感，震源断层2、断层3引发海啸在该地区的首波和最大波高接近2.5 m，而震源断层1、断层6（断层编号依据2006年美国地质调查局资料）引发海啸在该地区的首波和最大波高只在1.0 m上下。Abstract: The coastal area of South China overlooks the Manila trench across the South China Sea. Once a tsunami triggered by the earthquake in the trench that has frequent seismic activity, the tsunami energy may spread northward to the coast of China and influence the region. Based on the focal parameters of Manila trench which is divided into six segments, COMCOT (cornell multi-grid coupled tsunami model) is adopted to simulate six potential tsunamis. The six earthquakes are close to Mw9.0, and the epicenter of them distributes from north to south across the trench. By analyzing the height and arrival time of tsunami wave which propagates to the coastal area of South China, the results indicate that, after 1.9 h when the tsunami occurs, the lead wave firstly propagates to the eastern coast of Hainan Island. It costs 3.6 h to get to Qiongzhou Strait and the south of Fujian Province, and the maximum wave can reach 1.5-2.0 m. 2.6 h later, the region of Hong Kong, Macao and Swabue of Guangdong Province can detect the lead waves, and the maximum wave in this region is about 2.5 m. 8.9 h later, the wave travels to the Beibu Gulf, and the wave amplitude is less than 30 cm. The sensitivity analysis shows that the sensitivity to source location is different in different coastal regions of the South China. In general, the most sensitive area is the region of Hong Kong, Macao and most parts of Guangdong Province. If the tsunami is triggered by Fault2 and Fault3, the maximum wave height is near 2.5 m. However, in the case of Fault1 and Fault6, the maximum wave height is about 1 m.
-
Keywords:
- south China coastal area /
- tsunami /
- COMCOT /
- numerical simulation /
- manila trench /
- fault
-
数据融合最初仅仅局限于军事领域，是针对异质传感器产生的海量数据进行一致性估计的方法[1-3]。20世纪90年代后，数据融合的概念被引入中国[4]。文献[5]将数据融合的结构分为3层，即数据层、特征层和决策层，在不同层进行数据融合会得到不同的结果，其差异体现在融合精度和处理效率等方面，此结构也作为基本的数据融合结构应用于其他领域。
在大地测量领域，文献[6]从函数模型误差和随机模型误差两方面详细讨论了不同的融合结果，具有较大的理论意义，相关的数据融合应用也较为广泛，体现在GPS与INS组合导航[7]、GPS与InSAR[8]、GPS不同系统间的融合[9]、遥感图像及地图数据融合[9, 10]以及3S集成[11]等方面。在工业测量中，相关的数据融合理论方法并不常见，文献[12]介绍了在盾构自动导向中综合利用全站仪和倾斜仪数据联合平差解算盾构姿态的方法，体现了采用多源异质传感器进行数据融合的优势。
在mm级精度的工业测量中，是采用全站仪进行特定点的三维坐标的采集，全站仪是测距仪(获取斜距)和经纬仪(获取角度)的集成，其获取的坐标实际是测距数据和测角数据的融合结果；而在船舶、飞机制造等工业测量领域，由于一些大型构件无法一站式获取全部待测点位三维信息，需采用移站测量(moving station method, MSM)的方法[13]，通过在不同测站上重复测量公共点，将不同测站的坐标系统一起来。但是，随着移站次数的增加，测量误差会不断积累，影响整体测量的精度。
本文以工业测量整体平差(overall adjustment, OA)为例，讨论了两种不同数据融合方法，即基于原始观测值(水平角、竖直角、斜距)的平差模型和基于坐标观测值(直角坐标系中的三维坐标)的平差模型，从解算精度、多余观测值的敏感性以及随机模型偏差对结果的影响等3方面，详细对比了两者的区别并和传统的移站测量结果进行比较，得到在工业测量中融合策略选取的相关结论。
1 工业测量整体平差函数模型
1.1 基于原始观测值的模型
在传统的工业测量中，采用全站仪获取数据是使用最多的方法之一。在每站测量前经过整平后，全站仪的原始观测值为测站点到待测点的水平方向、竖直角和斜距，这3个观测值互相独立，并在全站仪出厂时会给出角度和距离测量的标称精度。通常所有的待测点无法在一个测站全部观测，故采用移站测量的方法将不同测站的坐标系通过重复测量公共点进行统一，一般都以起始测站的坐标系作为全局坐标系。
在基于原始观测值的整体平差函数模型中，设第j测站在全局坐标系中的坐标为$ \left( {{{\hat{X}}}_{\text{stat}-j}}, \ \ {{{\hat{Y}}}_{\text{stat}-j}}, \ \ {{{\hat{Z}}}_{\text{stat}-j}} \right) $，测量了n个全局坐标系中的待测点$ \left( {{{\hat{X}}}_{i}}, \ \ {{{\hat{Y}}}_{i}}, \ \ {{{\hat{Z}}}_{i}} \right) $，其中i=1, 2, …, n，则有：
$$ \left\{ \begin{array}{l} {L_{\rm{s}}} + {V_{\rm{s}}} = \sqrt {{{\left( {{{\hat X}_i} - {{\hat X}_{{\rm{stat}} - j}}} \right)}^2} + {{\left( {{{\hat Y}_i} - {{\hat Y}_{{\rm{stat}} - j}}} \right)}^2} + {{\left( {{{\hat Z}_i} - {{\hat Z}_{{\rm{stat}} - j}}} \right)}^2}} \\ {L_{\rm{H}}} + {V_{\rm{H}}} = \arctan \frac{{{{\hat Y}_i} - {{\hat Y}_{{\rm{stat}} - j}}}}{{{{\hat X}_i} - {{\hat X}_{{\rm{stat}} - j}}}} - {{\hat A}_{{\rm{azi}} - j}}\\ {L_{\rm{V}}} + {V_{\rm{V}}} = \arctan \frac{{{{\hat Z}_i} - {{\hat Z}_{{\rm{stat}} - j}}}}{{\sqrt {{{\left( {{{\hat X}_i} - {{\hat X}_{{\rm{stat}} - j}}} \right)}^2} + {{\left( {{{\hat Y}_i} - {{\hat Y}_{{\rm{stat}} - j}}} \right)}^2}} }} \end{array} \right. $$ (1) 式中，Ls、LH、LV分别为斜距、水平角、竖直角观测值；对应Vs、VH、VV为观测值的改正数；$ {{\hat{A}}_{\text{azi}-j}} $为第j测站在全局坐标系中水平角初始零方向参数。将式(1)线性化后有：
$$ \left\{ \begin{array}{l} {V_{\rm{s}}} = {p_1}\delta {{\hat X}_{{\rm{stat}} - j}} + {p_2}\delta {{\hat Y}_{{\rm{stat}} - j}} + {p_3}\delta {{\hat Z}_{{\rm{stat}} - j}} - {p_1}\delta {{\hat X}_i} + {p_2}\delta {{\hat Y}_i} + {p_3}\delta {{\hat Z}_i} - {l_s}\\ {V_H} = - {p_4}\delta {{\hat X}_{{\rm{stat}} - j}} + {p_5}\delta {{\hat Y}_{{\rm{stat}} - j}} + {p_4}\delta {{\hat X}_i} - {p_5}\delta {{\hat Y}_i} - \delta {{\hat A}_{{\rm{az}} - j}} - {l_H}\\ {V_V} = {p_6}\delta {{\hat X}_{{\rm{stat}} - j}} + {p_7}\delta {{\hat Y}_{{\rm{stat}} - j}} + {p_8}\delta {{\hat Z}_{{\rm{stat}} - j}} - {p_6}\delta {{\hat X}_i} + {p_7}\delta {{\hat Y}_i} + {p_8}\delta {{\hat Z}_i} - {l_V} \end{array} \right. $$ (2) $$ \left\{ \begin{array}{l} {p_1} = \frac{{{X_{{\rm{stat}} - {j_0}}} + {X_{{i_0}}}}}{{{S_{{i_0}}}}}\\ {p_2} = \frac{{{Y_{{\rm{stat}} - {j_0}}} + {Y_{{i_0}}}}}{{{S_{{i_0}}}}}\\ {p_3} = \frac{{{Z_{{\rm{stat}} - {j_0}}} + {Z_{{i_0}}}}}{{{S_{{i_0}}}}}\\ {p_4} = \frac{{{Y_{{\rm{stat}} - {j_0}}} + {Y_{{i_0}}}}}{{D_{{i_0}}^2}}\\ {p_5} = \frac{{{X_{{\rm{stat}} - {j_0}}} + {X_{{i_0}}}}}{{D_{{i_0}}^2}}\\ {p_6} = \frac{{\left( {{Z_{{\rm{stat}} - {j_0}}} - {Z_{{i_0}}}} \right) \cdot \left( {{X_{{\rm{stat}} - {j_0}}} - {X_{{i_0}}}} \right)}}{{{D_{{i_0}}} \cdot \left( {D_{{i_0}}^2 + {{\left( {{Z_{{\rm{stat}} - {j_0}}} - {Z_{{i_0}}}} \right)}^2}} \right)}}\\ {p_7} = = \frac{{\left( {{Z_{{\rm{stat}} - {j_0}}} - {Z_{{i_0}}}} \right) \cdot \left( {{Y_{{\rm{stat}} - {j_0}}} - {Y_{{i_0}}}} \right)}}{{{D_{{i_0}}} \cdot \left( {D_{{i_0}}^2 + {{\left( {{Z_{{\rm{stat}} - {j_0}}} - {Z_{{i_0}}}} \right)}^2}} \right)}}\\ {p_8} = \frac{{ - D_{{i_0}}^2}}{{{D_{{i_0}}} \cdot \left( {D_{{i_0}}^2 + {{\left( {{Z_{{\rm{stat}} - {j_0}}} - {Z_{{i_0}}}} \right)}^2}} \right)}}\\ {l_{\rm{s}}} = {L_{\rm{s}}} - {S_{{i_0}}};\\ {l_H} = {L_H} + {A_0} - {A_1}\\ {l_V} = {L_V} - \arctan \frac{{{Z_i} - {Z_{{\rm{stat}} - {j_0}}}}}{{{D_{{i_0}}}}} \end{array} \right. $$ (3) 式中，Xstat-j0、Ystat-j0、Zstat-j0为在全局坐标系中第j测站坐标的近似值；Xi0、Yi0、Zi0为在全局坐标系中第i目标点的坐标近似值；A0为初始零方向近似值；A1为通过测站和目标点近似值解算的初始零方向，两者的范围都为0°~360°；Si0和Di0分别为通过测站点和目标点近似值求得的点间距及其水平投影。由于式(1)为非线性函数，求解通过最小二乘法迭代，这些近似值在不断迭代的过程中经求解的改正数调整，直到求解的改正数小于某一给定的阈值δ。
在某次测量任务中，假设共测了s个测站，目标点一共有t个，以第一站的坐标系作为全局坐标系，则需要解得参数共有(s－1)×3+t×3+(s－1)个。需要注意的是，在t个目标点中某些点是可以在多个测站中重复观测的，也就是在移站测量中的公共点。参数解为：
$$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{\hat X}} = }\\ {\left( {\underbrace {\begin{array}{*{20}{c}} {\delta {{\hat X}_{{\rm{stat}} - 2}}}&{\delta {{\hat Y}_{{\rm{stat}} - 2}}}&{\delta {{\hat Z}_{{\rm{stat}} - 2}}}& \cdots &{\delta {{\hat X}_{{\rm{stat}} - s}}}&{\delta {{\hat Y}_{{\rm{stat}} - s}}}&{\delta {{\hat Z}_{{\rm{stat}} - s}}} \end{array}}_{s - 1个测站坐标}} \right.}\\ {\underbrace {\begin{array}{*{20}{c}} {\delta {{\hat X}_1}}&{\delta {{\hat Y}_1}}&{\delta {{\hat Z}_1}}& \cdots &{\delta {{\hat X}_t}}&{\delta {{\hat Y}_t}}&{\delta {{\hat Z}_t}} \end{array}}_{t\;个目标点坐标}}\\ {{{\left. {\underbrace {\begin{array}{*{20}{c}} {\delta {{\hat a}_{{\rm{azi}} - 2}}}& \cdots &{\delta {{\hat a}_{{\rm{azi}} - s}}} \end{array}}_{s - 1个测站初始零方向参数}} \right)}^{\rm{T}}}} \end{array} $$ (4) 由于第一个测站的测站坐标和初始零方向都为0，因此该测站列立的误差方程仅有目标点的坐标参数。
观测方程的随机模型由仪器标定精度给定，测距精度σs=(a+b×Ls), a为固定误差，b为比例误差, Ls为测站点到目标点的斜距；由于一般测量目标点的角度值都是半测回测一次的结果，而标称精度一般为一测回方向的精度，所以实际的测角精度σangle为给定测角精度的$ \sqrt{2} $倍，在空间上带来的误差影响为σangle×Ls。此外，距离观测与角度观测之间是独立的，因此基于原始观测值的整体平差权阵为对角阵。对于一次观测目标点的3个原始观测值权阵为：
$$ {\mathit{\boldsymbol{P}}_1} = {\rm{diag}}\left[ {\begin{array}{*{20}{c}} {\frac{{\sigma _0^2}}{{\sigma _s^2}}}&{\frac{{\sigma _0^2}}{{{{\left( {{\sigma _{{\rm{angle}}}} \times {L_s}} \right)}^2}}}}&{\frac{{\sigma _0^2}}{{{{\left( {{\sigma _{{\rm{angle}}}} \times {L_s}} \right)}^2}}}} \end{array}} \right] $$ (5) 式中，σ0为给定的验前单位权中误差，一般设定为1 mm。验后单位权中误差由式(6)求得：
$$ \hat \sigma _0^2 = \frac{{{\mathit{\boldsymbol{V}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{V}}}}{{n - \left( {s - 1} \right) \times 4 - t \times 3}} $$ (6) 式中，n为总的观测值个数; V为式(2)中的残差向量，则求解参数$ \mathit{\boldsymbol{\hat{X}}} $的精度：
$$ {\mathit{\boldsymbol{D}}_{\hat X\hat X}} = \hat \sigma _0^2{\mathit{\boldsymbol{Q}}_{\hat X\hat X}} = \hat \sigma _0^2{\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{B}}} \right)^{ - 1}} $$ (7) 式中，矩阵B为式(2)中的系数阵。
1.2 基于坐标观测值的模型
基于坐标观测值的整体平差中，在第j测站，观测i号点的函数模型为：
$$ \left( {\begin{array}{*{20}{c}} {{L_X}}\\ {{L_Y}}\\ {{L_Z}} \end{array}} \right) + \left( {\begin{array}{*{20}{c}} {{V_X}}\\ {{V_Y}}\\ {{V_Z}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {\Delta {{\hat X}_j}}\\ {\Delta {{\hat Y}_j}}\\ {\Delta {{\hat Z}_j}} \end{array}} \right) + \left( {\begin{array}{*{20}{c}} {\cos {{\hat \gamma }_j}}&{ - \sin {{\hat \gamma }_j}}&0\\ {\sin {{\hat \gamma }_j}}&{\cos {{\hat \gamma }_j}}&0\\ 0&0&1 \end{array}} \right)\left( {\begin{array}{*{20}{c}} {{{\hat X}_i}}\\ {{{\hat Y}_i}}\\ {{{\hat Z}_i}} \end{array}} \right) $$ (8) 式中，$ \Delta {{\hat{X}}_{j}} $、$ \Delta {{\hat{Y}}_{j}} $、$ \Delta {{\hat{Z}}_{j}} $、$ {\hat{\gamma }}_{j} $为第j测站相对于全局坐标系的平移参数和绕Z轴的旋转角；$ {{\hat{X}}_{i}} $、$ {{\hat{Y}}_{i}} $、$ {{\hat{Z}}_{i}} $为i号待测点在全局坐标系中的三维坐标值。
式(8)中，求解站与站之间转换参数采用的观测值为坐标，而坐标值是由笛卡尔坐标系中原始观测的斜距、水平角和竖直角推算求得的：
$$ \left\{ \begin{array}{l} {L_X} = {L_{\rm{s}}} \cdot \cos {L_H} \cdot \sin {L_V}\\ {L_Y} = {L_{\rm{s}}} \cdot \sin {L_H} \cdot \sin {L_V}\\ {L_Z} = {L_{\rm{s}}} \cdot \cos {L_V} \end{array} \right. $$ (9) 可以看出，虽然原始观测值之间是独立的，但是坐标观测值是相关的，由误差传播可知，坐标观测值的方差阵为：
$$ \begin{array}{*{20}{c}} {{\mathit{\boldsymbol{D}}_{XYZ}} = \mathit{\boldsymbol{A}} \cdot {\rm{diag}}}\\ {\left[ {\begin{array}{*{20}{c}} {\sigma _s^2}&{{{\left( {{\sigma _{{\rm{angle}}}} \times {L_s}} \right)}^2}}&{{{\left( {{\sigma _{{\rm{angle}}}} \times {L_s}} \right)}^2}} \end{array}} \right] \cdot {\mathit{\boldsymbol{A}}^{\rm{T}}}} \end{array} $$ (10) 式中，
$$ \mathit{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} {\cos {L_H}\sin {L_V}}&{\sin {L_H}\sin {L_V}}&{\cos {L_H}\cos {L_V}}\\ {\sin {L_H}\sin {L_V}}&{\cos {L_H}\sin {L_V}}&{\sin {L_H}\cos {L_V}}\\ {\cos {L_V}}&0&{\sin {L_V}} \end{array}} \right] $$ (11) 仍设定验前单位权中误差为σ0，则对应权阵为：
$$ {\mathit{\boldsymbol{P}}_2} = {\left( {\frac{{{\mathit{\boldsymbol{D}}_{XYZ}}}}{{\sigma _0^2}}} \right)^{ - 1}} $$ (12) 与式(1)一样，式(8)也是非线性的函数，需进行线性化，可以参考大角度七参数三维基准转换模型[14]。需要注意的是，由于以第一测站的坐标系作为全局坐标系，其误差方程更为简单：
$$ \left( {\begin{array}{*{20}{c}} {{L_X}}\\ {{L_Y}}\\ {{L_Z}} \end{array}} \right) + \left( {\begin{array}{*{20}{c}} {{V_X}}\\ {{V_Y}}\\ {{V_Z}} \end{array}} \right) = \left( {\begin{array}{*{20}{c}} {{{\hat X}_i}}\\ {{{\hat Y}_i}}\\ {{{\hat Z}_i}} \end{array}} \right) $$ (13) 传统移站测量方法中，站与站之间通过观测公共点将坐标系统一并传递下去。实际公共点在两套坐标系的观测值中都含有误差，而实际上一般仅采用最小二乘方法处理，忽略了系数阵的误差，因此随着移站的增多，误差累积越来越大。
对比§1.1和§1.2可以看出，两者采用不同的函数模型，但最终都能通过整体平差求解出待测点在第一测站坐标系下的坐标，而且没有传统移站测量中出现的误差累积的问题；此外，两者的系数阵中都不存在含观测误差的观测值，因此通过传统最小二乘算法和准确的随机模型即可获得准确的参数估计，而避免采用整体最小二乘等较复杂的方法来提高解算精度。
2 算例与分析
为了对比研究两种数据融合模型的优劣，笔者利用索佳Net05高精度全站仪，围绕某旗杆墩架设了4个测站，并在墩上的4个面观测了24个贴有反射片的待测点，在每个测站上都可以观测到两个墩面共12个目标点，测站1观测1~12号点，测站2观测7~18号点，测站3观测13~24号点，测站4观测19~24号点和1~6号点，如图 1所示。
对于索佳Net05，当反射物为反射片时，测距标称精度为σs = 0.5+1× ppm mm，测角精度为0.5″，则半测回某方向的精度为σangle=0.5×$ \sqrt{2} $≈0.7″，采用§1.1中的融合方法及式(7)的精度评定公式，以第一测站坐标系作为全局坐标系，得到该坐标系下除测站1的其他3个测站的测站坐标及初始零方向参数和24个待测点的坐标，解算结果见表 1、表 2。同理可以得到按照§1.2的融合方法的解算结果，见表 3、表 4。
表 1 基于原始观测值的整体平差结果Table 1. Result of Overall Adjustment Based on Original Data类型 编号 X/m Y/m Z/m $ \hat{\sigma } $X/mm $ \hat{\sigma } $Y/mm $ \hat{\sigma } $Z/mm 待测点 1 -1.847 0 -11.186 5 0.764 1 0.2 0.2 0.1 2 -1.483 9 -10.830 7 0.779 2 0.2 0.1 0.1 3 -1.064 8 -10.431 0 0.796 8 0.1 0.2 0.1 4 -1.872 2 -11.153 5 0.230 3 0.2 0.2 0.1 5 -1.521 4 -10.806 9 0.196 2 0.2 0.1 0.1 6 -1.098 7 -10.396 3 0.173 7 0.1 0.2 0.1 7 -0.786 9 -10.499 8 0.823 1 0.1 0.2 0.1 8 -0.454 7 -10.842 2 0.828 7 0.2 0.2 0.2 9 -0.108 5 -11.190 8 0.835 5 0.2 0.2 0.2 10 -0.768 9 -10.449 1 0.260 0 0.1 0.2 0.1 11 -0.417 1 -10.811 3 0.272 4 0.2 0.2 0.1 12 -0.082 2 -11.157 0 0.318 2 0.2 0.2 0.2 13 -0.093 4 -11.504 2 0.759 5 0.2 0.2 0.2 14 -0.439 7 -11.837 9 0.779 4 0.3 0.2 0.2 15 -0.826 0 -12.211 0 0.767 9 0.3 0.2 0.2 16 -0.054 3 -11.531 9 0.179 0 0.2 0.2 0.2 17 -0.442 7 -11.904 3 0.191 9 0.3 0.2 0.2 18 -0.798 3 -12.253 4 0.168 8 0.3 0.2 0.2 19 -1.116 3 -12.211 5 0.733 2 0.3 0.2 0.1 20 -1.468 5 -11.854 8 0.745 2 0.3 0.2 0.1 21 -1.845 4 -11.474 2 0.735 1 0.2 0.2 0.1 22 -1.150 9 -12.234 0 0.235 4 0.3 0.2 0.1 23 -1.493 8 -11.888 1 0.222 7 0.3 0.2 0.1 24 -1.877 1 -11.502 8 0.192 6 0.2 0.2 0.1 2 11.670 7 -12.423 4 -0.018 8 0.3 2.5 0.1 测站 3 -1.140 5 -23.007 8 -0.013 3 2.1 0.1 0.1 4 -11.095 8 -11.317 5 0.039 3 0.1 1.8 0.1 表 2 测站其他参数平差结果Table 2. Result of Other Parament编号 初始零方向
/rad解算精度/(″) 验后单位权中误差/mm 2 -2.922 588 41.6 0.5 3 0.741 673 37.4 4 5.704 257 37.9 表 3 基于坐标观测值的整体平差结果Table 3. Result of Overall Adjustment Based on Coordinate Data类型 编号 X/m Y/m Z/m $ \hat{\sigma } $X/mm $ \hat{\sigma } $Y/mm $ \hat{\sigma } $Z/mm 待测点 1 -1.847 1 -11.186 5 0.764 1 0.1 0.4 0.0 2 -1.483 8 -10.831 0 0.779 3 0.1 0.3 0.1 3 -1.064 9 -10.431 4 0.796 8 0.1 0.3 0.1 4 -1.872 2 -11.153 6 0.230 3 0.1 0.4 0.0 5 -1.521 4 -10.807 0 0.196 2 0.1 0.3 0.0 6 -1.098 7 -10.396 7 0.173 7 0.1 0.3 0.0 7 -0.786 8 -10.499 2 0.823 0 0.1 0.3 0.1 8 -0.455 0 -10.841 7 0.828 7 0.1 0.3 0.1 9 -0.108 4 -11.190 4 0.835 5 0.1 0.4 0.1 10 -0.768 7 -10.448 7 0.260 0 0.1 0.3 0.1 11 -0.417 4 -10.811 0 0.272 4 0.1 0.3 0.1 12 -0.082 0 -11.156 6 0.318 2 0.1 0.4 0.1 13 -0.093 3 -11.504 0 0.759 5 0.5 0.4 0.1 14 -0.439 8 -11.837 5 0.779 5 0.5 0.3 0.1 15 -0.825 9 -12.210 5 0.767 9 0.6 0.3 0.1 16 -0.054 5 -11.531 7 0.179 1 0.5 0.4 0.1 17 -0.442 8 -11.903 8 0.191 9 0.5 0.3 0.1 18 -0.798 1 -12.252 8 0.168 8 0.6 0.3 0.1 19 -1.116 1 -12.211 9 0.733 2 0.6 0.3 0.1 20 -1.468 5 -11.855 1 0.745 2 0.5 0.3 0.1 21 -1.845 6 -11.474 2 0.735 1 0.5 0.4 0.1 22 -1.150 8 -12.234 4 0.235 4 0.6 0.3 0.1 23 -1.494 0 -11.888 4 0.222 7 0.5 0.3 0.1 24 -1.877 4 -11.502 9 0.192 6 0.5 0.4 0.1 2 8.697 1 -14.663 0 0.018 8 5.2 1.5 0.0 测站 3 16.385 9 16.191 1 0.013 3 3.7 4.0 0.0 4 3.090 5 15.541 9 -0.039 3 4.4 2.2 0.0 表 4 测站其他参数平差结果Table 4. Result of Other Parament编号 绕Z轴旋转角/rad 解算精度/(″) 验后单位权中误差/mm 2 2.922 953 100.6 3 -0.742 009 98.3 1.7 4 0.579 411 91.0 从表 1、表 2可以看出，待测点的点位解算精度可达0.3 mm左右，且精度分布均匀。从表 3、表 4可以看出，点位解算精度基本在0.5 mm左右，且测站的定位参数解算精度也不如表 1、表 2的结果高。虽然从数学角度看，§1.2的融合方法中，观测值是由§1.1的原始观测值等价推导得到的，随机模型也是从§1.1的误差模型通过误差传播得来的，但是融合结果却不是完全等价的。
按照传统移站测量的思路，采用§1.1的函数模型进行分站平差，即第一测站观测的1~12点都为全局坐标系的坐标；第二测站观测7~12点来解算测站参数，根据解算得到的测站参数和13~18点的原始观测值解算出13~18点的全局坐标系的坐标；同理，第三测站也可求解出19~24点的全局坐标系坐标。此方法称为基于原始观测值的移站测量法。同理采用§1.2的函数模型进行传统移站测量的方法称为基于坐标观测值的移站测量法。
通过不同的方法得到24个待测点的全局坐标，有如下4种方案：①基于原始观测值的整体平差法；②基于坐标观测值的整体平差法；③基于原始观测值的移站测量法；④基于坐标观测值的移站测量法。分别对比表 1、表 2和表 3、表 4，可以看出方案①的参数解算精度较高。因此，以方案①解算的待测点坐标为基准，可以得到方案①与方案②、方案①与方案③以及方案②与方案④的坐标点位偏差，如图 2所示。
可以看出，方案①与方案②的各点点位偏差较均匀，范围在0.1~0.6 mm之间；而在传统的MSM中，无论是基于原始观测值或是坐标观测值，其解算的点位偏差随着搬站的增加而递增。
图 2中MSM采用6个公共点进行搬站，在经过2次搬站后测得的19~24号点与整体平差的点位偏差可达3~4 mm左右，而实际测量中一般不会采用这么多公共点进行移站测量。如若在24个待测点中选取12个点，如选取1、3、5、…、19、21、23号待测点，每两个测站之间只有3个公共点，即减少了多余观测数，可以得到两种OA结果与对应传统MSM的点位偏差，如图 3所示。
![]()
图 3 3个公共点时整体平差与移站测量的点位偏差Figure 3. Coordinate Differences Between OA and MSM with 3 Common Points从图 2和图 3可以看出，基于原始观测值的OA和MSM受多余观测值数量的影响较小，而基于坐标观测值的OA和MSM受多余观测数的影响非常明显。当搬站公共点仅有3个的时候，经两次搬站后测量的目标点与整体平差点位偏差可达将近3 cm，因此基于坐标观测值的MSM的点位偏差更加难以控制。
为了进一步比较当随机模型与实测数据存在偏差，即存在随机模型误差时，方案①和方案②的解算稳定性，在所有斜距原始观测值分别添加中误差为1 mm和5 mm的白噪声(标称精度仅为0.5 mm)，而随机模型仍然采用仪器的标称精度计算，可以得到两种方案在相同额外白噪声影响下解算的坐标偏差，如图 4所示。
![]()
图 4 添加不同白噪声对两种解算点位偏差影响Figure 4. Influence on the Coordinate Deviation when Adding Different White Noise方案①在外加1 mm斜距误差后解算的24个待测点位与不加误差时偏差在0.2~1.5 mm，在外加5 mm斜距误差后解算的24个待测点位与不加误差时偏差在1.1~5.9 mm；同样的条件，方案②的偏差分别在0~0.8 mm和1.1~2.9 mm。从以上偏差的范围和稳定性都可以看出，当随机模型含有误差时，方案②参数解算稳定性更高。
3 结语
1) 当随机模型的误差不明显，即仪器实际测量精度和标称精度相差不大时，采用基于原始观测值的整体平差融合算法来解算待定点坐标的方法精度最高。在传统的移站测量解算中，采用基于原始观测值的方法也比基于坐标观测值的结果要稳定。
2) 当多余观测数较少，即公共点个数较少时，若采用传统移站测量法解算点位坐标，不宜采用基于坐标观测值的融合方法，采用基于原始观测值的融合方法会更可靠。
3) 当随机模型含有明显误差时，基于坐标观测值的整体平差法的融合结果要比基于原始观测值的整体平差融合结果稳定可靠。
-
![]()
图 7 18个验潮站记录的6个断层单元引发的海啸首波波高和到时
Figure 7. Wave Heights of Head-Wave Recorded at 18 Tidal Stations
![]()
图 8 18个验潮站记录的6个断层单元引发的海啸首波到时
Figure 8. Arrival Time of Head-Wave Recorded at 18 Tidal Stations
![]()
图 9 18个验潮站记录到6个断层单元引发海啸的最大海啸波
Figure 9. Wave Heights of Maximum-Wave Recorded at 18 Tidal Stations
断层 经度 纬度 长度/km 宽度/km 位移/m 深度/km 走向角 倾角 滑移角 1 120.5°E 20.2°N 160 200 20 40 354° 10° 90° 2 119.8°E 18.7°N 180 200 20 40 22° 20° 90° 3 119.3°E 17.0°N 240 200 20 40 2° 28° 90° 4 119.2°E 15.1°N 170 200 20 40 356° 20° 90° 5 119.6°E 13.7°N 140 200 20 40 344° 22° 90° 6 120.5°E 12.9°N 100 200 20 40 331° 26° 90° 
下载: 导出CSV
表 2 模拟区域二级网格设置参数
Table 2 Parameters Configured for Two-level Grids in Simulated Domain
网格 分辨率 范围 行列数 控制方程 曼宁系数/(m-1/3·s) Grid1 2′ 10°N~30°N, 105°E~125°E 751×601 线性 缺省 Grid2 Layer21
Layer220.5′
0.5′17N°~22°N, 106.5°E~111°E
21°N~24.5°N, 111.5°E~118.5°E600×540
420×840非线性
非线性0.026
0.026
下载: 导出CSV
表 3 今村-饭田海啸强度分级表
Table 3 Imamura-Iida Intensity Scale of Tsunami
海啸强度/m 波高范围 灾害程度说明 海啸能量/J -1 < 50 cm 验潮站方能观测到海啸痕迹，通常无灾害 2.5×1018 0 1 m左右 渔船、水产设施有损坏；100 km长海岸范围内波高50~80 cm，可能造成小灾害 1.3×1019 1 2~3 m 海滨低地房屋有损坏，船舶流失；200 km长海岸范围内波高约1 m 6.4×1019 2 4~6 m 多数房屋浸水、流失，人畜溺死；300 km长海岸范围内波高约1.5 m 3.2×1020 3 10 m左右 局部达到20 m；破坏区长达400 km，波高约2.4 m，造成显著破坏 1.6×1021 4 20 m左右 局部超过30 m；500 km海岸范围内波高约4 m，造成巨大破坏和灾害损失 8.0×1021
下载: 导出CSV
-
[1] 王锋, 刘昌森, 章振铨.中国古籍中的地震海啸记录[J].中国地震, 2005, 21(3):437-443 http://www.doc88.com/p-988346240743.html Wang Feng, Liu Changsen, Zhang Zhenquan. Earthquake Tsunami Record in Chinese Ancient Books[J]. Earthquake Research in China, 2005, 21(3):437-443 http://www.doc88.com/p-988346240743.html
[2] 张俊香, 李平日, 黄光庆.新奥尔良飓风灾难与华南沿海台风暴潮[J].热带地理, 2006, 26(3):218-222 http://kns.cnki.net/KCMS/detail/detail.aspx?filename=rddd200603004&dbname=CJFD&dbcode=CJFQ Zhang Junxiang, Li Pingri, Huang Guangqing. Hurricane Disasters in New Orleans of USA and Typhoon Storm Surge[J]. Tropical Geography, 2006, 26(3):218-222 http://kns.cnki.net/KCMS/detail/detail.aspx?filename=rddd200603004&dbname=CJFD&dbcode=CJFQ
[3] 杨马陵, 魏柏林.南海海域地震海啸潜在危险的探析[J].灾害学, 2005, 20(3):41-47 http://kns.cnki.net/KCMS/detail/detail.aspx?filename=zhxu200503008&dbname=CJFD&dbcode=CJFQ Yang Maling, Wei Bailin. The Potential Seismic Tsunami Risk in South China Sea and Its Surrounding Region[J]. Journal of Catastrophology, 2005, 20(3):41-47 http://kns.cnki.net/KCMS/detail/detail.aspx?filename=zhxu200503008&dbname=CJFD&dbcode=CJFQ
[4] Wu T R, Huang H C. Modeling Tsunami Hazards from Manila Trench to Taiwan[J]. Journal of Asian Earth Sciences, 2009, 36(1):21-28 doi: 10.1016/j.jseaes.2008.12.006
[5] 王培涛, 于福江, 赵联大, 等. 2011年3月11日日本地震海啸越洋传播及对中国影响的数值分析[J].地球物理学报, 2012, 55(9):3088-3096 doi: 10.6038/j.issn.0001-5733.2012.09.026 Wang Peitao, Yu Fujiang, Zhao Lianda, et al. Numerical Analysis of Tsunami Propagating Generated by the Japan Mw9.0 Earthquake on Mar.11 in 2011 and its Impact on China Coasts[J]. Chinese Journal of Geophysics, 2012, 55(9):3088-3096 doi: 10.6038/j.issn.0001-5733.2012.09.026
[6] UNESCO-IOC. IUGG/IOC Time Project:Numerical Method of Tsunami Simulation with the Leap-frog Scheme[OL]. http://www.jodc.go.jp/info/ioc_doc/Manual/122367eb.pdf, 1997
[7] Titov V V, Gonzalez F I. Implementation and Testing of the Method of Splitting Tsunami (MOST) Model[OL]. http://docs.lib.noa.gov/noa_documents/OAR/PMEL/TM_ERL_PMEL/TM_ERL_PMEL_112.pdf, 1997
[8] Liu P L F. et al. Computer Programs for Tsunami Propagation and Inundation[OL]. http://ceeserver.cee.cornell.edu/pll-group/doc/COMCOT_User_Manual_v1_7.pdf, 1998
[9] George D L. Introduction Tsunami Claw User's Guide[OL]. http://faculty.washington.edu/rj1/pubs/icm06/TsunamiClawDoc.pdf, 2006
[10] 于福江, 叶琳, 王喜年.1994年发生在台湾海峡的一次地震海啸的数值模拟[J].海洋学报, 2001, 23(6):32-39 http://kns.cnki.net/KCMS/detail/detail.aspx?filename=seac200106003&dbname=CJFD&dbcode=CJFQ Yu Fujiang, Ye Lin, Wang Xinian. The Simulation of Tsunami Happened in the Taiwan Strait in 1994[J]. Acta Oceanlogica Sinica, 2001, 23(6):32-39 http://kns.cnki.net/KCMS/detail/detail.aspx?filename=seac200106003&dbname=CJFD&dbcode=CJFQ
[11] 李平日, 方国祥, 黄光庆.海平面上升对珠江三角洲经济建设的可能影响及对策[J].地理学报, 1993, 48(6):527-534 http://www.cqvip.com/QK/90059X/199306/1008802.html Li Pingri, Fang Guoxiang, Huang Guangqing. Impacts on Sea Level Rising on the Economic Development of Zhujiang Delta and Countermeasures[J]. Acta Geographica Sinica, 1993, 48(6):527-534 http://www.cqvip.com/QK/90059X/199306/1008802.html
[12] 陈志豪, 李家彪, 吴自银, 等.马尼拉海沟几何形态特征的构造演化意义[J].海洋地质与第四纪, 2009, 29(2):59-65 http://www.cqvip.com/QK/96122X/200902/1000652960.html Chen Zhihao, Li Jiabiao, Wu Ziyin, et al. Tectonic Evolution Implication of Geometry Shape Characteristics for Manila Trench[J]. Marine Geology & Quaternary Geology, 2009, 29(2):59-65 http://www.cqvip.com/QK/96122X/200902/1000652960.html
[13] 尚继宏. 马尼拉海沟中北段俯冲带特征对比及区域构造动力学研究[D]. 北京: 中国科学院大学, 2008 http://cdmd.cnki.com.cn/Article/CDMD-80068-2008117177.htm Shang Jihong. Tectonic Dynamics Research and Subducting Characteristics Comparison Between Middle and Northern Part Manila Subducting Belt[D]. Beijing:University of Chinese Academy of Sciences, 2008 http://cdmd.cnki.com.cn/Article/CDMD-80068-2008117177.htm
[14] 臧绍先, 宁杰远.西太平洋俯冲带的研究及其动力学意义[J].地球物理学报, 1996, 39(2):188-201 http://kns.cnki.net/KCMS/detail/detail.aspx?filename=dqwx199602005&dbname=CJFD&dbcode=CJFQ Zang Shaoxian, Ning Jieyuan. Study on the Subduction Zone in Western Pacific and its Implication for the Geodynamics[J]. Acta Geophysica Sinica, 1996, 39(2):188-201 http://kns.cnki.net/KCMS/detail/detail.aspx?filename=dqwx199602005&dbname=CJFD&dbcode=CJFQ
[15] 潘文亮, 王盛安. COMCOT数值模式的介绍和应用[J].海洋预报, 2009, 26(3):45-52 doi: 10.11737/j.issn.1003-0239.2009.03.006 Pan Wenliang, Wang Sheng'an. Introduction and Application of COMCOT Model[J]. Marine Forecasts, 2009, 26(3):45-52 doi: 10.11737/j.issn.1003-0239.2009.03.006
[16] Liu P L F, Wang X M, Andrew J. Tsunami Hazard and Early Warning System in South China Sea[J]. Journal of Asian Earth Sciences, 2009, 36(1):2-12 doi: 10.1016/j.jseaes.2008.12.010
[17] Nguyen P H, Bui Q C, Vu P H, et al. Scenario-Based Tsunami Hazard Assessment for the Coast of Vietnam from the Manila Trench Source[J]. Physics of Earth and Planetary Interiors, 2014, 236:95-108 doi: 10.1016/j.pepi.2014.07.003
[18] Megawati K, Shaw F, Sieh K, et al. Tsunami Hazard from the Subduction Megathrust of the South China Sea. Part Ⅰ:Source Characterization and the Resulting Tsunami[J]. Journal of Asian Earth Sciences, 2009, 36(1):13-20 doi: 10.1016/j.jseaes.2008.11.012
[19] 陈传绪, 吴时国, 赵昌垒.马尼拉海沟北段俯冲带输入板块的不均一性[J].地球物理学报, 2014, 57(12):4063-4073 doi: 10.6038/cjg20141218 Chen Chuanxu, Wu Shiguo, Zhao Changlei. Incoming Plate Variation Along the Northern Manila Trench[J]. Chinese Journal of Geophysics, 2014, 57(12):4063-4073 doi: 10.6038/cjg20141218
[20] 王培涛, 于福江, 范婷婷, 等.海啸波传播的线性核非线性特征及近海陆架效应影响的数值研究[J].海洋学报, 2014, 36(5):18-29 doi: 10.3969/j.issn.0253-4193.2014.05.003 Wang Peitao, Yu Fujiang, Fan Tingting, et al. Numerical Study on the Linear/nonlinear Characteristics and the Impacts of Continental Shelf Effects of the Tsunami Waves Propagating[J]. Acta Oceanologica Sinica (in Chinese), 2014, 36(5):18-29 doi: 10.3969/j.issn.0253-4193.2014.05.003
[21] Iida K. Magnitude, Energy and Generation Mechanisms of Tsunamis and a Catalogue of Earthquakes Associated with Tsunamis[C]. 10th Pacific Science Congress, Honolulu, 1963 http://www.researchgate.net/publication/284501461_Magnitude_energy_and_generation_mechanisms_of_tsunamis_and_a_catalogue_of_earthquakes_associated_with_tsunamis
-
期刊类型引用(4)
1. 朱杰，郑加柱，陈红华，杨静，胡平昌，陆敏燕. 结合POI数据的南京市商业中心识别与集聚特征研究. 现代测绘. 2022(06): 34-39 . 
百度学术
2. 金澄，安晓亚，陈占龙，马啸川. 矢量居民地多边形多级图划分聚类方法. 武汉大学学报(信息科学版). 2021(01): 19-29 . 百度学术
3. 张铭龙，何贞铭. 基于因子分析法的城市商业中心抽取研究. 地理空间信息. 2021(08): 58-60+64+5 . 百度学术
4. 李卫东，张铭龙，段金龙. 基于POI数据的南京市空间格局定量研究. 世界地理研究. 2020(02): 317-326 . 百度学术
其他类型引用(3)
计量
- 文章访问数:
- HTML全文浏览量:
- PDF下载量:
- 被引次数: 7
