地形图的绘制
海底地形图的绘制模型
华东师范大学 江** 常* 黄*
海底地形图的绘制模型
摘要:
该“海底地形图的绘制模型”首先通过各个角度的假设,建立一个较为理想的海洋测绘环境;考虑该地形图能够表示为在一个三维笛卡尔坐标内的函数,而“插值双三次样条函数”能够较好地对图形进行模拟,我们对函数的形式做出一些改进,以满足海底地形图的实际要求;为了使绘制的图形更贴近与现实,我们将待测的海洋区域以网格方式划分为若干小区域,通过对每个小区域的地形图的绘制,并将所有的小区域内的图形连接,最后生成整个待测海洋区域的海底地形图。
一. 问题重述
海洋测绘船利用声纳绘制海底的地形图。测绘船上的声纳向海底发射声脉冲,随后接收从海底反射的脉冲。发射的范围为与指向海底的铅垂线夹角从2°—30°之间。船只以2米/秒的速度行进,声脉冲在海水中传播的速度约为1500米/秒。
试建立绘制海底地形图的数学模型,并对绘制海底地形图的方法提出具体建议。
二. 模型假设
1. 造成海面的不平坦因素有很多,比如潮汐,风浪等,假设这些因素的影响在
测量过程中间可以忽略。
2. 讨论区域的海底曲面是光滑的;更确切些,可以认为曲面的一阶、二阶导数
是连续的。因为我们可以认为讨论区域为浅水海域,由于长期的海水水流作用,形成的是以砾石或沙为主要组成部分的海底,不存在珊瑚礁、水底峡谷、山脊等不可预料的突变地形。为了方便分析误差,假设海底临近的两点之间的与水平面的夹角在15°±15°范围内,且反射回接收装置的声波线与铅垂线之间的夹角在16°±14°范围内。
3. 声波在传播途中受海水介质不均匀分布和海面、海底的影响和制约,会产生
折射、散射、反射和干涉,会产生声线弯曲、信号起伏和畸变,造成传播途径的改变,以及出现声阴区,严重影响声纳的作用距离和测量精度[1]。假设海水介质非常均匀,造成声波在海水中传播方向不稳定的因素可以忽略。 4. 由于海水的浓度,温度,密度,压力等等因素起伏变化不定,声音在其中的
传播速度肯定也在起伏变化,由于这种变化没有固定的规律可循,假设声波在海水中以一个稳定的速度传播。
5. 假设回波接受装置接受到的回波是可辨别的,即能够确认是发射装置刚刚发
射的声波的回波,不会出现辨别不出回波信号的问题。
6. 假设声波在到达海底前引起的回波都是不明显的,即假设收到的回波信号是
遇到海底地面返回的。
7. 假设用于测量的声纳设备本身是精准的,至少在测量过程中不会因为声纳设
备本身的原因而导致采集的数据出错。
8. 由于测量条件的限制 ,假设测量误差在小于10%的范围内都是允许的。
三. 问题分析 相关技术背景: 1 声纳传播方式
题目中要求建立利用声纳绘制海底地形图的数学模型,所以首先有必要了解一下声纳的工作原理。
声波是目前已知的唯一能在海水中远程传播的波,声纳(SONAR )是利
用声波在水下的传播特性,通过电声转换装置和信号处理,完成对水下物体进行探测和定位识别的方法及所用设备的总称[2]。当我们通过声纳向水中发射声波脉冲时,声波能够以一定的速度在水中沿直线传播,并能从水底反射回波。为了得到某个目标位置的深度值,我们可以根据声波的这一特性,通过计算发射声波脉冲和回波到达的时间差来求得目标位置的深度。这是利用声纳测量海底地形的基本原理。
为了更好的理解利用声纳测量海底地形的原理和过程,我们需要了解一下声波的传播过程。波在均匀各向同性介质中传播时,波面及波前的形状不变,波线(波的传播方向成为波线)也保持为直线,沿途不会改变波的传播方向。但是当遇到障碍物或是从一种介质传播到另一种介质中时,波面的形状和波线方向将发生改变。通过惠更斯原理可以知道,在介质中,波传到的各点不论在同一波前或不同波前上,都可看作是发射子波的波源[3]。因此,测绘船发射的声脉冲是以一个球形波阵面的形式向海底传播的,而且在遇到海底障碍返回时,又是以球形波阵面的形式返回。
2 曲面造型方法
海底曲面可由插值双三次样条函数去模拟,前提是测出一些已知点的深度。关
于插值双三次样条函数说明如下:[5]
P150
定义:设在uw 平面上的矩形区域R:[a,b]×[c,d]上给定一个矩形网络分割△=△u ×△w, 其中
△u : a =u0
=b,
△w: c =w0
=d.
凡在R 上满足下述两条件的一个函数x(u,w)称为双三次样条函数。
●在每个子矩形Rij :[Ui -1,U i ] ×[Wi -1,W i ](i=1,2„n;j=1,2,„,m )上,x(u,w)关于u 和w 都是三次多项式函数,即 x(u,w) =
e, f =0
∑
3
ij B ef (u -u i -1) e (w -w j -1) f ; (1.1)
∂∂+βx (u , w )
● 在整个R 上,函数x(u,w)的偏导数,(α, β=0,1,2)是连续的,或
∂u αu β
4
简记为x(u,w) ∈C 2(R ) 。
此外,在给定一数组{x ij }(i=0,1,2„n; j=0,1,2,„,m )之后,假如x(u,w)在节点处再满足插值条件,即: ● x(u,w)=
x ij (i=0,1,2„n; j=0,1,2,„,m )(1.2)
则称x(u,w)为插值双三次样条函数。
给出以下记号:
(1) 矩形R 的四条边界上的所有节点处的一阶法向偏导数 (1.4)
'
p αj =x U (u α, w j ) (j=0,1,2,„,m ; α=0,n).
' q i β=x w (u i , w β) (i=0,1,2,„,n; β=0,m)
(2) 矩形R 的四个顶点处的二阶混合偏导数 S αβ
由存在唯一性定理[5] P152可知:给定uw 平面上的矩形R 、它的一种矩形分割△和(m+3)(n+3)个常数
' '
=x uw (u α, w β) (α=0,n;β =0,m). (1.5)
x ij , p αj ,, q i β, S αβ (i=0,1,…n; j=0,1,…,m; α=0,n;
β=0,m)
之后,恰好存在一个插值双三次样条函数x(u,w),使它是以(1.2)为插值条件,以(1.4)
和(1.5)为边界条件的。
当我们应用上述定理到单个子矩形Rij 时,便可知道:只要给定了Rij 的四个定点上的函数值,两个方向的一阶偏导数以及二阶混合偏导数等16个数据,也就是给定了四阶方阵
i -1, j
⎡⎤
x i , j ⎢x i , j -1⎥q i , j -1q ij
⎢⎥
[C]ij =⎢p p S S ⎥
i -1, j -1i -1, j i -1, j -1i -1, j ⎢⎥
p ij S i , j -1S ij ⎣p i , j -1⎦
x i -1, j -1x i -1, j
q i -1, j -1
q
那么,Rij 上的双三次函数x(u,w)就被唯一地决定了。矩阵[C]ij 中各元素的意义
分别是
' '
x i , j =x (u i , w j ); p ij = xu (u i , w j ), q ij = xw (u i , w j ), S ij =x uw (u i , w j ) .
' '
且易得出Rij 上的双三次函数一定可唯一地写成
x(u,w)=[(u-ui -1)][A(hi )][C]ij [A(gj )][(w-wj -1)], (u,w) ∈Rij , ( 1.7 )
式中记号[(u-u i -1)]代表关于变元(u-u i -1) 的1×4阶矩阵:
T
[(u-ui -1)]=[(u-ui -1) 3, (u-ui -1) 2, (u-ui -1),1],
[A(h i )]代表4×4阶矩阵:
h i =u i -u i -1(i=1,2,...,n ) g j =w j -w j -1(j=1, 2, …,m)
⎡2
⎢h 3⎢i
3⎢A (h i ) =-2⎢h i ⎢0⎢⎢⎣1
2
-3h i 3h i 2001h i 22-h i 10
1⎤h i 2⎥⎥1⎥-h i ⎥ ⎥0⎥0⎥⎦
符号T 代表矩阵的转置。
这样一来,双三次样条函数x(u,w)在子矩形Rij 上的表示式(1.7)完全决定于矩阵[C]ij 。[C]ij 的左上角四个元素由插值条件(1.2)直接给出,其余诸p,q,S 共12个元素则需要按照边界条件(1.4)和(1.5),并求解连续性方程获得。
首先,固定w=wj , 则函数x(u, wj ) ∈S (u ,△u ). 由三次样条函数的m 次连续性方程[5] P34,我们有
λi p i -1, j +2p ij +μi p i +1, j =3[λi
x ij -x i -1, j
h i
+μi
x i +1, j -x ij
h i +1
]
(i=1,2,„n-1), (1.9)
式中λi =
h i h i +1
μ
h i +h i +1 , i =h i +h i +1
.
在边界条件 (1.4) 式中,(1.9) 的边界切向 p 0j 和p nj 是已知的,由此解出 p ij (i=0,1,„n) 。然后取j=0,1,„m ,我们便求得R 上所有节点(u i ,w j )处沿u 方向的一阶偏导数
p ij (i=0,1,„,n;j=0,1,„m).
对称地,当固定u= ui 时,则函数x(ui ,w) ∈S (u ,△w ). 于是得出相应的m 连续性方程
λq
*
j i , j -1+2q ij +μq *
j i , j +1
=3[λ
*j
x ij -x i , j -1
g j
+μ
*j
x i , j +1-x ij
g j +1
]
(j=1,2,„ m-1), (1.10) 式中记
λ
*j =
g j +1g j +g j +1
,
μ
*j =
g j g j +g j +1
。其边界切向q i 0, 和q im 也是从
(1.4)获得的,由此我们解出q ij (j=0,1,„m). 这样,R 上所有节点(u i ,w j )处沿w 方向的一阶偏导数q ij (i=0,1,„m )就全部确定了。
为了计算所有节点处的二阶混合偏导数S ij ,我们首先固定u= ui , 则函数
'
(u i ,w ) ∈S (w ,△w )。它在各节点w j 处的函数值x U (u i ,w j )= x U
'
p ij
已从(1.9)求得。因此我们得到关于变元w 的三次样条函数的m 连续性方程
λS i , j -1+2S ij +μS *j *
j i , j +1
=3[λ
*j
p ij -p i , j -1
g j
+μ
*j
p i , j +1-p ij
g j +1
]
(j=1,2,„,m-1) (1.11)
假如还给定了连续性方程(1.11)的边界切向S i0和S im 的话,那末就解得S ij (j=0,1,„m). 这样R 上所有节点(u i ,w j )处的二阶混合偏导数(i=0,1,„,n;j=0,1,„m) 也就全部确定了。
其实,(1.11)的边界切向S i 0和S im 是这样求得的:沿R 的两条边界w=w0与
S ij
w=wm , 函数 (u, w0) 与x W (u, wm ) ∈S (u; △u ). 它们在各节点ui 处的函数
'
' 值x W (ui , w0)=qi 0和x W ( ui , wm )=qim 可从(1.10)获得。于是我们有关于
'
变元u 的m 连续性方程
λi S i -1, β+2S i β+μi S i +1, β=3[λi
(i=1,2,„,n-1), β
q i β-q i -1, β
h i
+μi
q i +1, β-q i β
h i +1
]
=0, m . (1.12)
它们的边界条件由(1.5)式提供。由此解出S i 0和S im (i=0,1,„n). 至此,矩阵(1.6)中所有各元均已求得,因此R 上的双三次样条函数x(u,m)也就被确定了,它们在各子矩形R ij 上的表达式就是(1.7)。
余下的问题是:参数平面(u,w)上的定义域R 及其分割△应该如何确定呢?通常的方法是,适当固定某行j 0, 而取u 方向的节点间距h i ,(i=1,2,...,n).再适当固定某列i 0,取w 方向的节点间距g j (j=1,2,...,m).若令定义域R=[0,H]×[0,G],这里的
H=
n
∑h , G=∑g
i i =1
m
j
.
j =1
这样,我们就确定了定义域R 及其分割△。这种取法实际上是构造了两个方向上的累加弦长参数分割。累加弦长参数三次样条曲线具有良好的光顺性,因此现在构造的样条曲面的网络曲线将是光顺的。运用这种方法的前提,就是型值点x(u,w)之间的间隔比较均匀。倘若一侧间距很大,而另一侧间距很小,象扇形分布那样,我们将不会收到好的效果。前所讲的所谓适当固定某行j 0,是指这一行上各点x(i, j0) 的间隔大致是一个平均值。
因此我们作以下假设:
讨论区域的海底曲面是光滑的;更确切些,可以认为曲面的一阶、二阶导数是连续的。因为我们可以认为讨论区域为浅水海域,由于长期的海水水流作用,形成的是以砾石或沙为主要组成部分的海底,不存在珊瑚礁、水底峡谷、山脊等不可预料的突变地形。
至此,所有问题都解决了;由测绘船测出规则点的深度值,则利用以上插值双三次样条函数的求解(规定边界上的一阶偏导数和二阶混合导数为0)可得出海底曲面;如果测的是随机点的深度值,则可由它来估算规则点的深度值,再求解;方法如下:
由随机点的深度值推测出分割节点上的水深值,我们采用在地球科学用得较多的IDW (inverse distance weighting)方法。
IDW 方法是地球科学中常用的计算机插值技术,最早可见于60年代的地质轮廓描述和矿产推测的有关文献中。IDW 的基本思想是估计已知数据点对推测点的
影响,假定推测点的值为周围已知点的值的线性组合。有前假定,每个已知点对推测点的影响反比于它们的距离的某种方式,这儿我们用距离的p 次方。设已知点(u i , w i ) 的值为X i , 则 X=其中C i
∑C X
i i =1
N i =1
N
i
,
=K /D i
p
,1/K =
p
1/D ∑i ,D i 为点(u,w)和点(u i , w i )之间的距离,
这儿K 起着规范化的作用。
对于IDW 推测出的未知点的值,显然在已知点的值的范围之内,因此,只有靠近已知点的附近才会出现极值情况。上式中的p 反映了未知点对已知点的偏倚情况,p 越大,表示未知点的值受近距离点的值的影响越大(相对于远距离点)。 例如,我们取p=3,依照上式就可以根据已知节点深度值计算出所有未知格点的水深值(已知随机点的深度值用船去测),然后用前面介绍的双三次样条插值方法(规定边界上的一阶偏导数和二阶混合导数为0)就可求出海底拟合曲面。
四. 符号定义
V swave 声波在海水中的传播速度
V s h i p 测绘船的航速
T s w a v e 声波从发射到接收所经过的时间
Z o c e a n 海底某一点(x,y )处的水深
∆X s h i 在Tswave 时间内测绘船走过的位移 p
S swave 在Tswave 时间内声波走过的路程 S rwave 反射声波走过的路程
α 海底临近的两点之间的倾斜角度
θ 反射声波线与铅垂线之间的夹角
δZ byocean 由于海底地形起伏引起的深度测量误差 δZ byship 由于轮船航行速度引起的深度测量误差 D average 被测海域的平均深度 L side 划分后得到的小区域的边长
D range 一次声纳发射可以覆盖的区域的直径长度
五.模型建立 概述:
在现有的假设和技术背景下,我们现建立如下测绘海底地形的方案: 首 先将被测绘的海底区域以网格方式划分为若干小区域,然后分别
测量网格上每个交叉点处的海底深度;通过双三次样条函数的计算方法,得出每个小区域的曲面图形;最后通过曲面整合的方式,得到整个海底地形图。
误差分析:
1 海底地形引起的误差
由于海底地形是呈波动起伏状的,所以向海底传播的声脉冲的球形波阵面最前沿不会同时到达海底并反射,而是这个波阵面上的一点或几点首先遇到海底地面并反射回波,这个过程如图1.1所示。
图1.1
从图中可以看出测绘船回波接受装置接受到的回波信号将不是位于轮
船铅垂线正下方的A 点反射回的回波信号,而是位于声脉冲覆盖范围内的
另外一点B 点的回波。而且B 点是声脉冲发射出来以后波阵面到达的第一
个反射点,也就是整个声脉冲覆盖区域内与发射装置直线距离最近的海底表
面的一点。这样测量的结果自然会与理论值有一个偏差δZ byocean ,这是由于
海底地形起伏引起的深度测量误差,我们可以利用全微分计算公式估算一下
这个误差造成的影响是否在我们设定的误差范围之内。
图1.2
根据图1.2所示,可以得到
⨯e c o θs +S r w a v ⨯e s i n θ⨯t a αn Z o c e a =n S r w a v
=S rwave ⨯(cos θ+sin θ⨯tan α) (2.1)
假设∆θ,∆α是实际测量时θ和α的偏差,则全增量的绝对值∆Zocean
就是利用上述公式所产生的误差由于∆θ和∆α都很小,我们可以用dZ ocean
来近似地代替∆Zocean ,这样就得到Z ocean 的误差为
∆Z ocean ≈dZ ocean
= ∂Z ocean ∂Z ∆θ+ocean ∆α ∂θ∂α
∂Z ocean ∂Z ocean +δα ≤∂θ∂α
2()=S ⨯cos θtan α-sin θδθ+sin θsec αδα rwave []
(2.2)
将θ=16°,α=15°,δθ= 14°,δα=15°代入式2.2,可以得到Z ocean
的绝对误差约为
δZ byocean =S rwave ⨯[(cos 16︒tan 15︒-sin 16︒)⨯14︒+sin 16︒sec 215︒⨯15︒]
S rwave ≈0. 0773
由此可以得到的相对误差约为 δZ byocean
Z ocean 0. 0773S rwave =⨯100% S rwave ⨯cos 16︒+sin 16︒tan 15︒≈0. 0773⨯100% 1. 0351
≈7. 5%
显然这一误差在我们的允许范围内,所以在实际测量过程中我们可以用
S rwave 的值来代替Z ocean 的值。
2 船速引起的误差
接着我们需要考虑到船速对测量结果的影响。
图1.3
如上图1.3所示,船在向正下方的A 点发出声脉冲时处在测绘方向上的
B 点,此后经过时间T swave 船有了一段位移
∆X ship
S swave =V ship ⨯T swave (2.3) 在C 点,船收到了回波信号。声波所走过的路程 =V swave ⨯T swave =Z ocean +S rwave (2.4)
所以实际测到的A 点深度值为 11(Z ocean +S rwave ) (2.5) S = swave 22
由三角关系可得
S rwave =(Z ocean )2+∆X ship 2 (2.6)
1
2 所以由于船运动所引起的绝对误差为 δZ byship =S swave -Z ocean
1
=(Z ocean +S rwave )-Z ocean 2
=1(S rwave -Z ocean ) (2.7) 2
由于Z ocean 未知,此绝对误差不好估算,我们可以通过估算相应的相对
误差来看其影响。
δZ byship
Z ocecan 1(S rwave -Z ocean )= Z ocean
=1⎛ 2⎝Z ocean 2+∆X ship 2
Z ocean ⎫-Z ocean ⎪⎭
1⎛ ⎝= Z ocean 2+V ship ⨯T swave 2-Z ocean ⎫⎪⎭
Z ocean (2.8)
由于V ship = 2m/s, V swave = 1500m/s,而海洋的平均深度在几十米到几
千米的范围内[4],因此测量时声波在海水中传播到反射回来所花的时间很
短,一般为几秒钟,这段时间里穿走过的路程很短,因此可以认为
T swave
将(2.9)代入(2.8)可得 2⨯Z ocean ≈V swave (2.9) δZ byship
Z ocean ⎛1 2 ≈⎝(Z ocean )2⎫⎛2⨯Z ocean ⎫⎪⎪+ V ⨯-Z ocean ⎪ ship ⎪V ⎪swave ⎝⎭⎭2Z ocean ⨯100%
1⎛
2 ⎝ =(Z ocean )22⎫2⨯Z ocean ⎫⎛+ 2⨯⎪-Z ocean ⎪⎪1500⎭⎝⎭
Z ocean
2⨯100% ⎫1⎛⎛4⎫ =+ -1⎪⨯100% ⎪⎪2 ⎝1500⎭⎝⎭
≈0. 00018%
这个误差远远小于所允许的误差范围,所以我们基本可以忽略船速对测
量结果的影响。
海底区域划分
最后我们需要考虑如何对一片海域划分,以得到所需的数据。很明显,
划分时既不能将区域划分的太大,这样得到的数据太少,不足以反映出实际
的海底地形特征;同时划分也不能太小,这样虽然会更真实地反映出实际的
海底地形特征,但是会造成数据量过多,影响计算的效率。从这种考虑角度
出发,我们认为划分的前提是利用较少的测量次数尽可能的覆盖整个区域,
从而可以从一定程度上确保采集到的数据能够比较完整的反映整个测量区
域的海底地形分布。
考虑要测量的区域是一片如图1.3所示的正方形区域 ABCD ,我们现在
讨论如何将该区域进行划分,即等分的划分为(m+1) (n+1)个小正方形。
图1.4
如果能够确定每个划分后的小区域的边长L side ,则划分就可以得到。
如图1.5所示,若划分的小区域的大小满足如下条件:当测绘船分别在该小
区域的对角线上的两个网格交叉点进行测量时,两处的声纳测量范围刚好相
互外切。由于划分的小区域都处于声纳的测量范围,因此得到的数据能够较
真实的反映该区域的海底地形,而且又不会产生划分过细的情况。
图1.5
由于声脉冲发射的范围为与指向海底的铅垂线夹角从2°--30°之间,
因此可以借助被测海域的平均深度D average ,估算出一次声纳发射可以覆盖
的范围,以此来确定L side 。
D range =2⨯D average ⨯tan 30︒=2
L side =D range ⨯sin 45︒=2D average (2.10) 3326D average ⨯=D average (2.11) 323
因此在划分小区域之前,还要进行一项工作,即估测北侧海域的平均深
度,这可以通过使用上述方法随机测量几组数据来估测。
模型测试:
模型已基本建立,现对模型进行模拟测试,以确定模型的可行性。
A 假设:
假设测绘的海域大小为4km * 4km = 16平方公里
假设该海域平均深度为1200m (以南海的平均深度为例)
按照已分析的海域划分方式(根据公式1.11),该海域应被划分为约1km *1km
的小区域,共4 * 4 = 16块
B 计算
根据假设的海域平均深度1200m ,在1200 ± 1200 * 30% 以内,随机产生5 * 5
个交叉点数据:
根据公式(1.9)(1.10)(1.11)(1.12)„„,所计算出的矩阵p q s为:
p =
z(u,w) =
(-[***********]3/[**************]0*(u-1)^3+[***********]1/12
[1**********]120*(u-1)^2-6183/5)*(w-1)^3+(629309/280*(u-1)^3-886377/280
*(u-1)^2+6183/5)*(w-1)^2+(-171867/140*(u-1)^3+47721/28*(u-1)^2-5643/5
)*(w-1)-57115/56*(u-1)^3+409167/280*(u-1)^2-108
该小区域的海底地形图为:
地形图尚未模拟出来
通过将4 * 4个小区域的海底地形图的连接,生成的整个待测海域的地形图为:
地形图尚未模拟出来
六.模型分析及改进
双三次样条拟合是很有用的一种方法,它具有局部性和连续光滑性;但是
IDW 算出的结点水深处于已知数据点的凸包之中,它还不是“保山脊”的,
即相似已知点与不相似已知点具有相同的影响。另一缺点是模型无法给出误
差的解析估计。
直接用差分法国用IDW 方法等求出足够精细分割的节点上的水深值,采用计
算机绘图的“以值代曲”原理用直线联结相邻叠点来描述曲面,计算量会比
较大,但能更多地反映已知数据点的影响。
参考文献:
[1] 声纳技术
[2] 声纳技术
[3] 惠更斯原理
[4] 挑战地球深渊—海沟
[5] 苏步青,刘鼎元 计算几何 上海科学技术出版社 1981
附录:
测试用例中各个子矩形上的曲面表达式:(坐标轴为u 和w)
[
(-[***********]3/[**************]0*(u-1)^3+[***********]1/
[**************]0*(u-1)^2-6183/5)*(w-1)^3+(629309/280*(u-1)^3-886377/28
0*(u-1)^2+6183/5)*(w-1)^2+(-171867/140*(u-1)^3+47721/28*(u-1)^2-5643/5)*
(w-1)-57115/56*(u-1)^3+409167/280*(u-1)^2-108, (-5014601/3920*(u-1)^3+7801413/3920*(u-1)^2-11244/5)*(w-2)^3+(310463/56
0*(u-1)^3-464127/560*(u-1)^2+11244/5)*(w-2)^2+(-[**************]29/1759
[1**********]*(u-1)^3+[***********]/[**************]0*(u-1)^2-9197
/5)*(w-2)-1849/16*(u-1)^3+15621/80*(u-1)^2-2047/5, ([***********]9/[**************]2*(u-1)^3-[***********]99/
[**************]20*(u-1)^2-1171/2)*(w-3)^3+(-56641/28*(u-1)^3+476853/14
0*(u-1)^2+1171/2)*(w-3)^2+(2022723/3920*(u-1)^3-755211/784*(u-1)^2-867)
*(w-3)+437951/560*(u-1)^3-143523/112*(u-1)^2+563/2, (374533/140*(u-1)^3-518397/140*(u-1)^2-1613/5)*(w-4)^3+(-374533/140*(u-1
)^3+518397/140*(u-1)^2+1613/5)*(w-4)^2+([***********]29/9851624
184872960*(u-1)^3-[***********]67/[**************]20*(u-1)^2-597
9/10)*(w-4)+206869/280*(u-1)^3-266901/280*(u-1)^2+2753/10]
[ ([***********]3/[**************]0*(u-2)^3-1419772
[1**********]1/[**************]*(u-2)^2+[**************]83/[**************]-[***********]/[**************]0*u)*(w-1)^3+(-831179/280*(u-2)^3+553443/140*(u-2)^2-70583/140+115173/280*u)*(w-1)^2+(46661/28*(u-2)^3-30333/14*(u-2)^2-3621/35-38391/140*u)*(w-1)+364569/280*(u-2)^3-250113/140*(u-2)^2+85067/140-38391/280*u,
(4355587/3920*(u-2)^3-3155247/1960*(u-2)^2-714653/392+559023/3920*u)*(w-2)^3+(-378613/560*(u-2)^3+233631/280*(u-2)^2+549697/280+627/112*u)*(w-2)^2+(-[***********]/[**************]0*(u-2)^3+[***********]/[**************]*(u-2)^2-[***********]1/[**************]-[**************]43/[**************]0*u)*(w-2)+8311/80*(u-2)^3-6057/40*(u-2)^2-3339/8+3507/80*u,
(-[***********]/[**************]2*(u-2)^3+[***********]51/[**************]20*(u-2)^2-[***********]3/[**************]20-[***********]3/[**************]0*u)*(w-3)^3+(38729/28*(u-2)^3-186381/70*(u-2)^2+16859/35+104091/140*u)*(w-3)^2+(-1297473/3920*(u-2)^3+322401/392*(u-2)^2-218409/392-1483941/3920*u)*(w-3)-345613/560*(u-2)^3+298119/280*(u-2)^2+12073/56-121377/560*u,
(-72355/28*(u-2)^3+113973/35*(u-2)^2-181319/70+17361/28*u)*(w-4)^3+(72355/28*(u-2)^3-113973/35*(u-2)^2+181319/70-17361/28*u)*(w-4)^2+(-[***********]3/[**************]0*(u-2)^3+[***********]67/[**************]20*(u-2)^2-[***********]57/[**************]20+[**************]83/[**************]*u)*(w-4)-227291/280*(u-2)^3+127083/140*(u-2)^2-78279/140+17361/56*u]
[ (-[***********]61/[**************]0*(u-3)^
3+[***********]/[**************]*(u-3)^2-[***********]/[**************]+[***********]/[**************]*u)*(w-1)^3+(1067823/280*(u-3)^3-39321/8*(u-3)^2+243459/70-20574/35*u)*(w-1)^2+(-62553/28*(u-3)^3+56283/20*(u-3)^2-182137/70+13716/35*u)*(w-1)-442293/280*(u-3)^3+84039/40*(u-3)^2-30661/35+6858/35*u,
(-4216851/3920*(u-3)^3+29607/16*(u-3)^2-2609689/980+251199/980*u)*(w-2)^3+(58641/112*(u-3)^3-95511/80*(u-3)^2+447997/140-9909/28*u)*(w-2)^2+([***********]7/[**************]0*(u-3)^3-[***********]3/[**************]0*(u-3)^2-[***********]/[**************]+[**************]31/[**************]*u)*(w-2)+3489/80*(u-3)^3+12819/80*(u-3)^2-9827/20+1053/20*u,
([***********]/[**************]0*(u-3)^3-[**************]15/[1**********]456*(u-3)^2-[***********]3/[**************]0+[***********]/[**************]0*u)*(w-3)^3+(-127683/140*(u-3)^3+29739/20*(u-
3)^2+191043/70-30249/70*u)*(w-3)^2+(-579119/3920*(u-3)^3-3489/80*(u-3)^2-1981089/980+53583/196*u)*(w-3)+284261/560*(u-3)^3-62943/80*(u-3)^2-23889/140+1713/28*u,
(261323/140*(u-3)^3-53271/20*(u-3)^2+126051/70-21684/35*u)*(w-4)^3+(-261323/140*(u-3)^3+53271/20*(u-3)^2-126051/70+21684/35*u)*(w-4)^2+([***********]/[**************]*(u-3)^3-[**************]43/[1**********]
456*(u-3)^2-[***********]/[**************]0-[***********]/[**************]0*u)*(w-4)+183847/280*(u-3)^3-36669/40*(u-3)^2+19547/14-10842/35*u]
[ ([***********]79/[**************]0*(u-4)^3-[**************]8
3887/[**************]0*(u-4)^2+[***********]1/[**************]0-[***********]89/[**************]0*u)*(w-1)^3+(-847257/280*(u-4)^3+563841/140*(u-4)^2-284643/70+286407/280*u)*(w-1)^2+(278107/140*(u-4)^3-184713/70*(u-4)^2+159137/70-95469/140*u)*(w-1)+291043/280*(u-4)^3-38883/28*(u-4)^2+62753/35-95469/280*u,
(-304399/3920*(u-4)^3-487119/1960*(u-4)^2-370735/98+2861673/3920*u)*(w-
2)^3+(77881/560*(u-4)^3+105519/280*(u-4)^2+811693/140-655719/560*u)*(w-2)^2+(-[***********]/[**************]2*(u-4)^3+[***********]/[**************]*(u-4)^2-[***********]3/[**************]0+[***********]7/[**************]2*u)*(w-2)-28963/80*(u-4)^3+11643/40*(u-4)^2-8371/4+40317/80*u,
(-[***********]/[**************]*(u-4)^3+[**************]03/[1**********]040*(u-4)^2-[**************]01/[1**********]520+[***********]/[**************]*u)*(w-3)^3+(125651/140*(u-4)^3-43719/35*(u-4)^2+82347/35-27201/140*u)*(w-3)^2+(-268063/784*(u-4)^3+1257141/1960*(u-4)^2-8933/98-1007619/3920*u)*(w-3)-55347/112*(u-4)^3+206091/280*(u-4)^2-6911/28+5841/560*u,
(-355889/140*(u-4)^3+279057/70*(u-4)^2-3043/35-48561/140*u)*(w-4)^3+(355889/140*(u-4)^3-279057/70*(u-4)^2+3043/35+48561/140*u)*(w-4)^2+(-[***********]/[**************]*(u-4)^3+[**************]23/[1**********]040*(u-4)^2-[**************]49/[1**********]520+[**************]63/[**************]*u)*(w-4)-245737/280*(u-4)^3+196443/140*(u-4)^2+41351/70-48561/280*u]