Regularisation Technique for a Distributed Parameter Identification Problem

This paper examined the problem of ill-posedness of solution in identifying parameters from a given groundwater flow model. The solution approach to the problem was attempted by the method of Parameter Transformation coupled with Tikhonov Regularisation with and without Truncation which has not been explored. Convergence of the method was assessed by the L-Curve criterion. Numerical examples were presented to illustrate the efficiency of the proposed Regularisation Technique. Tikhonov Regularisation with Truncation turns to give a more realistic solution estimates when examined numerically, compared to that of Regularisation without Truncation.


Introduction
The idea behind Regularisation is to construct a family of continuous maps from data space to model space, such that the parameter of interest can be approximately recovered from noisy data.For a small noise level, the parameter of interest can be approximated very close to the true parameter value as possible.Yeh (2015) noted in the context of groundwater flow parameter estimation that the parameter estimation problem is often ill-posed and beset by instability and non-uniqueness, particularly if one seeks parameters distributed in space-time.Various techniques for alleviating the ill-posedness of parameter identification problems have been proposed by incorporating prior information into the objective function.Gomez-Hernandez et al. (1997) and Doherty (2003) made a submission in their studies that it will be very difficult to estimate field parameters in groundwater modelling without any complexity.A common approach to this complexity is the use of the Regularisation, which adds a penalty function or term to convert the constrained problem to an unconstrained problem.Different Regularisation techniques have been proposed by many researchers.Some of the Regularisation techniques explored for constraining the distributed parameter identification problems are Pilot Point Regularisation (Doherty, 2003;Kowalski et al., 2004 andAlcolea et al., 2006), Levenuberg-Marquardt Regularisation Scheme (Hanke, 1997), Preference-Base Regularisation (Neuman, 1973), Tikhonov Regularisation (Tikhonov, 1963).Parameter Choice rules (Nguyen et al., 2015), and Total Variation Regularisation (Chan & Tai, 2003;Vogel, 2002).Most of the regularisation schemes mentioned above are Optimisation schemes.
Other statistical related regularisation models have also been explored.Most of the statistical related regularisation models can be categorized into: Smoothing penalty functions, State-space Modeling in Time Series and Spatio-Temporal data, and Bayesian and Maximum Likelihood Modeling.Examples of studies where statistical approaches have been applied in smoothing penalties are the works of Wood (2004), and Chouldechova and Hastie (2015).Both studies used Generalised Additive Models in fitting and smoothing sparse data from field measurements to estimate model parameters of interest.The work of Harvey (2001), Demel and Du (2015), and Pebesma (2012) are typical examples of studies where state-space modeling has been explored.In recent times, an area that have gained much attention is Bayesian modeling.Zhu et al. (2014), Li andGoel (2006), andKneib (2008) are some of the studies that have looked at the problem from Bayesian Regularisation perspective.Most of the Bayesian regularisation approaches focuses on controlling penalties arising from fitted models to avoid overfitting of data by introducing priors to smooth uncertain variables and parameters.
the parameter estimation problem (Demissie et al., 2014;Nguyen et al., 2015).In view of this, the study aims at developing a parameter dependent Regularisation model using the basic idea proposed by Tikhonov to help stabilise the ill-posedness generated in the solution for the groundwater flow model.

The Regularisation Method
Consider a distributed parameter system with observation of the solution as in Equation (1) where  represents the parameter to be estimated, (), the parameter dependent operator, u the state variable, and C the state-to-observation map.If the parameters in the model are represented by an abstract space, such that  denotes the parameter space which contains the parameter . denote the state space which contains the vector  ⊂ , and  denotes the observation space which contains the observed data   .For uniformity sake, assume that all the three spaces described above are Hilbert Spaces.Suppose the parameter to observation map (:  → ) can be formulated as The inexact measured data can be modelled as () where  denotes noisy data or the inadequacies of the model.The problem of estimating  given the observed data   in Equation ( 1) is expressed as a constrained least squares problem by minimising the function in Equation ( 4) The constrained minimisation problem in Equation ( 4) is transformed by the PTM as 2 2 min ( ) where () is the coefficient matrix that depends on the state variable .Formulating the unconstrained least squares minimisation function () by assuming that the forward problem for solving for  in Equation ( 1) is well-posed, one obtains Equation ( 6) Incorporating a Regularisation term to address the problem of data insufficiency (or ill-posedness) yields the unconstrained regularised least square minimisation function in Equation ( 7) The term () reg Jq in Equation ( 7) is the Regularisation functional and  is a Regularisation parameter.The Regularisation function is defined as , where  0 is an initial guess of the solution, the matrix  is typically either the identity matrix   or a ( − ) ×  discrete approximation of the  − ℎ derivative operator.The side constraint () The transformed minimisation problem in Equation ( 9) is equivalent to minimising the Equation (10), with the state to observation map C taken or assumed as unity.
Solving the minimisation problem in Equation ( 10), gives: Taking the gradient of the function with respect to the parameter  yields:

