CN106597546A - Method for detecting river bottom surface sediment physical properties - Google Patents

Abstract

The present invention provides a method for detecting river bottom surface sediment physical properties and belongs to the acoustics and geophysics field. The method includes the following steps that: sub-bottom profile acoustic data are extracted; sub-bottom profile data are analyzed; an attenuation coefficient alpha is calculated; a porosity is searched according to a porosity-reflection coefficient curve; a permeability kappa is calculated; a sediment velocity c is calculated; a permeability-attenuation coefficient curve is calculated, a permeability is searched according to the attenuation coefficient alpha; after a plurality of times of iterative computation, when the difference of a permeability beta<n> and a current iteration value beta<n-1> is small, namely, (|beta<n-1>-beta<n>|/beta<n> is smaller than 3%, iteration is terminated; a sediment density rho is calculated, and an attenuation coefficient R(f) is calculated; and the permeability kappa, the sediment velocity c, the attenuation coefficient R(f), the attenuation coefficient alpha and the sediment density rho are outputted. With the method adopted, a variety of sediment physical property parameters can be obtained; deposited sediment information is enriched; the river bottom surface sediment physical property parameters can be obtained just through using a sub-bottom profiler, and therefore, manpower and material resources can be saved, and work efficiency can be improved; more accurate sediment parameter information can be obtained; and data support can be provided for research such as river dynamics research.

Description

A kind of method of detection river bed top layer silt physical property

Technical field

The invention belongs to acoustics and geophysics field, and in particular to a kind of method of detection river bed top layer silt physical property.

Background technology

At present, the acquisition of river bed alluvial top layer silt physical parameter (density, permeability, porosity etc.) is mainly taken by machinery Sample is obtained, and the perturbation of mechanical sampling causes sample tests to produce deviation.Additionally, mechanical sampling spends a large amount of manpower things Power, it is less efficient.Acoustic method is also commonly used to detect river sediment characteristic, presently mainly judges that silt becomes silted up according to sonar image Product structure, the extraction to silt physical parameter still lacks effective ways.

The present invention is the new method that a kind of utilization sub-bottom profiler detects river bed top layer silt parameter.Based on pore media Acoustic theory, by analyzing sub-bottom profiler signal is received, and acoustic reflection coefficient and the decay for extracting sub-bottom profiler is special Levy, obtain the parameters such as river sediment density, velocity of wave, porosity and permeability.The present invention does not need mechanical sampling test result to make It is foundation, receiving signal merely with sub-bottom profiler can obtain related silt parameter.Compared with mechanical sampling, with quick Efficiently the characteristics of.

The content of the invention

In sum, in order to overcome the deficiencies in the prior art, the invention provides a kind of detection river bed top layer silt physical property Method, involved formula and parameter be as follows：

K is wave number,For angular frequency, f is frequency, and c is pore media velocity of wave, ρ_fAnd ρ_sRespectively pore-fluid density and Sand grain density, and for known constant.ρ is silt body density.η represents pore water viscosity, and κ is permeability, K_γIt is silt The bulk moduluses of grain, μ is the modulus of shearing of silt framework, K_bAnd K_fThe volume of silt skeleton modulus and pore-fluid is represented respectively Modulus, the deviation that F is produced for pore-fluid increases with frequency.The velocity of wave c of water_ω=1500kg/m³, K_r=3.2 × 10¹⁰Pa, K_f =2.395 × 10⁹Pa, η=0.001Kg/m*s.

ρ=β ρ_f+(1-β)ρ_s (2)

τ_x=K₀τ_z (13)

τ_y=K₀τ_z (14)

τ_z=(1- β) (ρ_s-ρ_f)gz (15)

The span of δ is 0.1～0.3, K₀=0.5, g are acceleration of gravity, and z is the depth below water-sediment interface Degree, δ_b=0.1, δ_s=0.1.

F is determined by (16)-(19) formula：

Attenuation quotient α is determined by (20)-(24) formula：

A=C²-HM (20)

Attenuation quotient α=Im { k₁Im represents imaginary part (24)

(β, κ f) are determined reflection R by (25) formula：

Wherein, B₁=ρ_fω², B₂=Hk₁ ²-Ck₁ ²G₁, B₃=Hk₂ ²-Ck₂ ²G₂C₂=Mk₁ ²G₁-Ck₁ ², C₃=Mk₂ ²G₂-Ck₂ ²；

2.1 is as follows from the computational methods of shallow seismic profile extracting data reflection R (f)：Transmitting and reception signal, from Sub-bottom profiler receives the reflection coefficient that water-top layer sediment interface is extracted in signal.Transmission signal S (f)=E (f) H_s F (), received spectrum can be expressed as：

r₁For transducer to water-bed distance, E (f) for source signal frequency spectrum, H_s(f) for Sonar system transfer function, E^* F () is the complex conjugate function of E (f), R (f) is reflection coefficient.The secondary wave signal of change reflection coefficient received using sonar, two Spectrum of the subwave after compression is represented by：

r₂For the water-bed distance of transducer distance and depth of water sum.Water-bed reflected signal reflects through air-water interface, Then received by transducer by underwater reflection, reception signal is secondary wave.Reflection coefficient is -1.Formula (26) is obtained divided by formula (27) To reflection coefficient：

2.2 relaxation time t_rWith the method for building up of reflected signal frequency displacement Δ f relation curves

The computational methods of signal transient frequency：

S (t) is signal,For Hilbert conversion, E (t) be function envelope, φ (t) be phase place, φ (t)=- jlog_e [z (t)/E (t)], the instantaneous frequency of s (t) is

The signal of sub-bottom profiler transmitting is linear FM signal S (f), is calculated by the computational methods of signal transient frequency Instantaneous frequency f of S (f)₁。

The transfer function of silt:

Wherein f is frequency, and x is propagation distance of the sound wave in silt,V is propagation of the sound wave in silt Speed.Only consider frequency shift effect of the sound wave in communication process, reflected signal is expressed as B'(f)=S (f) H_a(f).Take x=1m, By the computational methods of signal transient frequency, numerical computations obtain B'(f) instantaneous frequency f₂, Δ f=f₁-f₂.Obtain t_rWith frequency displacement The relation curve of Δ f；

2.3 methods for calculating attenuation quotient α：

According to 2.2, from A'(f) middle extraction river-top layer silt reflected signal instantaneous frequency f₁With top layer silt-bottom mud Instantaneous frequency f of husky interface reflected signal₂, Δ f=f₁-f₂, by t_rRelaxation time t is obtained with the relation curve of frequency displacement Δ f_r,

By formula (30), attenuation quotient of the sound wave in silt：

α=k'f² (31)

2.4 methods that porosity is searched according to porosity-reflection coefficient curve：

(unit is phi, φ=- log to porosity β with median particle diameter φ₂d_mm, d_mmFor the particle diameter of sand grain) experience it is public Formula

β=0.208+0.0943 φ -0.00334 φ² (32)

The empirical equation of permeability κ and porosity β, i.e. Kozeny-Carman formula

D is silt median particle diameter, and K is that empirical is estimated, when granule is circular tube shaped, K=2；When granule is spherical, K =5,

By formula (25) and formula (33), the relation curve template of median particle diameter d, porosity β, frequency f and reflection coefficient is obtained.

A kind of method of detection river bed top layer silt physical property, its step is as follows:

(1) shallow seismic profile sonic data is extracted, by 2.1 acoustic reflection coefficient R (f) is calculated；

(2) by 2.2, shallow seismic profile data are analyzed, calculates frequency displacement Δ f, find relaxation time t_r；

(3) FM signal mid frequency is taken by 2.3, f, calculates attenuation quotient α；

(4) by 2.4, porosity is searched according to porosity-reflection coefficient curve, convolution (32), (33) calculate permeability κ, here κ is only as initial value；

(5) silt velocity of wave c is calculated by formula (1) and formula (33)；

(6) permeability-attenuation quotient curve is calculated, according to the attenuation quotient α of step (2), inquires about permeability；

(7) repeat step (4)-(7), calculate, when porosity β through successive ignition_nWith previous iteration value β_n-1Difference is less When (| β_n-1-β_n|/β_n<3%), iteration is stopped；

(8) silt density p is calculated by formula (2), reflection R (f) is calculated by formula (25)；

(9) permeability κ, silt velocity of wave c, reflection R (f), attenuation quotient α, silt density p are exported.

Beneficial effect：

(1) the method can obtain various silt physical parameters, enrich shoal materials information.

(2) the method can obtain underwater sediment physical parameter merely with shallow bottom section plotter, use manpower and material resources sparingly, and improve Work efficiency.

(3) the method can obtain more accurate silt parameter information, and for the research such as river dynamics data supporting is provided.

Description of the drawings

Fig. 1 is sampling process schematic diagram of the present invention；

Fig. 2 is calculation flow chart of the present invention；

Fig. 3 is relation curve (f=3700Hz) figure of reflection R (f) and porosity β；

Fig. 4 is unit apart from frequency displacement Δ f and relaxation time t_rRelation curve；

Fig. 5 is that sub-bottom profiler receives 47 μ s of signal sampling interval.

Specific embodiment

With reference to embodiment, technical scheme is further described in detail.

Embodiment 1

The method of a kind of detection river bed top layer silt physical property, by taking Xiaolangdi reservoir area as an example：

The measuring point of 36 section of Xiaolangdi reservoir area 1, sub-bottom profiler transmitted signal bandwidth 500-7000Hz, Shi Kuan 20ms.The frequency displacement Δ f at the upper and lower interface of top layer silt be 88Hz, the frequency displacement Δ f=18Hz of unit distance, attenuation quotient α= 0.04, reflection R=0.17.Through calculating, permeability κ=4.46 × 10^-13, porosity β=67.66%, silt velocity of wave c =1370m/s, silt density p=1549.4kg/m³。

Embodiment 2

The method of a kind of detection river bed top layer silt physical property, by taking the section of Reservoir Area of Sanmenxia 4 as an example：

The measuring point of 4 section of Reservoir Area of Sanmenxia 1, sub-bottom profiler transmitted signal bandwidth 500-7000Hz, when width 20ms.It is single Position distance frequency displacement Δ f be 71Hz, attenuation quotient α=0.094, reflection R=0.16.Through calculating, permeability κ=3.89 ×10^-11, porosity β=0., silt velocity of wave c=1503m/s, silt density p=1412kg/m³。

Embodiment 3

The method of a kind of detection river bed top layer silt physical property, by taking the section of Reservoir Area of Sanmenxia 8 as an example：

The measuring point of 8 section of Reservoir Area of Sanmenxia 1, sub-bottom profiler transmitted signal bandwidth 500-7000Hz, when width 20ms.It is single The frequency displacement Δ f=104Hz of position distance, attenuation quotient α=0.136, reflection R=0.34.Through calculating, permeability κ= 6.49×10^-12, porosity β=0.46, silt velocity of wave c=1623m/s, silt density p=1921kg/m³。

Although above with a general description of the specific embodiments the present invention is described in detail, On the basis of the present invention, it can be made some modifications or improvements, this will be apparent to those skilled in the art.Cause This, without departing from theon the basis of the spirit of the present invention these modifications or improvements, belong to the scope of protection of present invention.

Claims (2)

1. the method for a kind of detection river bed top layer silt physical property, it is characterised in that：Involved formula and parameter is as follows：

K is wave number,For angular frequency, f is frequency, and c is pore media velocity of wave, ρ_fAnd ρ_sRespectively pore-fluid density and silt Granule density, and for known constant, ρ is silt body density, and η represents pore water viscosity, and κ is permeability, K_γIt is the body of sand grain Product module amount, m is the modulus of shearing of silt framework, K_bAnd K_fThe bulk moduluses of silt skeleton modulus and pore-fluid, F are represented respectively The deviation produced for pore-fluid increases with frequency, the velocity of wave c of water_ω=1500kg/m³, K_r=3.2 × 10¹⁰Pa, K_f= 2.395×10⁹Pa, η=0.001Kg/m*s；

$k=\frac{\stackrel{\&OverBar;}{\mathrm{\&omega;}}}{c},\stackrel{\&OverBar;}{\mathrm{\&omega;}}=2\mathrm{\&pi;}f---\left(1\right)$

ρ=β ρ_f+(1-β)ρ_s (2)

$H=\frac{{(K_r-K_b)}^2}{D-K_b}+K_b+\frac{4\mathrm{\&mu;}}{3}---\left(3\right)$

$C=\frac{K_r(K_r-K_b)}{D-K_b}---\left(4\right)$

$M=\frac{K_r^2}{D-K_b}---\left(5\right)$

$D=K_r\&lsqb;1+\mathrm{\&beta;}(\frac{K_r}{K_f}-1)\&rsqb;---\left(6\right)$

$K_b=K_{br}+{\mathrm{jK}}_{bi}=K_{br}(1+j\frac{{\mathrm{\&delta;}}_b}{\mathrm{\&pi;}})---\left(7\right)$

$\mathrm{\&mu;}={\mathrm{\&mu;}}_r(1+j\frac{{\mathrm{\&delta;}}_s}{\mathrm{\&pi;}})---\left(8\right)$

$K_{br}=\frac{2{\mathrm{\&mu;}}_r(1+\mathrm{\&delta;})}{3(1-2\mathrm{\&delta;})}---\left(9\right)$

${\mathrm{\&mu;}}_r=1.835\&times;{10}^5{e}^{-1.12}\sqrt{{\mathrm{\&tau;}}_a\left(z_s\right)}---\left(10\right)$

$e=\frac{\mathrm{\&beta;}}{1-\mathrm{\&beta;}}---\left(11\right)$

${\mathrm{\&tau;}}_a=\frac{1}{3}({\mathrm{\&tau;}}_x+{\mathrm{\&tau;}}_y+{\mathrm{\&tau;}}_z)---\left(12\right)$

τ_x=K₀τ_z (13)

τ_y=K₀τ_z (14)

τ_z=(1- β) (ρ_s-ρ_f)gz (15)

The span of δ is 0.1～0.3, K₀=0.5, g are acceleration of gravity, and z is water-sediment interface depth below, δ_b =0.1, δ_s=0.1；

F is determined by (16)-(19) formula：

$F\left(\mathrm{\&epsiv;}\right)=\frac{\frac{\mathrm{\&epsiv;}}{4}T\left(\mathrm{\&epsiv;}\right)}{1-\frac{2j}{\mathrm{\&epsiv;}}T\left(\mathrm{\&epsiv;}\right)}---\left(16\right)$

$T=\frac{-\sqrt{j}{J}_1\left(\mathrm{\&epsiv;}\sqrt{j}\right)}{J_0\left(\mathrm{\&epsiv;}\sqrt{j}\right)}---\left(17\right)$

$\mathrm{\&epsiv;}=a_p\sqrt{\frac{{\mathrm{\&omega;\&rho;}}_{\mathrm{\&omega;}}}{\mathrm{\&eta;}}}---\left(18\right)$

$a_p=\sqrt{\frac{8\&times;1.35\mathrm{\&kappa;}}{\mathrm{\&beta;}}}---\left(19\right)$

Attenuation quotient α is determined by (20)-(24) formula：

A=C²-HM (20)

$b={\mathrm{hm\&omega;}}^2+{\mathrm{\&rho;\&omega;}}^{2}M-jH\frac{\mathrm{\&omega;}F\mathrm{\&eta;}}{\mathrm{\&kappa;}}-{\mathrm{C\&rho;}}_f{\mathrm{\&omega;}}^2-{\mathrm{C\&rho;}}_f{\mathrm{\&omega;}}^2---\left(21\right)$

$c=-{\mathrm{\&rho;m\&omega;}}^4+j\frac{{\mathrm{\&rho;F\&eta;\&omega;}}^3}{\mathrm{\&kappa;}}+{\mathrm{\&rho;}}_f{\mathrm{\&omega;}}^4---\left(22\right)$

$k_1=\sqrt{\frac{-b+\sqrt{b^2-4ac}}{2a}},k_2=\sqrt{\frac{-b-\sqrt{b^2-4ac}}{2a}}---\left(23\right)$

Attenuation quotient α=Im { k₁Im represents imaginary part (24)

(β, κ f) are determined reflection R by (25) formula：

$R(\mathrm{\&beta;},\mathrm{\&kappa;},f)=\frac{(B_2-A_2{B}_1)(B_2{C}_3-C_2{B}_3)+(B_2{A}_3-A_2{B}_3)(B_1{B}_2+C_2{B}_1)}{(B_2+A_2{B}_1)(B_2{C}_3-C_2{B}_3)-(B_2{A}_3-A_2{B}_3)(B_1{B}_2+C_2{B}_1)}---\left(25\right)$

Wherein, B₁=ρ_fω², B₂=Hk₁ ²-Ck₁ ²G₁, B₃=Hk₂ ²-Ck₂ ²G₂C₂=Mk₁ ²G₁-Ck₁ ², C₃=Mk₂ ²G₂-Ck₂ ²；

2.1 is as follows from the computational methods of shallow seismic profile extracting data reflection R (f)：Transmitting and signal is received, from shallow Layer section plotter receives the reflection coefficient that water-top layer sediment interface is extracted in signal, transmission signal S (f)=E (f) H_sF (), connects Receiving signal spectrum can be expressed as：

$A^{\&prime;}\left(f\right)=\frac{1}{2r_1}E\left(f\right)H_s\left(f\right)E^{*}\left(f\right)R\left(f\right)---\left(26\right)$

r₁For transducer to water-bed distance, E (f) for source signal frequency spectrum, H_s(f) for Sonar system transfer function, E^*(f) For the complex conjugate function of E (f), R (f) is reflection coefficient.The secondary wave signal of change reflection coefficient received using sonar, secondary wave Spectrum after compression is represented by：

$C^{\&prime;}\left(f\right)=\frac{-1}{2r_2}E\left(f\right)H_s\left(f\right)E^{*}\left(f\right)R^2\left(f\right)---\left(27\right)$

r₂It is distance and depth of water sum of the transducer apart from the bottom, water-bed reflected signal reflects through air-water interface, then Received by transducer by underwater reflection, reception signal is secondary wave.Reflection coefficient is -1.Formula (26) obtains anti-divided by formula (27) Penetrate coefficient：

$R\left(f\right)=-\frac{r_2{C}^{\&prime;}\left(f\right)}{r_1{A}^{\&prime;}\left(f\right)}---\left(28\right);$

2.2 relaxation time t_rWith the method for building up of reflected signal frequency displacement Δ f relation curves

The computational methods of signal transient frequency：

$z\left(t\right)=s\left(t\right)+j\stackrel{\&CenterDot;}{s}\left(t\right)=E\left(t\right)e^{j\mathrm{\&phi;}\left(t\right)}---\left(29\right)$

S (t) is signal,For Hilbert conversion, E (t) be function envelope, φ (t) be phase place, φ (t)=- j log_e[z (t)/E (t)], the instantaneous frequency of s (t) is

The signal of sub-bottom profiler transmitting is linear FM signal S (f), and by the computational methods of signal transient frequency S (f) is calculated Instantaneous frequency f₁,

The transfer function of silt:

$H_a\left(f\right)=e^{-k^{\&prime;}{f}^{2}x}---\left(30\right)$

Wherein f is frequency, and x is propagation distance of the sound wave in silt,V is spread speed of the sound wave in silt, Only consider frequency shift effect of the sound wave in communication process, reflected signal is expressed as B'(f)=S (f) H_aF (), takes x=1m, by believing The computational methods of number instantaneous frequency, numerical computations obtain B'(f) instantaneous frequency f₂, Δ f=f₁-f₂.Obtain t_rWith frequency displacement Δ f Relation curve；

2.3 methods for calculating attenuation quotient α：

According to 2.2, from A'(f) middle extraction river-top layer silt reflected signal instantaneous frequency f₁Hand over top layer silt-bottom sediment Instantaneous frequency f of interface reflected signal₂, Δ f=f₁-f₂, by t_rRelaxation time t is obtained with the relation curve of frequency displacement Δ f_r,

By formula (30), attenuation quotient of the sound wave in silt：

α=k'f² (31)

2.4 methods that porosity is searched according to porosity-reflection coefficient curve：

(unit is phi, φ=- log to porosity β with median particle diameter φ₂d_mm, d_mmFor the particle diameter of sand grain) empirical equation

β=0.208+0.0943 φ -0.00334 φ² (32)

The empirical equation of permeability κ and porosity β, i.e. Kozeny-Carman formula

$\mathrm{\&kappa;}=\frac{d^2{\mathrm{\&beta;}}^3}{36K{(1-\mathrm{\&beta;})}^2}---\left(33\right)$

D is silt median particle diameter, and K is that empirical is estimated, when granule is circular tube shaped, K=2；When granule is spherical, K=5, By formula (25) and formula (33), the relation curve template of median particle diameter d, porosity β, frequency f and reflection coefficient is obtained.

2. a kind of method of detection river bed top layer silt physical property, its step is as follows:

(1) shallow seismic profile sonic data is extracted, by 2.1 acoustic reflection coefficient R (f) is calculated；

(2) by 2.2, shallow seismic profile data are analyzed, calculates frequency displacement Δ f, find relaxation time t_r；

(3) FM signal mid frequency is taken by 2.3, f, calculates attenuation quotient α；

(4) by 2.4, porosity is searched according to porosity-reflection coefficient curve, convolution (32), (33) calculate permeability κ, this In κ only as initial value；

(5) silt velocity of wave c is calculated by formula (1) and formula (33)；

(6) permeability-attenuation quotient curve is calculated, according to the attenuation quotient α of step (2), inquires about permeability；

(7) repeat step (4)-(7), calculate, when porosity β through successive ignition_nWith previous iteration value β_n-1When difference is less (| β_n-1-β_n|/β_n＜ 3%), stop iteration；

(8) silt density p is calculated by formula (2), reflection R (f) is calculated by formula (25)；

(9) permeability κ, silt velocity of wave c, reflection R (f), attenuation quotient α, silt density p are exported.