Ground penetrating radar(GPR)method,with advantages of high resolution,high efficiency,continuity and nondestructive,has become one of the most important geophysical exploration methods in the field ofshallow structure exploration,especially playing an important role in quality nondestructive testing in shallowengineering such as highway engineering,construction engineering(Saarenketo and Scullion,2000; Maierhofer,2003; Alani et al,2013; Huang and Zhang,2013). China’s highway structural layer usually adopts asphaltconcrete surface plus semi-rigid base. When GPR is applied in thickness detection and disease(such as cracks,oid,differential settlement,looseness)investigation of highway structural layer,the material of highway structural layer is usually simplified to homogeneous medium or layered homogeneous medium(Lu et al.,2007).But actual highway pavement material(such as asphalt concrete,cement stabilized macadam)is multiphaseheterogeneous mixture of a variety of substances in accordance with certain volume fraction(Cai,2008). Takeasphalt concrete as an example,it is composed of aggregate,asphalt mortar and air with different volume fraction. Aggregate,asphalt mortar and air usually have different sizes,shapes,dielectric properties,and distributer and omly in space(Wu,2009; Ding et al.,2012),and with obvious multiphase,discrete,r and om medium properties. Porous medium widespread in nature is usually multiphase discrete r and om medium(Cai and Yu,2011).When high-frequency electromagnetic wave propagates in multiphase discrete r and om medium,scattering willoccur,resulting in a large number of scattering wave arrivals which lead to a corresponding r and om feature inreceiving waveforms(Jiang et al,2013; Dai and Wang,2013). In particular,road GPR antenna is with relativelyhigh dominant frequency,relatively short detection wavelength(about 3 to 5 times the maximum particle size),and relatively high scattering,which amplifies the wave field changes caused by r and om distribution of dielectric properties in medium space and will significantly affect the propagation of high-frequency electromagneticwaves. Thus,for high-frequency electromagnetic waves,the physical parameters of asphalt concrete changedramatically,the electromagnetic scattering is strong,it is unable to be simplified to homogeneous medium,therefore,a multiphase discrete r and om medium model more in line with the actual situation should be established for numerical simulation of GPR wave propagation,to provide a theoretical basis for improving GPRresolution and quantitative interpretation of medium properties of engineering materials,which has importantscientific significance.

The authors adopted core drilling method to obtain asphalt concrete core samples,applied vector networkanalyzer to measure permittivity of the main composition,obtained spatial distribution image of permittivity,studied r and om statistical distribution characteristics of permittivity in the core sample space; then,bycombining volume fraction of all compositions of asphalt concrete,proposed modeling algorithm of multiphasediscrete r and om medium,constructed multiphase discrete r and om medium model more consistent with theactual situation of asphalt concrete,and conducted forward simulation and comparative analysis of GPR wavefield characteristics in homogeneous medium model,continuous r and om medium model and multiphase discreter and om medium model.

Theory multiphase discrete r and om medium can be described with r and om medium model theory. Inr and om medium model theory,the physical parameter of the medium is treated as a spatial r and om variable and its distribution characteristic is described based on statistical method. R and om medium model can be dividedinto stationary and non-stationary types. Stationary r and om medium model is characterized by statistics ofautocorrelation function autocorrelation length,autocorrelation angle,mean,st and ard deviation(Xu et al.,2007). Gaussian(Ergintav and Canitez,1997; Zhao,etc.,2013),exponential(Ikelle et al,1993.) and mixed(Xi and Yao,2002; Jiang et al,2013.)ellipsoidal autocorrelation functions are widely used to construct r and ommedium model. In the r and om medium,autocorrelation length describes the average scale of inhomogeneousanomaly in different directions,autocorrelation angle describes the trend of inhomogeneous anomaly,me and escribes average characterization of r and om medium,st and ard deviation describes dispersion degree of r and ommedium.

Stationary r and om medium model assumes that the statistical characteristic of medium is a constant,which has many limitations and cannot fully describe complex,non-stationary actual medium. Xi and Yao(2005)proposed a new method of constructing non-stationary r and om medium model via spatial variable autocorrelation function(local autocorrelation function). This method constructs various forms of non-stationaryr and om medium models by selecting local autocorrelation length in horizontal and vertical directions. Studieshave shown that: non-stationary r and om medium model can more flexibly describe various complex,heteroge-neous actual media. Expressions of Gaussian,exponential and mixed two-dimensional non-stationary ellipsoidalautocorrelation functions are as follows: Gaussian ellipsoidal autocorrelation function:

