自然界的许多现象, 包括海水表面波 ^[1–6]、一般流体湍流 ^[7–12]或气体湍流^[13,14], 以及等离子体边缘湍流^[15–20]等, 都可以用一个非线性系统来模拟. 在这种结构系统中, 对给定输入的响应由线性和非线性机制决定. 这种非线性系统包含线性和二次非线性输入-输出关系. 即此系统由一组源(“输入”)信号和响应(“输出”)信号来进行描述. 于是其可被看作一个“黑盒”. 用一个包含线性元素、二次非线性元素和更高阶的非线性元素的网络来近似此“黑盒”, 就可能通过测试来了解这个系统的动力学. 最简化的这种网络包含一个信号输入和信号输出系统, 且是纯线性的, 因为可以假定忽略系统的非线性贡献. 但是在许多结构系统中非线性却是不可忽略的^[1]. 针对模型系统的输入与输出信号谱, 可将模型系统分为两类: 一是高斯型模型, 输入信号谱是高斯信号; 二是非高斯型模型, 输入信号谱是非高斯的. 对高斯型和非高斯型非线性系统的讨论很多, 已有一些系统性的研究^[3,21–28]. 湍动激励-响应型系统是一类特殊的非高斯型非线性系统. 系统的输入是湍流波动谱信号. 湍流波动信号泛指一般湍流和向湍流过渡(transit to turbulence)的流体的状态信号, 可以是扰动位移信号$ X\left(t\right) $、扰动电势信号$ \varphi \left(t\right) $或者扰动电场信号$ E\left(t\right) $等. 这些信号的测量值构成湍流数据. 一般在一个系统中可以得到具有随机信号特征的大量湍流数据, 基于谱估计方法来进行处理. 谱估计分析是利用大量特征数据来估计随机变量的参数, 包括随机变量的数学期望(系综平均)、矩、方差和相关系数等. 根据不同的湍流类型, 湍动激励-响应型系统又可分为输入信号为近高斯信号谱的二次系统和输入是完全非高斯型信号谱的系统. 后者的代表是边缘等离子体湍流. 前者是常见湍动激励型系统, 其特征是作为输入的湍流波动信号功率谱分布接近于高斯分布. 对于这种湍流, 已有成熟的模拟手段, 易于进行仿真研究. 例如, 根据湍流的特点, Hasegawa和Mimma^[29]给出了一组关系来模拟充分发展湍流. 充分发展湍流是指已经发展成为稳态(stationary)的湍流. 这种湍流的特征是湍流能量在湍流内部传递、流动, 与外界基本上没有能量交换. 利用高斯色噪声作初始信号, 可以仿真近高斯输入型充分发展湍流.

出于计算效率和结果有效性的考虑, 不同种类模型适宜应用不同的算法求解. 本文将专门研究输入信号功率谱为近高斯分布型的特定湍动激励型系统, 给出一个建模方案, 并基于互相关技术, 根据所测量的输入信号和输出信号来确定这个线性系统模型. 然后设计合适的算法求解系统参数, 求解线性传递函数和二次非线性传递函数. 应用此非线性网络, 可以估计随机激励的非线性响应. 针对基于双谱分析技术的三波耦合系统建模, 人们已经提出了3种求解算法, 分别是Ritz法^[23–25]、Kim方法^[26]和完备迭代法^[27]. 其中, Shen等^[27]首先引入的完备迭代法是Ritz法的改进版, 主要特点是保留了四阶矩, 保留了所有有用的谱信息, 从而计算精度更高. 本文将考虑应用完备迭代法作为湍动激励型模型的求解算法, 并在仿真和应用分析中验证其有效性.

船只通常在与入射波谱相关的频率下发生相当轻微的振荡, 但当停泊在随机海洋中时, 可以发生非常大幅度的低频振荡. 这种现象通常被称为低频“漂移振荡”. 这些漂移振荡支配着系泊船舶的运动. 文献[1, 2]研究了系泊船舶对入射海浪场的漂移响应, 并分别基于互谱分析技术进行建模, 然后分别采用回归分析法进行求解. 这类现象的特点就是作为输入的湍流波动信号功率谱接近高斯分布. 上述工作的不足之处在于他们对所建模型缺乏系统性的描述, 而且求解算法效率不高. 本文将从模拟分析这一问题入手, 通过完善海波激励-系泊船舶漂移响应现象的建模分析, 给出了一种通用的基于双谱分析技术建模的系统方法, 提出并应用完备迭代法^[27]进行模型求解. 谱分析技术普遍适用于大数据量随机过程的数值分析, 是与湍流数据分析相关研究的主要手段. 在案例计算分析中, 仿真了海波激励-系泊船舶响应现象和一种充分发展湍流问题, 并进行了模型验算. 本文研究可以重现文献[1, 2]的结果, 但模型更系统、简明, 求解算法效率更高. 明确了基于完全迭代算法的程序应用到包括特定湍动激励-响应型二次非线性系统时特别准确. 另外, 在模拟海波-船舶系统和其他两个系统时, 本文在输入与输出信号谱计算中都加了汉宁(Hanning)窗^[30], 提高了计算精度和置信度. 最后, 讨论了模型和算法的误差来源及相应改进计算精度的措施.

