Rock physics

#Rock physics| 来源: 网络整理| 查看: 265

Rock physics aims to characterize rock properties based on the behavior of seismic waves propagating through them. This requires consideration of how the composition of a rock dictates its stress-strain relationship and thus seismic response. The effect of pore fluids is of particular interest due to its applicability to the hydrocarbon industry. In a standard seismic interpretation workflow rock physics is used to relate impedance and elastic parameters derived from seismic data to specific rock properties. This constrains what seismic data is physically capable of resolving and the non-uniqueness associated with a specific interpretation.

Contents 1 Mechanics 1.1 Hooke's Law 1.2 The Wave Equation 2 Elastic Moduli 2.1 Young's Modulus 2.2 Shear Modulus 2.3 Bulk Modulus 2.4 Poisson's Ratio 2.5 Moduli Relations 3 Effective Moduli 3.1 Voigt Bounds 3.2 Reuss Bounds 3.3 Voigt-Reuss-Hill average 3.4 Hashin-Shtrikman Bounds 3.5 Modified Hashin-Shtrikman Bounds 3.6 Kuster and Toksöz formulation 3.7 Self-consistent approximations (SCA) 3.8 Differential Effective Medium Model (DEM) 4 Fluid Effects 4.1 Gassmann's Equation 4.2 Biot's Equation 5 Cross Plots 6 Anisotropy 7 Important papers 8 References 9 See also Mechanics Hooke's Law

A material is linearly elastic if its stress-strain relationship is governed by Hooke's Law,

σ i j = c i j k l ε k l {\displaystyle \sigma _{ij}=c_{ijkl}\varepsilon _{kl}}

where σ i j {\displaystyle \sigma _{ij}} is the second order stress tensor, ε k l {\displaystyle \varepsilon _{kl}} is the second order strain tensor, and c i j k l {\displaystyle c_{ijkl}} is the fourth order elastic tensor, defining a material's resistance to strain due to a stress in each direction.

The fourth order elasticity or stiffness tensor requires 34 terms to specify the stress strain relationship for an arbitrary homogeneous material. Conservation of angular momentum requires the stress tensor to be symmetric, as any asymmetries would imply rotations, which are not described by stress. Stress describes deformation only. This also requires the strain tensor to be symmetric, implying that the elasticity tensor is symmetric, reducing the number of required components to 36. Energy conservation requires thermal energy and strain energy to balance, further reducing the number of required components to 21.

If one makes the simplifying assumption that the material is isotropic, that is, the strain experienced by a material is independent of the direction in which the stress is applied, the stiffness tensor requires only two parameters, Lamé's constants, λ and μ. Hooke's Law can then be expressed as

σ i j = λ δ i j ε k k + 2 μ ε i j {\displaystyle \sigma _{ij}=\lambda \delta _{ij}\varepsilon _{kk}+2\mu \varepsilon _{ij}} The Wave Equation

Equations describing the waves resulting from deformation of elastic materials can be derived by applying Hooke's Law to Cauchy's Equation of motion, given by:

ρ ∂ 2 u i ∂ t 2 = ∂ σ j i ∂ x j {\displaystyle \rho {\frac {\partial ^{2}u_{i}}{\partial t^{2}}}={\frac {\partial \sigma _{ji}}{\partial x_{j}}}}

where ρ is the density of the material and u_{i} is the displacement of the material, related to the infinitesimal strain tensor by

ε i j = 1 2 ( ∂ u i ∂ x j + ∂ u j ∂ x i ) {\displaystyle \varepsilon _{ij}={\tfrac {1}{2}}({\frac {\partial u_{i}}{\partial x_{j}}}+{\frac {\partial u_{j}}{\partial x_{i}}})}

Hooke's Law is used to describe the stress tensor in Cauchy's Equation along with the above definition of the infinitesimal strain tensor to give an equation for material displacement as a function of time and space. Applying the divergence to this equation describes the rate at which volumetric perturbations propagate through the material giving the P-wave equation. For isotropic materials, this simplifies to

∂ 2 Θ ∂ t 2 = λ + 2 μ ρ ▽ 2 Θ {\displaystyle {\frac {\partial ^{2}\Theta }{\partial t^{2}}}={\frac {\lambda +2\mu }{\rho }}\bigtriangledown ^{2}\Theta }