Exponential ellipsoidal autocorrelation function:

Mixed ellipsoidal autocorrelation function:

The above two-dimensional non-stationary ellipsoidal autocorrelation functions describe slow changeof autocorrelation functions R(x',y',x₁,y₁)on a largescale with the spatial coordinates(x',y'). In twodimensional non-stationary r and om medium,the relationship between global coordinates and local coordinates are shown in Fig. 1.

(Fig.1)

Fig.1 Relationship between global coordinate and local coordinate in two-dimensional non-stationary random medium

In the non-stationary r and om medium model,local autocorrelation length a = a(x',y'),b = b(x',y'),and direction angle θ = θ(x',y')change slowly withspatial coordinates(x',y'). a = a(x',y'),b = b(x',y') and θ = θ(x',y')respectively denotes autocorrelationlength and autocorrelation angle of r and om mediumat point(x',y')in x,y directions. Compare the aboveformulas,and it can be seen that mixed ellipsoidalautocorrelation function not only incorporates features of exponential(r = 1) and Gaussian(r = 0)ellipsoidalautocorrelation functions,but also can further exp and the type and scope to simulate actual r and om mediumthrough roughness factor r(Guo et al.,2012).

As medium physical parameters described in stationary and non-stationary r and om medium model havespatial continuous r and om variation,the model is also called continuous r and om medium model. Not only thephysical parameters of multiphase discrete r and om medium show multiphase,discrete and r and om distributioncharacteristic in space,but also the compositions meet certain volume fraction. Therefore,the above continuousr and om medium model cannot fully describe the multiphase discrete r and om medium. The modeling algorithmof r and om medium model based on spatial r and om distribution statistical characteristics and volume fraction ofeach composition of asphalt concrete should be improved to build a multiphase discrete r and om medium modelwhich is more in line with the actual situation of asphalt concrete.

Permittivity is a basic parameter to analyze the interaction of electromagnetic wave and medium. Thestatistical characteristic of spatial r and om distribution of permittivity in asphalt concrete is the basis for reconstruction of asphalt concrete medium model. The authors drilled middle surface course core sample of asphalt concrete of a highway in western Henan Province,transversely cut it and obtained a two-dimensional crosssection of 7 cm×7 cm as shown in Fig. 2a. With E5071C vector network analyzer(Wu et al.,2011),the authors measured relative permittivity of aggregate,asphalt at frequency b and of 0.5~2 GHz several times,thenadopted the mean as relative permittivity of aggregate and asphalt,respectively ε_r = 8.2,ε_r = 2.8. As there isrelatively few asphalt mortar in asphalt concrete,the distribution is scattered,it is unable to directly measureits permittivity. Therefore,according to volume fraction and permittivity of asphalt and mineral powder inasphalt mortar,this paper adopts complex refractive index model(Birchak et al.,1974)which estimates thatthe effective permittivity ε_r = 5.44.

(Fig.2)

Fig.2 Transverse profile of asphalt concrete (a) and its permittivity in spatial distribution image (b)

Based on different gray levels of aggregate,asphalt mortar and air in the asphalt concrete profile,separateaggregate,asphalt mortar and air with image threshold segmentation technique,and assign correspondingpermittivity to obtain spatial distribution image ε_r(x,y)of asphalt concrete permittivity,as shown in Fig. 2bbelow.

Based on spatial distribution image ε_r(x,y)of permittivity,volume fraction of composition can be calculated(for instance,aggregate ratio P_g=77.86%,asphalt mortar ratio P_j=21.70%,porosity P_k=0.43%). Also,statistical characteristics of spatial r and om distribution such as mean of relative permittivity ε_mv=7.5694,st and ard deviation σ=1.2175 can be obtained.

Autocorrelation function is an important parameter to reflect the spatial distribution characteristic ofr and om medium and determine the type of r and om medium( Klimeš,2002; Liu et al.,2007). According toWinner-Khintchine theorem,it can be known that the power spectrum of permittivity spatial distribution image and its autocorrelation function are a Fourier transform pair. Therefore,the power spectrum of permittivityspatial distribution image can be calculated first,then conduct two dimensional inverse Fourier transformof it to obtain its autocorrelation function image,then conduct normalization and binarization and extractcharacteristic parameters such as autocorrelation length,autocorrelation angle of autocorrelation function.