Incorporating Prior Orders
Normally the penalty function imposed as side constraint to the minimisation Equation ( 13) to recover q depends on the extent of error within the solution or the field data considered.For Order Zero Regularisation, substituting n L I  into Equation ( 13) gives: In the case where the component in a solution oscillates with moderate amplitudes, order zero Regularisation may not be ideal since it dampens components that are low in magnitudes.To dampen such moderate components or error term, a penalty term that is not too low for rapid change in solution is required.This penalty term results in another form of Regularisation called "Order One Regularisation".For this reason, a moderate penalty term is added to Equation ( 14).
The objective function of order one is thus formulated as in Equation ( 15) The Equation ( 15) is minimised by a solution   of the form where is an ( − 1) ×  first derivative operator.The solution to Equation ( 16) is given by Equation ( 17) On the other hand, if the penalty term is much stronger than that of order one, a penalty term that seeks to recover the parameter to give a much better solution is used.This type of Regularisation term is an Order Two Regularisation and is an extension of order one.It is based on minimising the function defined in Equation ( 18) as subject to the solution system of the form 2 2 2 2 0 ( ) where is an ( − 2) ×  second derivative operator.Equation ( 19) has solution of the form given by Equation ( 20) The estimated solutions of the parameter  in Equations ( 17) and ( 20) are termed as Regularisation of Order One and Two respectively.They are used to dampen moderate and large error components.

Assumptions of the Method
In estimating flow parameters for a distributed parameter system using Regularisation techniques, certain assumptions need to be factored into the equation to achieve a desired solution.From the unconstrained regularised minimisation function in Equation ( 8), with noisy data   = (   ) +   , ∀  ≥ 1, the following assumptions were factored into the model: 1 The true solution for the unknown parameter   should lie or be defined in a deterministic dimensional space. ii.
The mean square convergence of the error within the data should be negligible if not zero as the number of parameters to be estimated increases indefinitely.Thus, if

( ) , n true
Error e q q   then 0 iii.

L-Curve Criterion
The L-curve criterion for selecting the Regularisation parameter depends on the trade-off between over-Regularisation and under-Regularisation at the corner of the curve.The characterisation of the corner of the L-Curve can be determined either by inspection or by estimating the point of maximum curvature.The point of maximum curvature for the distributed parameter system in Equation ( 1) according to Hansen and O'Leary, (1993) can be formulated as: where () is the solution norm, () is the residual norm, and '( ) S  is a derivative operator defined as . The parameter  is called the filter factor and  the Regularisation parameter.In this paper, both techniques are explored to determine the Optimal Regularisation parameter.

Numerical Problem 1
The numerical problem 1 is based on the steady-state diffusion equation given by Equation ( 22) The interval chosen for the model example was [0,1], with homogeneous Dirichlet boundary conditions  as in Equation ( 23) (0) = 0, and (1) = 1 (23) The source term or forcing function for the model problem was taken as Dirac delta () function, as defined in Equation (24) The parameter of interest is to determine the unknown diffusion coefficient q(x), given a measurement of the solution u(x), called the observed data (  ), and the source or sink term .The solution approach for solving the problem after discretising Equation ( 22) is by the Parameter Transformation Method (PTM) coupled with Regularisation.A detailed solution approach for the problem using the PTM is discussed in Acquah et al. (2018).For this study, the focus is how to apply Regularisation techniques to recover the unknown estimated diffusion coefficient.The fitness plot for the true and estimated parameters using the PTM with least squares is shown in Figure 2.

