一种简单高效的迭代翻译近似方法，用于模拟平稳的非高斯随机向量过程

本文提出了一种改进的迭代翻译近似法，通过引入Hermite多项式展开降低计算复杂度，并直接计算最近正定半定矩阵解决非高斯随机向量过程的协方差不兼容问题，验证了该方法在计算效率（提升数百倍）和适用性（支持复数值协方差矩阵）上的优势。

　　

郭朝明|张楠|孙启凯|顾庆阳

北京交通大学土木工程学院，中国北京100044

摘要

基于Cholesky分解的迭代平移近似方法（ITAM）由于其高计算精度和解决不兼容性的强大能力，被广泛用于模拟平稳非高斯随机向量过程。然而，其局限性包括计算效率低、由于Cholesky分解导致的鲁棒性差、结果不唯一，以及无法处理复数值的交叉谱密度矩阵（CSDM），这些限制了其在工程应用中的进一步发展。为了解决这些问题，本研究提出了一种简单高效的ITAM方法，该方法结合了Hermite多项式展开来减少计算开销，并直接计算最接近的正半定高斯CSDM以解决不兼容性问题。通过四个数值示例验证了所提方法的有效性和优越性。结果表明，所提方法的计算效率提高了数百倍，其精度与现有方法相当或更好，并且适用于复数值CSDM。

引言

平稳非高斯随机过程的模拟是土木工程中的一个关键课题，在可靠性分析和不确定性量化中有广泛应用[1]、[2]、[3]、[4]、[5]。在过去的几十年中，特别是谱表示方法（SRM）由于具有高精度和高效性而变得成熟[6]、[7]、[8]、[9]。与高斯过程不同，模拟非高斯随机过程更具挑战性，因为样本统计量必须同时满足规定的功率谱密度函数（PSDF）（或相关函数）和规定的边际分布函数（或更高阶矩）[10]。

为了解决这一挑战，已经开发了两种主要方法。第一种方法将传统的SRM扩展到更高阶，从而能够直接模拟非高斯过程。然而，由于缺乏多谱函数，这种方法的适用性受到限制[11]、[12]、[13]、[14]、[15]。第二种方法基于规定的非高斯PSDF和边际分布函数（或更高阶矩），通过相关畸变来确定一个有效的高斯PSDF（称为潜在的高斯PSDF），然后使用SRM生成潜在的高斯样本，再利用边际分布函数（或更高阶矩）将其映射为非高斯样本。这种方法的关键在于在畸变相关性的同时解决不兼容性问题[16]。

这种相关畸变方法的一个主要实现是多项式类型方法，它基于更高阶矩建立显式表达式来完成相关畸变[17]、[18]、[19]、[20]、[21]。虽然这种方法高效且准确，但由于对更高阶矩范围的严格限制，其适用性受到限制。另一种重要的相关畸变实现是基于平移过程理论的迭代平移近似方法（ITAM）。ITAM使用等概率平移，迭代地确定从非高斯PSDF到潜在高斯PSDF的畸变。这种方法因其高精度、广泛的适用性和确定极值分布的能力而受到青睐[16]。

早期的ITAM主要针对随机标量过程，但它们的计算精度不高，并且每次迭代都需要通过SRM生成样本，导致计算效率低下[22]、[23]、[24]。无样本ITAM的引入通过消除迭代中的SRM提高了精度并显著提高了计算效率[10]、[25]。后来，针对标量过程的无样本ITAM自然扩展到了向量过程——基于Cholesky分解的ITAM被提出，该方法使用Cholesky分解来确保潜在高斯交叉谱密度矩阵（CSDM）的有效性，这是一个里程碑[26]。

然而，尽管仍然不需要样本，基于Cholesky分解的ITAM的计算效率已不再符合预期。因为与随机标量过程不同，高维向量过程需要更多的平移过程计算，而且这些过程中的相关畸变会产生显著的计算开销。此外，Cholesky分解容易因数值误差而失败，降低了方法的鲁棒性。最终结果的不唯一性也复杂化了实际应用。最后，基于Cholesky分解的ITAM无法处理复数值CSDM，进一步限制了其在工程中的应用。

尽管基于Cholesky分解的ITAM后来被扩展到多维情况，并考虑了波浪效应和随机波的模拟，但上述问题仍未得到解决[27]、[28]、[29]。最近，基于特征值分解的ITAM用特征值分解取代了Cholesky分解，但它没有提高计算效率，并且数值稳定性较差[30]。

显然，提高基于Cholesky分解的ITAM的计算效率并避免Cholesky分解引起的问题是一个迫切的需求。基于此，提出了一种简单高效的ITAM方法用于模拟非高斯随机向量过程。本文的其余部分组织如下：第2节介绍平移向量过程理论。第3节回顾基于Cholesky分解的方法。第4节介绍所提出的方法。第5节通过四个数值示例验证所提方法的有效性和优越性。第6节总结本文。

节选

平移向量过程理论

考虑一个均值为零、标准差为单位且具有m个变量的平稳高斯随机向量过程Y(t)=?[Y₁(t), Y₂(t), …, Y_m(t)]^T，其CSDM表示为S_G(ω)，其交叉相关矩阵（CCM）表示为R_G(τ)：$R_{Gij}\left(\tau \right)=E\left[Y_i\left(t\right)Y_j(t+\tau )\right]$

R_G(τ)可以通过对S_G(ω)应用Wiener-Khinchin变换得到：$R_{Gij}\left(\tau \right)=E\left[Y_i\left(t\right)Y_j(t+\tau )\right]$

$R_{Gij}\left(\tau \right)={\int }_{?\infty }^{+\infty }$${\rho }_{ij}\left(\tau \right){\sigma }_{Gi}\sigma Gj$

i, j?=?1, 2, …, m。

R_G(τ)可以通过对S_G(ω)应用Wiener-Khinchin变换得到：$R_{Gij}\left(\tau \right)={\int }_{?\infty }^{+\infty }$

ρij(τ)σGiσGj