The C-ISR method for facies inversion

Abstract

In the present Doctoral thesis, we develop and present the Cosimulated Iterative Spatial Resampling (C-ISR) method for stochastic solution of the inverse problem of hydrofacies characterization in a groundwater flow system and we establish the theoretical and mathematical background of the method. The development of the method stems from the need to characterize the geological formations before any other geostatistical estimation of hydrological parameters, due to the complexity of soil types as natural entities. The spatial distribution of hydrofacies could be estimated first by classical geostatistical simulations, such as Boolean methods, Sequential indicator simulation (SIS) and truncated Gaussian simulation (TGS) or Plurigaussian simulation (PGS), then they are populated with heterogeneous hydraulic and transport parameters. However, the solution of the inverse problem is considered the most efficient practice to model the structure of hydrofacies while the physical law governing the groundwater flow system is taken into account. The parameters of the physical law, such as hydrofacies distribution, are defined by optimizing the response of the system while solving the physical law (forward problem), compared to the observations of a physical variable such as the hydraulic head. In most cases, the inverse problems are ill-posed, which means they have not a unique solution, the solutions do not depend continuously on data, or a solution does not exist. Moreover, the solutions of the problem may be affected significantly by small changes in the observations and cause large variations in the parameters estimation, making the system ill-conditioned. So, the limited number of observations and the nature of system may impose difficulties in solving an inverse problem. Thus, the inverse problem is a classical optimization problem when the system of equations is linear, where the best-unbiased estimator can be found by the minimization of a least squares error criterion. In ill-conditioned or ill-posed problems, the objective is not only to find an unbiased estimator, but also one with a stable behavior. In weakly linear or non-linear systems, the objective function is linearized and iterative gradient methods are applied to find the estimator. However, the iterative process of optimization in those cases is often a demanding task with high computational cost, while gradient methods can get stuck at local optima. Therefore, the use of stochastic methods allowing the sampling of model space, such as the Markov chain Monte Carlo (McMC) methods, is often preferable. McMC methods are often used under the Bayesian perspective in solving the inverse problem. In Bayesian inference, the prior information is updated iteratively as new members are added to the chain. The members of the chain constitute the final a-posteriori distribution of the parameters. The McMC instead of Monte Carlo methods preferentially visit the model space where the posterior density is high, even if the dimension of the model space is large. The implementation of McMC requires the definition of a transition kernel from an accepted model to a new model, a criterion to accept a member in the chain and a criterion to interrupt the chain. In this work, we adopt the iterative spatial resampling (ISR) technique as the transition kernel, either when the objective is to sample the posterior distribution of parameters or to reach an optimal solution. In sampling the posterior, the candidate models must be independent of the last accepted member of the chain, while the criterion of Metropolis-Hastings or Mosegaard and Tarantola (1995) could be used to accept a model to the chain. In the case of optimization, we adopt the implementation of independent Markov chains to reach an optimal solution. A candidate model is accepted in a chain when its likelihood is better than the last member of the chain, while the interruption of the chain is stochastic with the probability to interrupt becoming higher when the likelihood of a candidate model is high. This way, bias with the opposite direction is introduced to the samples of the posterior. Thus, an approach to posterior distribution is achieved. This mechanism reaches an optimal solution avoiding a large number of geostatistical simulations and forward problem solutions. The novelty of the C-ISR algorithm is the iterative use of cosimulation between hydrofacies and the reference data as an auxiliary variable, in order to gradually improve the path to the optimal solution within a constantly improving Markov chain. Cosimulation for modeling a discrete variable such as the hydrofacies distribution has not been applied in inversion yet, due to the nonlinearity relation between the response variable and unknown parameters. In most cases, the response variable is used as indirect data to evaluate the prior models and drive the search path. The C-ISR method exploits the available information on the relationship between hydrofacies and the physical variable, to produce valid realizations by using cosimulation. More specifically, our method is based on the approach that, instantaneously, the Normal scores transformed hydrological measurements can be correlated with the Gaussian variable of hydrofacies, through a linear coregionalization model, although the underground flow problem is not linear. The approach is used repeatedly within a Markov chain, while the members of the chain result from an iterative spatial resampling transition kernel. In this case, apart from their indirect use for inversion, groundwater pressure measurements are used directly, in order to evaluate the prior geological model of the subsurface. This, results in a narrower and more informed prior distribution, due to the support of the reference variable. The effectiveness of our method is demonstrated by an example application on a synthetic underdetermined inverse problem in aquifer characterization. The results show that the C-ISR method is faster and more accurate as compared to plain ISR.

Στην παρούσα διδακτορική διατριβή αναπτύσσουμε και παρουσιάζουμε τη μέθοδο Cosimulated Iterative Spatial Resampling (C-ISR) για την στοχαστική επίλυση του αντίστροφου προβλήματος το οποίο σχετίζεται με τον χαρακτηρισμό των υδρογεωλογικών φάσεων του υπεδάφους και θέτουμε το θεωρητικό και μαθηματικό υπόβαθρο της μεθόδου. Η ανάπτυξη της μεθόδου βασίστηκε στην ανάγκη του χαρακτηρισμού των υδρογεωλογικών φάσεων του υπεδάφους πριν από τον προσδιορισμό/εκτίμηση άλλων υδρογεωλογικών παραμέτρων, όπως της υδραυλικής αγωγιμότητας, του πορώδους και της ζώνης κορεσμού. Η πολυπλοκότητα της γεωλογικής διαδικασίας, οι φυσικές και χημικές αντιδράσεις, όπως επίσης η δυσκολία αναπαράστασης της υδραυλικής αγωγιμότητας ως συνεχούς μεταβλητής (Matheron, 1967; Emsellem and De Marsily, 1971), καθιστούν απαραίτητη τη μοντελοποίηση των υδρογεωλογικών φάσεων εκ των προτέρων. Επιπλέον, ο χαρακτηρισμός των πετρωμάτων είναι πάγια πρακτική στις γεωλογικές μελέτες σε σχέση με άλλες παραμέτρους όπως η αγωγιμότητα και αυτό οδηγεί σε πιο ακριβή προσδιορισμό της χωρικής κατανομής των φάσεων, αφού υπάρχει μεγαλύτερη διαθεσιμότητα δεδομένων. Η χωρική κατανομή των υδρογεωλογικών φάσεων μπορεί να εκτιμηθεί με κλασικές μεθόδους γεωστατιστικής προσομοίωσης όπως οι Boolean μέθοδοι, η Sequential indicator simulation (SIS) και η truncated Gaussian simulation (TGS) ή Plurigaussian simulation (PGS). Στη συνέχεια άλλες υδρογεωλογικές παράμετροι μπορούν να προσδιοριστούν. Η προσομοίωση της χωρικής κατανομής των υδρογεωλογικών φάσεων μπορεί να γίνει πιο ακριβής, λαμβάνοντας υπόψη βοηθητικές μεταβλητές μέσω από κοινού προσομοίωσης, όπως για παράδειγμα την υδραυλική πίεση, καθώς συχνά υπάρχουν διαθέσιμα δεδομένα της πίεσης σε υδρογεωλογικές λεκάνες ενδιαφέροντος. Ωστόσο, η από κοινού προσομοίωση απαιτεί τη γνώση της συσχέτισης μεταξύ των φάσεων και της βοηθητικής μεταβλητής, ενώ δεν αξιολογείται η πιθανοφάνεια των παραμέτρων, καταλήγοντας έτσι σε μη αποδεκτές λύσεις. Κατά κανόνα, παρά το πλήθος των εκάστοτε διαθέσιμων δειγμάτων από γεωτρήσεις, ο αριθμός αυτός δεν είναι ποτέ αρκετός για την επίτευξη της επιθυμητής ακρίβειας στην απεικόνιση του υπεδάφους. Η πιο σύγχρονη πρακτική αντιμετώπιση του προβλήματος αυτού, είναι η επίλυση του φυσικού νόμου που διέπει το φαινόμενο μεταφοράς ροής (ευθύ πρόβλημα) υποθέτοντας ότι οι παράμετροι του συστήματος είναι γνωστές, ενώ στη συνέχεια γίνεται η ρύθμισή τους. Η πρακτική αυτή είναι γνωστή ως επίλυση του Αντίστροφου Προβλήματος, καθώς οι υδρογεωλογικές φάσεις αντιμετωπίζονται πλέον ως παράμετροι ενός συστήματος και η φυσική μεταβλητή (π.χ. υδραυλική πίεση) ως τυχαία μεταβλητή. Στο αντίστροφο πρόβλημα, οι μετρήσεις της φυσικής μεταβλητής και η απόκριση του φυσικού νόμου χρησιμοποιούνται για την κατανόηση της συμπεριφοράς των παραμέτρων και του προσδιορισμού τους. Αν και η πρακτική του προσδιορισμού των φάσεων είναι πιο αποτελεσματική επιλύοντας το αντίστροφο πρόβλημα καθώς χρησιμοποιείται μία επιπλέον πληροφορία, δηλαδή ο φυσικός νόμος, μία σειρά από άλλα προβλήματα δημιουργούνται στα οποία οι ερευνητές καλούνται να δώσουν απαντήσεις. Πιο συγκεκριμένα, τα αντίστροφα προβλήματα είναι συνήθως ασθενώς/κακώς τεθειμένα (ill-posed), δηλαδή επιδέχονται από καμία λύση έως άπειρες λύσεις. Επίσης, τα αντίστροφα προβλήματα μπορεί να είναι μη επιλύσιμα εξαιτίας της υπολογιστικής ακρίβειας και τότε χαρακτηρίζονται ως ασταθή συστήματα (ill-conditioned systems). Τα προβλήματα αυτά δημιουργούνται κυρίως λόγω του περιορισμένου αριθμού δειγμάτων της φυσικής μεταβλητής αλλά και της φύσης του συστήματος. Στην περίπτωση που το σύστημα εξισώσεων είναι γραμμικό, το αντίστροφο πρόβλημα είναι ένα κλασσικό πρόβλημα βελτιστοποίησης, στο οποίο ο βέλτιστος αμερόληπτος εκτιμητής των παραμέτρων μπορεί να δοθεί από το κριτήριο των ελαχίστων τετραγώνων. Στην περίπτωση γραμμικών ασθενών ή ασταθών συστημάτων, πέραν από την αναζήτηση ενός αμερόληπτου εκτιμητή, αντικείμενο είναι η εύρεση ενός εκτιμητή ο οποίος μπορεί να είναι μεροληπτικός αλλά να δίνει πιο δυνατές ή σταθερές λύσεις στο σύστημα των εξισώσεων. Στην περίπτωση (ασθενώς) μη γραμμικών συστημάτων, η αντικειμενική συνάρτηση γραμμικοποιείται και η διαδικασία εύρεσης του βέλτιστου εκτιμητή γίνεται επαναληπτικά με τους αλγόριθμους κλίσης. Η επαναληπτική διαδικασία εύρεσης τους βέλτιστου εκτιμητή στις περιπτώσεις ασθενώς μη γραμμικών συστημάτων ή τελείως μη γραμμικών συστημάτων συχνά είναι μία επίπονη διαδικασία με μεγάλο υπολογιστικό κόστος (μεγάλος αριθμός παραμέτρων), ενώ είναι πολύ πιθανό το αποτέλεσμα του βέλτιστου εκτιμητή να αναφέρεται σε ένα τοπικό βέλτιστο της αντικειμενικής συνάρτησης. Για τον λόγο αυτό η χρήση στοχαστικών μεθόδων που επιτρέπουν το συστηματικό δειγματισμό του πεδίου τιμών των παραμέτρων, όπως οι μέθοδοι Markov chain Monte Carlo, είναι προτιμότερη, ξεπερνώντας τα προβλήματα των αιτιοκρατικών ή άλλων στοχαστικών μεθόδων. Στην αντιστροφή του προβλήματος οι μέθοδοι McMC χρησιμοποιούνται συνήθως σε συνδυασμό με τη στατιστική κατά Bayes, δηλαδή η εκ προοιμίου (a-priori) κατανομή των παραμέτρων αναθεωρείται επαναληπτικά όσο νέα μέλη προστίθενται στην αλυσίδα. Τα μέλη της αλυσίδας αποτελούν δείγματα της μεταγενέστερης (a-posteriori) κατανομής των παραμέτρων. Οι μέθοδοι McMC έχουν την ικανότητα να «επισκέπτονται» το χώρο των παραμέτρων όπου η πυκνότητα της μεταγενέστερης κατανομής είναι μεγάλη ακόμα και όταν ο χώρος των παραμέτρων είναι πολλών διαστάσεων. Η εφαρμογή των McMC απαιτεί τον ορισμό ενός πυρήνα μετάβασης από το ένα μέλος της αλυσίδας σε ένα υποψήφιο μέλος, ενός κριτηρίου επιλογής των μελών της αλυσίδας και ενός κριτηρίου για να διακοπεί η αλυσίδα. Ο καθορισμός των παραπάνω κριτηρίων πρέπει να είναι τέτοιος ώστε τα μέλη της αλυσίδας να είναι ανεξάρτητα μεταξύ τους και να αποτελούν δείγματα της a-posteriori κατανομής. Η διακοπή της αλυσίδας συνήθως εξαρτάται από τον επιθυμητό αριθμό των δειγμάτων και τον διαθέσιμο χρόνο για την υλοποίηση της αλυσίδας. Στη παρούσα διατριβή, τα υποψήφια μέλη της αλυσίδας δημιουργούνται με γεωστατιστική προσομοίωση και με τυχαίο χωρικό δειγματισμό του αμέσως προηγούμενου μέλους της αλυσίδας και των a-priori δεδομένων, δηλαδή ένα προσαρμοσμένο αλγόριθμο Gibbs Sampling. Όταν ο στόχος της έρευνας είναι ο ορισμός της πρότερης κατανομής, τότε ένα υποψήφιο μέλος της αλυσίδας κρίνεται βάσει του κριτηρίου που χρησιμοποιείται στον αλγόριθμο Metropolis-Hastings ή του κριτηρίου που προτείνεται από τους Mosegaard and Tarantola (1995). Τα υποψήφια μέλη της αλυσίδας πρέπει να κρίνονται αφού έχει επιτευχθεί η ανεξαρτησία τους από το προηγούμενο μέλος της αλυσίδας. Στην περίπτωση που ο στόχος είναι η εύρεση ενός βέλτιστου εκτιμητή, υιοθετούμε την υλοποίηση ανεξάρτητων Μαρκοβιανών αλυσίδων, ενώ κάθε μέλος της αλυσίδας γίνεται αποδεκτό εφόσον η πιθανοφάνειά του είναι μεγαλύτερη από το προηγούμενο μέλος της αλυσίδας. Αυτό οδηγεί στην επιλογή μεροληπτικού δείγματος της ύστερης κατανομής, αφού τα τελευταία μέλη κάθε αλυσίδας που αποτελούν την κατανομή έχουν προκύψει με υψηλή πιθανοφάνεια. Η διακοπή κάθε αλυσίδας γίνεται με στοχαστικό κριτήριο, στο οποίο η πιθανότητα διακοπής της αλυσίδας αυξάνεται όσο η πιθανοφάνεια των μελών μίας αλυσίδας αυξάνεται. Με αυτόν τον τρόπο προκύπτει ένα μεροληπτικό δείγμα της a-posteriori κατανομής αλλά με αντίθετη κατεύθυνση της μεροληψίας που εισάγεται με την επιλογή των μελών κάθε αλυσίδας. Έτσι, επιτυγχάνεται μία προσέγγιση της ύστερης κατανομής. Αυτός ο μηχανισμός οδηγεί στην εύρεση του βέλτιστου εκτιμητή αποφεύγοντας έναν μεγάλο αριθμό γεωστατιστικών υλοποιήσεων και της επίλυσης του ευθέως προβλήματος. Σε προηγούμενες μελέτες οι οποίες χρησιμοποιούν τις μεθόδους McMC, δεν λαμβάνεται υπόψη η συσχέτιση της φυσικής μεταβλητής και των παραμέτρων για την παραγωγή των γεωστατιστικών υλοποιήσεων. Σε αυτή την περίπτωση, οι μετρήσεις της πίεσης του υπόγειου νερού χρησιμοποιούνται μόνον εμμέσως, για την αξιολόγηση των πρότερων γεωλογικών μοντέλων του υπεδάφους. Ο περιορισμός αυτός οφείλεται στην μη γραμμική συσχέτιση μεταξύ υδραυλικών μετρήσεων και παραμέτρων του υπεδάφους. Ως εκ τούτου, οι ευρέως γνωστές μέθοδοι των από κοινού kriging και από κοινού προσομοίωσης που βασίζονται στην γραμμική εκτίμηση, δεν μπορούν να χρησιμοποιηθούν απευθείας, γιατί καταλήγουν σε μη αποδεκτές λύσεις. Αντίθετα, η μέθοδος C-ISR εκμεταλλεύεται την πληροφορία των μετρήσεων της φυσικής μεταβλητής για την παραγωγή των γεωστατιστικών υλοποιήσεων. Πιο συγκεκριμένα, η μέθοδος στηρίζεται στην προσέγγιση ότι, στιγμιαία, οι μετασχηματισμένες σε Normal Scores υδρογεωλογικές μετρήσεις, μπορούν να συσχετιστούν με την Γκαουσιανή μεταβλητή των υδρογεωλογικών φάσεων, μέσω γραμμικού μοντέλου συμμεταβλητότητας, παρότι το πρόβλημα της υπόγειας ροής δεν είναι γραμμικό. Η προσέγγιση αυτή γίνεται επαναληπτικά εντός μιας Μαρκοβιανής αλυσίδας, της οποίας τα μέλη προκύπτουν χρησιμοποιώντας ως πυρήνα μετάβασης έναν επαναληπτικό χωρικό δειγματισμό του προηγούμενου μέλους. Στην περίπτωση αυτή, οι μετρήσεις της πίεσης του υπόγειου νερού, εκτός από την έμμεση χρήση τους για την επίλυση του αντίστροφου προβλήματος, χρησιμοποιούνται και άμεσα, για την υλοποίηση των πρότερων γεωλογικών μοντέλων του υπεδάφους. Με τον τρόπο αυτό προκύπτει μια στενότερη και πιο ενημερωμένη πρότερη κατανομή, λόγω της συνδρομής της μεταβλητή αναφοράς.