V p = λ + 2 μ ρ {\displaystyle V_{p}={\sqrt {\frac {\lambda +2\mu }{\rho }}}}

where Θ is the dilatation, Θ = ▽ ∙ u → {\displaystyle \Theta =\bigtriangledown \bullet {\vec {u}}} . Similarly, applying the curl operator describes the rate at which shape perturbations propagate through a material giving the S-wave equation. For isotropic materials this gives

∂ 2 ω ∂ t 2 = μ ρ ▽ 2 ω {\displaystyle {\frac {\partial ^{2}\omega }{\partial t^{2}}}={\frac {\mu }{\rho }}\bigtriangledown ^{2}\omega } V s = μ ρ {\displaystyle V_{s}={\sqrt {\frac {\mu }{\rho }}}}

where ω is the rotation vector, ω = ▽ × u → {\displaystyle \omega =\bigtriangledown \times {\vec {u}}} .

It is important to note that these velocities assume a homogeneous, isotropic, perfectly elastic material. All energy associated with deformation is conserved. No energy is lost to friction etc., so velocities do not depend on wave frequency. These assumptions may not always be accurate, in practice.

Elastic Moduli

The following definitions assume a linearly elastic, homogenous isotropic material.

Young's Modulus

Young's modulus, E (Pa), describes a material's ability to stretch or compress by quantifying how resistant a material is to extensional deformation.

σ i i = E ε i i {\displaystyle \sigma _{ii}=E\varepsilon _{ii}}

In the hydrocarbon industry, Young's modulus is used to design hydraulic fracturing mechanisms. Rigid rocks with a high value of Young's modulus will result in narrower fractures than rocks with a lower rigidity.

Shear Modulus

The shear modulus, µ (Pa) or G describes the rigidity of a material, its resistance to shear strain/angular distortions.

σ i j = 2 μ ε i j {\displaystyle \sigma _{ij}=2\mu \varepsilon _{ij}} for i≠j

Note that fluids have a shear modulus of zero as they cannot withstand any shear stress.

Bulk Modulus

The bulk modulus, K (Pa) of a material describes its ability to withstand volumetric stress,

σ i i = 3 K Θ {\displaystyle \sigma _{ii}=3K\Theta }

The inverse of the bulk modulus is defined as the compressibility, β

K = 1 β {\displaystyle K={\frac {1}{\beta }}}

The bulk modulus is related to the Lamé parameters by

K = λ + 2 3 μ {\displaystyle K=\lambda +{\tfrac {2}{3}}\mu }

When compared to the isotropic form of Hooke's Law, this expresses that volumetric deformation depends on normal and shear stiffnesses. The bulk modulus quantifies its total volumetric stiffness. The Lamé parameter λ can therefore be thought of as a material's resistance to exclusively normal stresses, whereas µ quantifies the material's resistance to shape deformations, that can result from both normal and non-normal stress regimes.

Poisson's Ratio

Poisson's ratio, ν defines the transverse strain that results from longitudinal strain. That is, when a material is stretched in one direction, Poisson's ratio defines how much it will contract in the other directions. For an isotropic material deformed in the x direction,

ν = − d ε y d ε x = − d ε z d ε x {\displaystyle \nu =-{\frac {d\varepsilon _{y}}{d\varepsilon _{x}}}=-{\frac {d\varepsilon _{z}}{d\varepsilon _{x}}}}

The value of Poisson's Ratio varies from -1 to 0.5. A material is perfectly incompressible if its Poisson's ratio is 0.5, such as rubber. Materials that do not deform transversally due to longitudinal strain have a Poisson ratio of 0, such as cork. Materials with negative Poisson's ratios expand in transverse directions when stretched longitudinally due to their hinge-like molecular structures and are known as auxetic materials. A material is described as an ideal Poisson solid when λ=µ giving ν=0.25. This is roughly true for many rocks.

Poisson's ratio can also be expressed in terms of the materials P and S wave velocities,

ν = ( V p V s ) 2 − 2 2 ( V p V s ) 2 − 2 {\displaystyle \nu ={\frac {({\tfrac {V_{p}}{V_{s}}})^{2}-2}{2({\tfrac {V_{p}}{V_{s}}})^{2}-2}}}

Because µ=0 for fluids, Vs=0, making Poisson's ratio very sensitive to changes in fluid content.

