一种适用于台风区域的沙波运移预测方法与流程

本发明涉及海洋工程领域，特别是涉及一种预测台风影响区域的沙波运移时，综合非线性谱分析和波浪输沙理论的方法。

背景技术：

海底蕴藏着丰富的矿物资源和油气资源。随着陆地上的可供利用的资源日渐枯竭，人类开始更多地开发和利用海洋资源。近年来海洋油气业的快速发展，海洋石油日益成为我国原油增量的主要来源。海底沙波是一种很常见的近岸海底地貌，其垂向剖面呈起伏状。在海洋水动力作用下，沙波的形态和位置都在随时间(年)变化，属于不稳定的地貌类型。沙波运移可能造成海底管道(线)的悬空，甚至发生断裂、海洋平台垮塌、航道淤积等重大危害。准确预测沙波运移能够减少经济损失，并服务于海洋工程结构物的选址设计。

目前常用的沙波运移预测技术主要有以下几种：

现场实测：采用多波束测量方法对沙波区进行多次水深调查后比较历次测量结果，确定沙波运移方向和运移速率。很明显这种方法耗资巨大，而且调查时间受气象条件和海况的制约，其间隔一般都在一年以上，根本无法实时监测和预测沙波运移。在可预计的未来，观测的时间和空间分辨率及长度仍然十分有限。

经验公式计算：研究者普遍认为，一般情况下近底床流动是沙波运移的主要动力，沙波运移速率主要取决于流速大小，其运移方向应和流速方向一致。相关研究提出了一系列不同的沙波运移速率计算公式。这些公式都是依据推移质输沙率、沙波波高等参数来计算沙波运移速率，其中推移质输沙律由底床剪切流速或平均流速计算而得。这些公式一般为工程界和力学界所采用，早期用于预测河道中的沙波运移速度。后来也有应用于海洋中沙波运移的模拟，但对于往复潮流影响下的沙波运移并不适用。

动力学计算：其从沙质底床和海流相互作用的角度研究强潮流区的沙波的形成和运移。目前国际上常用的近海和区域海洋数值模型包括pom、ecom、roms、fvcom等，它们被广泛应用于模拟海洋潮汐、流场和温盐场等。然而方法存在计算消耗大的特点。

谱分析方法：该方法将沙波看成一系列余弦波的叠加，采用二阶谱分析并考虑沙波的非线性，并将沙波分成自由波和束缚波，其中束缚波为非线性项，描述非对称(尖峰平底或锯齿形)沙波。但是单纯谱分析方法的预测结果不稳定，而且不适用于台风影响区域。

技术实现要素：

本发明的目的是要提供一种预测台风影响区域的沙波运移时，综合非线性谱分析和波浪输沙理论的方法。

特别地，本发明提供一种适用于台风区域的沙波运移预测方法，包括如下步骤：

步骤100，设置一个以年为单位的时间段，取海底底床同一位置首尾两个时段的沙波图像；

步骤200，利用谱分析方法，将沙波看成一系列余弦波的叠加，通过参数估计得到自由波傅里叶变换得到的波谱、自由波的线速度以及二阶传递函数,计算波谱能量平均速度并获取该时间段内沙波的总运移；

步骤300，首先计算该时间段内单次台风中合成浪的波高和周期，再利用非线性输沙公式计算其引起的泥沙运移，然后叠加计算所有台风引起的泥沙运移，以得到该时间段内所有台风引起的台风运移；

步骤400，利用总运移减去台风运移以得到无台风因素下的常规运移，用常规运移除以该时间段即得到沙波在单位时间内的常规运移速度；

步骤500，根据波谱、二阶传递函数和常规运移速度即可预测目的年份的沙波剖面。

在本发明的一个实施方式中，所述步骤200中利用谱分析方法过程如下：

将沙波分成自由波和束缚波，先获取沙波不运移状态下的沙波形状：

b(x)＝a(x)+βa²(x)(1)

其中a(x)是一系列余弦波的叠加形成的输入，b(x)是输出，β是二阶传递函数；li是特征波长；对式(1)两边同时做傅里叶变换得：

b(k)是b(x)的傅里叶变换，a(k)是a(x)的傅里叶变换，对上式进行对角项简化：

变形得到：