(1)Power spectrum of spatial distribution image ε_r(x,y)of permittivity

Power spectrum of spatial distribution image ε_r(x,y)of permittivity can be calculated with the followingformula(Balboa and Grzywacz,2003; Amin and Subbalakshmi,2007).

In the formula(4)

where k_x,k_y are spatial frequency of spatial distribution image of permittivity in x,y directions respectively,M,N are sampling points of spatial distribution image of permittivity in x,y,directions respectively.

(2)Autocorrelation function R_εε(x,y) corresponding to ε_r(x,y)

According to Winner-Khintchine theorem,autocorrelation function R_εε(x,y) of spatial distribution imageε_r(x,y) of permittivity is the two-dimensional inverse Fourier transform of power spectrum Γ (k_x,k_y). That is

Autocorrelation function corresponding to spatial distribution image ε_r(x,y)of permittivity is shown in Fig. 3.

(Fig.3)

Fig.3 Auto correlation function image

(3)Normalized autocorrelation function

To st and ardize the degree of correlation of permittivity ε_r(x,y)at different spatial locations,conductlinear normalization of autocorrelation function R_εε(x,y),so that R_εε(x,y) is in the range of 0 to 1,i.e.0 ≤ |R_εε(x,y)| ≤ 1. The autocorrelation function image after linear normalization is shown in Fig. 4.

(Fig.4)

Fig.4 Auto correlation function image after normalization

(4)Extract characteristic parameters of autocorrelation function

For linearly normalized autocorrelation function,its autocorrelation length is the distance between 1 and e^-1(Zhang and Sundararajan,2006). Therefore,to facilitate the extraction of autocorrelation length a,b ofautocorrelation function in x,y directions,autocorrelation angle θ and other characteristic parameters,assignthe autocorrelation function value greater than or equal to e^-1 as 1,while that smaller than e^-1 as zero(Fig. 5).

(Fig.5)

Fig.5 Auto correlation function image after binary-conversion (figure b is the area bounded by dotted line in figure a)

Based on the above statistics and estimation method,multiphase discrete r and om medium model parameters such as spatial r and om distribution statistical characteristics of permittivity spatial distribution image(Fig. 2b),autocorrelation function characteristic parameters,and volume fraction of composition can be obtained,as shown in Table 1.

Table 1()

Table 1 Model parameters of multiphase discrete random medium

Table 1 Model parameters of multiphase discrete random medium

The shape of the autocorrelation function in Fig. 5 is similar to ellipse. We have calculated autocorrelationfunctions of a lot of asphalt concrete sections,and all their shapes are similar to ellipse. Therefore,we believethat autocorrelation function of asphalt concrete approximates ellipsoidal autocorrelation function. Ellipsoidalautocorrelation functions used to construct r and om medium model mainly include Gaussian,exponential and mixed type. Considering that mixed ellipsoidal autocorrelation function can further exp and type and scopeto simulate actual r and om medium through roughness factor r,this paper builds multiphase discrete r and ommedium model based on mixed ellipsoidal autocorrelation function.

Based on the r and om distribution statistical characteristic quantity of permittivity in the space of asphaltconcrete and the autocorrelation function characteristic parameter of permittivity,and by combining the volumefraction of components of asphalt concrete,we improve the modeling algorithm of continuous r and om medium,establish the modeling algorithm of multiphase discrete r and om medium which can accurately describe thespatial r and om distribution statistical characteristic of electromagnetic parameters of engineering materials and the quantitative constraints of volume fraction of composition. Specific modeling method of multiphase discreter and om medium model is as follows:

(1)According to the estimation results of autocorrelation function characteristic parameters such as autocorrelation length a,b,and autocorrelation angle θ in permittivity spatial distribution,select an appropriateroughness factor r,conduct two-dimensional Fourier transform of mixed ellipsoidal autocorrelation functionR_M(x,y),and obtain power spectrum R_M(k_x,k_y) of r and om process.

(2)Generate r and om information. Generate uniformly distributed independent two-dimensional r and omfield Φ(k_x,k_y) in interval [0,2π] with r and om number generator.

