(19)中华 人民共和国 国家知识产权局 (12)发明 专利申请 (10)申请公布号 (43)申请公布日 (21)申请 号 202111497671.X (22)申请日 2021.12.09 (71)申请人 河南理工大 学 地址 454000 河南省焦作市高新区世纪大 道2001号 (72)发明人 任安康　徐克科　 (74)专利代理 机构 郑州立格知识产权代理有限 公司 41126 代理人 崔卫琴 (51)Int.Cl. G06F 30/20(2020.01) G06F 17/16(2006.01) G01V 1/30(2006.01) G01V 1/36(2006.01) G06F 119/10(2020.01) (54)发明名称 一种有色噪声中精确估计震后松弛时间的 方法 (57)摘要 本发明公开了一种有色噪声中精确估计震 后松弛时间的方法， 包括以下步骤： 选择白噪声 和闪烁噪声的组合， 使用最大似然技术估计震前 GNSS时间序列， 获取震前信号； 删除地震后GNSS 时间序列中的地震前信号； 构建震后GNSS时间序 列的随机模型。 根据实际需要设定阻尼因子λ； 使用线性最小二乘法计算地震后 变形的振幅Al、 Ae和χ2(X)； 计算法线矩阵N和L， 使用线性最小 二乘法计算Al、 Ae和χ2(X+δX)； 当χ2(X+δX)迭 代次数减少到λ可忽略的程度时， 停止迭代。 本 发明既不要求初值的准确性， 又能解决有色噪声 的干扰， 且算法的自适应性显著提高， 估计松弛 时间的结果基本不受初值影 响， 进而有效地规避 信号提取所携带的大量 误差。 权利要求书2页 说明书7页 附图1页 CN 114169171 A 2022.03.11 CN 114169171 A 1.一种有色噪声中精确估计震后松弛时间的方法， 其特 征在于： 包括以下步骤： 1)选择白噪声和闪烁噪声 的组合， 使用最大似然技术估计震前GNSS时间序列， 从而获 取震前信号； 2)删除地震后GNS S时间序列中的地震前信号； 3)利用震前噪声参数构建震后GNSS时间序列的随机模型； 其中， 震前噪声参数为步骤 1)中利用震前GNS S时间序列计算出的不同噪声类型之间的比例； 4)根据实际需要设定阻尼因子 λ； 5)给出拟合参数集的松弛时间初始猜测X， 使用线性最小二乘法计算地震后变形的对 数型震后衰减函数振幅Al、 指数型震后衰减函数振幅Ae和改正前松弛时间卡方 χ2(X)； 6)计算法线矩阵N和L， 使用线性最小二乘法计算Al、 Ae和改正后松弛时间卡方χ2(X+δ X)； N为正规矩阵， 向量L＝ATC‑1V， C为白噪声和闪烁噪声组合模型， A为设计矩阵， V残差向 量； 7)如果 χ2(X+δ X)≥χ2(X)， λ增加10倍， 则修改正 规矩阵并返回步骤6)； 8)如果χ2(X+δX)＜χ2(X)， 使用更新向量X+δX， λ减少10倍， 修改正规矩阵， 然后返回步骤 6)； 9)当χ2(X+δ X)迭代次数减少到 λ可忽略的程度时， 停止迭代； 10)最终估计的参数向量X， 即为震后松弛时间。 2.根据权利要求1所述的有色噪声中精确估计震后松弛时间的方法， 其特征在于： 步骤 1)中地震前信号包括台站位置初值、 GNS S站速率、 阶跃项、 周期信号和偏移。 3.根据权利要求1所述的有色噪声中精确估计震后松弛时间的方法， 其特征在于： 步骤 3)随机模型计算如下： GNSS时间序列的统计模型描述 为白噪声和闪烁噪声的组合， 其形式如下： C＝σ2(cos2(φ)I+si n2(φ)Jk)     (1) 其中， 验后方差σ 由残差直接计算， 噪声比例参数φ使用极大似然估计计算得到， I为单 位矩阵， Jk为有色噪声对应的协方差矩阵； Jk通过转换矩阵T转换 得到， 具体方法是 Jk＝TTT        (2) 而且 其中， ψn为转换矩阵T中的元素， ΔT 为采样间隔， k为噪声管谱指数， Γ( ·)为伽玛函数， n！ 为n的阶乘；权　利　要　求　书 1/2 页 2 CN 114169171 A 2震后形变模型用下述公式描述： Al,Ae, τl， τe通过拟合一条震后位移曲线得到， 并使用迭代非线性最小二乘法来估计这 些参数； 其中， δXpsd(ti)为ti时刻的震后形变位移量， τl和 τe分别为对数型、 指数型震后 松弛 时间， T为 地震发生时刻。 4.根据权利要求3所述的有色噪声中精确估计震后松弛时间的方法， 其特征在于： 步骤 5)计算地震后变形振幅具体为： 为了获得弛豫时间的整体解， 估计松弛时间后， 地震后变形的振幅用线性最小二乘法 估计， 因此仅对松弛时间应用泰勒级数， 将这个系列的扩展 写成： 利用公式(6)中给 出的参数和 测量值之间的关系， 形成一组观测方程： AX＝V, 其中， A为关于τl和 τe的函数， X＝[dτl dτe]T是参数的更新向量， V是观测 到的和计算出 的地震后位移之间的差值； 基于最小二乘法创建法方程： NX＝L     (7) 其中， 正规矩阵为 N＝ATC‑1A， 向量为L＝ATC‑1V， C为白噪声和闪烁噪声组合模型。 5.根据权利要求4所述的有色噪声中精确估计震后松弛时间的方法， 其特征在于： 步骤 6)具体包括： 利用Levenberg ‑Marquardt方法解决非线性参数估计， 由于法向矩阵仅给出特定点的 趋势， 而不给 出斜率延伸的程度， 用以下规则定义 N的新矩阵N'＝[n ′ij]替换正规矩阵： 其中， λ非常大， 矩阵N'被迫接近对角， 并且由于正规矩阵的对角元素被放大， 因此逼近 解的步长被减小； 另一方面， 当 λ接 近零时， N'将接 近真实法线矩阵。权　利　要　求　书 2/2 页 3 CN 114169171 A 3