b(k)＝h⁽¹⁾(k)a(k)+j⁽²⁾(k)f{a²(x)}(5)

h⁽¹⁾(k)是线性一阶传递函数，j⁽²⁾(k)是平方二阶传递函数，f{·}表示傅里叶变换，而f^-1{·}表示逆傅里叶变换；

再加入时间t来获取沙波运移时海底底床演变的自由波表示：

令h⁽¹⁾(k)＝1，j⁽²⁾(k)＝β(k)，束缚波利用二阶传递函数β(k)进行描述，将自由波和束缚波叠加得到输出波形为：

b(x,t)＝a(x,t)+f^-1{β(k)f{a²(x,t)}}(7)

至此，非线性沙波通过自由波a(x)或其傅里叶变换得到的波谱a(k)、自由波的线速度v(k)以及二阶传递函数β(k)实现描述。

在本发明的一个实施方式中，所述步骤200中，通过参数估计得到自由波傅里叶变换得到的波谱、自由波的线速度以及二阶传递函数的过程如下；

步骤210，将初始波形的谱分为自由波和束缚波，初始迭代时，假定自由波位于谱峰波数附近，而束缚波位于二倍的谱峰波数附近，得到自由波波谱a(k)；

步骤211，通过式(6)和式(7)计算得到初始波形的输出波形b(x,t)，将实际测量得到的波形图设为目标图形z(x,t)，将目标图形z(x,t)减去输出波形b(x,t)的差值看成是一系列余弦波的叠加并计算出来，将差值对应得到的波谱叠加到初始波形的波谱a(k)上；

步骤212，得到新的波谱a(k)后，再次计算，又得到新的输出波形b(x,t)，根据精度需要，反复迭代，直至达到目的后收敛。

得到波谱a(k)、自由波的线速度v(k)以及二阶传递函数β(k)之后，波谱能量平均速度通过下式计算：

式中u为波谱能量平均速度，v(k)和a(k)分别代表波数为k的自由波对应的速度和频谱值大小；

计算时段内的沙波总运移ls＝utt；tt为初始时刻与末时刻的时间间隔。

在本发明的一个实施方式中，所述步骤300中，计算该时间段内单次台风中合成浪的波高和周期的过程如下：

采用holland气压模型计算台风风场，描述以台风中心为起点的任意剖面的气压分布，计算为：

式中，r为距离；p为距离台风中心r处的气压；p0为台风中心最低气压；δp＝pn-p0为台风中心气压降，pn为台风外围环境气压；a和b为台风形状参数；

设rmax为最大风速半径，rmax＝a^1/b，不考虑科氏力，通过地转风方程得到台风的风速剖面：

式中，vr为距离台风中心r处的台风环流风速，vrmax为环流最大风速；

采用陈奇礼经验公式先计算涌浪和总浪的波高，根据叠加原理得到风浪的波高，进而求出合成浪的波高及周期：

式中hr为合成浪波高，tr为混合浪周期；hrf为风浪波高，trf为风浪周期；hru为涌浪波高，tru为涌浪周期。

在本发明的一个实施方式中，所述步骤300中，计算单次台风引起的泥沙运移过程如下：

利用watanabe公式计算单次台风中一个周期内波浪的净输沙率：

式中qb,net是以体积计的单宽净输沙率，d为泥沙粒径；ω为泥沙沉速，ψ为shields数，其中临界shields数ψc＝0.11，τbm为底部剪应力最大值，ρ为水的密度，ρs为泥沙密度，g为重力加速度，fw为底部摩擦系数；底部水平轨迹速度最大值h为波浪的波高，t为波浪的周期，k为波数，h为水深；

底部摩擦系数fw通过与雷诺数和相对糙率的关系式计算，过程如下：

式中雷诺数底部水质点运动振幅ks为底床粗糙高度，取ks＝2d90；

再通过rubin公式即可计算出单次台风中沙波的运移速率：

式中ug为沙波运移速率，hs为沙波高度，γ为沉积物容重，qs为底沙输运率。

在本发明的一个实施方式中，所述步骤400中，获取常规运移速度vc的公式如下：

式中，tt为沙波运移时间，ls为时段内总运移，lt为时段内的台风引起的沙波运移。