(3)Generate r and om power spectrum function. According to the theory of r and om process,the productof r and om process power spectrum R_M(k_x,k_y) and two-dimensional r and om field Φ(k_x,k_y) is r and om powerspectrum function.

(4)Conduct two-dimensional inverse Fourier transform of r and om power spectrum function and obtainr and om disturbance.

(5)St and ardize mean and st and ard deviation of r and om disturbance,obtain continuous two-dimensionalr and om medium model with R_M(k_x,k_y) as autocorrelation function and with specified mean and st and ard deviation.

(6)Determine local size R of parameters of multiphase discrete r and om medium model. Suppose aggregatevolume fraction of asphalt concrete is P_g and air volume fraction is P_k,then the volume fraction of asphaltmortar is P_j = 1 - P_g - P_k,the relative permittivity of aggregate,asphalt mortar and air are ε₁=8.2,ε₂=5.44,ε₃=1. Take R × R as local region segmentation model,enter the local region.

(7)According to the model values,corrode grid points in the local region from big to small in turn. Assignε1 to value of grid point,and judge whether the local region meets aggregate volume fraction P_g. If so,go toStep(8). Otherwise,repeat Step(7).

(8)According to the model values,corrode grid points in the local region from small to big in turn. Assignε₃ to value of grid point,and judge whether the local region meets aggregate volume fraction P_k. If not,repeatStep(8). Otherwise,go to the next local region,and perform(7)–(8)until all local regions meet aggregate ratioP_g and porosity P_k.

(9)After circulation of the entire model,assign ε₂ to all grid points with model values other than ε1 orε₃. Thus,a multiphase discrete r and om medium model with self-organized structur,with aggregate volumefraction P_g,air volume fraction P_k and asphalt mortar volume fraction P_j.

Take spatial distribution statistical characteristic quantity(mean,st and ard deviation)of asphalt concretepermittivity and autocorrelation function characteristic parameter(correlation length autocorrelation angle)inTable 1 as model parameters,based on mixed ellipsoidal autocorrelation function with roughness factor r=2.5 and according to Steps(1)to(5)in the above modeling method,establish continuous r and om medium model,asshown in Fig. 6. On the basis of Fig. 6,combine volume fraction of each composition,and according to Steps(6)–(9)in the above modeling method,establish multiphase discrete r and om medium model which describes spatialr and om distribution statistical characteristic of asphalt concrete permittivity and also conforms to volumefraction quantitative constraint,as shown in Fig. 7. Fig. 7 better reflects organizational structure characteristicsof asphalt concrete,that is,aggregate plays a role of skeleton while asphalt mortar plays a role of filling.

(Fig.6)

Fig.6 Continuous random medium model

(Fig.7)

Fig.7 Multiphase discrete random medium model

Take spatial distribution statistical characteristic quantity of permittivity,fraction of all composition and autocorrelation function characteristic parameter in Table 1 as model parameters,change the roughness factorr in mixed ellipsoidal autocorrelation function in order,establish multiphase discrete r and om medium model ofroughness factor r = 0,0.1,0.2,· · ·,5. Part of the multiphase discrete r and om medium model is shown in Fig. 8.

(Fig.8)

(a) Roughness factor r = 1; (b) Roughness factor r = 2; (c) Roughness factor r = 3; (d) Roughness factor r = 4. Fig.8 Impact of roughness factor on multiphase discrete random medium models

As can be seen from Fig. 8,when other model parameters are constant,the smaller roughness factor ofmixed ellipsoidal autocorrelation function is,more dispersed the generated model value of multiphase discreter and om medium is; the bigger roughness factor is,the better aggregation effect of model value of multiphasediscrete r and om medium is. As roughness factor becomes bigger and bigger,aggregation effect of its generatedmultiphase discrete r and om medium model becomes better,which means larger aggregate particle size. In orderto compare the degree of similarity between the above multiphase discrete r and om medium models and spatialdistribution image(Fig. 2b)of asphalt concrete permittivity,this paper estimated autocorrelation function and characteristic parameter of the above multiphase discrete r and om medium models,and calculated correlationcoefficient between them and autocorrelation function of permittivity spatial distribution image. Change curveof correlation coefficient of roughness factor r is shown in Fig. 9.