Moduli Relations Conversion formulas Homogeneous isotropic linear elastic materials have their elastic properties uniquely determined by any two moduli among these; thus, given any two, any other of the elastic moduli can be calculated according to these formulas. K = {\displaystyle K=\,} E = {\displaystyle E=\,} λ = {\displaystyle \lambda =\,} μ = {\displaystyle \mu =\,} ν = {\displaystyle \nu =\,} M = {\displaystyle M=\,} Notes ( K , E ) {\displaystyle (K,\,E)} K {\displaystyle K} E {\displaystyle E} 3 K ( 3 K − E ) 9 K − E {\displaystyle {\tfrac {3K(3K-E)}{9K-E}}} 3 K E 9 K − E {\displaystyle {\tfrac {3KE}{9K-E}}} 3 K − E 6 K {\displaystyle {\tfrac {3K-E}{6K}}} 3 K ( 3 K + E ) 9 K − E {\displaystyle {\tfrac {3K(3K+E)}{9K-E}}} ( K , λ ) {\displaystyle (K,\,\lambda )} K {\displaystyle K} 9 K ( K − λ ) 3 K − λ {\displaystyle {\tfrac {9K(K-\lambda )}{3K-\lambda }}} λ {\displaystyle \lambda } 3 ( K − λ ) 2 {\displaystyle {\tfrac {3(K-\lambda )}{2}}} λ 3 K − λ {\displaystyle {\tfrac {\lambda }{3K-\lambda }}} 3 K − 2 λ {\displaystyle 3K-2\lambda \,} ( K , μ ) {\displaystyle (K,\,\mu )} K {\displaystyle K} 9 K μ 3 K + μ {\displaystyle {\tfrac {9K\mu }{3K+\mu }}} K − 2 μ 3 {\displaystyle K-{\tfrac {2\mu }{3}}} μ {\displaystyle \mu } 3 K − 2 μ 2 ( 3 K + μ ) {\displaystyle {\tfrac {3K-2\mu }{2(3K+\mu )}}} K + 4 μ 3 {\displaystyle K+{\tfrac {4\mu }{3}}} ( K , ν ) {\displaystyle (K,\,\nu )} K {\displaystyle K} 3 K ( 1 − 2 ν ) {\displaystyle 3K(1-2\nu )\,} 3 K ν 1 + ν {\displaystyle {\tfrac {3K\nu }{1+\nu }}} 3 K ( 1 − 2 ν ) 2 ( 1 + ν ) {\displaystyle {\tfrac {3K(1-2\nu )}{2(1+\nu )}}} ν {\displaystyle \nu } 3 K ( 1 − ν ) 1 + ν {\displaystyle {\tfrac {3K(1-\nu )}{1+\nu }}} ( K , M ) {\displaystyle (K,\,M)} K {\displaystyle K} 9 K ( M − K ) 3 K + M {\displaystyle {\tfrac {9K(M-K)}{3K+M}}} 3 K − M 2 {\displaystyle {\tfrac {3K-M}{2}}} 3 ( M − K ) 4 {\displaystyle {\tfrac {3(M-K)}{4}}} 3 K − M 3 K + M {\displaystyle {\tfrac {3K-M}{3K+M}}} M {\displaystyle M} ( E , λ ) {\displaystyle (E,\,\lambda )} E + 3 λ + R 6 {\displaystyle {\tfrac {E+3\lambda +R}{6}}} E {\displaystyle E} λ {\displaystyle \lambda } E − 3 λ + R 4 {\displaystyle {\tfrac {E-3\lambda +R}{4}}} 2 λ E + λ + R {\displaystyle {\tfrac {2\lambda }{E+\lambda +R}}} E − λ + R 2 {\displaystyle {\tfrac {E-\lambda +R}{2}}} R = E 2 + 9 λ 2 + 2 E λ {\displaystyle R={\sqrt {E^{2}+9\lambda ^{2}+2E\lambda }}} ( E , μ ) {\displaystyle (E,\,\mu )} E μ 3 ( 3 μ − E ) {\displaystyle {\tfrac {E\mu }{3(3\mu -E)}}} E {\displaystyle E} μ ( E − 2 μ ) 3 μ − E {\displaystyle {\tfrac {\mu (E-2\mu )}{3\mu -E}}} μ {\displaystyle \mu } E 2 μ − 1 {\displaystyle {\tfrac {E}{2\mu }}-1} μ ( 4 μ − E ) 3 μ − E {\displaystyle {\tfrac {\mu (4\mu -E)}{3\mu -E}}} ( E , ν ) {\displaystyle (E,\,\nu )} E 3 ( 1 − 2 ν ) {\displaystyle {\tfrac {E}{3(1-2\nu )}}} E {\displaystyle E} E ν ( 1 + ν ) ( 1 − 2 ν ) {\displaystyle {\tfrac {E\nu }{(1+\nu )(1-2\nu )}}} E 2 ( 1 + ν ) {\displaystyle {\tfrac {E}{2(1+\nu )}}} ν {\displaystyle \nu } E ( 1 − ν ) ( 1 + ν ) ( 1 − 2 ν ) {\displaystyle {\tfrac {E(1-\nu )}{(1+\nu )(1-2\nu )}}} ( E , M ) {\displaystyle (E,\,M)} 3 M − E + S 6 {\displaystyle {\tfrac {3M-E+S}{6}}} E {\displaystyle E} M − E + S 4 {\displaystyle {\tfrac {M-E+S}{4}}} 3 M + E − S 8 {\displaystyle {\tfrac {3M+E-S}{8}}} E − M + S 4 M {\displaystyle {\tfrac {E-M+S}{4M}}} M {\displaystyle M}

