CN114139363B - Method for simulating space variation non-stationary earthquake motion time course - Google Patents

Method for simulating space variation non-stationary earthquake motion time course Download PDF

Info

Publication number
CN114139363B
CN114139363B CN202111410406.3A CN202111410406A CN114139363B CN 114139363 B CN114139363 B CN 114139363B CN 202111410406 A CN202111410406 A CN 202111410406A CN 114139363 B CN114139363 B CN 114139363B
Authority
CN
China
Prior art keywords
frequency
time
points
point
simulation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202111410406.3A
Other languages
Chinese (zh)
Other versions
CN114139363A (en
Inventor
赵宁
李小龙
李钰岩
陆哲明
陈晓伟
徐志龙
宋鑫
张耀文
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Sichuan Agricultural University
Original Assignee
Sichuan Agricultural University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Sichuan Agricultural University filed Critical Sichuan Agricultural University
Priority to CN202111410406.3A priority Critical patent/CN114139363B/en
Publication of CN114139363A publication Critical patent/CN114139363A/en
Application granted granted Critical
Publication of CN114139363B publication Critical patent/CN114139363B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/14—Fourier, Walsh or analogous domain transformations, e.g. Laplace, Hilbert, Karhunen-Loeve, transforms
    • G06F17/141—Discrete Fourier transforms
    • G06F17/142—Fast Fourier transforms, e.g. using a Cooley-Tukey type algorithm
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/02—Reliability analysis or reliability optimisation; Failure analysis, e.g. worst case scenario performance, failure mode and effects analysis [FMEA]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Computational Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Computing Systems (AREA)
  • Discrete Mathematics (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

The invention discloses a method for simulating a spatial variation non-stationary earthquake motion time course, which relates to the technical field of civil engineering earthquake-proof design and has the technical scheme key points that: s1, determining frequency interpolation points based on given seismic motion coherence function gamma (omega) and power spectral density function S (omega, t)
Figure DDA0003373540690000011
And time-frequency interpolation point
Figure DDA0003373540690000012
S2, calculating frequency interpolation point
Figure DDA0003373540690000013
Delayed coherence matrix for n-point seismic field to be simulated
Figure DDA0003373540690000014
Computing time-frequency interpolation points
Figure DDA0003373540690000015
Self-evolving power spectrum at n earthquake motion simulation points
Figure DDA0003373540690000016
(j ═ 1,2, …, n); s3, pair
Figure DDA0003373540690000017
Performing Cholesky decomposition to obtain
Figure DDA0003373540690000018
Will be provided with
Figure DDA0003373540690000019
(j ═ 1,2, …, n) decomposed into principal coordinates
Figure DDA00033735406900000110
And feature vectors
Figure DDA00033735406900000111
S4, is prepared from
Figure DDA00033735406900000112
And
Figure DDA00033735406900000113
interpolation results in B (omega),
Figure DDA00033735406900000114
And
Figure DDA00033735406900000115
and S5, generating the spatial multi-point seismic motion samples by using a random simulation formula based on a spectrum representation method and utilizing Fast Fourier Transform (FFT). The invention can avoid the problems of large power spectrum matrix decomposition and time frequency function decoupling calculation amount and large computer storage space requirement caused by a plurality of simulation points, time points and frequency points in the complete non-stable earthquake motion simulation, and meets the precision requirement of the simulation.

Description

Method for simulating space variation non-stationary earthquake motion time course
Technical Field
The invention relates to the technical field of civil engineering earthquake-resistant design, in particular to a method for simulating a spatial variation non-stationary earthquake motion time course.
Background
According to the requirement of earthquake-proof specification, when a large-span structure is designed for earthquake resistance, nonlinear power time-course analysis under the action of space variation earthquake motion is required. However, the actual measurement seismic record is often insufficient, so that a large amount of seismic motion time courses reflecting the real seismic motion characteristics need to be obtained through a numerical simulation mode.
At present, the simulation aiming at the spatial variation non-stationary earthquake motion time course is mainly based on a non-stationary power spectrum model and a coherent function model for representing earthquake motion characteristics, and a spectrum representation method is adopted for simulation. In order to improve the simulation efficiency of the spatial multi-point seismic motion field, the method generally adopts the means of intrinsic orthogonal decomposition (POD) and the like to decouple the completely non-stationary power spectrum function, and then introduces Fast Fourier Transform (FFT). However, for the case that the number of analog points is large and the number of time discrete points and frequency discrete points is large, the decoupling calculation requires a large amount of computer memory and time. This results in inefficient simulation of spatially variant non-stationary seismic motion fields.
Therefore, it is necessary to provide a more efficient method for modeling the spectral representation of the spatially variant non-stationary seismic motion field to solve the above problems.
Disclosure of Invention
The invention aims to provide a method for simulating a spatial variation non-stationary earthquake motion time course aiming at the technical problems.
The technical purpose of the invention is realized by the following technical scheme: a method for simulating a spatial variation non-stationary seismic motion time course specifically comprises the following steps:
s1, determining frequency interpolation points based on given seismic motion coherence function gamma (omega) and power spectral density function S (omega, t)
Figure BDA0003373540670000021
And time-frequency interpolation point
Figure BDA0003373540670000022
S2, calculating frequency interpolation point
Figure BDA0003373540670000023
Delayed coherence matrix for n-point seismic field to be simulated
Figure BDA0003373540670000024
Computing time-frequency interpolation points
Figure BDA0003373540670000025
Self-evolving power spectrum at n earthquake motion simulation points
Figure BDA0003373540670000026
S3, pair
Figure BDA0003373540670000027
Performing Cholesky decomposition to obtain
Figure BDA00033735406700000222
By intrinsic orthogonal decomposition (POD)
Figure BDA0003373540670000028
Decomposed into a small number of principal coordinates
Figure BDA0003373540670000029
And feature vectors
Figure BDA00033735406700000210
S4, is prepared from
Figure BDA00033735406700000211
And
Figure BDA00033735406700000212
interpolation results in B (omega),
Figure BDA00033735406700000213
And
Figure BDA00033735406700000214
and S5, generating the spatial multi-point seismic motion samples by using a random simulation formula based on a spectrum representation method and utilizing Fast Fourier Transform (FFT).
Further, the specific method of step S1 is:
(1) determining frequency interpolation points
Figure BDA00033735406700000215
Selecting uniformly distributed frequency interpolation points:
Figure BDA00033735406700000216
wherein, ω is1And ωuRespectively a first frequency point and a last frequency point,
Figure BDA00033735406700000217
the number of frequency interpolation points is satisfied
Figure BDA00033735406700000218
N is the total number of frequency discrete points;
(2) determining time-frequency interpolation points
Figure BDA00033735406700000219
Selecting time-frequency interpolation points with uniform distribution
Figure BDA00033735406700000220
The determination method is as follows:
Figure BDA00033735406700000221
Figure BDA0003373540670000031
in the formula, T0Is the total duration of the power spectrum;
Figure BDA0003373540670000032
and
Figure BDA0003373540670000033
the number of time-frequency interpolation points along the frequency direction and the time direction respectively; and is
Figure BDA0003373540670000034
M is the total number of time discrete points.
Further, the specific method of step S2 is:
calculating a delay coherence matrix by substituting the value of the frequency interpolation point into a coherence function:
Figure BDA0003373540670000035
then, substituting the time and frequency values of the time-frequency interpolation points by adopting the same mode to obtain the self-evolution power spectrum of each earthquake motion simulation point
Figure BDA0003373540670000036
Further, the specific method of step S3 is:
(1) for is to
Figure BDA0003373540670000037
Performing Cholesky decomposition to obtain
Figure BDA0003373540670000038
And will be
Figure BDA0003373540670000039
The decomposition is as follows:
Figure BDA00033735406700000310
in the formula, the superscript T represents the transpose of a matrix or vector,
Figure BDA00033735406700000311
is the lower triangular matrix, expressed as:
Figure BDA00033735406700000312
(2) each is decomposed by POD
Figure BDA00033735406700000313
Obtain corresponding main coordinates
Figure BDA00033735406700000314
And feature vectors
Figure BDA00033735406700000315
(3) At all time-frequency interpolation points
Figure BDA00033735406700000316
The values are expressed as a constant matrix, i.e.:
Figure BDA0003373540670000041
each column of the constant matrix is regarded as a column vector
Figure BDA0003373540670000042
A column vector decomposed by the following eigenvectors:
Figure BDA0003373540670000043
finding an optimal set of orthonormal bases
Figure BDA0003373540670000044
The projection of the column vector will be maximized thereon, where ΦqIs the q-th feature vector; lambda [ alpha ]qIs the qth eigenvalue; r is this
Figure BDA0003373540670000045
The correlation matrix of the column vectors can be calculated by:
Figure BDA0003373540670000046
in obtaining
Figure BDA0003373540670000047
After each feature vector, the projection of each column vector, i.e., the principal coordinate, is determined by the following equation:
Figure BDA0003373540670000048
wherein, aqIs the q projection vector;
(4) and performing descending order recombination on the characteristic values, reserving a low-order characteristic value containing more energy, and approximately expressing L as:
Figure BDA0003373540670000049
wherein N isΦThe number of low order terms selected to meet the accuracy requirement. Based on the above formula, each
Figure BDA00033735406700000410
Using its principal coordinates
Figure BDA00033735406700000411
And feature vectors
Figure BDA00033735406700000412
Expressed as:
Figure BDA00033735406700000413
wherein the content of the first and second substances,
Figure BDA0003373540670000051
and
Figure BDA0003373540670000052
are respectively a vector aqAnd phiqA corresponding discrete function.
Further, the specific method of step S4 is:
at the point of finding the interpolation
Figure BDA0003373540670000053
And
Figure BDA0003373540670000054
then B (omega) at other time-frequency coordinates,
Figure BDA0003373540670000055
And
Figure BDA0003373540670000056
and (3) performing approximate calculation by adopting an interpolation technology, and approximately reconstructing results at the original time and frequency coordinates according to an interpolation result, namely:
Figure BDA0003373540670000057
further, the specific method of step S5 is:
(1) based on a spectral representation method, the simulation formula of any random process is as follows:
Figure BDA0003373540670000058
wherein t is a time coordinate; omegalIs the l-th frequency point; Δ ω is the frequency interval; e is an exponential function; i is an imaginary unit; phi is aklIs a set of random phase angles, a deterministic quantity for a particular sample; re represents the real part of the acquired complex number; hjk(ω, t) is the element of the decomposed spectral matrix, calculated by:
Figure BDA0003373540670000059
wherein, thetajk(ω) is HjkThe phase angle of (ω, t), which is equal to the coherence function γ for seismic oscillationsjk(ω) a phase angle;
(2) substituting equations (13) and (15) for equation (14) yields the following simulation equation:
Figure BDA00033735406700000510
and by exchanging the summation order:
Figure BDA0003373540670000061
then n is performed on the formula (17)ΦAnd performing a secondary FFT operation, namely generating seismic motion samples of n points in space.
In conclusion, the invention has the following beneficial effects:
1. the method simulates the complete non-stationary earthquake motion, only needs Cholesky decomposition at a few frequency interpolation points, and the results of the rest frequency points are obtained approximately through interpolation, thereby greatly simplifying the calculation of the spectral matrix decomposition. In addition, POD decoupling is only carried out on the power spectrum functions at a few time-frequency interpolation points, and decoupling results of other time-frequency points are obtained approximately through interpolation, so that the calculation amount of time-frequency function decoupling is greatly reduced. Finally, the number of FFT operations is further reduced. The efficiency of the whole earthquake dynamic simulation process is greatly improved under the condition of keeping enough precision.
2. The method can provide the input earthquake dynamic load time course for the time course analysis of the earthquake-resistant design of large-span structures, such as a transmission tower line system, an oil pipeline and the like.
Drawings
FIG. 1 is a flow chart in an embodiment of the invention;
FIG. 2 is a graph of a correlation function validation graph for simulating seismic oscillation in an embodiment of the present invention;
Detailed Description
The present invention is described in further detail below with reference to fig. 1.
Example (b): a method for simulating a spatially variant non-stationary seismic motion time course is shown in FIG. 1, and specifically comprises the following steps:
s1, determining frequency interpolation points based on given seismic motion coherence function gamma (omega) and power spectral density function S (omega, t)
Figure BDA0003373540670000062
And time-frequency interpolation point
Figure BDA0003373540670000063
The specific method comprises the following steps:
(1) determining frequency interpolation points
Figure BDA0003373540670000071
Selecting uniformly distributed frequency interpolation points:
Figure BDA0003373540670000072
wherein, ω is1And omegauRespectively a first frequency point and a last frequency point,
Figure BDA0003373540670000073
the number of frequency interpolation points is satisfied
Figure BDA0003373540670000074
N is the total number of frequency discrete points;
(2) determining time-frequency interpolation points
Figure BDA0003373540670000075
Selecting evenly distributed time-frequency interpolation points
Figure BDA0003373540670000076
The determination is as follows:
Figure BDA0003373540670000077
Figure BDA0003373540670000078
in the formula, T0Is the total duration of the power spectrum;
Figure BDA0003373540670000079
and
Figure BDA00033735406700000719
the number of time-frequency interpolation points along the frequency direction and the time direction respectively; and is
Figure BDA00033735406700000710
M is the total number of time discrete points.
S2, calculating frequency interpolation point
Figure BDA00033735406700000711
Delayed coherence matrix for n-point seismic field to be simulated
Figure BDA00033735406700000712
Computing time-frequency interpolation points
Figure BDA00033735406700000713
Self-evolving power spectrum at n earthquake motion simulation points
Figure BDA00033735406700000714
The specific method comprises the following steps:
calculating a delay coherence matrix by substituting the value of the frequency interpolation point into a coherence function:
Figure BDA00033735406700000715
then, the same as described above is usedSubstituting the mode into the time and frequency values of the time-frequency interpolation point to obtain the self-evolution power spectrum of each earthquake motion simulation point
Figure BDA00033735406700000716
S3, pair
Figure BDA00033735406700000717
Performing Cholesky decomposition to obtain
Figure BDA00033735406700000718
By intrinsic orthogonal decomposition (POD)
Figure BDA0003373540670000081
Decomposed into a small number of principal coordinates
Figure BDA0003373540670000082
And feature vectors
Figure BDA0003373540670000083
The specific method comprises the following steps:
(1) to pair
Figure BDA0003373540670000084
Performing Cholesky decomposition to obtain
Figure BDA0003373540670000085
And will be
Figure BDA0003373540670000086
The decomposition is as follows:
Figure BDA0003373540670000087
in the formula, the superscript T represents the transpose of a matrix or vector,
Figure BDA0003373540670000088
is the lower triangular matrix, expressed as:
Figure BDA0003373540670000089
(2) each is decomposed by POD
Figure BDA00033735406700000810
Obtain corresponding main coordinates
Figure BDA00033735406700000811
And feature vectors
Figure BDA00033735406700000812
(3) At all time-frequency interpolation points
Figure BDA00033735406700000813
The values are expressed as a constant matrix, i.e.:
Figure BDA00033735406700000814
each column of the constant matrix is regarded as a column vector
Figure BDA00033735406700000815
A column vector decomposed by the following eigenvectors:
Figure BDA00033735406700000816
finding an optimal set of orthonormal bases
Figure BDA00033735406700000817
The projection of the column vector will be maximized thereon, where ΦqIs the q-th feature vector; lambda [ alpha ]qIs the qth eigenvalue; r is this
Figure BDA00033735406700000818
A correlation matrix of column vectors, which can be expressed byCalculating by the formula:
Figure BDA00033735406700000819
in obtaining
Figure BDA0003373540670000091
After each feature vector, the projection of each column vector, i.e., the principal coordinate, is determined by the following equation:
Figure BDA0003373540670000092
wherein, aqIs the q projection vector;
(4) and (3) performing descending reorganization on the characteristic values, reserving low-order characteristic values containing more energy, wherein L is approximately expressed as:
Figure BDA0003373540670000093
wherein, NΦThe number of low order terms selected to meet the accuracy requirement. Whereby each one
Figure BDA0003373540670000094
Using its principal coordinates
Figure BDA0003373540670000095
And feature vectors
Figure BDA0003373540670000096
Expressed as:
Figure BDA0003373540670000097
wherein, the first and the second end of the pipe are connected with each other,
Figure BDA0003373540670000098
and
Figure BDA0003373540670000099
are respectively a vector aqAnd phiqA corresponding discrete function.
S4, is prepared from
Figure BDA00033735406700000910
And
Figure BDA00033735406700000911
interpolation results in B (omega),
Figure BDA00033735406700000912
And
Figure BDA00033735406700000913
the specific method comprises the following steps:
at the point of finding the interpolation
Figure BDA00033735406700000914
And
Figure BDA00033735406700000915
b (omega) at other time-frequency coordinates,
Figure BDA00033735406700000916
And
Figure BDA00033735406700000917
and (3) performing approximate calculation by adopting an interpolation technology, and approximately reconstructing results at the original time and frequency coordinates according to an interpolation result, namely:
Figure BDA00033735406700000918
s5, generating space multipoint earthquake motion samples by using a random simulation formula based on a spectrum representation method and utilizing Fast Fourier Transform (FFT), wherein the method specifically comprises the following steps:
(1) based on a spectral representation method, the simulation formula of any random process is as follows:
Figure BDA0003373540670000101
wherein t is a time coordinate; omegalIs the l-th frequency point; Δ ω is the frequency interval; e is an exponential function; i is an imaginary unit; phi is aklIs a set of random phase angles, a deterministic quantity for a particular sample; re represents the real part of the acquired complex number; hjk(ω, t) is the element of the decomposed spectral matrix, calculated by:
Figure BDA0003373540670000102
wherein, thetajk(ω) is HjkThe phase angle of (ω, t), which is equal to the coherence function γ for seismic oscillationsjk(ω) a phase angle;
(2) substituting equations (13) and (15) for equation (14) yields the following simulation equation:
Figure BDA0003373540670000103
and by exchanging the summation order:
Figure BDA0003373540670000104
then n is performed on the formula (17)ΦAnd performing a secondary FFT operation, namely generating seismic motion samples of n points in space.
The method of the invention is further illustrated by the following example of simulating a horizontally distributed non-stationary seismic motion field:
(1) the simulation points are 50 points evenly distributed along the 2000m horizontal straight line direction.
(2) Assuming the self-evolving power spectral density function at each point as:
Figure BDA0003373540670000105
wherein the zero-mean stationary process has the same power spectral density, described as Kanai-Tajimi spectra:
Figure BDA0003373540670000111
wherein S is0＝0.1cm2/s3(ii) a For K1(f) Take f g115/(2 pi) Hz; for K2(f) Take fg2＝5/(2π)Hz；ζg0.25. The general time modulation function is:
Figure BDA0003373540670000112
wherein c is 0.125 and d is 0.25.
Figure BDA0003373540670000113
Wherein, epsilon is 16s and mu is 4 s.
The time invariant coherence function is considered to be:
Figure BDA0003373540670000114
Figure BDA0003373540670000115
θ(f)＝θ0[1+(f/f0)b]-1/2
wherein v isjkRepresents the distance between two points; v is the propagation speed of the wave, and V is 2000 m/s; α is 0.147; β ═ 0.736; theta0＝5210m/s；f01.09 Hz; and b is 2.67.
(3) Determining uniformly distributed frequency interpolation points based on the defined power spectral density function and coherence function of seismic evolution
Figure BDA0003373540670000116
Figure BDA0003373540670000117
Wherein, ω is1And omegauRespectively a first frequency point and a last frequency point, respectively taken as 0.0767rad/s and 50 pi rad/s,
Figure BDA0003373540670000118
the number of the interpolation points of the frequency of the coherent function is 48, and the condition is satisfied
Figure BDA0003373540670000119
N is the total number of frequency discrete points, and N is 2048;
determining time-frequency interpolation points
Figure BDA0003373540670000121
The coordinates of each point are as follows:
Figure BDA0003373540670000122
Figure BDA0003373540670000123
wherein, ω is1Take 0.0767rad/s, omegauTake 50 π rad/s, T0The sample was taken for 40.96s,
Figure BDA0003373540670000124
for the number of frequency interpolation points of the evolution power spectral density function, 64 is taken,
Figure BDA00033735406700001221
and 48, taking the number of the time interpolation points of the evolution power spectral density function.
(4) Calculating frequency interpolation points
Figure BDA0003373540670000125
Delayed coherence matrix for n-point seismic field to be simulated
Figure BDA0003373540670000126
Computing time-frequency interpolation points
Figure BDA0003373540670000127
Self-evolving power spectrum at n earthquake motion simulation points
Figure BDA0003373540670000128
The specific method comprises the following steps:
calculating a delay coherence matrix by substituting the value of the frequency interpolation point into a coherence function:
Figure BDA0003373540670000129
then, substituting the time and frequency values of the time-frequency interpolation points by adopting the same mode to obtain the self-evolution power spectrum of each earthquake motion simulation point
Figure BDA00033735406700001210
(5) To pair
Figure BDA00033735406700001211
Performing Cholesky decomposition to obtain
Figure BDA00033735406700001212
By intrinsic orthogonal decomposition (POD)
Figure BDA00033735406700001213
Decomposed into a small number of principal coordinates
Figure BDA00033735406700001214
And feature vectors
Figure BDA00033735406700001215
The specific method comprises：
To pair
Figure BDA00033735406700001216
Performing Cholesky decomposition to obtain
Figure BDA00033735406700001217
And will be
Figure BDA00033735406700001218
The decomposition is as follows:
Figure BDA00033735406700001219
in the formula, the superscript T represents the transpose of a matrix or vector,
Figure BDA00033735406700001220
is the lower triangular matrix, expressed as:
Figure BDA0003373540670000131
each is decomposed by POD
Figure BDA0003373540670000132
Obtain corresponding main coordinates
Figure BDA0003373540670000133
And feature vectors
Figure BDA0003373540670000134
At all time-frequency interpolation points
Figure BDA0003373540670000135
The values are expressed as a constant matrix, i.e.:
Figure BDA0003373540670000136
each column of the constant matrix is regarded as a column vector
Figure BDA0003373540670000137
A column vector decomposed by the following eigenvectors:
Figure BDA0003373540670000138
finding an optimal set of orthonormal bases
Figure BDA0003373540670000139
The projection of the column vector will be maximized thereon, where ΦqIs the q-th feature vector; lambda [ alpha ]qIs the qth eigenvalue; r is this
Figure BDA00033735406700001310
The correlation matrix of the column vectors can be calculated by:
Figure BDA00033735406700001311
in obtaining
Figure BDA00033735406700001312
After each feature vector, the projection of each column vector, i.e., the principal coordinate, is determined by the following equation:
Figure BDA00033735406700001313
wherein, aqIs the q projection vector;
and performing descending order recombination on the characteristic values, reserving a low-order characteristic value containing more energy, and approximately expressing L as:
Figure BDA0003373540670000141
wherein N isΦTaking N as the number of the low-order items selected to meet the precision requirementΦ4. Whereby each one
Figure BDA0003373540670000142
Using its principal coordinates
Figure BDA0003373540670000143
And feature vectors
Figure BDA0003373540670000144
Expressed as:
Figure BDA0003373540670000145
wherein, the first and the second end of the pipe are connected with each other,
Figure BDA0003373540670000146
and
Figure BDA0003373540670000147
are respectively a vector aqAnd phiqA corresponding discrete function.
(6) By
Figure BDA0003373540670000148
And
Figure BDA0003373540670000149
interpolation results in B (omega),
Figure BDA00033735406700001410
And
Figure BDA00033735406700001411
the specific method comprises the following steps:
at the point of finding the interpolation
Figure BDA00033735406700001412
And
Figure BDA00033735406700001413
b (omega) at other time-frequency coordinates,
Figure BDA00033735406700001414
And
Figure BDA00033735406700001415
and (3) performing approximate calculation by adopting an interpolation technology, and approximately reconstructing results at the original time and frequency coordinates according to an interpolation result, namely:
Figure BDA00033735406700001416
(7) a random simulation formula based on a spectrum representation method is used for generating space multipoint earthquake motion samples by utilizing Fast Fourier Transform (FFT), and the specific method is as follows:
based on a spectral representation method, the simulation formula of any random process is as follows:
Figure BDA00033735406700001417
where N is the number of frequency steps, 2048 is taken, and Δ ω is ωuN is the frequency increment, ωlL Δ ω is a frequency point coordinate, t is a time coordinate; omegalIs the l-th frequency point; e represents an exponential function; i is an imaginary unit; phi is a unit ofklIs a set of random phase angles, at 0,2 pi]Uniformly distributed in the range, and the specific sample is a determined amount; re represents the real part of the acquired complex number; hjk(ω, t) is the element of the decomposed spectral matrix, calculated by:
Figure BDA0003373540670000151
wherein, thetajk(ω) is HjkThe phase angle of (ω, t), which is equal to the coherence function γ for seismic oscillationsjk(ω) a phase angle;
the following simulation formula of the multipoint earthquake motion is obtained from the above steps:
Figure BDA0003373540670000152
and by exchanging the summation order:
Figure BDA0003373540670000153
then proceed the above formula to nNΦAnd a secondary FFT operation, namely generating seismic motion samples.
Through the simulation method, the seismic motion time course of 50 points in a seismic field can be simulated, the autocorrelation function and the cross-correlation function of 1000 samples are calculated, and then the autocorrelation function and the cross-correlation function are compared with a given target value. The verification results of the autocorrelation function and the cross-correlation function of the simulated seismic motion time course are respectively shown in fig. 2, and it can be known that the simulated value is consistent with the target value, so that the rationality of the simulation method can be explained. Further, comparing the simulation efficiency of the method of the present invention with that of the conventional method, the result shows that the method of the present invention only needs 4.2 seconds for simulating one sample, while the conventional method needs 149.9 seconds, which shows that the simulation efficiency is greatly improved by using the method of the present invention.
In the embodiment of the invention, the method is adopted to simulate the complete non-stationary earthquake motion, Cholesky decomposition is only needed to be carried out at a few frequency interpolation points, and the results of the rest frequency points are approximately obtained through interpolation, so that the calculation of the spectral matrix decomposition is greatly simplified. In addition, POD decoupling is only carried out on the power spectrum functions at a few time-frequency interpolation points, and decoupling results of other time-frequency points are obtained approximately through interpolation, so that the calculation amount of time-frequency function decoupling is greatly reduced. Finally, the number of FFT operations is further reduced. Under the condition of keeping enough precision, the efficiency of the whole earthquake dynamic simulation process is greatly improved. The method can provide the input earthquake dynamic load time course for the time course analysis of the earthquake-resistant design of large-span structures, such as a transmission tower line system, an oil pipeline and the like.
The present embodiment is only for explaining the present invention, and it is not limited to the present invention, and those skilled in the art can make modifications of the present embodiment without inventive contribution as needed after reading the present specification, but all of them are protected by patent law within the scope of the claims of the present invention.

Claims (6)

1. A method for simulating a spatial variation non-stationary seismic motion time course is characterized by comprising the following steps: the method specifically comprises the following steps:
s1, determining a frequency interpolation point based on the given seismic motion coherence function gamma (omega) and power spectral density function S (omega, t)
Figure FDA0003373540660000011
And time-frequency interpolation point
Figure FDA0003373540660000012
S2, calculating frequency interpolation point
Figure FDA0003373540660000013
Delayed coherence matrix for n-point seismic field to be simulated
Figure FDA0003373540660000014
Computing time-frequency interpolation points
Figure FDA0003373540660000015
Self-evolving power spectrum at n seismographic simulation points
Figure FDA0003373540660000016
S3, pair
Figure FDA0003373540660000017
Performing Cholesky decomposition to obtain
Figure FDA0003373540660000018
By intrinsic orthogonal decomposition (POD)
Figure FDA0003373540660000019
Decomposed into a small number of principal coordinates
Figure FDA00033735406600000110
And feature vectors
Figure FDA00033735406600000111
S4, is prepared from
Figure FDA00033735406600000112
And
Figure FDA00033735406600000113
interpolation results in B (omega),
Figure FDA00033735406600000114
And
Figure FDA00033735406600000115
and S5, generating the spatial multi-point seismic motion samples by using a random simulation formula based on a spectrum representation method and utilizing Fast Fourier Transform (FFT).
2. The method as claimed in claim 1, wherein the method comprises the following steps: the specific method of step S1 is:
(1) determining frequency interpolation points
Figure FDA00033735406600000116
Selecting uniformly distributed frequency interpolation points:
Figure FDA00033735406600000117
wherein, ω is1And ωuRespectively, the first frequency point and the last frequencyThe point(s) is (are) such that,
Figure FDA00033735406600000118
the number of frequency interpolation points is satisfied
Figure FDA00033735406600000119
N is the total number of frequency discrete points;
(2) determining time-frequency interpolation points
Figure FDA00033735406600000120
Selecting evenly distributed time-frequency interpolation points
Figure FDA00033735406600000121
The determination method is as follows:
Figure FDA0003373540660000021
Figure FDA0003373540660000022
in the formula, T0Is the total duration of the power spectrum;
Figure FDA0003373540660000023
and
Figure FDA0003373540660000024
the number of time-frequency interpolation points along the frequency direction and the time direction respectively; and is
Figure FDA0003373540660000025
M is the total number of time discrete points.
3. The method as claimed in claim 1, wherein the method comprises the following steps: the specific method of step S2 is:
calculating a delay coherence matrix by substituting the value of the frequency interpolation point into a coherence function:
Figure FDA0003373540660000026
then, substituting the time and frequency values of the time-frequency interpolation points by adopting the same mode to obtain the self-evolution power spectrum of each earthquake motion simulation point
Figure FDA0003373540660000027
4. The method as claimed in claim 1, wherein the method comprises the following steps: the specific method of step S3 is:
(1) to pair
Figure FDA0003373540660000028
Performing Cholesky decomposition to obtain
Figure FDA0003373540660000029
And will be
Figure FDA00033735406600000210
The decomposition is as follows:
Figure FDA00033735406600000211
in the formula, the superscript T represents the transpose of a matrix or vector,
Figure FDA00033735406600000212
is the lower triangular matrix, expressed as:
Figure FDA00033735406600000213
(2) each is decomposed by POD
Figure FDA00033735406600000214
Obtain corresponding main coordinates
Figure FDA00033735406600000215
And feature vectors
Figure FDA00033735406600000216
(3) At all time-frequency interpolation points
Figure FDA0003373540660000031
The values are expressed as a constant matrix, i.e.:
Figure FDA0003373540660000032
each column of the constant matrix is regarded as a column vector
Figure FDA0003373540660000033
A column vector decomposed by the following eigenvectors:
Figure FDA0003373540660000034
finding an optimal set of orthonormal bases
Figure FDA0003373540660000035
The projection of the column vector will be maximized thereon, where ΦqIs the q-th feature vector; lambda [ alpha ]qIs the qth eigenvalue; r is this
Figure FDA0003373540660000036
The correlation matrix of the column vectors can be calculated by:
Figure FDA0003373540660000037
in obtaining
Figure FDA0003373540660000038
After each feature vector, the projection of each column vector, i.e., the principal coordinate, is determined by the following equation:
Figure FDA0003373540660000039
wherein, aqIs the q projection vector;
(4) and performing descending order recombination on the characteristic values, reserving a low-order characteristic value containing more energy, and approximately expressing L as:
Figure FDA00033735406600000310
wherein N isΦThe number of low-order terms selected to meet the accuracy requirement; based on the above formula, each
Figure FDA0003373540660000041
By its principal coordinate
Figure FDA0003373540660000042
And feature vectors
Figure FDA0003373540660000043
Expressed as:
Figure FDA0003373540660000044
wherein, the first and the second end of the pipe are connected with each other,
Figure FDA0003373540660000045
and
Figure FDA0003373540660000046
are respectively a vector aqAnd phiqA corresponding discrete function.
5. The method as claimed in claim 1, wherein the method for simulating the time-course of the spatially variant non-stationary seismic activity comprises: the specific method of step S4 is:
at the point of finding the interpolation
Figure FDA0003373540660000047
And
Figure FDA0003373540660000048
then B (omega) at other time-frequency coordinates,
Figure FDA0003373540660000049
And
Figure FDA00033735406600000410
and (3) performing approximate calculation by adopting an interpolation technology, and approximately reconstructing results at the original time and frequency coordinates according to an interpolation result, namely:
Figure FDA00033735406600000411
6. the method as claimed in claim 1, wherein the method for simulating the time-course of the spatially variant non-stationary seismic activity comprises: the specific method of step S5 is:
(1) based on a spectral representation method, the simulation formula of any random process is as follows:
Figure FDA00033735406600000412
wherein t is a time coordinate; omegalIs the l-th frequency point; Δ ω is the frequency interval; e is an exponential function; i is an imaginary unit; phi is aklIs a set of random phase angles, a deterministic quantity for a particular sample; re represents the real part of the acquired complex number; hjk(ω, t) is the element of the decomposed spectral matrix, calculated by:
Figure FDA00033735406600000413
wherein, thetajk(ω) is HjkThe phase angle of (ω, t), which is equal to the coherence function γ for seismic oscillationsjk(ω) a phase angle;
(2) substituting equations (13) and (15) for equation (14) yields the following simulation equation:
Figure FDA0003373540660000051
and by exchanging the summation order:
Figure FDA0003373540660000052
then n is performed on the formula (17)ΦAnd performing a secondary FFT operation, namely generating seismic motion samples of n points in space.
CN202111410406.3A 2021-11-25 2021-11-25 Method for simulating space variation non-stationary earthquake motion time course Active CN114139363B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111410406.3A CN114139363B (en) 2021-11-25 2021-11-25 Method for simulating space variation non-stationary earthquake motion time course

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111410406.3A CN114139363B (en) 2021-11-25 2021-11-25 Method for simulating space variation non-stationary earthquake motion time course

Publications (2)

Publication Number Publication Date
CN114139363A CN114139363A (en) 2022-03-04
CN114139363B true CN114139363B (en) 2022-07-01

Family

ID=80391607

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111410406.3A Active CN114139363B (en) 2021-11-25 2021-11-25 Method for simulating space variation non-stationary earthquake motion time course

Country Status (1)

Country Link
CN (1) CN114139363B (en)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115526125B (en) * 2022-09-22 2023-09-19 四川农业大学 Numerical cutoff-based efficient simulation method for random wind speed field

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106682277A (en) * 2016-12-06 2017-05-17 西南交通大学 Fast simulation method for non-stationary random process
CN107480325A (en) * 2017-07-03 2017-12-15 河海大学 The non-stationary non-gaussian earthquake motion time history analogy method of spatial variability
CN107657127A (en) * 2017-10-11 2018-02-02 河海大学 One-dimensional multi-variate random process efficient analogy method based on energy interpolations such as frequency domains
CN110059286A (en) * 2019-03-07 2019-07-26 重庆大学 A kind of structure non stationary response efficient analysis method based on FFT

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106682277A (en) * 2016-12-06 2017-05-17 西南交通大学 Fast simulation method for non-stationary random process
CN107480325A (en) * 2017-07-03 2017-12-15 河海大学 The non-stationary non-gaussian earthquake motion time history analogy method of spatial variability
CN107657127A (en) * 2017-10-11 2018-02-02 河海大学 One-dimensional multi-variate random process efficient analogy method based on energy interpolations such as frequency domains
CN110059286A (en) * 2019-03-07 2019-07-26 重庆大学 A kind of structure non stationary response efficient analysis method based on FFT

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
Efficient Nonstationary Stochastic Response Analysis for Linear and Nonlinear Structures by FFT;Ning Zhao 等;《Journal of Engineering》;20190220;全文 *
Efficient simulation of fully non-stationary random wind field;Tianyou Tao 等;《Mechanical Systems and Signal Processing》;20200911;全文 *
地震动随机场的POD降维表达;刘章军等;《中国科学:技术科学》;20181114(第05期);全文 *

Also Published As

Publication number Publication date
CN114139363A (en) 2022-03-04

Similar Documents

Publication Publication Date Title
Egholm et al. Modeling the flow of glaciers in steep terrains: The integrated second‐order shallow ice approximation (iSOSIA)
Masri et al. A nonparametric identification technique for nonlinear dynamic problems
Fukumori A partitioned Kalman filter and smoother
US8334803B1 (en) Method for simulating noisy radar target echoes
CN110083920B (en) Analysis method for random response of non-proportional damping structure under earthquake action
Zhang et al. Assimilation of ice motion observations and comparisons with submarine ice thickness data
CN113806845B (en) POD interpolation-based non-stationary wind field efficient simulation method
CN114139363B (en) Method for simulating space variation non-stationary earthquake motion time course
Li et al. An efficient parallel Krylov-Schur method for eigen-analysis of large-scale power systems
Ohtani et al. Fast computation of quasi-dynamic earthquake cycle simulation with hierarchical matrices
Berman et al. Improvement of analytical dynamic models using modal test data
Troullinos et al. Coherency and model reduction: State space point of view
CN110059286A (en) A kind of structure non stationary response efficient analysis method based on FFT
Özgökmen et al. Assimilation of drifter observations in primitive equation models of midlatitude ocean circulation
Wang et al. Response surface method using grey relational analysis for decision making in weapon system selection
Juhlin et al. Optimal sensor placement for localizing structured signal sources
Wyatt Development and assessment of a nonlinear wave prediction methodology for surface vessels
Gershgorin et al. A nonlinear test model for filtering slow-fast systems
Cerv et al. Wave motion in a thick cylindrical rod undergoing longitudinal impact
Chen et al. Data-driven arbitrary polynomial chaos expansion on uncertainty quantification for real-time hybrid simulation under stochastic ground motions
Auton et al. Investigation of procedures for automatic resonance extraction from noisy transient electromagnetics data
Cooper et al. Non-intrusive polynomial chaos for efficient uncertainty analysis in parametric roll simulations
Gebbie Subduction in an eddy-resolving state estimate of the northeast Atlantic Ocean
CN115526125B (en) Numerical cutoff-based efficient simulation method for random wind speed field
Tuan et al. A fuzzy finite element algorithm based on response surface method for free vibration analysis of structure

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant