- Email: [email protected]
A simple method for estimating space-dependent parameters in soil science
Proceedings of the 18th World Congress The International Federation of Automatic Control Milano (Italy) August 28 - September 2, 2011
A simple method for estimating space-dependent parameters in soil science Philippe Neveux Universit´e d’Avignon et des Pays de Vaucluse UMR UAPV-INRA 1114 Environnement M´editerran´een et Mod´elisation des Agro-Hydrosyst`emes (EMMAH) Facult´e des sciences 84029 Avignon cedex 1 - France email: [email protected] Abstract: The problem of parameter estimation for a distributed parameter system (DPS) is treated. The DPS under consideration appears in Soil Science where it models a measurement method of the soil water content (namely Time Domain Reflectometry (TDR)). The model consists in a set of partial differential equation (PDE) known as the telegraphist equation. The measurement method provides the current at only one position. This increases the complexity of the parameter estimation problem. The solution proposed takes advantage of a property of the measured signal. Compared to existing solutions, the proposed technique has a very reduced computational burden. Examples show its ability to estimate space-dependent parameters with a reasonable confidence interval and a reduced computational burden compared to existing solutions. Keywords: Distributed Parameter System, Parameter identification, Telegraphist Equation, Soil Science. 1. INTRODUCTION Soil science is an area that develops complex models. They can consist of sets of Partial Derivative Equations (PDE) and/or Fractional Derivative Equations. The objective is to represent phenomenon precisely enough in order to understand the behavior of soils under various conditions. Consequently, these models require a parametric estimation in order to be tractable for study and prediction. Depending on the model structure and the available measurements, this task can be very complex. Among Soil Science areas, Time Domain Reflectometry (TDR) is known as a valuable indirect technique for the measurement of soil water content. Indirect because TDR permits to obtain an apparent relative permittivity of the soil εr . Once the apparent εr is known, the determination of the apparent water content is done using the relation between the soil apparent relative permittivity εr and the soil water content (Topp [1980]). Soil water content permits to get information on the soil status. This information is required to model various phenomenon. Hence, an accurate knowledge of the water content is necessary.
equation). The parameter estimation of this model is a very complex task. As a matter of fact, the estimation is based on the fitting of the model and the measurement made at only one position. Furthermore, the discretization of this distributed parameter system results in a model with a large number of interpolation nodes that entails a heavy computational burden. The parameters to be estimated are directly related to the value of εr at the interpolation nodes. Hence, the number of unknown parameters is large as one needs at least 500 interpolation nodes in the Finite Element Method to simulate the TDR intensity record. In order to reduce the size of the problem, we can consider two cases for the representation of the profile of εr in the soil:
Though efficient, this technique only brings a general information on the soil water content. In order to gain more detailed information, the community of Soil Science has developed several algorithms in order to estimate the soil water content at various locations in the soil from the TDR measurements (Greco [2006], Heimovaara et al. [2004], Oswald et al. [2003], Schlaeger [2005]).
(1) The permittivity is considered as a stair function (i.e. the soil is considered has a stratified medium) (Oswald et al. [2003]). In Oswald et al. [2003], the question is the optimal choice for the number of stair steps, their localization and of course, the magnitude of the stair steps.This problem is optimized by means of Genetic Algorithm (GA) search algorithm. (2) The permittivity profile can be regarded as a polynomial function (i.e. the soil is considered as a continuum) if the profile is known to be monotone (Greco [2006]). If such an information is not available, a series of interpolation points together with a linear interpolation between consecutive points, permits to approximate the permittivity profile (Greco and Guida [2008]). In both cases, a GA algorithm is used to optimize the estimate according to a fitness function.
In the time domain analysis, the TDR line is well described by a set of Partial Derivative Equations (the telegraphist 978-3-902661-93-7/11/$20.00 © 2011 IFAC
2925
10.3182/20110828-6-IT-1002.02746
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
Though interesting, these approaches lead to time consuming techniques due to extensive computations. In comparison, the measurement process takes approximately 100 nano-seconds while the optimization by means of GA requires several days. In order to obtain a faster estimate of the permittivity, we propose a simple method that aims at taking advantage of all the information present in the recorded signal.
the measured current at z = 0. Generally, parameter estimation or inverse problems for systems described by PDE require the knowledge of the state variables at various positions z (Banks and Kunisch [1989], Coca and Billings [2000], Lundstedt and He [1996], Seinfeld [1969]). Unfortunately, in TDR, this is not possible. Hence, the problem is not easy to solve. This explains why search algorithms are usually employed.
2. STATEMENT OF THE PROBLEM `a TDR technique consists in at least two rods that are ”plugged” in the soil (see figure 1). At the inlet of the rods (i.e. z = 0), we impose a step voltage u(t) and we measure the variation of the current at the inlet of the line (i.e. i(z = 0, t)). The fluctuations in the magnitude of the current are essentially due to the fluctuations of the permittivity εr along the rods. The theory of transmission line is applied to TDR (Magnusson et al. [1992]). As a matter of fact, the TDR probe can be considered as a transmission line with an inlet circuit (the impedance of the cable tester device) and outlet circuit (representing the loss of EM energy at the end of the line). The dimensionless model (the telegraphist equations) corresponding to this wave propagation phenomenon is given below (Bolvin et al. [2005]):
where:
u(t) -
`a
0 v(z, t)
Soil
?i(z, t)
εr (z) 1 z ?
RT Fig. 1. Schematic view of the TDR line. 3. MAIN RESULT
µr ∂t i = −∂z v − ai εr ∂t v = −∂z i − bv
- i(z, t) and v(z, t) are the intensity and the voltage at time t and position z ∈ [0, 1], - εr (z) is the relative permittivity of the soil, - µr ≃ 1 is the relative magnetic permeability of the soil, - a is related to the resistivity of the rods and b to the permittivity of the soil through the relation (Noborio [2003]): √ b ≡ b0 εr In the following, the contribution of the parameter a will be neglected as it is proposed in (Schlaeger [2005]), - u(t) is a deterministic input signal. It is a step function with specified raise time, - the boundary conditions: v(z = 0, t) = u(t) v(z = 1, t) RT - RT is without loss of generality representing the resistivity of the soil at the end of the line, - ∂x f is the partial derivative of function f with respect to the variable x. i(z = 1, t) =
The design of the estimation technique will be done in three steps. Firstly, an interesting property will be explored when the loss parameter b0 is set to zero. It will permit to get the basis of the estimation technique. Then, the effect of parameter b0 will be explored and the results obtained will allow us to achieve the design of the estimation technique. 3.1 First Step: the lossless case To understand the behavior of a model, one can develop a methodology based on sensitivity analysis. But in our situation, the idea is to take advantage of a property of this model when the parameter b0 is null, i.e., the lossless case. As a matter of fact, it is known that for an homogeneous medium (i.e. εr = constant), the magnitude of the signal √ i(z = 0, t) equates εr . Consequently, for inhomogeneous media, we have to verify whether this property still holds. A permittivity profile has been generated and the corresponding TDR signal has been simulated (see figure 2). In order to verify the property, the following procedure has been applied:
The objective of this work is to find a strategy that permits to obtain a fast estimate of parameter b0 and permittivity εr (z) from a single measurement i(z = 0, t) with a reduced computational cost.
Step 1- select the magnitude of the wavefront → is (t), Step 2- square the selected signal → i2s (t) ≡ εˆ0r (t), Step 3- project εˆ0r (t) onto z-space according to the timepropagation constraint: Zz p t(z) = 2 εˆ0r (ζ) dζ
The difficulty of the problem lies in the fact that the only available information for parameter estimation is
to get the estimated profile εˆ0r (z). In the following, this projection operation will be denote π(ˆ ε0r (t)).
2926
0
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
3.2 Second Step: the lossy case
16
Permittivity
14
In section 3.1, an estimation of the permittivity profile εr (z) has been obtained within an interesting confidence interval when the parameter b0 is null. In this section, we will evaluate the influence of the parameter b0 on the estimation strategy described in section 3.1. Consequently, we will consider that b0 varies from 0 up to 0.40. This range of variation has been validated experimentally.
12 10 8 6 0
0.1
0.2
0.3
0.4
0.5 Depth
0.6
0.7
0.8
0.9
1
The TDR line has been simulated for 5 different values of b0 (see figure 4). The strategy described in section 3.1 has been applied and the results have been plotted in figure 5. Clearly, it appears that as b0 increases:
5 4 Current
3 2 1
(1) the distance between εr and its estimate εˆ0r increases. (2) the distance between the two levels seems to decrease.
0 −1 −2 0
1
2
3
4
5
6
Globally, we can easily admit that there exists a transformation between the permittivity profile εr (z) and the estimate obtained for any value of b0 . Furthermore, it appears that parameter b0 plays a key-role in this transformation, i.e.:
7
Time
Fig. 2. Synthetic profile and its simulated TDR signal.
εˆΓr (z) = φ(b0 , εr (t), z, Γ) where Γ is a set of tuning parameters.
Permittivity
15
From plots in figure 5, we can consider that the transformation can be rewritten as follows:
10
5 0
0.1
0.2
0.3
0.4
0.5 Depth
0.6
0.7
0.8
0.9
1
0.1
0.2
0.3
0.4
0.5 Depth
0.6
0.7
0.8
0.9
1
εˆΓr (z) = π εˆ0r (t) + δ Γ (t)
Furthermore, according to Remark 1, we can see that flat sections in figure 5 enlarge for z ≈ 1 when b0 increases. In other words, it suggests that the time-propagation τ of the overall TDR line (i.e. the time between t = 0 and the total reflection) plays also a role in the determination of an acceptable estimate. Consequently, we should impose an additional constraint in the estimation algorithm.
Relative error (%)
20 10 0 −10
5
Fig. 3. Estimated profile estimation error.
εˆ0r (z)
4
(dashed) and its relative 3
From figure 3, it appears that the property appears to be valid within a confidence interval. Due to the complexity of the problem, estimating the permittivity profile εr (z) within a confidence interval of approximately ±15% is very interesting espescially if it requires a short time of calculation. Consequently, we consider that this strategy is very interesting and should be confronted to real cases (i.e. b0 6= 0). Remark 1. On the estimated signal, flat sections appear for z ≈ 0 and z ≈ 1. This is due to the projection in Step 3, according to the time-propagation constraint. This is absolutely normal as the time-propagation originates at t = 0 while our technique samples the magnitude of the wavefront in the time interval [t1 ; t2 ]. Hence, we get the projection in the space interval [z1 ; z2 ] with (z1 ; z2 ) 6= (0; 1). Hence, for convenience, we set εˆr ([0; z1 ]) = εˆr (z = z1 ) and εˆr ([z2 ; 1]) = εˆr (z = z2 )
Current
−20 0
2 0 0.1 0.2 0.3 0.4
1
0
−1
−2 0
1
2
3
4
5
6
7
Time
Fig. 4. Simulated TDR signal for various values of b0 . 3.3 Third Step: The Estimation Technique The objective of this section is the estimation of the spacevarying parameter εr (z) and the parameter b0 from the
2927
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
|τ − τ Γ | ≤ 1, 5% (3) τ S8- Optimize J with an appropriate optimization algorithm
20
The complexity of function δ Γ (t) and the number p of tuning parameters is of crucial importance in our context. Indeed, our objective is to reduce drastically the computational burden compared to GA based estimation techniques. Consequently, it is critical to choose δ Γ (t) such that it requires less than 2 tuning parameters. If it exceeds 2 parameters, due to the poor amount of available information on system states, we will be forced to use GA technique. Remark 2. According to standard notation, the estimation problem can be written as follows:
15
εr 0 0.1 0.2 0.3 0.4
10
5 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
min J
1
q∈Q
Fig. 5. Estimates of εr applying the procedure in section 3.1 for different values of b0 .
with - q ∈ ℜp+1 such that q = (b0 , γ1 , · · · , γp ) - Q is the set of acceptable values for q, i.e., values such that constraint (3) is satisfied.
measurement of the current at the inlet, say i(0, t). In the light of the previous sections, we know that the permittivity profile εr (z) and parameter b0 can be estimated as follows: S1- select the first wavefront of the TDR trace, say is (t), S2- construct the corresponding time interval ∆t = [t1 ; t2 ], S3- square is (t) and obtain εˆ0r (t), S4- using ∆t, construct the normalized time tn on the interval ∆tn = [0; 1], S5- set εˆr (tn ) = εˆr (t), S6- consider that: εˆΓr (z) = π εˆr (tn ) + δ Γ (tn ) (1) where: · the function δ Γ (tn ) is defined by the user, · Γ ∈ ℜp is the vector of p tuning parameters γi appearing in δ Γ (tn ). S7- Define the optimization criterion: J=
kis (t) − ˆis (t)k22 kis (t)k22
(2)
where ˆis (t) is the estimate of is (t) the intensity during the first wavefront, for a given pair (b0 , Γ). In practice, the data acquisition system provides the propagation time τ corresponding to the soil response. We will use this information in our optimization loop. The propagation time τ is given by: Z1 p εr (z)dz τ =2 0
Due to the discretization of the numerical problem, this quantity can be well approximated by: p τ = 2 εr (z)
Hence, for a given Γ, we calculate τ Γ using (1). The time propagation constraint is such that:
4. EXAMPLE 4.1 Estimation of a synthetic profile In the present paper, we desire to obtain a fast estimate of εr and b0 . Consequently, we propose a simple for δ(tn ), i.e.: δ Γ (tn ) = −γ1 tn εˆr (tn ) where γ1 is a tuning parameter to be optimized according to (2). Consequently, we will have to estimate 2 parameters which can be easily realized by a grid (e.g. by sampling the parameter space (b0 , γ1 )). Of course, this simple form is not the optimal form but permits to obtain a fast approximation of the actual εr (z). We have chosen to decouple δ Γ (tn ) from b0 in order to speed up the optimization operation as there are only two parameters to optimize. Of course, one can set δ Γ (tn ) = −b0 γ1′ tn εˆr (tn ) and obtain similar results. We have simulated the TDR line with the synthetic profile in section 3 and the parameter b0 = 0.37. A grid search has been performed in order to find the optimal parameters, say (ˆb0 , γˆ1 ), in the range ([0; 0.50], [−1; 0]) with a discretization step of 0, 01 for each parameter. Potentially, we should realize 5000 runs, but, due to the propagation time constraint, we will have about 350 runs. Though efficient, this optimization approach is not optimal. A gradient based optimization is under study. Remark 3. In practice, the parameter b0 can be approximated employing the relation given in (Noborio [2003]). This could drastically reduce the grid range in b0 and speed up the optimization process. Furthermore, it could permit to increase the complexity of δ Γ (tn ) with an additional tuning parameter. The proposed technique is compared to the algorithm based on GA presented in (Greco and Guida [2008]). As
2928
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
16
−3
x 10
0.4
Permittivity
14
2.4
0.39
12
2.2 0.38
10
2 0.37
8
1.8 0.36
6
b
0
0
0.1
0.2
0.3
0.4
0.5 Depth
0.6
0.7
0.8
0.9
1
20
Relative Error
10
1.6
0.35
1.4
0.34
1.2
0.33
1 0.8
0.32
0
0.6
0.31
−10
0.4 0.3 0.5
−20
0.51
0.52
0.53
0.54
−30
0.55 γ1
0.56
0.57
0.58
0.59
0.6
−40 −50 0
0.1
0.2
0.3
0.4
0.5 Depth
0.6
0.7
0.8
0.9
Fig. 7. Synthetic example - Contour plot of criterion (2) for the proposed method.
1
Fig. 6. Synthetic example - Estimated profile by the proposed technique (dashed) vs Greco and Guida [2008] (dash-dotted). Method
b0
γ1
J · 10−4
Our approach Greco and Guida [2008]
0.33 0.35
0.51
3.03 20.0
Date 11 : 15am 4 : 15pm
Proposed technique b0 γ1 J · 10−2 0.32 0.06
0.07 0.19
3.86 2.96
b0 0.34 0.09
Greco J · 10−2 7.83 5.48
Table 2. Experimental validation - Results of the estimation.
Table 1. Synthetic profile: Results of the optimization
4.2 Experimental validation
proposed by the authors, 12 interpolation points have been equally spaced on [0; 1] with a linear interpolation between two successive points. The GA has been launched with a population of 200 individuals. The convergence has been obtained after 207 generations which corresponds to 41400 runs. First of all, in terms of computational burden, our approach clearly reduces the load (350 runs vs 41400 runs). Concerning the quality of the estimates, figures 6 and 7 together with table 1 permit to make the decision. Plots in figure 6 show that the general feature of the profile has been captured by both techniques. But if we consider the relative estimation error, it is clear that the confidence interval of our approach is smaller. Criterion (2) for the proposed approach has been plotted (figure 7). It can be seen that the minimum obtained by the grid search lies on the limit of the acceptable domain for γ1 . As a matter of fact, the constraint on the propagation time (3) defines the set of acceptable values of γ1 . Relaxing the constraint (for example 1.75% instead of 1.5%) may give a better minimum but controlling visually the inlet current would show us that the solution is not coherent. In both cases, the estimated value of b0 is close to the real value (see table 1). Due to the complexity of the problem, both algorithms tend to obtain a compromise between the general feature of the permittivity profile and the effect of parameter b0 on the magnitude of the current. The compromise leads a reasonable solution thanks to the time-propagation constraint. Finally, on this synthetic problem, the proposed approach gives interesting results compared to a GA based technique with a very reduced numerical burden.
The previous subsection has shown the ability of the proposed method to handle a numerically generated problem, i.e., the soil properties were corresponding to the hypotheses made to build up the transmission line model. In practice, these hypotheses are not systematically met. Hence, in order to test the robustness of the proposed method to practical situations, the following experimental process has been developed. The experimental apparatus consists in a cylinder (diameter 30cm, height 40cm) filled with sand (fine/very fine 97%, coarse 1.3%, silt 1.7%). The bottom of the cylinder is in contact with a suction table connected to a Mariotte bottle. The soil was equipped with 6 three-rods TDR probes. A TDR probe is positioned vertically; it provides the signal used for parameter estimation. The 5 other probes are positioned horizontally. Considering that the medium is homogeneous in a layer, these 5 probes give an estimation of the permittivity at their positions (i.e. 10cm, 18cm, 20cm, 22cm and 30cm from soil surface). TDR waveforms are generated with a Campbell TDR100 system connected to a SDMX50 multiplexer and a CR1000 data acquisition system. We consider the drying of the sand starting from semi-wet status, i.e., the top of the column is dry while the bottom is wet. The experiment starts at Ts = 11 : 15am and ends at Te = 4 : 15pm. The numerical conditions for the test are similar to the previous section. The results are given in figure 8 and table 2. From the plots in figure 8, it can be seen that in both cases, the proposed method provides estimates that fit correctly
2929
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
provided comparable results with a very reduced numerical cost.
25
Permittivity
20 15
REFERENCES
10
H.T. Banks and K. Kunisch. Estimation techniques for distributed parameter systems. Birkh¨ auser: BostonBasel-Berlin, 1989. H. Bolvin, A. Chambarel, and Ph. Neveux. A numerical tool for transmission lines. In Proc. of International Conference on Computational Science, pages 271–278, Atlanta, USA, 2005. D. Coca and S.A. Billings. Direct parameter identification of distributed parameter systems. Int. J. of Systems Science, 31:11–17, 2000. R. Greco. Soil water content inverse profiling from single tdr waveforms. J. of Hydrology, 317:325–339, 2006. R. Greco and A. Guida. Field measurements of topsoil profiles by vertical tdr probes. J. of Hydrology, 348: 442–451, 2008. T.J. Heimovaara, J.A. Huismann, J.A. Vrugt, and W. Bouten. Obtaining the spatial distribution of water content along a tdr probe using the scem-ua bayesian inverse modelling scheme. Vadose Zone Journal, 3:1128– 1145, 2004. J. Lundstedt and S. He. A time-domain optimization yechnique for the simultaneous reconstruction of the characteristic impedance, resistance and conductance of a transmission line. J. of Electromagnetic waves and application, 10:581–601, 1996. P.G. Magnusson, G.C. Alexander, and V.K. Tripathi. Transmission Lines and Waves Propagation. 3rd Edition CRC: Boca Raton, 1992. K. Noborio. Measurement of soil water content and electrical conductivity by time domain reflectometry: a review. Computers and electronics in agriculture, 31: 213–237, 2003. B. Oswald, H.R. Benedickter, W. Bachtold, and H. Fluhler. Spatially resolved water content profiles from inverted time domain reflectometry signals. Water Resources Research, 39:15.1–15.19, 2003. S. Schlaeger. A fast tdr-inversion technique for the reconstruction of spatial soil moisture content. Hydrology and Earth system sciences discussions, 2:971–1009, 2005. J.H. Seinfeld. Identification of parameters in partial diferential equations. Chemical engineering science, 24: 65–74, 1969. G.C. Topp. Electromagnetic determination of soil water content: measurements in coaxial transmission lines. Water ressources research, 16:574–582, 1980.
5 0 0
0.1
0.2
0.3
0.4
0.5 Depth
0.6
0.7
0.8
0.9
1
0.1
0.2
0.3
0.4
0.5 Depth
0.6
0.7
0.8
0.9
1
12
Permittivity
10 8 6 4 2 0 0
Fig. 8. Experimental validation - Estimated profile by the proposed technique (solid) vs Greco and Guida [2008] (dashed). the experimental points. Compared to the quality of the estimate of Greco’s technique, our approach seems to be quite equivalent for a very reduced computational burden. From table 2, it can be seen that the optimal values of criterion (2) are quite equivalent. Greco and Guida [2008] gives lower values but the visual impression in figure 8 shows our approach behaves as ”regression estimator”. The variation of b0 from 11 : 15am to 4 : 15pm is consistent. When b0 is close to zero, it corresponds to a dry soil which the case at 4 : 15pm. Finally, we can conclude that the proposed approach permits to obtain consistent estimations of model parameters. 5. CONCLUSIONS The parameter estimation problem for a particular DPS has been treated. This DPS appears in Soil Science when soil water content is indirectly measured by TDR. Indeed, TDR signal is directly related to soil apparent permittivity. The estimation of the permittivity profile of soils is necessary in order to gain information on soil status. In order to develop an estimation technique requiring a reduced numerical burden, an interesting property of this model has been explored. A numerical study has permitted to understand the influence of model parameters onto the measured signal and to elaborate an estimation technique. The technique consists in selecting a part of the measured signal. On the other hand, a nonlinear transformation with a reduced number of tuning parameters is defined by the user. Applying this transformation to the selected part of the measured signal, provides the estimate of soil permittivity profile. The convergence to physically tractable solutions is ensured by inserting constraints into the optimization of the tuning parameters. A grid optimization of model parameters on synthetic and experimental data has shown the efficiency of the proposed technique. Compared to techniques based on GA, it has
2930
A simple method for estimating space-dependent parameters in soil science Philippe Neveux Universit´e d’Avignon et des Pays de Vaucluse UMR UAPV-INRA 1114 Environnement M´editerran´een et Mod´elisation des Agro-Hydrosyst`emes (EMMAH) Facult´e des sciences 84029 Avignon cedex 1 - France email: [email protected] Abstract: The problem of parameter estimation for a distributed parameter system (DPS) is treated. The DPS under consideration appears in Soil Science where it models a measurement method of the soil water content (namely Time Domain Reflectometry (TDR)). The model consists in a set of partial differential equation (PDE) known as the telegraphist equation. The measurement method provides the current at only one position. This increases the complexity of the parameter estimation problem. The solution proposed takes advantage of a property of the measured signal. Compared to existing solutions, the proposed technique has a very reduced computational burden. Examples show its ability to estimate space-dependent parameters with a reasonable confidence interval and a reduced computational burden compared to existing solutions. Keywords: Distributed Parameter System, Parameter identification, Telegraphist Equation, Soil Science. 1. INTRODUCTION Soil science is an area that develops complex models. They can consist of sets of Partial Derivative Equations (PDE) and/or Fractional Derivative Equations. The objective is to represent phenomenon precisely enough in order to understand the behavior of soils under various conditions. Consequently, these models require a parametric estimation in order to be tractable for study and prediction. Depending on the model structure and the available measurements, this task can be very complex. Among Soil Science areas, Time Domain Reflectometry (TDR) is known as a valuable indirect technique for the measurement of soil water content. Indirect because TDR permits to obtain an apparent relative permittivity of the soil εr . Once the apparent εr is known, the determination of the apparent water content is done using the relation between the soil apparent relative permittivity εr and the soil water content (Topp [1980]). Soil water content permits to get information on the soil status. This information is required to model various phenomenon. Hence, an accurate knowledge of the water content is necessary.
equation). The parameter estimation of this model is a very complex task. As a matter of fact, the estimation is based on the fitting of the model and the measurement made at only one position. Furthermore, the discretization of this distributed parameter system results in a model with a large number of interpolation nodes that entails a heavy computational burden. The parameters to be estimated are directly related to the value of εr at the interpolation nodes. Hence, the number of unknown parameters is large as one needs at least 500 interpolation nodes in the Finite Element Method to simulate the TDR intensity record. In order to reduce the size of the problem, we can consider two cases for the representation of the profile of εr in the soil:
Though efficient, this technique only brings a general information on the soil water content. In order to gain more detailed information, the community of Soil Science has developed several algorithms in order to estimate the soil water content at various locations in the soil from the TDR measurements (Greco [2006], Heimovaara et al. [2004], Oswald et al. [2003], Schlaeger [2005]).
(1) The permittivity is considered as a stair function (i.e. the soil is considered has a stratified medium) (Oswald et al. [2003]). In Oswald et al. [2003], the question is the optimal choice for the number of stair steps, their localization and of course, the magnitude of the stair steps.This problem is optimized by means of Genetic Algorithm (GA) search algorithm. (2) The permittivity profile can be regarded as a polynomial function (i.e. the soil is considered as a continuum) if the profile is known to be monotone (Greco [2006]). If such an information is not available, a series of interpolation points together with a linear interpolation between consecutive points, permits to approximate the permittivity profile (Greco and Guida [2008]). In both cases, a GA algorithm is used to optimize the estimate according to a fitness function.
In the time domain analysis, the TDR line is well described by a set of Partial Derivative Equations (the telegraphist 978-3-902661-93-7/11/$20.00 © 2011 IFAC
2925
10.3182/20110828-6-IT-1002.02746
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
Though interesting, these approaches lead to time consuming techniques due to extensive computations. In comparison, the measurement process takes approximately 100 nano-seconds while the optimization by means of GA requires several days. In order to obtain a faster estimate of the permittivity, we propose a simple method that aims at taking advantage of all the information present in the recorded signal.
the measured current at z = 0. Generally, parameter estimation or inverse problems for systems described by PDE require the knowledge of the state variables at various positions z (Banks and Kunisch [1989], Coca and Billings [2000], Lundstedt and He [1996], Seinfeld [1969]). Unfortunately, in TDR, this is not possible. Hence, the problem is not easy to solve. This explains why search algorithms are usually employed.
2. STATEMENT OF THE PROBLEM `a TDR technique consists in at least two rods that are ”plugged” in the soil (see figure 1). At the inlet of the rods (i.e. z = 0), we impose a step voltage u(t) and we measure the variation of the current at the inlet of the line (i.e. i(z = 0, t)). The fluctuations in the magnitude of the current are essentially due to the fluctuations of the permittivity εr along the rods. The theory of transmission line is applied to TDR (Magnusson et al. [1992]). As a matter of fact, the TDR probe can be considered as a transmission line with an inlet circuit (the impedance of the cable tester device) and outlet circuit (representing the loss of EM energy at the end of the line). The dimensionless model (the telegraphist equations) corresponding to this wave propagation phenomenon is given below (Bolvin et al. [2005]):
where:
u(t) -
`a
0 v(z, t)
Soil
?i(z, t)
εr (z) 1 z ?
RT Fig. 1. Schematic view of the TDR line. 3. MAIN RESULT
µr ∂t i = −∂z v − ai εr ∂t v = −∂z i − bv
- i(z, t) and v(z, t) are the intensity and the voltage at time t and position z ∈ [0, 1], - εr (z) is the relative permittivity of the soil, - µr ≃ 1 is the relative magnetic permeability of the soil, - a is related to the resistivity of the rods and b to the permittivity of the soil through the relation (Noborio [2003]): √ b ≡ b0 εr In the following, the contribution of the parameter a will be neglected as it is proposed in (Schlaeger [2005]), - u(t) is a deterministic input signal. It is a step function with specified raise time, - the boundary conditions: v(z = 0, t) = u(t) v(z = 1, t) RT - RT is without loss of generality representing the resistivity of the soil at the end of the line, - ∂x f is the partial derivative of function f with respect to the variable x. i(z = 1, t) =
The design of the estimation technique will be done in three steps. Firstly, an interesting property will be explored when the loss parameter b0 is set to zero. It will permit to get the basis of the estimation technique. Then, the effect of parameter b0 will be explored and the results obtained will allow us to achieve the design of the estimation technique. 3.1 First Step: the lossless case To understand the behavior of a model, one can develop a methodology based on sensitivity analysis. But in our situation, the idea is to take advantage of a property of this model when the parameter b0 is null, i.e., the lossless case. As a matter of fact, it is known that for an homogeneous medium (i.e. εr = constant), the magnitude of the signal √ i(z = 0, t) equates εr . Consequently, for inhomogeneous media, we have to verify whether this property still holds. A permittivity profile has been generated and the corresponding TDR signal has been simulated (see figure 2). In order to verify the property, the following procedure has been applied:
The objective of this work is to find a strategy that permits to obtain a fast estimate of parameter b0 and permittivity εr (z) from a single measurement i(z = 0, t) with a reduced computational cost.
Step 1- select the magnitude of the wavefront → is (t), Step 2- square the selected signal → i2s (t) ≡ εˆ0r (t), Step 3- project εˆ0r (t) onto z-space according to the timepropagation constraint: Zz p t(z) = 2 εˆ0r (ζ) dζ
The difficulty of the problem lies in the fact that the only available information for parameter estimation is
to get the estimated profile εˆ0r (z). In the following, this projection operation will be denote π(ˆ ε0r (t)).
2926
0
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
3.2 Second Step: the lossy case
16
Permittivity
14
In section 3.1, an estimation of the permittivity profile εr (z) has been obtained within an interesting confidence interval when the parameter b0 is null. In this section, we will evaluate the influence of the parameter b0 on the estimation strategy described in section 3.1. Consequently, we will consider that b0 varies from 0 up to 0.40. This range of variation has been validated experimentally.
12 10 8 6 0
0.1
0.2
0.3
0.4
0.5 Depth
0.6
0.7
0.8
0.9
1
The TDR line has been simulated for 5 different values of b0 (see figure 4). The strategy described in section 3.1 has been applied and the results have been plotted in figure 5. Clearly, it appears that as b0 increases:
5 4 Current
3 2 1
(1) the distance between εr and its estimate εˆ0r increases. (2) the distance between the two levels seems to decrease.
0 −1 −2 0
1
2
3
4
5
6
Globally, we can easily admit that there exists a transformation between the permittivity profile εr (z) and the estimate obtained for any value of b0 . Furthermore, it appears that parameter b0 plays a key-role in this transformation, i.e.:
7
Time
Fig. 2. Synthetic profile and its simulated TDR signal.
εˆΓr (z) = φ(b0 , εr (t), z, Γ) where Γ is a set of tuning parameters.
Permittivity
15
From plots in figure 5, we can consider that the transformation can be rewritten as follows:
10
5 0
0.1
0.2
0.3
0.4
0.5 Depth
0.6
0.7
0.8
0.9
1
0.1
0.2
0.3
0.4
0.5 Depth
0.6
0.7
0.8
0.9
1
εˆΓr (z) = π εˆ0r (t) + δ Γ (t)
Furthermore, according to Remark 1, we can see that flat sections in figure 5 enlarge for z ≈ 1 when b0 increases. In other words, it suggests that the time-propagation τ of the overall TDR line (i.e. the time between t = 0 and the total reflection) plays also a role in the determination of an acceptable estimate. Consequently, we should impose an additional constraint in the estimation algorithm.
Relative error (%)
20 10 0 −10
5
Fig. 3. Estimated profile estimation error.
εˆ0r (z)
4
(dashed) and its relative 3
From figure 3, it appears that the property appears to be valid within a confidence interval. Due to the complexity of the problem, estimating the permittivity profile εr (z) within a confidence interval of approximately ±15% is very interesting espescially if it requires a short time of calculation. Consequently, we consider that this strategy is very interesting and should be confronted to real cases (i.e. b0 6= 0). Remark 1. On the estimated signal, flat sections appear for z ≈ 0 and z ≈ 1. This is due to the projection in Step 3, according to the time-propagation constraint. This is absolutely normal as the time-propagation originates at t = 0 while our technique samples the magnitude of the wavefront in the time interval [t1 ; t2 ]. Hence, we get the projection in the space interval [z1 ; z2 ] with (z1 ; z2 ) 6= (0; 1). Hence, for convenience, we set εˆr ([0; z1 ]) = εˆr (z = z1 ) and εˆr ([z2 ; 1]) = εˆr (z = z2 )
Current
−20 0
2 0 0.1 0.2 0.3 0.4
1
0
−1
−2 0
1
2
3
4
5
6
7
Time
Fig. 4. Simulated TDR signal for various values of b0 . 3.3 Third Step: The Estimation Technique The objective of this section is the estimation of the spacevarying parameter εr (z) and the parameter b0 from the
2927
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
|τ − τ Γ | ≤ 1, 5% (3) τ S8- Optimize J with an appropriate optimization algorithm
20
The complexity of function δ Γ (t) and the number p of tuning parameters is of crucial importance in our context. Indeed, our objective is to reduce drastically the computational burden compared to GA based estimation techniques. Consequently, it is critical to choose δ Γ (t) such that it requires less than 2 tuning parameters. If it exceeds 2 parameters, due to the poor amount of available information on system states, we will be forced to use GA technique. Remark 2. According to standard notation, the estimation problem can be written as follows:
15
εr 0 0.1 0.2 0.3 0.4
10
5 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
min J
1
q∈Q
Fig. 5. Estimates of εr applying the procedure in section 3.1 for different values of b0 .
with - q ∈ ℜp+1 such that q = (b0 , γ1 , · · · , γp ) - Q is the set of acceptable values for q, i.e., values such that constraint (3) is satisfied.
measurement of the current at the inlet, say i(0, t). In the light of the previous sections, we know that the permittivity profile εr (z) and parameter b0 can be estimated as follows: S1- select the first wavefront of the TDR trace, say is (t), S2- construct the corresponding time interval ∆t = [t1 ; t2 ], S3- square is (t) and obtain εˆ0r (t), S4- using ∆t, construct the normalized time tn on the interval ∆tn = [0; 1], S5- set εˆr (tn ) = εˆr (t), S6- consider that: εˆΓr (z) = π εˆr (tn ) + δ Γ (tn ) (1) where: · the function δ Γ (tn ) is defined by the user, · Γ ∈ ℜp is the vector of p tuning parameters γi appearing in δ Γ (tn ). S7- Define the optimization criterion: J=
kis (t) − ˆis (t)k22 kis (t)k22
(2)
where ˆis (t) is the estimate of is (t) the intensity during the first wavefront, for a given pair (b0 , Γ). In practice, the data acquisition system provides the propagation time τ corresponding to the soil response. We will use this information in our optimization loop. The propagation time τ is given by: Z1 p εr (z)dz τ =2 0
Due to the discretization of the numerical problem, this quantity can be well approximated by: p τ = 2 εr (z)
Hence, for a given Γ, we calculate τ Γ using (1). The time propagation constraint is such that:
4. EXAMPLE 4.1 Estimation of a synthetic profile In the present paper, we desire to obtain a fast estimate of εr and b0 . Consequently, we propose a simple for δ(tn ), i.e.: δ Γ (tn ) = −γ1 tn εˆr (tn ) where γ1 is a tuning parameter to be optimized according to (2). Consequently, we will have to estimate 2 parameters which can be easily realized by a grid (e.g. by sampling the parameter space (b0 , γ1 )). Of course, this simple form is not the optimal form but permits to obtain a fast approximation of the actual εr (z). We have chosen to decouple δ Γ (tn ) from b0 in order to speed up the optimization operation as there are only two parameters to optimize. Of course, one can set δ Γ (tn ) = −b0 γ1′ tn εˆr (tn ) and obtain similar results. We have simulated the TDR line with the synthetic profile in section 3 and the parameter b0 = 0.37. A grid search has been performed in order to find the optimal parameters, say (ˆb0 , γˆ1 ), in the range ([0; 0.50], [−1; 0]) with a discretization step of 0, 01 for each parameter. Potentially, we should realize 5000 runs, but, due to the propagation time constraint, we will have about 350 runs. Though efficient, this optimization approach is not optimal. A gradient based optimization is under study. Remark 3. In practice, the parameter b0 can be approximated employing the relation given in (Noborio [2003]). This could drastically reduce the grid range in b0 and speed up the optimization process. Furthermore, it could permit to increase the complexity of δ Γ (tn ) with an additional tuning parameter. The proposed technique is compared to the algorithm based on GA presented in (Greco and Guida [2008]). As
2928
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
16
−3
x 10
0.4
Permittivity
14
2.4
0.39
12
2.2 0.38
10
2 0.37
8
1.8 0.36
6
b
0
0
0.1
0.2
0.3
0.4
0.5 Depth
0.6
0.7
0.8
0.9
1
20
Relative Error
10
1.6
0.35
1.4
0.34
1.2
0.33
1 0.8
0.32
0
0.6
0.31
−10
0.4 0.3 0.5
−20
0.51
0.52
0.53
0.54
−30
0.55 γ1
0.56
0.57
0.58
0.59
0.6
−40 −50 0
0.1
0.2
0.3
0.4
0.5 Depth
0.6
0.7
0.8
0.9
Fig. 7. Synthetic example - Contour plot of criterion (2) for the proposed method.
1
Fig. 6. Synthetic example - Estimated profile by the proposed technique (dashed) vs Greco and Guida [2008] (dash-dotted). Method
b0
γ1
J · 10−4
Our approach Greco and Guida [2008]
0.33 0.35
0.51
3.03 20.0
Date 11 : 15am 4 : 15pm
Proposed technique b0 γ1 J · 10−2 0.32 0.06
0.07 0.19
3.86 2.96
b0 0.34 0.09
Greco J · 10−2 7.83 5.48
Table 2. Experimental validation - Results of the estimation.
Table 1. Synthetic profile: Results of the optimization
4.2 Experimental validation
proposed by the authors, 12 interpolation points have been equally spaced on [0; 1] with a linear interpolation between two successive points. The GA has been launched with a population of 200 individuals. The convergence has been obtained after 207 generations which corresponds to 41400 runs. First of all, in terms of computational burden, our approach clearly reduces the load (350 runs vs 41400 runs). Concerning the quality of the estimates, figures 6 and 7 together with table 1 permit to make the decision. Plots in figure 6 show that the general feature of the profile has been captured by both techniques. But if we consider the relative estimation error, it is clear that the confidence interval of our approach is smaller. Criterion (2) for the proposed approach has been plotted (figure 7). It can be seen that the minimum obtained by the grid search lies on the limit of the acceptable domain for γ1 . As a matter of fact, the constraint on the propagation time (3) defines the set of acceptable values of γ1 . Relaxing the constraint (for example 1.75% instead of 1.5%) may give a better minimum but controlling visually the inlet current would show us that the solution is not coherent. In both cases, the estimated value of b0 is close to the real value (see table 1). Due to the complexity of the problem, both algorithms tend to obtain a compromise between the general feature of the permittivity profile and the effect of parameter b0 on the magnitude of the current. The compromise leads a reasonable solution thanks to the time-propagation constraint. Finally, on this synthetic problem, the proposed approach gives interesting results compared to a GA based technique with a very reduced numerical burden.
The previous subsection has shown the ability of the proposed method to handle a numerically generated problem, i.e., the soil properties were corresponding to the hypotheses made to build up the transmission line model. In practice, these hypotheses are not systematically met. Hence, in order to test the robustness of the proposed method to practical situations, the following experimental process has been developed. The experimental apparatus consists in a cylinder (diameter 30cm, height 40cm) filled with sand (fine/very fine 97%, coarse 1.3%, silt 1.7%). The bottom of the cylinder is in contact with a suction table connected to a Mariotte bottle. The soil was equipped with 6 three-rods TDR probes. A TDR probe is positioned vertically; it provides the signal used for parameter estimation. The 5 other probes are positioned horizontally. Considering that the medium is homogeneous in a layer, these 5 probes give an estimation of the permittivity at their positions (i.e. 10cm, 18cm, 20cm, 22cm and 30cm from soil surface). TDR waveforms are generated with a Campbell TDR100 system connected to a SDMX50 multiplexer and a CR1000 data acquisition system. We consider the drying of the sand starting from semi-wet status, i.e., the top of the column is dry while the bottom is wet. The experiment starts at Ts = 11 : 15am and ends at Te = 4 : 15pm. The numerical conditions for the test are similar to the previous section. The results are given in figure 8 and table 2. From the plots in figure 8, it can be seen that in both cases, the proposed method provides estimates that fit correctly
2929
18th IFAC World Congress (IFAC'11) Milano (Italy) August 28 - September 2, 2011
provided comparable results with a very reduced numerical cost.
25
Permittivity
20 15
REFERENCES
10
H.T. Banks and K. Kunisch. Estimation techniques for distributed parameter systems. Birkh¨ auser: BostonBasel-Berlin, 1989. H. Bolvin, A. Chambarel, and Ph. Neveux. A numerical tool for transmission lines. In Proc. of International Conference on Computational Science, pages 271–278, Atlanta, USA, 2005. D. Coca and S.A. Billings. Direct parameter identification of distributed parameter systems. Int. J. of Systems Science, 31:11–17, 2000. R. Greco. Soil water content inverse profiling from single tdr waveforms. J. of Hydrology, 317:325–339, 2006. R. Greco and A. Guida. Field measurements of topsoil profiles by vertical tdr probes. J. of Hydrology, 348: 442–451, 2008. T.J. Heimovaara, J.A. Huismann, J.A. Vrugt, and W. Bouten. Obtaining the spatial distribution of water content along a tdr probe using the scem-ua bayesian inverse modelling scheme. Vadose Zone Journal, 3:1128– 1145, 2004. J. Lundstedt and S. He. A time-domain optimization yechnique for the simultaneous reconstruction of the characteristic impedance, resistance and conductance of a transmission line. J. of Electromagnetic waves and application, 10:581–601, 1996. P.G. Magnusson, G.C. Alexander, and V.K. Tripathi. Transmission Lines and Waves Propagation. 3rd Edition CRC: Boca Raton, 1992. K. Noborio. Measurement of soil water content and electrical conductivity by time domain reflectometry: a review. Computers and electronics in agriculture, 31: 213–237, 2003. B. Oswald, H.R. Benedickter, W. Bachtold, and H. Fluhler. Spatially resolved water content profiles from inverted time domain reflectometry signals. Water Resources Research, 39:15.1–15.19, 2003. S. Schlaeger. A fast tdr-inversion technique for the reconstruction of spatial soil moisture content. Hydrology and Earth system sciences discussions, 2:971–1009, 2005. J.H. Seinfeld. Identification of parameters in partial diferential equations. Chemical engineering science, 24: 65–74, 1969. G.C. Topp. Electromagnetic determination of soil water content: measurements in coaxial transmission lines. Water ressources research, 16:574–582, 1980.
5 0 0
0.1
0.2
0.3
0.4
0.5 Depth
0.6
0.7
0.8
0.9
1
0.1
0.2
0.3
0.4
0.5 Depth
0.6
0.7
0.8
0.9
1
12
Permittivity
10 8 6 4 2 0 0
Fig. 8. Experimental validation - Estimated profile by the proposed technique (solid) vs Greco and Guida [2008] (dashed). the experimental points. Compared to the quality of the estimate of Greco’s technique, our approach seems to be quite equivalent for a very reduced computational burden. From table 2, it can be seen that the optimal values of criterion (2) are quite equivalent. Greco and Guida [2008] gives lower values but the visual impression in figure 8 shows our approach behaves as ”regression estimator”. The variation of b0 from 11 : 15am to 4 : 15pm is consistent. When b0 is close to zero, it corresponds to a dry soil which the case at 4 : 15pm. Finally, we can conclude that the proposed approach permits to obtain consistent estimations of model parameters. 5. CONCLUSIONS The parameter estimation problem for a particular DPS has been treated. This DPS appears in Soil Science when soil water content is indirectly measured by TDR. Indeed, TDR signal is directly related to soil apparent permittivity. The estimation of the permittivity profile of soils is necessary in order to gain information on soil status. In order to develop an estimation technique requiring a reduced numerical burden, an interesting property of this model has been explored. A numerical study has permitted to understand the influence of model parameters onto the measured signal and to elaborate an estimation technique. The technique consists in selecting a part of the measured signal. On the other hand, a nonlinear transformation with a reduced number of tuning parameters is defined by the user. Applying this transformation to the selected part of the measured signal, provides the estimate of soil permittivity profile. The convergence to physically tractable solutions is ensured by inserting constraints into the optimization of the tuning parameters. A grid optimization of model parameters on synthetic and experimental data has shown the efficiency of the proposed technique. Compared to techniques based on GA, it has
2930