S = ± E 2 + 9 M 2 − 10 E M {\displaystyle S=\pm {\sqrt {E^{2}+9M^{2}-10EM}}}

There are two valid solutions. The plus sign leads to ν ≥ 0 {\displaystyle \nu \geq 0} . The minus sign leads to ν ≤ 0 {\displaystyle \nu \leq 0} .

( λ , μ ) {\displaystyle (\lambda ,\,\mu )} λ + 2 μ 3 {\displaystyle \lambda +{\tfrac {2\mu }{3}}} μ ( 3 λ + 2 μ ) λ + μ {\displaystyle {\tfrac {\mu (3\lambda +2\mu )}{\lambda +\mu }}} λ {\displaystyle \lambda } μ {\displaystyle \mu } λ 2 ( λ + μ ) {\displaystyle {\tfrac {\lambda }{2(\lambda +\mu )}}} λ + 2 μ {\displaystyle \lambda +2\mu \,} ( λ , ν ) {\displaystyle (\lambda ,\,\nu )} λ ( 1 + ν ) 3 ν {\displaystyle {\tfrac {\lambda (1+\nu )}{3\nu }}} λ ( 1 + ν ) ( 1 − 2 ν ) ν {\displaystyle {\tfrac {\lambda (1+\nu )(1-2\nu )}{\nu }}} λ {\displaystyle \lambda } λ ( 1 − 2 ν ) 2 ν {\displaystyle {\tfrac {\lambda (1-2\nu )}{2\nu }}} ν {\displaystyle \nu } λ ( 1 − ν ) ν {\displaystyle {\tfrac {\lambda (1-\nu )}{\nu }}} Cannot be used when ν = 0 ⇔ λ = 0 {\displaystyle \nu =0\Leftrightarrow \lambda =0} ( λ , M ) {\displaystyle (\lambda ,\,M)} M + 2 λ 3 {\displaystyle {\tfrac {M+2\lambda }{3}}} ( M − λ ) ( M + 2 λ ) M + λ {\displaystyle {\tfrac {(M-\lambda )(M+2\lambda )}{M+\lambda }}} λ {\displaystyle \lambda } M − λ 2 {\displaystyle {\tfrac {M-\lambda }{2}}} λ M + λ {\displaystyle {\tfrac {\lambda }{M+\lambda }}} M {\displaystyle M} ( μ , ν ) {\displaystyle (\mu ,\,\nu )} 2 μ ( 1 + ν ) 3 ( 1 − 2 ν ) {\displaystyle {\tfrac {2\mu (1+\nu )}{3(1-2\nu )}}} 2 μ ( 1 + ν ) {\displaystyle 2\mu (1+\nu )\,} 2 μ ν 1 − 2 ν {\displaystyle {\tfrac {2\mu \nu }{1-2\nu }}} μ {\displaystyle \mu } ν {\displaystyle \nu } 2 μ ( 1 − ν ) 1 − 2 ν {\displaystyle {\tfrac {2\mu (1-\nu )}{1-2\nu }}} ( μ , M ) {\displaystyle (\mu ,\,M)} M − 4 μ 3 {\displaystyle M-{\tfrac {4\mu }{3}}} μ ( 3 M − 4 μ ) M − μ {\displaystyle {\tfrac {\mu (3M-4\mu )}{M-\mu }}} M − 2 μ {\displaystyle M-2\mu \,} μ {\displaystyle \mu } M − 2 μ 2 M − 2 μ {\displaystyle {\tfrac {M-2\mu }{2M-2\mu }}} M {\displaystyle M} ( ν , M ) {\displaystyle (\nu ,\,M)} M ( 1 + ν ) 3 ( 1 − ν ) {\displaystyle {\tfrac {M(1+\nu )}{3(1-\nu )}}} M ( 1 + ν ) ( 1 − 2 ν ) 1 − ν {\displaystyle {\tfrac {M(1+\nu )(1-2\nu )}{1-\nu }}} M ν 1 − ν {\displaystyle {\tfrac {M\nu }{1-\nu }}} M ( 1 − 2 ν ) 2 ( 1 − ν ) {\displaystyle {\tfrac {M(1-2\nu )}{2(1-\nu )}}} ν {\displaystyle \nu } M {\displaystyle M}