如前所述, 近高斯分布型湍动激励系统的典型示例之一, 是系泊船舶随随机海浪漂移运动模型. 当驳舶系泊在随机海洋中时, 船舶在系泊系统的无阻尼固有频率处或其附近发生大幅度振荡. 这种现象被称为低频漂移振荡, 因为这种运动的频率远低于入射海浪的频率. 图1 ^[1]给出了表征这一现象的数据示例, 其中显示了随机海浪高度和系泊驳船摇摆响应的时间轨迹和功率谱. 图1清楚地表明, 系泊驳船的振荡频率(可以建模为具有大质量和小恢复刚度的动态系统)明显低于作为入射海浪的振荡频率. 海浪高度和驳船摇摆响应之间的互双谱在与低频驳船摇摆响应相同的不同频率组合处具有振幅峰值. 这意味着入射海浪的各种频谱分量通过二次非线性相互作用促使驳船作出低频摇摆响应. 因此, 在对随机海浪激励与系泊船舶相应低频漂移振荡之间的关系进行建模时, 必须考虑二次非线性. 考虑将互双谱分析应用于系泊船舶系统的漂移振荡数据分析, 以模拟从高频海浪输入到低频摇摆响应的非线性(二次)能量传递机制. 实验数据^[1]可以从系泊船舶在不规则波条件下的按比例(1∶48)定制模型的波盆(wave basin)试验中获得. 变化的参数包括波高、频率、系泊船舶航向和系泊刚度(mooring stiffness). 观察系泊船舶运动情况, 其结果应是控制参数的函数. 此处分析的数据适用于横浪中的系泊船舶. 图1给出了不规则海浪高度和系泊船舶摇摆响应的典型时间轨迹及其各自的自功率谱. 在此系统中, 首先, 需要进行线性谱分析, 以确定波高输入和系泊船舶摇摆响应之间的线性关系. 因为部分输出功率可能归因于波高输入和系泊船舶摇摆响应之间的线性关系. 还需要进行互谱分析, 因为部分输出功率可能来自二次关系作用的结果.

此外, 一类充分发展湍流的自功率谱具有近似高斯分布特征, 也可以视为近高斯输入型系统. 这种湍流的不同波数或频率处的波可以有自己的增长率, 能量在湍流系统内不同波数或频率之间传递. 此种充分发展湍流在波数空间里的典型的自功率谱和湍流信号如图2所示. 其中, 湍流输入信号谱与输出信号谱只有很小的变化. 模拟充分发展湍流系统的构建过程是: 以一个高斯色噪声作为初始输入, 连续5次通过满足湍流发展特征^[24]的封闭非线性系统, 最后得到全发展湍流的信号谱, 如图2(a)所示. 相应湍流模拟信号(部分)如图2(b)所示. 一般地, 分析充分发展湍流, 涉及分析湍流波增长率及系统内能量传递函数的大小. 对上述两类代表性的近高斯输入激励型非线性系统, 可以用双谱分析技术进行建模分析, 完成各自所需的模型描述(求模型参数)与系统分析工作.

在海波和充分发展湍流系统中, 已知输入与输出信号, 可以考虑输入信号经过线性和非线性变换过程, 生成输出信号. 这个系统可以看作一个“黑盒子”. 当系统由二次非线性作用主导, 高于二阶的非线性项如三阶、四阶、五阶等高阶非线性项可以忽略时, 可以用二次非线性系统来进行建模分析. 例如本文所用事例, 高于二阶的非线性可以忽略不计, 并且输入-输出关系可以通过高达二阶的Volterra级数来合理描述. 由于波或模式相互作用等现象在频域中最容易理解, 本工作的重点是在频域, 因此, 建模将根据传递函数而不是内核本身进行, 其中传递函数被定义为内核的傅里叶变换. 内核对应于系统的线性和二次脉冲响应. 关键思想是用线性和二阶传递函数的并行组合来建模未知的非线性系统. 因此, 我们的模型可以表示为

其中$ X\left({f}_{m}\right) $和$ Y\left({f}_{m}\right) $分别是输入(海浪)和输出(驳船摇摆响应)的N点离散傅里叶变换. 在(1)式中的$ {H}_{1}\left({f}_{m}\right) $称为线性传递函数, $ {H}_{2}\left({f}_{i}, {f}_{j}\right) $称为二次传递函数, 且$ {H}_{2}\left({f}_{j}, {f}_{i}\right)={H}_{2}\left({f}_{i}, {f}_{j}\right) $; $ {\varepsilon }_{m} $代表系统误差.

一般地, 当将频域空间换为波数空间来描述此系统会更加方便. 即简单地将频率$ {f}_{m} $, $ {f}_{i} $和$ {f}_{j} $用波数$ k, {k}_{1} $和$ {k}_{2} $代替, 而用$ {L}_{k} $和$ {Q}_{k}^{{k}_{1}, {k}_{2}} $取代$ {H}_{1}\left({f}_{m}\right) $和$ {H}_{2}\left({f}_{i}, {f}_{j}\right) $, 我们的模型可扩展为

即将整个系统可用输入信号$ {X}_{k} $与输出信号$ {Y}_{k} $来描述.

方程(2)分为线性元项、二次元项与误差项三部分, $ {L}_{k} $为线性传递函数, $ {Q}_{k}^{{k}_{1}, {k}_{2}} $为二次传递函数, 一般都是复杂的量. 通过实验测量可得的量是输入信号$ x\left(s\right) $与输出信号$ y\left(s\right) $, 其傅里叶变换分别为(2)式中的$ {X}_{k} $和$ {Y}_{k} $. 误差项$ {\varepsilon }_{k} $是与(2)式中前两项在统计上独立的一个过程的傅里叶变换, 它是测量中的固有噪声和不能用线性项与二次项描述的系统误差. 图3给出了方程(2)所描述的黑盒子(非线性系统)的示意模型.

在湍流内部能量传递的估测算法中, 最终的目标就是利用实验中所测得的输入信号$ x\left(s\right) $与输出信号$ y\left(s\right) $, 基于(4)式给出的二次非线性方程, 估计出线性传递函数$ {L}_{k} $与二次传递函数$ {Q}_{k}^{{k}_{1}, {k}_{2}} $, 从而准确描述该黑盒子系统, 并可根据新的输入信号$ x' (s) $计算出其对应的输出信号$ y'(s) $. 对于充分发展湍流系统, 考虑如下形式的三波耦合非线性方程^[25]:

其中, $ \phi \left(k, t\right) $是在波数$ k $处的密度波动$ \varphi \left(x, t\right) $的傅里叶分量, 定义为$ \varphi \left(x, t\right)=\displaystyle\sum\nolimits _{k}\phi \left(k, t\right){{\mathrm{e}}}^{{\mathrm{i}}kx} $; $ {\varLambda }_{k}^{L} $是线性耦合系数, 其中包含了增长率($ {\varLambda }_{k}^{L} $的实部)和平均色散($ {\varLambda }_{k}^{L} $的虚部); $ {\varLambda }_{k}^{Q} $是二次非线性耦合系数. (3)式描述对象是自然界的很多波动系统, 包括海洋浅水波及边缘等离子体全发展湍流系统. 其中, 四波耦合及五波、六波······耦合相对三波耦合可以忽略. 对于四波耦合不可忽略的非线性系统如海洋深水波^[31], 则必须考虑此项.

当使用有限差分方法处理^[19]后, (3)式可以转换为(2)式. 其中, 主要的变量替换为: $ {X}_{k} = \phi (k, t) $, $ {Y}_{k}=\phi \left(k, t+\tau \right) $, 以及

其中, $ \varTheta \left(k, t\right) $代表时刻$ t $处的相位. 确定传递函数$ {L}_{k} $与$ {Q}_{k}^{{k}_{1}, {k}_{2}} $后, 可以间接确定耦合系数, 进而确定湍流波增长率与非线性能量传递函数. 本文研究目标是提供一种算法技术, 以便能通过实验测量所得的波动信号$ x\left(s\right) $和$ y\left(s\right) $, 高效地估计线性传递函数与二次传递函数.

首先用$ {X}_{k} $的复共轭$ {X}_{k}^{*} $乘以(1)式, 并针对许多组统计上类似的实现(realization)作系综平均, 得到二阶矩方程:

接着使用$ {X}_{{k}_{1}'}^{*}{X}_{{k}_{2}'}^{*} $乘以(1)式并作系综平均, 得到包含四阶矩$\langle{X}_{{k}_{1}}{X}_{{k}_{2}}{X}_{{k}_{1}'}^{*}{X}_{{k}_{2}'}^{*}\rangle $的方程, 其中$ {k}_{1}'+{k}_{2}'=k, {\mathrm{且}}\left({k}_{1}', {k}_{2}'\right)\ne \left({k}_{1}, {k}_{2}\right) $. 具体方程如下:

在尖括号里的统计量是自功率谱$ {P}_{k}^{X}= \left\langle{{X}_{k}{X}_{k}^{*}}\right\rangle $、互功率谱$ {C}_{k}^{YX} = \left\langle{{Y}_{k}{X}_{k}^{*}}\right\rangle $、自双谱$ {B}_{{k}_{1, }{k}_{2}}^{XXX} = \langle{X}_{k} {X}_{{k}_{1}}^{*}{X}_{{k}_{2}}^{*}\rangle $、互双谱$ {B}_{{k}_{1, }{k}_{2}}^{YXX}=\left\langle{{Y}_{k}{X}_{{k}_{1}}^{*}{X}_{{k}_{2}}^{*}}\right\rangle $, 以及四阶矩$ \langle{{X}_{{k}_{1}}{X}_{{k}_{2}}{X}_{{k}_{1}'}^{*}{X}_{{k}_{2}'}^{*}}\rangle $.

为得到$ {L}_{k} $的计算式, 直接对(4)式进行移项:

类似地, 对(5)式进行移项, 可得$ {Q}_{k} $的计算式: 对$ i=1, {\mathrm{ }}2, \cdots , {\mathrm{k}}{\mathrm{n}}{\mathrm{u}}{\mathrm{m}} $ ($ {\mathrm{k}}{\mathrm{n}}{\mathrm{u}}{\mathrm{m}} $代表在表达式$ k= {k}_{1}+{k}_{2} $, 且$ {k}_{1}\geqslant {k}_{2} $中, k所对应的成对出现的$ {k}_{1}, {k}_{2} $的组数),

即:

有了(6)式和(8)式, 就可以通过迭代法求解$ {L}_{k} $与$ {Q}_{k}^{{k}_{1}, {k}_{2}} $的值. 这种方法可称为“完备迭代法”^[27], 因为迭代过程中保留了四阶矩, 从而保留了谱的全部信息. 基于(6)式和(8)式的求解算法如下. 一般地, 迭代法会从$ {L}_{k} $开始进行. 迭代开始时, 取初始值$ {Q}_{k}^{{k}_{1}\left(i\right), {k}_{2}\left(i\right)}=0, i={\mathrm{1, 2}}, \cdots , {\mathrm{k}}{\mathrm{n}}{\mathrm{u}}{\mathrm{m}} $. 作为初始的猜测, 直接忽略$ {Q}_{k}^{{k}_{1}, {k}_{2}} $项后, 得到线性传递函数以计算$ {L}_{k} $的初始值:

然后在迭代中, 只需要不断重复将$ {L}_{k} $代入求解$ {Q}_{k}^{{k}_{1}, {k}_{2}} $ ((8)式)、将$ {Q}_{k}^{{k}_{1}, {k}_{2}} $代入求解$ {L}_{k} $ ((6)式)的过程. 其中, (8)式或右边所使用的所有$ {Q}_{k}^{{k}_{1}\left(i\right), {k}_{2}\left(i\right)} $均为最近一次更新的结果. 另外, 对于一个给定的波数k, 所需的迭代次数是10次方量级, 最终要求前后两次迭代的结果差异小于千分之一. 由这种完备迭代法构建的程序, 在本文后面称之为“完备程序”.

对于高斯型, 方程(6)分子的第2项为0, (8)式等号右边第2项也为0. 这是因为^[23]斜度 (skewness)$ {\left(\langle{{x}^{3}(s)}\rangle/\langle{{x}^{2}(s)}\rangle\right)}^{3/2} $为0, 所以 $ \langle{X}_{k} {X}_{{k}_{1}}^{*} {X}_{{k}_{2}}^{*}\rangle = 0$. 这样, 可以直接解出$ {L}_{k} $. 高斯型系统可专门构建相对简单的计算传递函数的程序, 本文后面称之为“高斯输入型程序”.

在计算传递函数时, 可以通过计算相干系数来测试模型拟合的性能好坏. 相干系数定义如下:

线性相干性(linear coherency) $ {\gamma }_{{\mathrm{L}}}^{2}\left(k\right) $和二次相干性$ {\gamma }_{{\mathrm{Q}}}^{2}\left(k\right) $分别代表着线性传递函数与二次传递函数带来的输出功率份额. 线性传递和二次传递的关联所带来的输出功率成份由混合相干性$ {\gamma }_{{\mathrm{L}}{\mathrm{Q}}}^{2}\left(k\right) $表示. 注意, 仅当输入是非高斯型时, 才会出现项$ {\gamma }_{{\mathrm{L}}{\mathrm{Q}}}^{2}\left(k\right) $. 当输入为高斯型时, $ {\gamma }_{{\mathrm{L}}{\mathrm{Q}}}^{2}\left(k\right)\equiv 0 $. $ {\gamma }^{2}\left(k\right) $表示总相干性. 而$ {\gamma }_{{\mathrm{E}}}^{2}\left(k\right) $表示由于噪声等其他所有误差因素所造成的输出功率份额.

对两类近高斯激励系统进行仿真, 应用程序对模拟的非线性系统进行模拟计算, 检验模型和算法的有效性. 最后对模型算法应用于边缘等离子体湍流情况进行对比分析.

将基于波数空间的模型转换到频率空间, 回到用(1)式表述系统, 即可用本模型来模拟海波激励-船舶摇摆响应现象. 这里需要得到与图1输入-输出信号谱类似的、接近相同的仿真信号, 只需保证仿真系统与原系统属于同一类(近高斯输入型)非线性系统.

根据图1, 随机海波浪高测量数据的功率谱类(近)高斯分布. 海波与船泊摇动的采样频率统一取为1 Hz, 采样的数据一般可以分为$ 12 $个段, 每个段包含$ 128 $个样本. 据此利用向湍流过渡(transition to turbulence)数据的生成方法, 即“仿真数据生成子系统”采用向湍流过渡阶段的特征参数^[29]来构建, 以高斯色噪声为初始输入, 通过3次黑盒子^[23](仿真数据生成子系统), 来生成仿真输入信号谱. 通过3次黑盒子的目的是为高斯色噪声信号加入一些非线性因素. 然后, 对谱作傅里叶逆变换, 可得近似的输入(仿随机海波波高)信号. 取得实际摇摆位形信号的近似值, 作为输出(仿摇摆位移)信号, 并作快速傅里叶变换得输出信号谱, 结果如图4所示.

图4(a)—(d)分别表示输入信号谱密度、近似输入(随机海波波高仿真)信号、输出信号谱密度和近似输出(摇摆位形仿真)信号. 需要注意的是, 图4是通过模拟手段最容易获得的、与海浪激励船舶摇摆系统最接近的一个近高斯输入-响应非线性系统. 显然, 图4(d)所示船舶摇摆信号与图1(c)所示的测量信号基本一致, 但有稍许误差, 这就造成图4(c)所示船舶摇摆响应信号自功率谱与图1(d)稍有差别, 即在$ f=0.18 $Hz附近没有相干峰. 不过该相干峰幅度$ \sim O\left({10}^{-2}\right) $, 相对于最大谱功率($ \sim O\left({10}^{2}\right) $)来说非常小. 因此该仿真系统与真实系统在建模分析中定性上是一致的, 基于仿真系统的分析结果可以定性模拟真实系统的建模情况.

图4(a)所示海浪输入自功率谱图表明大多数入射海浪功率包含在0.13—0.21 Hz的频带内, 而图4(c)所示摇摆响应发生在约小于0.01 Hz的低频率. 为了测试信号的高斯性质, 计算入射海浪输入的互双谱. 高斯信号的互双谱应都为0^[2]. 图5(a), (b)显示了互双谱$ \left|\left\langle{{Y}_{i+j}{X}_{i}^{{\mathrm{*}}}{X}_{j}^{{\mathrm{*}}}}\right\rangle\right| $幅度的三维图和等高线图. 其中, $ {X}_{i}={X}_{{f}_{1}} $表示在频率$ {f}_{1} $处的输入谱. 同理, $ {X}_{j}={X}_{{f}_{2}} $表示在频率$ {f}_{2} $处的输入谱, $ {Y}_{i+j} $表示$ f={f}_{1}+{f}_{2} $处的输出谱. 在大约$ {f}_{1}= {f}_{2} = 0.17 $ Hz处出现的大峰值表明入射海浪中存在二次谐波分量, 这反过来导致非高斯统计.

图6所示为“完备”程序(即输入信号可以是 高斯、近高斯或非高斯信号)和高斯输入型程序(即输入信号只能为高斯信号)的计算结果. 其中 包含了线性传递函数$ {L}_{f} $的实部和虚部、线性传递函数的模$ \left|{L}_{f}\right| $、二次非线性传递函数模$ |{Q}_{f}^{{f}_{1}, {f}_{2}}| $ 的等高线图和三维结构图. 计算中的线性传递函数和二次传递函数收敛性如图7所示. 图示在迭代进行到大约第10步时, 计算结果已经达到收敛状态.

比较图6的绿线(对应非高斯型结果)及黑色点线(对应高斯型结果)可见, 采用高斯输入型程序计算的线性和非线性传递函数的结构和大小高度类似. 这是由于所使用的输入信号谱是近高斯型的. 注意到图6(c)绿线所示的线性传递函数的模与文献[1]中图9(a)及文献[2]中图6所示的结果都类似, 其在数量上小的差异并不反映模型的对错优劣, 而纯粹是由于本文生成的仿真数据较实测输入与输出信号更简单造成的. 同时, 注意到图6(c), (d)所示的、用完备程序计算的线性传递函数与二次传递函数形状分别与文献[2]中图5与图6类似, 差异较小, 证明仿真数据与计算结果比较符合实测系统.

图8显示完备程序计算的二次传递函数幅度$ |{Q}_{f}^{{f}_{1}, {f}_{2}}| $等位线图. 其中黑色线表示对角线$ {f}_{1}+ {f}_{2}=0 $. 红色线显示模型具有沿线 $ {f}_{1}+{f}_{2}=0.008 $ Hz的振幅峰值, 这是驳船摇摆响应峰值的频率. 这意味着低频下的大输出功率是由于二次波相互作用以及船舶系泊系统的线性响应的影响. 二次波相互作用使得入射海浪中差异等于0.008 Hz的所有成对频谱分量都有助于低频驳船摇摆响应.

为测试这种仿真的“拟合性能”, 利用(10)式—(14)式计算相干性. 所估计的总相干性与1相差值, 可看作是由于估计方法的系统误差以及由于统计量的估计值的方差而产生的. 图9(a)显示了应用完备程序计算的线性、二次和混和相干性对模型的总相干性的贡献. 结果表明, 线性相干性$ {\gamma }_{{\mathrm{L}}}^{2} $是小的, 但二次相干性$ {\gamma }_{{\mathrm{Q}}}^{2}\left(f\right) $相对很大. 这表明由二次传递函数带来的输出功率份额最多. 模型的总相干性$ {\gamma }^{2}\left(f\right) $接近于$ 1 $, 这表明该方法的总误差较小. 由此可推断, 此迭代方法能够对数据产生好的拟合. 另一方面, 应用高斯型输入类程序即考虑输入信号$ {X}_{f} $为高斯型时, 可以采用非迭代方法. 由于输入信号不是纯粹的高斯型, 结果导致了二次相干系数恒等于1, 暗示二次传递函数与真实值可以存在大的差异, 它导致总相干性的值大于1, 如图9(b)所示. 可见, 这里人为假定输入信号为高斯型的, 所引入的系统误差会很大. 因此, 近高斯激励-响应系统不能完全当作高斯系统进行处理. 特别是当输出信号谱分布接近于输入信号时(如下面将讨论的充分发展湍流问题), 计算结果对模型系统误差更加敏感.

文献[2]中图8也给出了基于实测数据的真实模型的谱相干系数. 其中, 文献[2]给出非高斯型的总相干系数值几乎等于1, 而高斯型的总相干系数普遍高于1, 其值介于1与3之间, 与本文图9模拟结果基本吻合. 文献[1]图8(a), (b)分别给出了线性传递函数的模$ \left|{L}_{f}\right| $和线性相干性(coherence)的平方值, 分别对应于本文图6(c)和图9(a)中$ {\gamma }_{{\mathrm{L}}}^{2} $. 本文计算的$ \left|{L}_{f}\right| $和$ {\gamma }_{{\mathrm{L}}}^{2} $与文献[1]的结果基本一致. 注意$ {\gamma }_{{\mathrm{L}}}^{2} $由$ \left|{L}_{f}\right| $通过(10)式进行定义. $ \left|{L}_{f}\right| $描述了计算结果的线性内容(线性传递函数), $ {\gamma }_{{\mathrm{L}}}^{2} $描述了线性项带给系统的贡献. 计算结果的合理性说明本文基于模拟数据所建立的模型与所用的模型求解算法对于这类近高斯输入型二次非线性系统的特性描述与系统分析是合适的.

当应用的是完备程序, 则利用激励方数据(输入信号谱), 通过“黑盒”一次, 可以得到输出信号谱, 如图10(a)所示. 以此输出谱作傅里叶逆变换得到摇摆位形, 如图10(b)所示. 当应用的是高斯型程序, 则用上述类似方法, 复验得到摇摆位形, 如图10(c)所示. 很明显, 图10(b), (c)所示摇摆位形之间有一定误差, 但都与用作正程序响应方数据的原摇摆信号谱有差距. 这是由于原信号变换中加了汉宁窗. 如果在对图10(b)所示输出信号谱作傅里叶逆变换时, 去除汉宁窗, 得到的摇摆信号就完全与初始信号相同, 如图10(d)所示. 这证明了我们的完备程序能正确计算随机海浪激励-船舶响应非线性系统的线性与非线性传递函数.

注意, 上述复验过程实际上可以看作是“黑盒子”的一次应用. 即根据测试得到的海浪波高信号来计算、预估系泊船舶摇摆位形.

在研究充分发展湍流内部能量传递问题中, 为了模拟充分发展湍流, 需要利用特定的、具有流体向湍流过渡阶段特征的线性与二次非线性传递函数. 即定义$ {L}_{k} $和$ {Q}_{k}^{{k}_{1}, {k}_{2}} $如下:

其中, $ k={k}_{1}+{k}_{2}, {k}_{1}\geqslant {k}_{2} $. $ \hat{k}={k}_{{\mathrm{N}}{\mathrm{y}}{\mathrm{q}}} $表示尼奎斯特波数.

如此定义的二次传递函数的振幅和形状是任意选择的, 其值与由Hasegawa-Mima方程^[29]所预测的振幅有相同的阶次, 因此是合理的. 线性传递函数按(15)式定义, 使得输入谱与输出谱在形状上是类似的, 而这正是研究稳态问题中所需要的. 为了方便起见, 考虑将(15)式和(16)式定义的波数空间问题转换到频率空间问题, 即以频率代替波数进行湍流的模拟. 实验上, 两个空间上分离的探针能够测量波动的时间变化. 我们不计算满足选择规则$ k = {k_1} + {k_2} $的不同波数的双谱, 而是对于频率成分$ f = {f_1} + {f_2} $的双谱进行计算. 设定采样频率为1 MHz, 模拟在空间相距$ \Delta x $的两个点连续采集100000组数据. 然后, 数据预处理中设定1280个实现, 每个实现包含128个样本数据. 数据的获取方式是, 将类似图3所示的5个相同的“黑盒子”类型连成一个序列^[23], 应用高斯色噪声信号作为第1个黑盒子的输入. 由于黑盒子的非线性属性, 其输出现在变成非高斯型, 第1个盒子的输出成为第2个黑盒子的输入. 依次类推. 仿真中, 利用第5个黑盒子的输入与输出作为模拟湍流系统最终的输入与输出, 如图11所示. 可见, 输入信号$ {X}_{k} $与输出信号$ {Y}_{k} $一样可以看作是近高斯型的.

图11所示充分发展湍流输入自功率谱图表明大多数入射湍流功率包含在160—200 kHz的频带内, 接下来计算入射湍流的自双谱. 图12(a), (b)显示了互双谱幅度$ \left|\left\langle{{Y}_{i+j}{X}_{i}^{{\mathrm{*}}}{X}_{j}^{{\mathrm{*}}}}\right\rangle\right| $的三维图和等高线图, 可以看出在较宽的频带内均出现的大峰值, 表明入射湍流中存在二次谐波分量, 即存在非高斯成分.

应用完备程序计算传递函数, 结果如图13所示. 若未做特殊说明, 表明使用的都是完备程序. 这里, 线性传递函数的实部和虚部分别如图13(a), (b)所示, 二次传递函数的模如图13(c)的等高线所示. 图13(d)表示二次传递函数模的三维结构. 可以看到, 线性传递函数和二次传递函数与定义此湍流系统时的预设值((15)式和(16)式)相当吻合. 这表明我们的程序可以正确计算近高斯输入激励全发展湍流系统的黑盒子参数, 如传递函数等.

图14给出了模型的线性相干性、二次相干性和混合相干性对总相干性的贡献. 线性相干性贡献在1左右, 而二次相干性贡献很小, 在0左右. 混合相干性在0附近. 总误差$ {\gamma }_{{\mathrm{E}}}^{2}\sim0 $, 说明计算准确率很高. 总相干性$ {\gamma }^{2} $(红色实线所示)基本等于1, 表明该方法的系统误差较小. 由此可推断, 此迭代方法能够对充分发展湍流数据有好的拟合. 以输入和输出功率谱及其所计算的线性和二次非线性传递函数为基础, 可以得到相位差、平均色散和增长率, 如图15所示.

最后, 根据本项工作计算的线性传递函数和二次传递函数, 以图11所示输入谱信号为输入数据集, 反算出输出信号, 对应于模拟的输出测试信号, 这些输出信号部分如图16所示. 注意模拟采样频率为1 MHz. 图中给出了对第2个探针在两个时刻开始测量的一小段湍流信号, 信号长度各为128个数据点, 信号时长0.128 ms. 其中第1个模拟起测时刻$ {t}_{0}=0.125\;{\mathrm{s}} $, 第2个起测时刻的间隔时间$ \Delta t=0.125\;{\mathrm{ }}{\mathrm{s}} $,