本发明基于统计、物理和实用三个角度，提出波谱能量平均的运移速度，预测结果稳定且准确。与已有的标量形式不同，计算台风运移过程作用矢量，并且对不同台风不同位置时的影响进行叠加，根据短时间一些台风影响大，但长时间范围内台风影响相互抵消的原因，计算常规波流相互作用引起的运移，并用其预测长时段波形。本文的模型能够准确拟合台风区的沙波波形，与已有其他方法相比，具有高效高精度的优点，可用于10～20年内的底床演变预测。

附图说明

图1是本发明一个实施方式的预测方法流程示意图；

图2是北部海湾地貌及研究区海底地形图及剖面位置示意图；

图3是图2中2004年和2007年的波形实测值和计算值示意图；其中初始波形为2004测量的波形，目标波形为2007测量的波形；

图4是本发明一个实施方式中台风时风浪和涌浪叠加的示意图；

图5是台风“达维”移动路径和holland模型计算的风场示意图；

图6是台风“达维”引起的研究区风流、涌浪及合成波浪的波高示意图；

图7是台风“达维”引起的研究区风流、涌浪及合成波浪的周期示意图；

图8是台风引起的沙波运移示意图；

图9是台风相互作用抵消示意图；

图10是本发明一个实施方式根据计算得到的参数预测2020年的沙波波形示意图。

具体实施方式

本方案是对一定年限间正常和台风天气数据进行分析，然后根据得到结果去预测指定年的沙波运移。处理过程如下：首先，根据采集的海底底床高度数据给定两个时段的初始波形和末时刻波形，如2004年至2007年；通过谱分析得到沙波的总位移，然后，考虑台风运移中位置的改变，计算单个台风的作用，将时段内所有台风影响叠加，便得到台风引起的总位移，计算过程中发现，虽然单个台风可能引起沙波大的位移，但是长时间范围内不同台风作用相互抵消；最后，将总位移减去台风总的台风位移，得到常规波流作用引起的位移，计算常规运移速度，即可利用该常规运移速度对未来的沙波波形进行预测。

如图1所示，本发明一个实施例公开的适用于台风区域的沙波运移预测方法一般包括如下步骤：

步骤100，设置一个以年为单位的时间段，取海底底床同一位置首尾两个时段的沙波图像；

这里的时间段可以是预测年之前的任何一个时间段，一般取3-5年即可。如图2所示，如取04-07年这个时间段，研究区2004年海底地形图及剖面位置(这里有六个剖面，即p1～p6)。选取wgs84大地坐标系，和墨卡托投影(utm)办法，将测量得到经纬高程坐标转换为xyz坐标，并导入surfer软件。在沙波明显的区域截取多个剖面，得到剖面的始末点的位置坐标和剖面上各点的底床高度。由于研究区域远离海岸，故而不存在泥沙的输入和输出，只考虑床面起伏的振荡部分，该部分认为是相应时间范围内的主要移动部分。调整剖面的平均高程为0，得到剖面各年的波形。

步骤200，利用谱分析方法，将沙波看成一系列余弦波的叠加，通过参数估计得到自由波傅里叶变换得到的波谱、自由波的线速度以及二阶传递函数；计算波谱能量平均速度并获取该时间段内沙波的总运移；

在将沙波看成一系列余弦波的叠加时，考虑沙波的非线性，将沙波分成自由波和束缚波。

(1)当不考虑沙波运移时，设沙波形状可表示为

b(x)＝a(x)+βa²(x)(1)

其中a(x)是输入，b(x)是输出，β是二阶传递函数。输入的a(x)是一系列余弦波的叠加，li是特征波长。对式(1)两边同时做傅里叶变换得：

b(k)是b(x)的傅里叶变换，a(k)是a(x)的傅里叶变换，估计双参数的二次传输函数h⁽²⁾(k1，k2)需要大量的时间，对上式进行对角项简化：

变形得到：

b(k)＝h⁽¹⁾(k)a(k)+j⁽²⁾(k)f{a²(x)}(5)

h⁽¹⁾(k)是线性(一阶)传递函数，j⁽²⁾(k)是平方(二阶)传递函数。f{·}表示傅里叶变换，而f^-1{·}表示逆傅里叶变换。

根据前述推导过程，可推导出沙波运移时的参数计算基础。