Implementation of the Tikhonov Regularisation for Example 1
Using the Regularisation weight of 16 15 14 0 1 [10 ,10 ,10 , ,10 ,10 ] 13), ( 17) and ( 20) are solved for the parameter   for order zero, one and two.Usually, a large λ favours a small solution norm at the expense of a small residual norm.Each chosen λ value gives a column regularised solution.The implementation of the Tikhonov Regularisation leads to several columns of regularised solutions for each order, having a convergent regularised solution.The convergent regularised solutions are not the Optimal solution to the problem.To determine the Optimal solution for computing the parameter of interest, the L-curve for Order Zero, One and Two are assessed to determine the Optimal Regularisation parameter and its corresponding Optimal solution.The simulated L-curve of Order Zero, One and Two for determining the Optimal Regularisation parameter and its corresponding Optimal solution are shown in Figure 1.
The solution curve for Order One depicts the L-shape nature required to estimate the Optimal Regularisation parameter, whilst that of Order Zero and Two for the test problem on the other hand, did not in any way depict an L-curve.This discrepancy in solution for can be attributed to either over-smoothing or under-smoothing.The Optimal Regularisation parameter corresponding to the Optimal solution was determined using the L-curve of Order One.

Determining the Optimal Regularisation Parameter
The choice of an Optimal Regularisation parameter () corresponding to the estimated parameter values for the numerical problem 1 was determined by the L-Curve of Order One.The vertical part of the Curve depicts over-Regularisation whilst the horizontal part under-Regularisation.The corner of the L-curve defines the point of maximum curvature.The Optimal Regularisation parameter is a minimiser of the estimation error norm.The plot of the solution norm against the residual norm for the Numerical Problem 1 shows a concentration of points at the corner of the L-curve.The values of the residual and the solution norms for the 1D test problem were computed and assessed for a range of values of the Regularisation parameter .The value corresponding to the Optimal Regularisation parameter was determined by inspection using the L-curve and the computed regularised solutions.
Table 1 shows the values of the residual norm error and the solution norm for a range of values of the Regularisation parameter .The points concentrated at the shape corner of the L-curve are bolded in red in Table 1 and have the parameter values  13 = 10 −4 ,  14 = 10 −3 and  15 = 10 −2 respectively.The Optimal Regularisation parameter value by inspection was  0 = 10 −3 .The value  0 = 10 −3 was obtained by observing the parameter with a minimal solution and residual error norm, as well as the ratio of solution to residual norm errors.The Optimal solution corresponding to the Optimal Regularisation parameter  0 was determined from the column regularised solutions computed.The Optimal solution corresponding to  0 = 10 −3 is indicated in Table 1.Table 1 shows the values of the residual norm error and the solution norm for a range of values of the Regularisation parameter  with weights of  = 10 −1 to  = 10 16 .
The Optimal solution is then used to estimate the parameter  of interest.The true and estimated parameter values are also shown in Table 1.Although, there are slight variations between the true and the estimated parameter values, the estimated parameter values approximate the true parameter values perfectly well.The model fitness plots of the regularised parameters for the Numerical Problem 1 is shown in Figure 2. The solution curve for the regularised parameters shows a perfect model fit between the estimated and the true parameter values compared to that of un-regularised.The relative norm-error between the estimated and true parameter values is 0.013619 , which is approximately zero signifies the stability of the estimated solution.was far less than 10 -3 , the value estimated by inspection using the L-curve method.In reality, one expects the () value computed to be greater than the value estimated by inspection to give room for improving upon the solution where necessary.In brief, the value obtained using Equation ( 21) signifies that the L-curve method over-determines the Optimal Regularisation parameter value by nearly two orders of magnitude.

Numerical Problem 2
The Numerical Problem 2 was adapted from Sun (1999).The problem is on a horizontal two-dimensional confined aquifer with boundaries AB, BC, CD, and AD.AD has a constant head boundary (  = 100 .), and the other heads are impervious, with length of  ̅̅̅̅ being 6500 m,  ̅̅̅̅ is 4500 m, and the aquifer is heterogeneous.In the model problem, the aquifer was divided into three zones (AEGD, EFHG, FBCH) with transmissivities  1 ,  2 ,  3 , respectively but their storage coefficient S being same.Also, it was assumed that there was a pumping well located at the third zone with constant pumping rate 2000  3 /, and at the beginning of the pumping, the head is constant everywhere in the aquifer with initial head data value   = 100 .Two scenarios were considered by the problem, the first was to assume values for the transmissivities and the storativity to estimate the head distribution by solving a defined forward problem generated out of the flow equations.The second scenario was to use the head distributions as an observed data to solve the inverse problem to recover the parameter of interest.A more detailed solution approach is discussed in Acquah et al. (2018).The problem is solved using the Parameter Transformation Method.The estimated parameters obtained are indicated in Table 2. Table 2 gives the true and estimated parameter values at selected time points.The estimated parameter values do not approximate to the true parameter values and has a minimal relative norm error value of 10.98314.

Implementation of Tikhonov Regularisation
Using the Regularisation weight of 16 15 14 0 1 [10 ,10 ,10 , ,10 ,10 ] , Equations ( 13), ( 17) and ( 20) are solved for the parameter   for the Orders Zero, One and Two respectively.The implementation of the Tikhonov Regularisation leads to several column regularised solutions for each order.The column regularised solutions for each Order has a convergent regularised solution.To assess Optimality, the L-curve plot of Order Zero, One and Two were determined.Figure 3 shows the L-Curve plot of Various Orders for Tikhonov Regularisation.Here, the full rank of the generated system (B(u)) in Problem 2, with its singular values was used to compute the regularised solutions.The solution curve for Order Two depicts the L-shape nature, but that of Order Zero and One show some irregularities on the solution curves of the L-curve plots.The irregularities on the solution curves of Order Zero and One can be attributed to the small singular values within the solution estimates.From Figure 3, it can be said that all the three orders show some resemblance to the L-shaped nature of the L-curve, but with distortions on the solution plots of order zero and one.In all three cases, the corner part needed to assess the Optimal Regularisation parameter with minimal error can be determined.
The Optimal Regularisation parameter for the Tikhonov case without truncation was observed on the L-curve plot of order two with a parameter value of  0 = 10 −11 .Table 3 gives the values of the Optimal regularised parameters using Tikhonov Regularisation of Order Two.Although, the estimated parameter values approximate the true parameter values to an extent, a norm-error value of 7.30020 true est qq  which is not approximaly zero, indicates that the Optimal regularised parameters still have some appreciable errors in the solution.For a solution curve that truly depict that of the L-Curve shape in all three cases, with no or very minimal error norm, some form of truncation needs be enforced.
where   ,   are the singular vectors and   the singular values.The regularised solution of the parameter q  for Order Zero, One and Two becomes respectively Equations ( 26), ( 27) and ( 28): The solution plots of the L-curve for Order Zero, One and Two for the regularised parameters are shown in

Results and Discussion
The paper has focused on the use of Regularisation techniques to recover estimated groundwater flow parameters.The Tikhonov Regularisation Method with and without TSVD, and the L-Curve Method were applied to estimate the flow parameters.Two numerical examples were considered and regularised to assess how the parameters of interest behaved.
The results show that for the studied one-dimensional test problem, the convergent regularised solutions for Order One and Two are accurate to about 15-digits of precision compared to the computed solutions of 10-digits.The Regularisation parameter corresponding to the Optimal solution (  ) for the 1D test problem was determined using L-curve of Order One with an Optimal parameter value of  0 = 10 −3 .
The solution curves for Order Zero and Two for the 1D test problem did not in any way depict an L-curve.This anomaly in solution was attributed to either over-smoothing or under-smoothing.The estimated parameters for the 1D test problem using the Optimal solution (  ) approximate perfectly well to the true parameter values with a relative norm error value of ‖  −   ‖ = 0.013619.The estimated maximum curvature value of  = 3.96365 × 10 −5 signifies that the L-curve method over-determines the Optimal Regularisation parameter of the 1D problem by nearly two orders of magnitude.For the two-dimensional test problem, it was evident from the study that the initial time set could not yield a solution even after Regularisation.When the time step was increased and solved using the PTM with least-squares, the solution became more stable with a norm error value of 10.9831 compared to the initial time step of 14.47566.Implementing Tikhonov Regularisation for Order Zero, One and Two yields a solution for the Optimal regularised parameter.The Optimal Regularisation parameter was determined using L-curve of Order Two with a norm-error value of 7.300020 . The estimated norm-error value of 7.300020, which is high and signifies that the estimated parameter values cannot be used for any prediction, though they approximate very well to the true parameter values.
Upon implementing Tikhonov Regularisation with TSVD to further reduce the error in solution, the L-curve plots of Order Zero, One and Two gave a smooth L-shaped nature needed to estimate the Optimal Regularisation parameter.The Optimal Regularisation parameter value was observed on the L-curve plot of Order Two with a parameter value of  0 = 10 −13 , and at a time step of  = 1.0 .The Optimal solution corresponding to the Optimal Regularisation parameter value of  0 = 10 −13 was determined from the column regularised solutions computed in Appendix A. The Optimal solution was used to estimate the parameter values.The estimated parameter values of the 2D test problem approximate very well to the true values with a minimal norm-error value of 0.0000043733.In conclusion, the application of the Tikhonov Regularisation Techniques together with Truncated Singular Value Decomposition in stabilising the computations of groundwater flow parameters from a discretised transformed system that depends continually on an observed data is achievable.