The stiffness matrix (6 by 6 in Voigt notation) in Hooke's law can be parametrized by only two components for homogeneous and isotropic materials. One may choose whichever pair one prefers among the elastic moduli given above. Some of the possible conversions are listed in the table.

Effective Moduli

Most geophysical theories make the simplifying assumption that inherently heterogenous Earth materials obtain some effective, averaged homogenous properties over the dimension and scale of measurement. It is therefore essential to have an understanding of how the components of a material are arranged to define its overall physical properties. The properties of the whole can be significantly different than the properties of its parts.

The process of defining the effective properties of a heterogenous material is referred to as homogenization. Given a material with small scale, rapidly varying properties, homogenization attempts to redefine the material over a larger volume such that it becomes homogenous in a statistical sense at a particular scale. These homogenized volumes are referred to as representative volume elements (RVE’s). The parameters that describe the behavior of these RVE’s are defined as its effective parameters. A material can be considered homogenous if its constituents are smaller than the wavelength of measurement.

Seismic studies therefore require effective elastic moduli to be defined. This depends on the elastic properties of the constituents and their geometric relationships. Because it is essentially impossible to obtain sufficient information to uniquely define the effective moduli of a natural substance, the theoretical maximum and minimum moduli, known as the Voigt and Reuss bounds are often used instead.[1].

More refined effective moduli can be defined by accounting for the arrangement of the material's components. In the case of a rock, the mineral matrix components are similar enough that a simple Voigt-Reuss-Hill average is sufficiently accurate. Pores, pore filling fluids, clays and cements must be considered separately as they can significantly change the effective moduli. There are two general methods of doing so, depending on the relative volume fraction of the component being considered. For large relative volumes a granular approach is usually taken, where the effective moduli is determined by considering how the contact forces on grains change the overall moduli. Such methods are used for high porosity, more unconsolidated to conventional sands. For low relative volumes, the components are considered as isolated deviations from a background medium, suitable for low porosity, unconventional isotropic rocks.

Voigt Bounds

The Voigt Bound is obtained by assuming the strain field remains constant throughout the material when subjected to an arbitrary average stress field. For a composite composed of N components, the effective Voigt moduli, Mv, are given by

M v = ∑ i N f i M i {\displaystyle M_{v}=\sum _{i}^{N}f_{i}M_{i}}

Where fi is the volume fraction of the ith component with a modulus Mi. It was proven by Hill [2] using energy conservation principles that this gives the maximum possible moduli of the composite.

Reuss Bounds

The Reuss Bound is obtained by assuming the stress field remains constant throughout the material in an arbitrary average strain field. Energy conservation principles imply that gives the minimum possible moduli for the material. For a composite composed of N components, the effective Reuss moduli, MR, are given by

1 M R = ∑ i N f i 1 M i {\displaystyle {\frac {1}{M_{R}}}=\sum _{i}^{N}f_{i}{\frac {1}{M_{i}}}}

Where fi is the volume fraction of the ith component with a modulus Mi.