(2)考虑沙波运移，加入时间t，底床演变的线性部分(自由波)表示为：

令h⁽¹⁾(k)＝1，j⁽²⁾(k)＝β(k)，非线性部分(束缚波)利用二阶传递函数β(k)进行描述，将自由波和束缚波叠加得到输出波形为：

b(x,t)＝a(x,t)+f^-1{β(k)f{a²(x,t)}}(7)

非线性沙波可通过自由波a(x)或其傅里叶变换得到的波谱a(k),自由波的线速度v(k)以及二阶传递函数β(k)进行描述。

其中的参数a(k)、v(k)以及β(k)的估计，利用如下的迭代法得到：

step1：将初始波形(如：t＝2004年时刻的波形)的谱分为自由波和束缚波，初始迭代时，假定自由波位于谱峰波数附近，而束缚波位于二倍的谱峰波数附近，得到自由波波谱a(k)；

step2：通过式(6)和式(7)计算得到b(x,t)；

step3：目标图形(如：t＝2007s时刻的波形)为z(x,t)，z(x,t)-b(x,t)即目标图形和计算图形的差值，将差值看成是一系列余弦波的叠加；这里目标图形是实际测量得到的波形图z(x,t)，计算波形是用式(7)计算得到的b(x,t)。

step4：将z(x,t)-b(x,t)剩余得这部分余弦波计算出来，再叠加到初始自由波组分上，(z(x,t)-b(x,t)是测量值和计算值的差值，这部分差值进行逆傅里叶变化，得到差值对应的波谱，然后再加到step1中的a(k)上面；

step5：得到新的自由波谱a(k)，再次计算，这样就得到了第二步的b(x,t)；

step6：根据精度需要，反复迭代，直至收敛。

以上结果得到的自由波线速度结果不稳定，从统计角度考虑，速度反应余弦波相位的改变，小振幅波速度意义不大；从物理角度考虑，沙波运移

以推移质平移为主，波形不会发生大的改变；从实践角度考虑，利用能量平均速度拟合误差差异很小，但是不同时段数据预测结果准确。基于以上三方面考虑，又振幅的平方代表了能量，引入波谱能量平均速度处理谱分析结果，如下：

式中u为波谱能量平均速度，v(k)和a(k)分别代表波数为k的自由波对应的速度和频谱值大小(即自由波振幅值)。计算时段内的沙波总运移ls＝utt；tt为初始时刻与末时刻的时间间隔。如图3所示，以其中的一个剖面p1为例：2004年和2007年的波形实测值和计算值如图，初始波形为2004测量的波形，目标波形为2007测量的波形。

步骤300，首先计算该时间段内单次台风中合成浪的波高和周期，再利用非线性输沙公式计算其引起的泥沙运移，然后叠加计算所有台风引起的泥沙运移，以得到该时间段内所有台风引起的台风运移；

这里以2005年最大的台风“达维”为例，先计算单个台风引起的位移。

(1)计算台风期间的风浪和涌浪的波高和周期，并计算合成浪的波高和周期。

台风期间波高较大，相对来说潮流作用小，故而计算台风作用时只考虑波浪作用下的沙波运移。考虑波浪非线性，计算波质点运动的不对称性所产生的净输沙率，选用watanabe公式，得到无量纲输沙率与无量纲剪应力的函数关系。泥沙沉速采用rubey公式计算。

采用holland气压模型计算台风风场，描述以台风中心为起点的任意剖面的气压分布，可写为：

式中，r为距离；p为距离台风中心r处的气压；p0为台风中心最低气压；δp＝pn-p0为台风中心气压降，pn为台风外围环境气压；a和b为台风形状参数，rmax为最大风速半径，rmax＝a^1/b。

不考虑科氏力，通过地转风方程，可得到台风的风速剖面：

式中，vr为距离台风中心r处的台风环流风速，vrmax为环流最大风速。

在圆对称风场中，风浪和涌浪的合成浪方向简化为如图4所示。

根据风(涌)浪波高和周期的陈奇礼经验公式，计算合成浪的波高和周期，由于风场利用圆对称近似，采用经验公式先计算涌浪和总浪的波高，根据叠加原理得到风浪的波高，进而求出风、涌浪的周期以及合成周期：

式中hr为合成浪波高，tr为混合浪周期；hrf为风浪波高，trf为风浪周期；hru为涌浪波高，tru为涌浪周期。

上述结果可参见图5所示的台风“达维”移动路径和holland模型计算的风场示意；图6所示的台风“达维”引起的研究区风流、涌浪及合成波浪的波高示意；图7所示的台风“达维”引起的研究区风流、涌浪及合成波浪的周期示意。

(2)计算单个台风引起的泥沙输运；用前述得到的合成浪波高和周期进行计算作用位移，然后换算到剖面的方向。

考虑波质点运动不对称性所产生的净输沙，利用watanabe公式计算波浪一个周期内的净输沙率：

式中qb,net是以体积计的单宽净输沙率，d为泥沙粒径，ω为泥沙沉速，用rubey公式计算，研究区内泥沙粒径为0.2mm，ψ为shields数，其中临界shields数ψc＝0.11，τbm为底部剪应力最大值，ρ为水的密度，ρs为泥沙密度，g为重力加速度，底部水平轨迹速度最大值h为波浪的波高，t为波浪的周期，k为波数，h为水深。

底部摩擦系数fw通过与雷诺数和相对糙率的关系式计算：

式中雷诺数底部水质点运动振幅ks为底床粗糙高度，取ks＝2d90。

利用rubin公式计算沙波的迁移速率：

式中ug为沙波运移速率，h为沙波高度，γ为沉积物容重，qs为底沙输运率(以重量计)。

台风在研究区域合成波浪方向即为台风对沙波作用方向，近似取研究区域中心经纬度作为作用点，计算单个台风对沙波剖面的作用位移矢量ft。对于不同的实测沙波剖面，台风作用的大小与其波高和水深的不同有关，台风对于沙波运移的作用为作用位移矢量ft在剖面方向上的投影。

(3)计算2004-2007年间多个台风引起的沙波运移lt。如图8所示：04-07年共计85个台风，先计算单个台风引起的作用运移矢量，然后将时段内所有台风造成的影响叠加。由图9的台风作用抵消示意图可知，从台风路径1和台风路径2过去的两个台风，分别产生了合成方向1和合成方向2的合成浪(即风浪+涌浪)，对于沙波运区域的作用方向相反。

对沙波产生运移影响的只有少数几个台风，这些台风风速大且运行路径接近研究区域，计算得到04-07年所有台风对p1的作用位移矢量ft为[-2.0,-0.7]m，表示台风的作用总方向为西南方。将ft投影到剖面方向，得到剖面沙波产生的位移，04-07年沙波运移-0.8m。并且可知，虽然单个台风可能引起沙波大的位移，但是长时间范围内不同台风作用相互抵消。故而可以在总位移中减去台风影响位移，得到常规波流作用引起的位移，以此做预测。

步骤400，利用总运移减去台风运移以得到无台风因素下的常规运移，用常规运移除以该时间段即得到沙波在单位时间内的常规运移速度；

常规运移速度式中，tt为沙波运移时间，ls为时段内总运移，lt为时段内的台风引起的沙波运移。

步骤500，根据波谱、二阶传递函数和常规运移速度即可预测目的年份的沙波剖面。

利用步骤200得到的波谱、二阶传递函数和和步骤400中得到的常规运移速度预测目的年份的沙波剖面。效果如图10所示。

本发明基于统计、物理和实用三个角度，提出波谱能量平均的运移速度，预测结果稳定且准确；与已有的标量形式不同，计算台风运移过程作用矢量，并且对不同台风不同位置时的影响进行叠加；发现短时间一些台风影响大，但长时间范围内台风影响相互抵消，计算常规波流相互作用引起的运移，并用其预测长时段波形。本文的模型能够准确拟合台风区的沙波波形，与已有其他方法相比，具有高效高精度的优点，可用于10～20年内的底床演变预测。

至此，本领域技术人员应认识到，虽然本文已详尽示出和描述了本发明的多个示例性实施例，但是，在不脱离本发明精神和范围的情况下，仍可根据本发明公开的内容直接确定或推导出符合本发明原理的许多其他变型或修改。因此，本发明的范围应被理解和认定为覆盖了所有这些其他变型或修改。