Conclusions and Future Work
The application of Tikhonov Regularisation techniques in stabilising the computations of the flow parameters by incorporating Order Zero, Order One and Order Two penalty functions that depend on squared norms derivative functions is possible.The L-curve criterion for determining an Optimal Regularisation parameter and its corresponding Optimal solution for the studied one-dimensional and two-dimensional test problems proved to be reliable and efficient in computing the Optimal solution.Also, the Modified Tikhonov Regularisation method with and without TSVD introduced in this study reduced the drastic effect of the small singular values on the solution.It was obvious from the study that the Tikhonov Regularisation method with TSVD gives a much better solution curve if appropriate rank-k approximation techniques are observed.
To conclude, the application of the Regularisation Methods in stabilising the computations of groundwater flow parameters from a discretised transformed system or PTM with least squares that depends continually on an observed data is feasible.Further research needs to be done to investigate if varied time step for computing inverse flow parameters can have a drastic effect on the flow parameters, and to investigate if geological information and other prior information can be directly incorporated into the inverse problems to increase stability and accuracy of the solution.

Figure 1 .
Figure 1.L-Curve of Various Orders for Numerical Problem 1

Figure 2 .
Figure 2. Fitness Plots for Regularised and Un-Regularised Parameter

Figure 3 .
Figure 3. L-Curve of Various Orders for Problem 3 Using Tikhonov7.2Implementation of Tikhonov Regularisation With TSVDTo obtain a smooth curve needed to determine the Optimal Regularisation parameter at the corner of the L-curve with minimal error estimate, some form of truncation needed to be enforced.Here, the proposed Tikhonov Regularisation with Truncated Singular Value Decomposition (TSVD) method is employed.Applying the rank-k approximation, the singular values for the generated matrix () in Numerical problem 2 with order  = 4 is determined.The estimated singular values are respectively 128.956679    , 2 8.226   Figure 4.After truncation, all the three orders depict a smooth L-shape nature.The truncation also reduced the range of values of the Regularisation parameters to values concentrated only at the corner of the curves.The corner of the L-curve in all three cases are evident and easy to identify.The Optimal Regularisation parameter values in all three cases were determined.The Optimal Regularisation parameter value with minimal error norm solution was observed at Order Two with an Optimal parameter value of  0 = 10 −13 .The Optimal solution corresponding to  0 = 10 −13 was then used to estimate the parameter  of interest.The Optimal regularised and un-regularised parameter values for Numerical Problem 2 are indicated in Table4.

Figure 4 .
Figure 4. L-Curve for various Orders Using Tikhonov with TSVD The Optimal regularised parameter values of the 2D test problem approximate very well to the true parameter values with a minimal relative norm-error value.The relative norm-error values of the regularised and unregularised parameters of 0.000043733 and 10.9831 respectively, signify that the regularised parameter values are more stable.In conclusion, Tikhonov Regularisation method with Truncated Singular Value Decomposition is efficient in reconstructing groundwater flow parameters with a very minimal error norm.

Table 1 .
Optimal Regularisation Parameter and its Corresponding Solution λ R(λ)=||A(q)  -f|| () = ‖‖ U at λ=10E-03To verify the Optimal Regularisation parameter value by inspection, the maximum curvature formula in Equation (21) was used to estimate the Optimal Regularisation value.The estimated maximum curvature value of

Table 2 .
True and Identified Parameters

Table 4 .
Optimal Regularised Parameters Using Tikhonov with TSVD