Voigt-Reuss-Hill average

A common simplistic approach to estimating the effective moduli is to calculate the average of the Voigt and Reuss bounds, referred to as the Voigt-Reuss-Hill (VRH) average.

M V R H = M V + M R 2 {\displaystyle M_{VRH}={\frac {M_{V}+M_{R}}{2}}}

This is relatively accurate when the moduli of the components have similar values and is commonly used to estimate the matrix moduli of multi-mineralic rocks.

Hashin-Shtrikman Bounds

Hashin and Shtrikman[3] use a more robust approach using the calculus of variations. They derive the effective elastic moduli that minimize the potential elastic energy of the composite. The strain energy of the homogenous, isotropic reference material, U0 is defined then altered by introducing an inclusion with different elastic properties, and the resultant change in the elastic strain energy, δU=U-U0 is considered. An arbitrary number of inclusions can be added as long as the assumptions required for homogenization remain valid. The change in elastic strain energy due to the inclusion give maximum bounds when the reference material is the phase with the highest moduli and minimum bounds when the reference material is the phase with the lowest moduli. This is analogous to the method used to derive Voigt and Reuss bounds, except that instead of averaging stress and strain fields, Hashin minimizes the strain energy. Obtaining expressions for the effective moduli requires and assumption of spherical symmetry. Thus Hashin-Shtrikman bounds assume spherical inclusions distributed symmetrically about the matrix. For a two phase medium the resultant effective moduli are given by

K H S = K 1 + f 2 ( K 2 − K 1 ) − 1 + f 1 ( K 1 + 4 3 μ 1 ) − 1 {\displaystyle K_{HS}=K_{1}+{\frac {f_{2}}{(K_{2}-K_{1})^{-1}+f_{1}(K_{1}+{\tfrac {4}{3}}\mu _{1})^{-1}}}} μ H S = μ 1 + f 2 ( μ 2 − μ 1 ) − 1 + 2 f 1 ( K 1 + 2 μ 1 ) 5 μ 1 ( K 1 + 4 3 μ 1 ) {\displaystyle \mu _{HS}=\mu _{1}+{\frac {f_{2}}{(\mu _{2}-\mu _{1})^{-1}+{\tfrac {2f_{1}(K_{1}+2\mu _{1})}{5\mu _{1}(K_{1}+{\tfrac {4}{3}}\mu _{1})}}}}}

Where the subscript 1 refers to the matrix material, and the subscript 2 refers to the inclusion material. The upper bound is obtained when the matrix has the maximum moduli, and the lower bound is obtained when the matrix has the minimum moduli. Note that if one of the components is a fluid, the lower bound for the shear modulus is zero, as the lower bound represents a suspension, which is effectively a fluid with a shear modulus of zero.

Modified Hashin-Shtrikman Bounds

The Modified Hashin-Shtrikman Bounds separates consolidated and unconsolidated elastic moduli trends according to the rock's critical porosity. Consider the rock's diagenetic history beginning with an unconsolidated, highly porous sediment that gradually gets compacted until a critical porosity, ϕc (0.36-0.4 for sandstones) is reached. After this critical porosity is reached, cementation takes over as the primary diagenetic process that decreases porosity causing a faster increase in elastic moduli.

Empirically, it has been found that the cementation trend can be modelled using the Upper Hashin-Shtrikman bound by replacing the porosity, or fractional fluid content, ϕ with ϕ ϕ c {\displaystyle {\tfrac {\phi }{\phi _{c}}}} , this is known as the Modified Upper Hashin-Shtrikman bound.

For porosities greater than the critical porosity, the compaction trend can be modelled using the Reuss bound, or the Modified Lower Hashin-Shtriman bound obtained by replacing the porosity ϕ with ϕ ϕ c {\displaystyle {\tfrac {\phi }{\phi _{c}}}} in the Lower Hashin-Shtrikman bound.

Kuster and Toksöz formulation

Kuster and Toksöz [4]derive expressions for effective moduli of a two phase isotropic material using first order scattering theory. This formulation emphasizes the effect of inclusion shape by the incorporation of an aspect ratio term,α.

It is assumed that the concentration, f, of inclusions is small, so that interactions between inclusions can be neglected. The concentration of inclusions can be taken to be small if f/α

【本文地址】

公司简介

联系我们