(Fig.9)

Fig.9 Change curve of correlation coefficient with roughness factors

As can be seen from Fig. 9,when roughness factor r=2.5,the correlation coefficient between autocorrelation function image of multiphase discrete r and om medium model and autocorrelation function image of asphalt concrete reaches the maximum,whichmeans the best correlation and the highest similarity.Therefore,this paper takes mixed ellipsoidal autocorrelation function with roughness factor r=2.5 as autocorrelation function expression of asphalt concretepermittivity.

According to statistical characteristic quantity of permittivity in asphalt concrete space,autocorrelationfunction characteristic parameter,estimation result of roughness factor of mixed ellipsoidal autocorrelationfunction,and based on modeling algorithm of multiphase discrete r and om medium,multiphase discrete r and ommedium with any volume fraction can be built. Next,take multiphase discrete r and om medium model withdifferent porosity as an example.

Porosity means the volume fraction of air in asphalt concrete,which is the most important factor affectingservice life and deformation property of pavement. For asphalt concrete pavement,the asphalt-aggregate ratiois relatively stable(Zhong et al.,2007). When the porosity increases,volume fraction of aggregate and asphaltmortar decreases in proportion. Take spatial distribution statistical characteristic quantity(mean,st and ard deviation)of asphalt concrete permittivity and autocorrelation function characteristic parameter(correlation length autocorrelation angle)in Table 1 asmodel parameters,based on mixed ellipsoidal autocorrelation function(roughness factor r=2.5),constructmultiphase discrete r and om medium model(Fig. 10)with porosity of 2%,4%,6%,8% respectively. Table 2 shows volume fraction of all compositions in multiphase discrete r and om medium with different porosity.

(Fig.10)

(a) Porosity=2%; (b) Porosity=4%; (c) Porosity=6%; (d) Porosity=8%. Fig.10 Multiphase discrete random medium models with different porosity

Table 2()

Table 2 Volume fractions of composition of multiphase discrete random medium models with different porosity

Table 2 Volume fractions of composition of multiphase discrete random medium models with different porosity

In order to study the propagation characteristics of GPR wave in multiphase discrete r and om medium,this paper conducted GPR forward simulation of multiphase discrete r and om medium model(porosity=4%) and homogeneous medium model based on FDTD method. The total number of model grid is 500×600,steplength of spatial grid is 0.1 cm,the upper part of the model is air medium with thickness of 10 cm; the centralarea of the model is multiphase discrete r and om medium or homogeneous medium,the thickness is 40 cm,thereis a vertical crack of 0.5 cm wide in the middle,the top and bottom of the crack is 15 cm from the interface ofupper and lower medium; the lower area of the model is homogeneous lossless medium(ε_r = 9)with thickness of 10 cm. Center frequency of GPR pulse excitation source is 1500 MHz. Fig. 11 is wave field snapshot of GPRwave in multiphase discrete r and om medium and homogeneous medium at the moment of 3 ns. As can beseen,in multiphase discrete r and om medium,GPR wave shows strong scattering,r and omly scattered waves aredisordered and mutually superimpose,which makes reflected waves at the top and bottom of the crack distort and discontinuous,difficult to be identified; while in homogeneous medium,GPR wave is not scattered withrelatively regular waveform and reflected wave at the top and bottom of the crack is clearly visible.

(Fig.11)

Fig.11 Wave field snapshots of multiphase discrete random medium model (a) and homogeneous medium model (b)

In order to facilitate comparative analysis,this paper put forward simulation results of multiphase discreter and om medium and homogeneous medium in the same cross-sectional profile as shown in Fig. 12. 0~50 cm inthe figure is for forward simulation profile of multiphase discrete r and om medium. The profile is with obviousr and om disturbance and “noise”,the hyperbolic reflected wave from the top and bottom of the crack is difficultto identify,the event of lower interface reflected wave is distorted,not continuous and with weak amplitude.This stems from GPR wave scattering in multiphase discrete r and om medium. Some scattered waves return tothe ground and are received by receiving antenna,which reduces signal to noise ratio and resolution of GPRecho. 51～100 cm in the figure is for forward simulation profile of homogeneous medium whose diffracted wave hyperbolic curve from the crack top and bottom isclearly visible,event of lower interface reflected wave isdistinct,straight and with relatively strong amplitude.