(6)式—(8)式所给定的估计传递函数的方法, 在将随机输入信号假定为高斯型分布时, 在数学上可以得到相当大的简化. 当输入信号$ x\left( s \right) $为高斯型时, 可以直接确定传递函数$ {L}_{k} $和$ {Q}_{k}^{{k}_{1}, {k}_{2}} $, 但这样的方法只能用于分析由高斯噪声源作为外部激励的系统^[22], 以及输入信号可以假定为高斯型的系统. 可是有许多系统, 比如湍动流体和等离子体, 却不允许对输入信号作这样的限定性假设.

另外, 对于湍流研究来说, 当我们假定四波耦合以及更高阶次耦合过程比三波耦合弱得多时, (1)式—(3)式是波-波耦合有关的方程中的最简单变形. 而当某些情况下, 三波非线性耦合被系统的色散属性所禁止时, 比如像对深水中的表面重力波^[31], 就必须考虑更高阶次项了.

影响模型计算误差的因素主要有下面几方面. 首先是物理模型造成的误差. 比如湍流建模受到三波相互作用占比程度的影响. 如果三波相互作用在一个物理过程中并不占支配地位, 则在建模时, 必须考虑四波甚至更多波之间的相互作用, 否则会带来大的系统误差. 另外, 有的系统计算需要基于特定的假设. 例如, 边缘等离子体的增长率和能量传递函数计算, 需要基于“稳态湍流”这一假定^[26], 相应的湍流模型必须适应这一要求.

其次是计算方法带来的误差. 例如, 在本文模型应用中, 一般要求在作傅里叶变换前对信号数据加汉宁窗或海明窗, 以减少频谱能量泄漏. 否则, 不加窗(即只加矩形窗)会造成小频谱能量被忽视, 计算结果欠精确.

其他影响精度的因素还包括测量数据的代表性, 以及样本长度和样本总段数必须达到最低要求. 对于海波激励系统来说, 每段样本长度取128、共取12—24段是足够的. 但对于湍流研究来说, 样本长度至少128以上, 总样本段数要求达到800及以上^[32], 方能得到比较精确的谱, 以保证计算精度.

需要说明的是, 在模拟充分发展湍流并应用本文的分析方法时, 不必局限于某类湍流, 只要保证湍流处于稳态(stationary)即可. 但不同形式的湍流会对计算结果(传递函数)产生影响. 因此, 需要针对不同的研究对象, 采取针对性措施, 包括附设额外条件, 以使分析对象限定在特定的物理范畴内, 避免计算结果失真或失去物理意义. 比如, 在研究边缘等离子体湍流的内部能量传递过程时, 需要限定系统输入与输出满足所谓的Kim条件^[26], 即: $\langle{Y}_{k}{Y}_{k}^{*}\rangle=\langle{X}_{k}{X}_{k}^{*}\rangle $. 这个条件限定了能量只在满足耦合条件的波之间传递. 对于系泊船舶与波动海波构成的系统, 海波不是充分发展湍流, 在研究船舶随海波的漂移时, 只要海波波动是随机的即可.

本文从仿真研究随机海浪激励系泊船舶漂移问题和一类充分发展湍流入手, 提出了一种通用的特定湍动激励-响应型二次非线性系统的建模方法. 作为输入的这种湍流的特点是其波动信号功率谱分布接近高斯分布. 研究基于Choi等^[1]和Kim等^[2]早年各自提出的随机海浪激励-系泊船舶漂移响应问题的建模方法, 本文改进、完善了这些方法, 对以随机海浪激励系泊船舶漂移响应和一类充分发展湍流为代表的非线性系统进行了扩展的、系统的建模分析, 首次应用完备迭代法进行了模型求解.

在仿真数据构建中, 我们都采用了具有向湍流过渡阶段特征的传递函数^[29]用作“仿真数据生成子系统”的参数. 然后应用算法(程序)对上述两个模拟系统进行了验证计算. 计算中重现了文献[1, 2]的结果, 模型的正确性和算法的有效性得到证明. 相干性分析表明, 随机海浪-船舶摇动系统二次相干性远大于线性相干性, 但近高斯输入型充分发展湍流的线性相干性更大. 研究表明, 用于求解该模型传递函数参数的“完备迭代法”适用于随机海波激励系泊船舶摇摆响应问题和充分发展湍流问题的计算. 研究结果为以近高斯输入谱为特征的特定湍动激励-响应类二次非线性系统的模型描述及模型求解提供了一种有效途径.