• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    Frequency-Domain Full-Waveform Inversion Based on Tunnel Space Seismic Data

    2022-02-13 09:18:56MingyuYuFeiChengJingpingLiuDichengPengZhijinTin
    Engineering 2022年11期

    Mingyu Yu, Fei Cheng*, Jingping Liu, Dicheng Peng, Zhijin Tin

    a Hubei Subsurface Multiscale Imaging Key Laboratory (SMIL), Institute of Geophysics and Geomatics, China University of Geosciences, Wuhan 430074, China

    b Hubei Key Laboratory of Marine Geological Resources, China University of Geosciences, Wuhan 430074, China

    Keywords:Tunnel seismic detection Full-waveform inversion Frequency select strategy Observation system

    A B S T R A C T Tunnel seismic detection methods are effective for obtaining the geological structure around the tunnel face,which is critical for safe construction and disaster mitigation in tunnel engineering.However,there is often a lack of accuracy in the acquired geological information and physical properties ahead of the tunnel face in the current tunnel seismic detection methods.Thus,we apply a frequency-domain acoustic full-waveform inversion(FWI)method to obtain high-resolution results for the tunnel structure.We discuss the influence of the frequency group selection strategy and the tunnel observation system settings regarding the inversion results and determine the structural imaging and physical property parameter inversion of abnormal geological bodies ahead of the tunnel face. Based on the conventional strategies of frequency-domain acoustic FWI, we propose a frequency group selection strategy that combines a low-frequency selection covering the vertical wavenumber and a high-frequency selection of antialiasing.This strategy can effectively obtain the spatial structure and physical parameters of the geology ahead of the tunnel face and improve the inversion resolution.In addition,by linearly increasing the side length of the tunnel observation system, we share the influence of the length of the two sides of the observation systems of different tunnels on the inversion results.We found out that the inversion results are the best when the side length is approximately five times the width of the tunnel face,and the influence of increasing the side observation length beyond this range on the inversion results can be ignored.Finally, based on this approach, we invert for the complex multi-stratum model, and an accurate structure and physical property parameters of the complex stratum ahead of the tunnel face are obtained,which verifies the feasibility of the proposed method.

    1. Introduction

    With the development of global economies and the increasing demand for modernization in construction, tunnels have become one of the optimal choices for traffic infrastructure construction in complex geological terrains. The geological conditions of the surrounding rock in a tunnel construction environment are complex, and various geological disasters, such as a collapse or water inrush, are often caused by unfavorable geological bodies, such as karsts and fracture zones. Because tunnel construction occurs in the geological body, there can be a substantial number of casualties and property losses during construction if the unfavorable geological bodies (belt) around the tunnel are not accurately predicted [1]. Therefore, it is critical to detect hidden geological hazards ahead of the tunnel face in advance to reduce the hidden dangers and ensure the safety of the construction site.

    Advanced tunnel detection is a type of geophysical technology used to detect hidden geological hazards ahead of the tunnel face using observation systems. Current geophysical methods used for advanced tunnel detection primarily include seismic, electromagnetic,electrical,and geological radar methods[2-5],among which the tunnel seismic detection method has become common for its long detection range and accurate for the prediction. This technique is based on the difference in the seismic wave velocity between abnormal geological bodies and the surrounding rock.During construction, as seismic waves are transmitted to the surrounding rock of the tunnel and the collected seismic data is processed, the distribution of abnormal geological bodies ahead of the tunnel face and the mechanical parameters of the rock are obtained for an early warning and guidance during tunnel construction. In the late 1970s, engineers in Germany and England explored the geological structure ahead of the tunnel working face by using the airy seismic phase of the channel wave. In the early 1990s, the Swiss Surveying Technology Co., Ltd. developed a set of advanced tunnel seismic prediction (TSP) systems. In the late 1990s, a US engineering company developed the true reflection tomography (TRT) technology, Zeng [6] and Inazaki et al. [7] proposed the vertical seismic profile method for tunnels, then developed into application by Zhao et al. [8] and Alimoradi et al. [9] at the beginning of 21st century.Based on the increase in tunnel construction, the currently available seismic advanced geological prediction methods include the TSP, horizontal seismic profiling(HSP),TRT,tunnel seismic tomography(TST),tunnel seismic while drilling(TSWD),tunnel geology prediction(TGP),and the negative apparent velocity method of seismic waves [10-15]. However,the observation system is limited to the tunnel environment, and the amount of data collected cannot meet the requirements of a high calculation accuracy for the tunnel geological body wave velocity.The imaging results obtained using the tunnel seismic detection method cannot accurately detect abnormal geological bodies. This method is only suitable for simple geological conditions[16].More accurate detection methods suitable for advanced tunnel geological prediction are necessary.

    The development of full-waveform inversion(FWI),which is an inversion technique that uses full wave field information to invert medium parameters in the seismic exploration field, has provided several opportunities. In the 1980s, researchers proposed a timedomain FWI based on the least-squares method and introduced this concept to the seismic exploration field [17-19]. Compared to the traditional inversion method that uses a single reflected wave or the first arrival wave data to obtain the attributed parameter imaging, FWI fully utilizes the full wave field information to achieve a higher resolution [20]. Therefore, this high-precision and high-resolution inversion method has been significantly praised for seismic wave field inversion and reconstruction and has gained increasing attention in the research and application of seismic velocity modeling [21,22]. Owing to the nonlinearity and cycle-skipping of FWI, the objective function has multiple local minima, which makes the inversion significantly dependent on the initial model [23-25]. To reduce the dependence of the inversion results on the initial model, early researchers proposed a multi-scale FWI method in the time domain that filters the seismic data to isolate frequencies [26]. In the 1990s, Pratt and Worthington [27] extended the theory of frequency-domain FWI.The inversion results of the low-frequency components in the frequency-domain FWI can be used as the initial model of the high-frequency components, which can directly achieve the effect of multi-scale inversion and reduce the dependence on the initial model [27-31]. Owing to this advantage, the frequency-domain FWI is widely used in seismic exploration.

    For advanced tunnel detection applications, Musayev et al.[32,33]first applied the full waveform inversion method in the frequency domain to tunnel and discussed whether the full waveform can successfully image the velocity field in a tunnel. Nguyen and Nestorovic′ [34,35] proposed a global optimization procedure for the FWI of two-dimensional (2D) tunnel seismic waves; they also used the elastic FWI enhanced with the parametric representation for locating the disturbance zones ahead of a tunnel face.Bharadwaj et al. [36] developed a seismic prediction system to enable imaging ahead of a tunnel-boring machine. Lamert et al.[37] proposed two flexible elastic time-domain FWI methods to predict the disturbance area ahead of a tunnel face.As a local case study, Zhang et al. [38] used the FWI method with ground penetrating radar (GPR) to distinguish unfavorable geological bodies within 20 m ahead of a tunnel face. Li [39] used the acoustic FWI method to predict the velocity interface of large geological bodies ahead of a tunnel face.Feng et al.[40]improved the reconstruction of tunnel lining defects using GPR profiles with FWI.

    Considering that the environment owns a limited tunnel seismic detection space and observation coverage, to improve the accuracy of tunnel seismic detection, we tested the frequencydomain acoustic FWI method for tunnel seismic detection. Using the abnormal low-speed body as example,we constructed a tunnel low-speed body model and its observation system based on the tunnel space and reconstructed the tunnel velocity model with frequency-domain acoustic FWI.By comparing different frequency group selection strategies of frequency-domain FWI, we analyzed the results of the frequency group selection and determined suitable options for the tunnel seismic method. Herein, we discuss the influence of the side length of the tunnel observation system on the inversion. Finally, we used a complex tunnel geological model to verify the effectiveness of the method and the parameter selection strategy.

    2. Frequency-domain acoustic FWI

    2.1. Theory of 2D frequency-domain acoustic FWI

    In an isotropic medium, the 2D acoustic wave equation in the frequency domain is expressed as follows [41]:

    where k(x, z)is the bulk modulus,ρ(x,z) is the density, d(x,z,ω) is the pressure field,(x,z)represents the 2D coordinates,ω is the frequency, and s(x, z,ω) is the source function.

    Because the pressure field d(x,z,ω)is linear with respect to the source s(x,z,ω), the discretized 2D acoustic wave equation can be simplified into the following large sparse linear equation:

    where A represents the impedance matrix of the frequency and medium properties. Considering that d(x, z, ω) and s(x, z, ω) are stored as vectors of dimension Nx× Nz, A(ω) is a finite difference operator matrix of (Nx,Nz) × (Nx,Nz), which can be solved using the lower-upper (LU) decomposition method. The wave equation(Eq. (1)) is discretized by the mixed-grid finite-difference (FD)method; in addition, the perfectly matched layer (PML) absorbing boundary condition is used to simulate the virtual boundary in the model and the free surface at the inner tunnel surface [42,43].

    We use the L2norm as our frequency-domain acoustic FWI’s objective function; it is given by the squared difference between the observed and calculated data, as follows:

    where Nω is the number of frequencies in the group, Nsrepresents the number of sources,Nrrepresents the number of receivers,dobs(s,r, ω) is the observed wave field, and dcal(s, r, ω) is the calculated wave field using the parameters of model m, H is conjugate transpose. E(m) is a function related to the model parameter m, and its gradient is calculated as follows:

    where JTis the transpose of the Jacobian matrix, which is derived from the partial derivatives of the wave field on the model parameters. Δd* is the complex conjugation of the wave field residuals,and Re is the real part of the complex number. The gradient of E(m) is calculated to obtain the perturbation of the model parameters in iterations to minimize the misfit between the observed and calculated wave fields [44]. The inversion error threshold and iteration number are set as the termination conditions for updating the model to ensure convergence at a practical cost. The model medium parameters satisfying the conditions are obtained when the inversion is terminated.

    2.2. Strategy for FWI frequency group selection in tunnel space

    2.2.1. Strategy for frequency group selection

    Frequency-domain FWI allows us to obtain large-scale information by first inverting low-frequency data corresponding to long wavelengths, then the information retrieved from middle- and high-frequency short wavelengths can depict detailed features. In the process of FWI, the low-frequency inversion results are used as the initial model for the subsequent middle- and highfrequency inversions, which can directly reach multi-scale inversion. Because the probability of non-convergence caused by the local extremum of the low-frequency data inversion is small,a relatively good initial model can be estimated,which can improve the convergence stability of the inversion process and accelerate the inversion convergence [45].

    In this manner,the selection of frequency groups for the FWI in the frequency domain is worth studying.Sirgue and Pratt[46]conducted a frequency selection method based on the continuity of vertical wavenumber coverage,in which the maximum wavenumber corresponding to the frequency value of the current stage should be equal to the minimum wavenumber corresponding to the frequency value of the next frequency stage, as follows:

    where fnrepresents the frequency value of the current stage, fn+1is the frequency value of the next frequency stage, and the vertical wavenumber kzis the vertical component of the wavenumber vector k. The range of kzdepends on the incident angle of the seismic wave and its relationship with the calculation results in the following equation:

    where α is the cosine of the maximum incident angle, hmaxis the maximum half offset of the current observation system, and z is the depth of the reflection layer.

    An extensive offset, high-density geophone distribution observation system, and horizontal reflecting strata are apparently required. However, the observation system of the tunnel space is limited,and the geological conditions of tunnels are often complex.This strategy can only obtain limited depth information ahead of the tunnel face. In addition, to satisfy the antialiasing condition,the sampling rate of the vertical wavenumber Δkzshould satisfy Δkz≤1/zmax[47,48]. Then, the relationship between the vertical wave number and background velocity in the calculation can be used to obtain the following formula:where c0is the velocity of the background, Δf represents the frequency value between two frequency stage intervals, and zmaxis the maximum imaging depth.

    This strategy for frequency sampling has been proven to have wider application possibilities, despite the relatively slow convergence rate of inversion due to the decrease in sampling points in the low-frequency region.

    2.2.2. Strategy for frequency group selection of tunnel FWI

    Tunnel seismic detection systems often have a limited offset of observation and its records own a large frequency distribution range with a high dominant frequency. Therefore, the frequency selection strategies in Section 2.2.1 are no longer applicable in the tunnel space. To improve the inversion resolution of tunnel FWI, we propose a method that combines the advantages of both as follows:

    In addition, we suggest using the background velocity vsfor the elastic FWI.

    Eq. (8) is a combined formula as shown in Fig. 1, in the lowfrequency range, it is used to obtain a more detailed lowfrequency group to ensure the large-scale inversion effect of the model. Eq. (7) is used as the judgment condition in Eq. (8), when the frequency group obtained cannot satisfy the antialiasing criterion, a constant frequency interval is chosen to ensure it.

    3. Discussion of parameters and calculation results

    3.1. Observation system and model design

    The tunnel seismic detection method arranges excitation holes in the tunnel face and the sidewalls with depths of 0.5 or 1.0 m.Seismic waves propagate to the wall rock of the tunnel through artificial excitation. When the impedance of these waves changes,part of the seismic waves will be reflected,and the other part continues to spread forward.The reflected seismic waves are recorded by the receivers and provide the seismic records. The processed records can predict geological changes ahead of the tunnel face and provide reliable geological data for tunnel construction[49,50]. Based on actual tunnel parameters, we built the tunnel low-speed anomalies model shown in Fig. 2(a), which is 200 m(X axis) × 30 m (Z axis), with a 100 m tunnel length and a 12 m tunnel face (width), along with three anomalies with a radius of 3 m located 14, 46, and 84 m ahead of the tunnel face. The tunnel space is full of air, the rock wall velocity is 4000 m·s-1, and the velocity of the anomalies is 3000 m·s-1.The observation system of the tunnel seismic detection was set as U-shaped, as shown in Fig. 2(b), including the tunnel face and both sidewalls. With a 1 m excitation hole depth, 57 shot points and 112 receivers are arranged from 51 to 101 m in the sidewalls and tunnel face with a 2 m shot interval length and a 1 m receiver interval length. In the forward modeling,the dominant frequency of the Ricker wavelet is 200 Hz; Figs. 2(c)-(e) present the real part of the monofrequency pressure field slices of 50,200,and 450 Hz.As indicated,when the velocity remains unchanged and the frequency increases,the wavelength of a single frequency slice in the frequency domain becomes shorter, the difference between the wavelength of the low-velocity anomaly body and that of the wall rock increases,and the details regarding the anomalous body position become more apparent.

    Fig. 1. Frequency group curve. In the beginning frequency range, the frequency interval Δf increases linearly with the frequency group number,when increasing to a threshold the frequency increases with equal Δf to satisfy the antialiasing condition.

    Fig.2. Tunnel model with observation system design and real parts of the pressure field in the frequency-domain.(a)Low speed anomalies model with a tunnel,where X is the length of the tunnel model and Z is the width. (b) Observation system of a tunnel seismic detection. (c)-(e) Real part of a mono-frequency slice for 50,200, and 450 Hz with shot point located in the middle of the tunnel face.

    With a sampling interval of 0.5 ms and a recording length of 160 ms, the acoustic wave pressure component records in the time-domain were obtained through the inverse Fourier transform from the frequency-domain.The shot point in Fig.3(a)is located at coordinates(51, 9), which is the beginning of the observation system on the side of the tunnel.The shot point in Fig.3(b)is located in the middle of the tunnel face with coordinates (102, 15). These records contain direct waves and diffracted waves caused by the interface of three low velocity anomalous bodies.

    3.2. Discussion of different frequency selection strategies

    Because low-frequency inversion provides long-wavelength information while high-frequency provides detailed wave propagation information, multi-scale FWI using the low-frequency inversion result as the initial model to invert for the higher frequency data helps the high-frequency inversion converge as well as describe the details [51]. Using the observation system shown in Fig.2(b),a frequency-domain acoustic FWI of the abnormal tunnel body model, shown in Fig. 2(a), was conducted. The forward simulation in the inversion process adopts a Ricker wavelet with a dominant frequency of 200 Hz. A homogeneous velocity model with a velocity of 3000 m·s-1was selected as the initial model.Based on the frequency spectrum of the Ricker wavelet, three frequencies(50,200,and 450 Hz)were selected from low to high.We obtained three root mean square error(RMSE)convergence curves and velocity models after the inversion, as shown in Fig. 4.

    Fig.3. Records of two shot points for the designed model.(a)Records of Fig.2(a)with a shot point located the side of the tunnel.(b)Records of Fig.2(a)with the shot point located at the middle of the tunnel face. These records contain direct waves and diffracted waves caused by the interface of three low velocity anomalous bodies.

    Based on Eqs. (6)-(8), three frequency groups between 10 and 500 Hz were obtained.To compare the results of the inversion frequency selection strategy,the frequency group obtained by Eq.(6)is named S,while Eq.(7)describes the W frequency group,and Eq.(8)describes the C frequency group.The wavelet spectrum and frequency group parameter curves are shown in Figs. 5(a) and (b),respectively.

    The three frequency groups shown in Fig.5(b)were used for the FWI with the same initial model shown in Fig. 6(a) and same iterations. The three velocity resolution results are shown in Fig. 6.Figs. 6(b)-(d) W, S, and C frequency groups, respectively. Compared to Fig. 4(d), the inversion respresent the velocity inversion result of theults shown in Figs.6(b)-(d)are significantly improved;however,the resolution results shown in Figs.6(b)-(d)remain different. First, we discuss the sensitivity of the velocity inversion result to the location of the abnormal bodies and the imaging accuracy of the inversion results to compare the effects of different frequency selection strategies. The inversion result of the W frequency group showed in Fig. 6(b) is accurate considering the location of the first shallow abnormal body.Additionally,its velocity profile curve has a relatively noticeable disturbance in the position range of the middle and deep abnormal bodies and its adjacent position shown in Fig.6(e),indicating that the inversion converges beyond the range of the abnormal bodies. The inversion result of group S is relatively accurate for the location of the shallow abnormal bodies shown in Fig.6(c);the velocity variations in Fig.6(e)are apparent in the location of the first abnormal body; however, are not as apparent in the location of the middle and deep abnormal bodies. The inversion result of the C frequency group is significantly close to reality in the shallow, middle, and deep abnormal body locations shown in Fig. 6(d), and the velocity curve of the middle profile shown in Fig. 6(e) presents a noticeable and accurate velocity disturbance in the position of the abnormal bodies.Overall, the inversion results obtained by the C frequency group have better convergence sensitivity in the position of the abnormal bodies, which demonstrates the contribution of the frequency group selection to the convergence stability of multi-scale FWI.

    The velocity profiles shown in Fig. 6(e) and Fig. 7 are comprehensively compared to discuss the accuracy of the inversion results.The inversion accuracy of the abnormal bodies of the three frequency groups is affected by the distance from the tunnel face,as shown by the comparison velocity profiles; the accuracy worsens with the distance, and the velocity inversion accuracies of the middle and far abnormal bodies are worse than those that are close.A comparison of the three frequency groups reveals that the velocity inversion accuracy of the S frequency group is relatively poor, whereas the velocity inversion accuracy of the W and C frequency groups is better.Their difference in the shallow section is minimal between the inversion results and the initial models,and the inversion result of the C frequency group reflects a higher resolution in the middle and deep sections.

    In summary, we demonstrated that the frequency-domain acoustic multi-scale FWI based on the tunnel seismic detection method can obtain satisfactory velocity inversion results, and the suitable frequency group selection method proposed in this study can obtain superior resolution results ahead of the tunnel face.

    Fig.4. Iteration of multi-frequency inversion.(a)Iteration curves for the inversion of three frequency groups;inversion result of frequencies of(b)50 Hz,(c)50 and 200 Hz,and (d) 50, 200, and 450 Hz. (a) demonstrates the convergence process of three frequencies in the iteration, and (b)-(d) describe the process of the multi-scale inversion results from low to high frequency, in which the inversion results can be observed from determining the location of the abnormal bodies.

    Fig.5. Parameters for acoustic FWI.(a)Ricker wavelet spectrum of the source wavelet;(b)parameters of three frequency groups.Based on Fig.5(a),the frequency range from 10 to 500 Hz were obtained.

    Fig.6. Initial model and inversion results of different frequency groups.(a)Initial model,(b)-(d)inversion results for models of W,S,C group frequencies obtained by Eqs.(7),(6),and(8)frequency selection strategies,respectively;(e)velocity profile curves in the grid axis Z=15 of the tunnel models.The comparison of inversion results reveals the difference in the accuracy of inversion results of different frequency group methods for abnormal bodies at different depths.

    Fig.7. Velocity curve comparison of the cross-sections of the three abnormal bodies based on inversion results in Figs.6(b)-(d).(a)Velocity curve of first abnormal body in grid axis X=114 of the tunnel models;(b)velocity curve of second abnormal body in grid axis X=147 of the tunnel models;(c)velocity curve of third abnormal body in grid axis X = 186 of the tunnel models.

    3.3. Discussion of tunnel model observation system

    To obtain more effective information regarding the geological body ahead of the tunnel face,we suggest using a U-shaped observation system for tunnel seismic detection.Referring to the propagation path of seismic waves in the tunnel space and the energy loss in the process of seismic wave propagation, we consider that the extended range of the observation system on both tunnel sides would no longer affect the inversion results after reaching a certain length.To confirm this,we designed five U-shaped tunnel observation system groups with side lengths of 10, 30, 50, 70, and 90 m,while all other parameters remained constant. The specific inversion results after applying these parameters to the frequencydomain acoustic FWI in the tunnel space are shown in Figs. 8 and 9.

    A comparison of Figs. 8(b)-(d) with Fig. 9(a) reveals that the longer the observation system, the more accurate the velocity inversion result. From Figs. 8(d)-(f), when the side length of the tunnel observation system changes beyond a certain range, its influence on the inversion results decreases; increasing the side length has a very weak influence on the inversion results. As demonstrated in the comparison of velocity curves across the tunnel cross-section in Fig. 9(b), the difference between the inversion results of the 70 m side length and that of the 90 m side length is significantly small,and the relative error between them is less than 0.4%.According to the comparison results,for the tunnel face width of this model, considering a built-in geophone and source depth,the optimal resolution results can be obtained by selecting a tunnel observation system with a side length of 70 m.

    Specifically,through multiple comparisons and simulations,we advise that within limits,the longer the arrangement side length of the tunnel observation system, the better the resolution of the inversion result. However, when the side length of the tunnel observation system is more than five times the width of the tunnel face,the influence of increasing the observation system side length on the inversion results is negligible. Overall, this conclusion can potentially aid theoretical studies regarding tunnel FWI and the actual construction of tunnel seismic detection.

    3.4. Complex model calculation and application

    To verify the effect of the parameters discussed in this study, a tunnel stratum model including abnormal bodies, low-velocity zones, and complex structures is designed for frequency-domain acoustic FWI based on tunnel seismic detection (Fig. 10(a)).

    Fig.8. Initial model and inversion results model of varying tunnel observation system side lengths.(a)Initial model,and(b)-(f)resulting velocity models with 10,30,50,70,and 90 m tunnel observation system side lengths.It can be seen that the velocity results of the same abnormal body obtained by inversion with different observation system and the differences of abnormal bodies with different depths in the inversion results.

    Fig. 9. Velocity curve of different tunnel observation system models in axis Z = 15 cross-section based on inversion results in Figs. 8(b)-(f). (a) Velocity results of the side lengths 10, 30, and 50 m, and (b) velocity results of side lengths 50, 70, and 90 m.

    The true model is designed based on limestone strata with wave velocities of 4500 to 5500 m·s-1, including low-speed thin-layer abnormal body strata with a velocity of 3000-4000 m·s-1and other complex geological conditions. The initial model for the forward calculation in the inversion uses the background velocity of the observation model (Fig. 10(b)). The three frequency selection strategies were used as the selection basis for the frequency group parameters to obtain a set of frequencies for the inversion. Using the tunnel observation system with a side length of 70 m as the tunnel observation system, the final velocity inversion results are shown in Figs. 10(c)-(e).

    The results of velocity inversion in Figs. 10(c)-(e) demonstrate that the frequency-domain acoustic FWI of the multi-scale strategy can inversely affect the velocity of the shallow sections of the complex tunnel geological conditions,and the closer it is to the tunnel face, the more accurate the results are. For the shallow part, the difference between the inversion results and the true model is minimal in all the three results, whereas the inversion results obtained by the frequency group selection strategy proposed in this study perform better in the middle and deep inversion,which indicates that the resolution in Fig. 10(e) is better than the results provided in Figs. 10(c) and (d).

    Comparing the profile velocities in Figs.10(e)and(b),Fig.11(d)indicates that the velocity obtained from the inversion results accurately corresponds to that of the true model, whereas there is a slight misfit at the far side from the tunnel face. The relative error between the inversion result and the real model ranged from 0.3% to 8% in different positions and increased with depth. The error in the near part was the smallest, and that in the furthest position was 6%.

    The complex geological model inversion results of the frequency-domain acoustic FWI based on tunnel seismic detection prove that the method can successfully invert the stratum information and obtain high-resolution results of complex geological information ahead of the tunnel face under suitable parameters.

    For tunnel exploration, all theoretical research and discussions are for practical applications; therefore, we used the TSP observation system with the parameters discussed in this study for a simple field FWI verification.The observation system for the field data acquisition is shown in Fig. 12.

    Because the data acquisition is conducted in the time domain,the frequency-domain acoustic FWI requires wave field separation and Fourier transform. Therefore, the X-component field records(Fig. 13(a)) of the TSP were used for filtering and wave field separation to obtain the P-wave component records (Fig. 13(b)) [52],and the direct wave velocity was used as the initial model(Fig.13(c)), as shown in Fig. 13(d).

    Fig.10. Comparison of tunnel complex model.(a)True model of tunnel complex model inversion,(b)initial model of tunnel complex model inversion,(c)-(e)inversion result models of W, S, C frequency groups obtained by Eqs. (7), (6), and (8) frequency selection strategies, respectively. For the shallow part, the difference between the inversion results and the true model is minimal in all the three results, whereas the inversion results in (e) performs better in the middle and deep inversion.

    Fig. 11. Velocity comparison of tunnel complex model using C frequency selection strategy. (a)-(c) Velocity model of the initial, true, and result models and (d) velocity profiles of tunnel models from (a)-(c) in the grid axis Z = 15.

    Fig.12. Tunnel seismic exploration field observation system with two receivers and 24 shot points.

    According to the inversion result velocity model in Fig. 13(d),a low-speed zone is located approximately 35-45 m away from the tunnel face with a velocity of 3600 m·s-1. In the subsequent excavation process, a weak layer with a velocity of 3700 m·s-1was encountered approximately 30 m ahead of the tunnel face, which proves that the inversion result using the parameters proposed in this study is relatively effective and practical,and will also help to provide a guide for future tunnel construction engineering.

    4. Conclusions

    We applied a frequency-domain acoustic FWI to the tunnel seismic detection method and determined the influence of the frequency group selection strategy and tunnel observation system settings. The specific analysis is summarized as follows:

    (1) The observation system of the TSP is limited by the tunnel space,and its records have a high dominant frequency with a large frequency distribution range,which makes the common frequency group selection strategies for FWI no longer applicable. Thus, we proposed a strategy for selecting the frequency group under the tunnel detection condition,which uses a frequency group selection strategy combining the low-frequency selection strategy covering the vertical wave number and the high-frequency selection strategy of antialiasing.In comparison,this strategy can obtain a better resolution result ahead of the tunnel face.

    Fig.13. Tunnel seismic records and the velocity models.(a,b)X-component field records and processed field records of receiver R1,(c)initial model for field data inversion,and (d) inversion result models of the C frequency group selection strategy.

    (2) Because tunnel construction areas are significantly limited,the U-shaped observation system is widely used in tunnel seismic detection theoretical simulations and practical applications.The Ushaped observation system includes the tunnel face and both sidewalls; however, there is no recommended side length for the FWI of the tunnel seismic detection method.Therefore,in this study,by linearly increasing the side length of the tunnel observation system within limits, if the equipment allows, longer side length arrangements of the tunnel observation system will produce a higher resolution of the inversion result. However, when the side length of the tunnel observation system is more than 5 times the width of the tunnel face,the influence of increasing the observation system side length on the inversion results is negligible.

    (3) The results of frequency-domain acoustic FWI of the tunnel seismic detection based on the complex geological model proved that the proposed method can successfully reverse the stratum information ahead of the tunnel face under the parameters discussed and obtain high-resolution results of complex geological information. The practical results of the TSP data further verified the effectiveness and practicality of the parameter selection strategy.This study can provide ideas and references for future theoretical research and practical applications.

    Compliance with ethics guidelines

    Mingyu Yu,Fei Cheng,Jiangping Liu,Daicheng Peng,and Zhijian Tian declare that they have no conflict of interest or financial conflicts to disclose.

    Acknowledgments

    This work was supported by the National Natural Science Foundation of China (41704146) and the Fundamental Research Funds for National Universities,China University of Geosciences(Wuhan)(CUGL180816).We thank the seiscope research group for developing the shared FWI algorithm.

    日韩熟女老妇一区二区性免费视频| 欧美日韩一级在线毛片| 免费av中文字幕在线| 欧美人与性动交α欧美软件| 国产深夜福利视频在线观看| 看免费成人av毛片| 国产精品一国产av| 亚洲精品国产色婷婷电影| 夜夜骑夜夜射夜夜干| 999精品在线视频| 亚洲成人国产一区在线观看 | 亚洲一级一片aⅴ在线观看| 国产免费视频播放在线视频| 建设人人有责人人尽责人人享有的| 十八禁高潮呻吟视频| 国产精品.久久久| 亚洲色图 男人天堂 中文字幕| 久久ye,这里只有精品| 女性被躁到高潮视频| 国产亚洲精品第一综合不卡| 国产伦人伦偷精品视频| 七月丁香在线播放| 2021少妇久久久久久久久久久| 精品人妻一区二区三区麻豆| av.在线天堂| 性少妇av在线| netflix在线观看网站| 18禁观看日本| 最近中文字幕2019免费版| 91aial.com中文字幕在线观看| 精品一区在线观看国产| 啦啦啦啦在线视频资源| 19禁男女啪啪无遮挡网站| 亚洲色图 男人天堂 中文字幕| 两个人看的免费小视频| 蜜桃国产av成人99| 久久精品人人爽人人爽视色| 日韩不卡一区二区三区视频在线| 日韩欧美精品免费久久| 亚洲第一av免费看| 天天躁夜夜躁狠狠久久av| 亚洲精品国产一区二区精华液| 五月天丁香电影| 在线观看一区二区三区激情| 国产精品久久久久久人妻精品电影 | 中国国产av一级| 99热全是精品| 国产探花极品一区二区| 69精品国产乱码久久久| 热re99久久国产66热| 在线亚洲精品国产二区图片欧美| 亚洲av电影在线进入| 少妇人妻 视频| 啦啦啦视频在线资源免费观看| 国产成人欧美| 成人亚洲精品一区在线观看| 七月丁香在线播放| 午夜激情av网站| 精品一品国产午夜福利视频| 欧美日韩成人在线一区二区| 啦啦啦中文免费视频观看日本| 亚洲国产欧美日韩在线播放| 免费高清在线观看日韩| 精品午夜福利在线看| 老司机深夜福利视频在线观看 | 熟妇人妻不卡中文字幕| 国产精品秋霞免费鲁丝片| 久久久久精品国产欧美久久久 | 啦啦啦在线观看免费高清www| 婷婷色综合www| 久久久国产欧美日韩av| a级毛片在线看网站| √禁漫天堂资源中文www| 我要看黄色一级片免费的| 极品少妇高潮喷水抽搐| 国产成人精品久久二区二区91 | 久久精品国产a三级三级三级| 国产高清不卡午夜福利| 国产亚洲欧美精品永久| 国产精品久久久久久精品电影小说| 精品亚洲成国产av| 51午夜福利影视在线观看| 久久久久国产一级毛片高清牌| 亚洲四区av| 天天影视国产精品| a级毛片在线看网站| av天堂久久9| 青春草国产在线视频| 欧美黑人欧美精品刺激| 七月丁香在线播放| 亚洲一区中文字幕在线| 欧美黄色片欧美黄色片| 日本色播在线视频| 婷婷色综合www| 久久狼人影院| 国产精品一国产av| 国产在视频线精品| 亚洲欧美成人精品一区二区| 日本wwww免费看| 女人爽到高潮嗷嗷叫在线视频| 在线观看人妻少妇| 这个男人来自地球电影免费观看 | 99精国产麻豆久久婷婷| 亚洲欧洲日产国产| 女人被躁到高潮嗷嗷叫费观| 欧美日韩综合久久久久久| 日韩中文字幕视频在线看片| 麻豆乱淫一区二区| 成人亚洲精品一区在线观看| 男女边吃奶边做爰视频| 97人妻天天添夜夜摸| 亚洲,欧美精品.| 只有这里有精品99| 亚洲欧美一区二区三区久久| 午夜免费观看性视频| 亚洲精品自拍成人| 91精品国产国语对白视频| 精品人妻一区二区三区麻豆| 最新的欧美精品一区二区| 午夜免费鲁丝| 国产又爽黄色视频| av免费观看日本| 久久韩国三级中文字幕| videos熟女内射| 满18在线观看网站| 欧美精品av麻豆av| 亚洲国产成人一精品久久久| 黄频高清免费视频| 黑人欧美特级aaaaaa片| 精品人妻一区二区三区麻豆| 免费高清在线观看视频在线观看| 亚洲欧美一区二区三区久久| 蜜桃在线观看..| 制服诱惑二区| 老司机亚洲免费影院| 另类亚洲欧美激情| 日韩制服骚丝袜av| 免费在线观看视频国产中文字幕亚洲 | 国产亚洲av片在线观看秒播厂| 国产日韩欧美亚洲二区| 两个人免费观看高清视频| 亚洲人成网站在线观看播放| 亚洲精品美女久久久久99蜜臀 | 九草在线视频观看| 久久免费观看电影| 1024视频免费在线观看| 亚洲熟女精品中文字幕| 又粗又硬又长又爽又黄的视频| 日韩欧美精品免费久久| 尾随美女入室| 少妇的丰满在线观看| 97精品久久久久久久久久精品| 国产精品二区激情视频| 狂野欧美激情性xxxx| 亚洲五月色婷婷综合| 国产欧美亚洲国产| 国产亚洲午夜精品一区二区久久| 韩国精品一区二区三区| 人成视频在线观看免费观看| 高清不卡的av网站| 久热这里只有精品99| 国产 一区精品| 丰满少妇做爰视频| 午夜福利乱码中文字幕| 精品久久久久久电影网| 国产亚洲一区二区精品| 欧美日韩一区二区视频在线观看视频在线| 两个人看的免费小视频| 搡老岳熟女国产| a级毛片在线看网站| av国产精品久久久久影院| 午夜福利视频在线观看免费| 日本色播在线视频| 国产精品免费视频内射| 国产高清不卡午夜福利| 青青草视频在线视频观看| 伊人久久国产一区二区| 国产黄色免费在线视频| 又黄又粗又硬又大视频| 不卡视频在线观看欧美| 2018国产大陆天天弄谢| av线在线观看网站| 精品福利永久在线观看| 亚洲精品乱久久久久久| netflix在线观看网站| 亚洲熟女毛片儿| 午夜福利一区二区在线看| 午夜福利免费观看在线| 人体艺术视频欧美日本| 色网站视频免费| 欧美精品人与动牲交sv欧美| 久久99精品国语久久久| 亚洲精品国产一区二区精华液| 一边摸一边抽搐一进一出视频| 午夜福利视频在线观看免费| 日韩熟女老妇一区二区性免费视频| 午夜福利在线免费观看网站| 两性夫妻黄色片| 精品午夜福利在线看| 久久精品国产综合久久久| 精品少妇一区二区三区视频日本电影 | 久久精品国产综合久久久| 青春草国产在线视频| 国产精品国产av在线观看| 亚洲国产欧美一区二区综合| 日韩一卡2卡3卡4卡2021年| 制服人妻中文乱码| 精品久久久精品久久久| 亚洲伊人色综图| 看十八女毛片水多多多| 亚洲天堂av无毛| 午夜久久久在线观看| 在线 av 中文字幕| 人人妻,人人澡人人爽秒播 | 考比视频在线观看| 超色免费av| 亚洲,一卡二卡三卡| 九九爱精品视频在线观看| 国产熟女欧美一区二区| 亚洲一区中文字幕在线| 午夜精品国产一区二区电影| 国产成人精品在线电影| 亚洲av男天堂| 9热在线视频观看99| 亚洲男人天堂网一区| 99热国产这里只有精品6| 亚洲精品国产av成人精品| 久久精品aⅴ一区二区三区四区| 午夜日本视频在线| 黑人猛操日本美女一级片| 免费黄色在线免费观看| 欧美黑人欧美精品刺激| 亚洲av国产av综合av卡| 久久久久久久精品精品| 午夜精品国产一区二区电影| 亚洲一码二码三码区别大吗| 9色porny在线观看| 亚洲成人免费av在线播放| 天天操日日干夜夜撸| bbb黄色大片| 国产在视频线精品| 亚洲视频免费观看视频| 1024视频免费在线观看| 精品午夜福利在线看| 丰满少妇做爰视频| 亚洲熟女毛片儿| av视频免费观看在线观看| av.在线天堂| 免费高清在线观看视频在线观看| 国产激情久久老熟女| 亚洲精品国产av成人精品| 男女无遮挡免费网站观看| 成年美女黄网站色视频大全免费| 99热网站在线观看| netflix在线观看网站| 18禁观看日本| 欧美最新免费一区二区三区| 国产麻豆69| xxxhd国产人妻xxx| 久久久久久久大尺度免费视频| 人人妻人人爽人人添夜夜欢视频| 91aial.com中文字幕在线观看| 最近中文字幕高清免费大全6| 国产精品久久久久久精品古装| 精品免费久久久久久久清纯 | 欧美日韩国产mv在线观看视频| 亚洲成人国产一区在线观看 | 777久久人妻少妇嫩草av网站| 久久女婷五月综合色啪小说| av在线播放精品| 少妇被粗大的猛进出69影院| 国产精品人妻久久久影院| 久久久久久久大尺度免费视频| 亚洲精品视频女| 亚洲综合精品二区| 最近中文字幕2019免费版| 国产99久久九九免费精品| 赤兔流量卡办理| 中文字幕色久视频| 日本91视频免费播放| 欧美日韩福利视频一区二区| 麻豆精品久久久久久蜜桃| 久久99一区二区三区| 热re99久久国产66热| 国产精品二区激情视频| 亚洲情色 制服丝袜| 精品视频人人做人人爽| 大片电影免费在线观看免费| 午夜福利视频在线观看免费| av片东京热男人的天堂| 亚洲伊人色综图| 69精品国产乱码久久久| 国产一级毛片在线| 一区在线观看完整版| 一边摸一边做爽爽视频免费| 亚洲精品视频女| 国产精品嫩草影院av在线观看| 亚洲欧美成人精品一区二区| 日本午夜av视频| 亚洲视频免费观看视频| av福利片在线| av女优亚洲男人天堂| 人人妻,人人澡人人爽秒播 | 国产老妇伦熟女老妇高清| 午夜精品国产一区二区电影| 波野结衣二区三区在线| 国产国语露脸激情在线看| 啦啦啦在线免费观看视频4| 欧美成人精品欧美一级黄| 亚洲成av片中文字幕在线观看| 天天操日日干夜夜撸| 欧美精品av麻豆av| 欧美日韩亚洲高清精品| 丝袜美腿诱惑在线| 夫妻午夜视频| 黄色视频在线播放观看不卡| 一区福利在线观看| 亚洲一卡2卡3卡4卡5卡精品中文| 超碰成人久久| kizo精华| 十分钟在线观看高清视频www| 女的被弄到高潮叫床怎么办| 一区在线观看完整版| 韩国高清视频一区二区三区| 街头女战士在线观看网站| 亚洲五月色婷婷综合| 国产精品 国内视频| 中文乱码字字幕精品一区二区三区| 国产高清不卡午夜福利| 一区二区三区乱码不卡18| 久久女婷五月综合色啪小说| netflix在线观看网站| 看十八女毛片水多多多| 国产免费一区二区三区四区乱码| 国产精品一区二区在线观看99| 天美传媒精品一区二区| 亚洲图色成人| 国产爽快片一区二区三区| 日韩中文字幕欧美一区二区 | 美女视频免费永久观看网站| 国产精品.久久久| 国产精品人妻久久久影院| 91精品国产国语对白视频| 99久久综合免费| 视频区图区小说| 美女扒开内裤让男人捅视频| 亚洲精品成人av观看孕妇| 国产成人精品久久二区二区91 | 一边摸一边做爽爽视频免费| 叶爱在线成人免费视频播放| 久久久精品94久久精品| 精品酒店卫生间| 亚洲av男天堂| 纯流量卡能插随身wifi吗| 国产男女超爽视频在线观看| 日韩av免费高清视频| 一边摸一边做爽爽视频免费| 男女免费视频国产| 啦啦啦在线免费观看视频4| 国产精品国产av在线观看| 亚洲av男天堂| 日韩中文字幕欧美一区二区 | 精品酒店卫生间| xxx大片免费视频| 免费黄网站久久成人精品| 久久久久久久久久久久大奶| 女人爽到高潮嗷嗷叫在线视频| 国产一区二区激情短视频 | av福利片在线| 国产乱来视频区| 男的添女的下面高潮视频| 超碰成人久久| 夫妻性生交免费视频一级片| 丰满迷人的少妇在线观看| 新久久久久国产一级毛片| 人人妻人人澡人人看| 国产精品一区二区精品视频观看| 欧美日韩一区二区视频在线观看视频在线| 国产成人a∨麻豆精品| 亚洲精品美女久久av网站| 久久ye,这里只有精品| 久久精品熟女亚洲av麻豆精品| 亚洲精品乱久久久久久| 久久精品国产a三级三级三级| 一区在线观看完整版| 国产 精品1| 免费观看av网站的网址| 嫩草影院入口| 19禁男女啪啪无遮挡网站| 又粗又硬又长又爽又黄的视频| 亚洲精品中文字幕在线视频| 国产精品久久久久久人妻精品电影 | 国产人伦9x9x在线观看| 久久久久久久国产电影| 91aial.com中文字幕在线观看| 国产亚洲av高清不卡| 黄频高清免费视频| 国产免费视频播放在线视频| 美女中出高潮动态图| 99国产综合亚洲精品| 热re99久久国产66热| 捣出白浆h1v1| 综合色丁香网| 一本一本久久a久久精品综合妖精| av一本久久久久| 日韩视频在线欧美| 免费高清在线观看视频在线观看| 卡戴珊不雅视频在线播放| 秋霞伦理黄片| 日本vs欧美在线观看视频| 黄色 视频免费看| 麻豆av在线久日| xxx大片免费视频| 满18在线观看网站| 久久 成人 亚洲| 免费在线观看完整版高清| av片东京热男人的天堂| 午夜精品国产一区二区电影| 久久精品亚洲av国产电影网| 欧美日韩精品网址| 纵有疾风起免费观看全集完整版| 国产免费一区二区三区四区乱码| 在线精品无人区一区二区三| 丰满乱子伦码专区| 欧美精品高潮呻吟av久久| 少妇人妻久久综合中文| 亚洲精品中文字幕在线视频| 这个男人来自地球电影免费观看 | 男人添女人高潮全过程视频| 国产成人91sexporn| 中文字幕制服av| 午夜av观看不卡| 一级黄片播放器| 欧美最新免费一区二区三区| 啦啦啦在线免费观看视频4| 人人妻人人爽人人添夜夜欢视频| 国产精品欧美亚洲77777| 久久午夜综合久久蜜桃| 丝袜在线中文字幕| 中文欧美无线码| 99热全是精品| 亚洲国产日韩一区二区| 久久午夜综合久久蜜桃| 亚洲精品久久成人aⅴ小说| 最新的欧美精品一区二区| 国产精品99久久99久久久不卡 | 中文字幕精品免费在线观看视频| 午夜福利免费观看在线| 日日摸夜夜添夜夜爱| 高清视频免费观看一区二区| 男女高潮啪啪啪动态图| h视频一区二区三区| 一本色道久久久久久精品综合| 亚洲少妇的诱惑av| 操出白浆在线播放| 999精品在线视频| 免费人妻精品一区二区三区视频| 国产黄色视频一区二区在线观看| 免费观看人在逋| 久久精品久久精品一区二区三区| 亚洲第一青青草原| 欧美国产精品一级二级三级| 欧美在线黄色| 99热网站在线观看| 亚洲专区中文字幕在线 | 丝袜脚勾引网站| 国产精品av久久久久免费| 少妇猛男粗大的猛烈进出视频| 日本黄色日本黄色录像| 搡老岳熟女国产| 亚洲欧美精品自产自拍| 亚洲综合精品二区| 久久久国产精品麻豆| 一二三四中文在线观看免费高清| 亚洲,欧美,日韩| 青春草亚洲视频在线观看| 久久久欧美国产精品| av福利片在线| 如何舔出高潮| 欧美精品一区二区免费开放| 亚洲欧美激情在线| 日韩大码丰满熟妇| 欧美日韩视频高清一区二区三区二| 精品一品国产午夜福利视频| 欧美日韩亚洲国产一区二区在线观看 | 精品人妻熟女毛片av久久网站| 制服人妻中文乱码| 亚洲视频免费观看视频| 丰满乱子伦码专区| h视频一区二区三区| 午夜久久久在线观看| 国产精品香港三级国产av潘金莲 | 国产精品女同一区二区软件| a级片在线免费高清观看视频| 国产精品免费大片| 男人添女人高潮全过程视频| 在线观看三级黄色| 国产xxxxx性猛交| 国产精品人妻久久久影院| 午夜影院在线不卡| av免费观看日本| 晚上一个人看的免费电影| 亚洲色图 男人天堂 中文字幕| 999精品在线视频| 日本色播在线视频| 国产精品蜜桃在线观看| 99香蕉大伊视频| 蜜桃国产av成人99| 日本午夜av视频| 纯流量卡能插随身wifi吗| 99香蕉大伊视频| 18禁观看日本| 操出白浆在线播放| 日韩中文字幕欧美一区二区 | 欧美日韩一级在线毛片| 午夜精品国产一区二区电影| 99精品久久久久人妻精品| 国产日韩欧美视频二区| 亚洲精品久久久久久婷婷小说| 久久精品久久久久久久性| 如日韩欧美国产精品一区二区三区| 亚洲精品中文字幕在线视频| av免费观看日本| 在线亚洲精品国产二区图片欧美| 欧美日韩一级在线毛片| 狠狠精品人妻久久久久久综合| 在线观看免费午夜福利视频| 天天躁夜夜躁狠狠躁躁| 国产一区二区在线观看av| 极品少妇高潮喷水抽搐| 女人精品久久久久毛片| 极品少妇高潮喷水抽搐| 国产成人欧美在线观看 | 亚洲综合精品二区| 在线看a的网站| 丝袜脚勾引网站| 国产精品国产三级专区第一集| 精品国产超薄肉色丝袜足j| 两性夫妻黄色片| 国产精品蜜桃在线观看| 亚洲欧洲国产日韩| 久久久久久人人人人人| 日韩制服骚丝袜av| 国产精品麻豆人妻色哟哟久久| 大片电影免费在线观看免费| 秋霞伦理黄片| 成人影院久久| 国产探花极品一区二区| 卡戴珊不雅视频在线播放| 国产av精品麻豆| 无限看片的www在线观看| 狠狠婷婷综合久久久久久88av| 亚洲精品av麻豆狂野| 99国产综合亚洲精品| 欧美国产精品va在线观看不卡| 午夜福利网站1000一区二区三区| 色婷婷久久久亚洲欧美| 日韩人妻精品一区2区三区| 麻豆精品久久久久久蜜桃| 国产一区二区三区av在线| 中文字幕亚洲精品专区| 亚洲伊人久久精品综合| 黑人巨大精品欧美一区二区蜜桃| 色吧在线观看| 午夜福利乱码中文字幕| 在线观看国产h片| 黄色视频不卡| 女人爽到高潮嗷嗷叫在线视频| 高清黄色对白视频在线免费看| 少妇 在线观看| 超色免费av| 伊人久久大香线蕉亚洲五| 日日撸夜夜添| 精品一品国产午夜福利视频| 精品第一国产精品| 大香蕉久久网| 久久精品国产综合久久久| 国产女主播在线喷水免费视频网站| 欧美成人午夜精品| 国产熟女午夜一区二区三区| 国产乱人偷精品视频| 国产男人的电影天堂91| 久久天堂一区二区三区四区| 国产成人av激情在线播放| 亚洲人成77777在线视频| 2021少妇久久久久久久久久久| 亚洲国产欧美一区二区综合| 国产精品久久久久久精品古装| 国产精品秋霞免费鲁丝片| 欧美97在线视频| 伦理电影免费视频| 亚洲精品美女久久av网站| 99九九在线精品视频| 国产精品一区二区在线观看99| 亚洲综合精品二区| 男女边摸边吃奶| 丝袜在线中文字幕| 日韩大片免费观看网站| 午夜免费观看性视频| 国产精品女同一区二区软件| 波野结衣二区三区在线| a 毛片基地| av免费观看日本| 一个人免费看片子| 搡老乐熟女国产| 亚洲成色77777| 下体分泌物呈黄色| 国产一区二区三区av在线| 伊人久久国产一区二区| 亚洲欧洲日产国产| 国产欧美亚洲国产| 精品久久久精品久久久| 天天躁夜夜躁狠狠久久av| 在线亚洲精品国产二区图片欧美| 亚洲一级一片aⅴ在线观看|