(Fig.12)

Fig.12 Forward simulation profiles of multiphase discrete random medium model and homogeneous medium model

Compared with continuous r and om mediummodel,multiphase discrete r and om medium model notonly describes multiphase,discrete and r and om distribution statistical characteristics(such as autocorrelation function,mean,st and ard deviation,etc.)ofasphalt concrete in space,but also further describesvolume fractions of all compositions,thus able to morecomprehensively and accurately describe the complexmultiphase discrete r and om medium of asphalt concrete. Take Table 2 as modeling parameters,constructcontinuous r and om medium model and correspondingmultiphase discrete r and om medium model,simulate and analyze propagation characteristics of GPR wavein continuous r and om medium and multiphase discrete r and om medium. Forward simulation parameter isidentical with those in 5.1,so only the central area of the model needs to be replaced with continuous r and ommedium model or multiphase discrete r and om medium model. Fig. 13 is wave field snapshot of GPR wave incontinuous discrete r and om medium and multiphase discrete r and om medium at the moment of 3 ns. Fig. 14 isforward simulation profile of continuous r and om medium and multiphase discrete r and om medium. 0~50 cmin the figure is for forward simulation profile of continuous r and om medium,while 51~100 cm in the figure isfor forward simulation profile of multiphase discrete r and om medium.

(Fig.13)

Fig.13 Wave field snapshots of continuous random medium model (a) and multiphase discrete random medium model (b)

(Fig.14)

Fig.14 Forward simulation profiles of continuous random medium model and multiphase discrete random medium model

As can be seen from Figs. 13,14:(1)GPR wave scattering occurs in both continuous r and om medium and multiphase discrete r and om medium,as scattered wave of r and om propagation mutually superimpose,thereoccurs r and om disturbance and “noise”,causing distortion of reflected wave from top and bottom end of the crack,and making the reflection event from mediuminterface offset and discontinuous;(2)Comparedwith multiphase discrete r and om medium,GPRwave scattering is stronger in continuous r and ommedium,with greater r and om disturbance,stronger“noise”,greater distortion of reflected wave ofanomalous body.

In GPR measured data interpretation of asphalt concrete,r and om disturbance caused by heterogeneity of multiphase discrete r and om medium isusually treated as “noise”,which virtually gives upa lot of potentially valuable information. Research and use of this information is believed to be able toimprove resolution capability of GPR so as to distinguish property information of materials and morehidden disease information.

This paper regards asphalt concrete as multiphase discrete r and om medium. Through core sample,westudied spatial r and om distribution statistical characteristic of permittivity,estimated its autocorrelation function and characteristic parameters,proposed multiphase discrete r and om medium modeling algorithm of quantitative constraint,constructed multiphase discrete r and om medium model with different porosity,conductedforward intermediate and comparative analysis of wave field characteristic of GPR wave in homogenous model,continuous r and om medium model and multiphase discrete r and om medium model,and reached the followingconclusions:

(1)Asphalt concrete permittivity shows multiphase,discrete,r and om distribution characteristic in space,a typical multiphase discrete r and om medium. Its autocorrelation function approximates ellipsoidal autocorrelation function;

(2)This paper presents a modeling algorithm for multiphase discrete r and om medium,the modeling resultshows that the multiphase discrete r and om medium model not only describes the statistical characteristic ofasphalt concrete permittivity in space,but also coincides with volume fraction of compositions. It can morecomprehensively and accurately describe the complex heterogeneous mixture medium of asphalt concrete. Also,the model provides new ways and means for describing other similar materials or media;

(3)Strong scattering occurs in propagation of GPR wave in multiphase discrete r and om medium,r and omscattered waves in disordered propagation interfere and superimpose with each other,as a result,the receivedwaveform shows corresponding r and om disturbance characteristic and “noise”. Studying the relations betweenr and om disturbance characteristic of GPR echo and multiphase discrete r and om medium model parameter willprovide reference and help in quantitative evaluation of attribute parameter of multiphase discrete r and ommedium.

This work was supported by the High School Key Scientific Research Project of Henan province(15A170005) and the Doctoral Fund of Henan Institute of Engineering(D2014008). In addition,I extend my heartfelt thanksto the two reviewers for their careful review and constructive modification suggestions for this paper.