The C-ISR method for facies inversion

Στην παρούσα διδακτορική διατριβή αναπτύσσουμε και παρουσιάζουμε τη μέθοδο 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 υδρογεωλογικές μετρήσεις, μπορούν να συσχετιστούν με την Γκαουσιανή μεταβλητή των υδρογεωλογικών φάσεων, μέσω γραμμικού μοντέλου συμμεταβλητότητας, παρότι το πρόβλημα της υπόγειας ροής δεν είναι γραμμικό. Η προσέγγιση αυτή γίνεται επαναληπτικά εντός μιας Μαρκοβιανής αλυσίδας, της οποίας τα μέλη προκύπτουν χρησιμοποιώντας ως πυρήνα μετάβασης έναν επαναληπτικό χωρικό δειγματισμό του προηγούμενου μέλους. Στην περίπτωση αυτή, οι μετρήσεις της πίεσης του υπόγειου νερού, εκτός από την έμμεση χρήση τους για την επίλυση του αντίστροφου προβλήματος, χρησιμοποιούνται και άμεσα, για την υλοποίηση των πρότερων γεωλογικών μοντέλων του υπεδάφους. Με τον τρόπο αυτό προκύπτει μια στενότερη και πιο ενημερωμένη πρότερη κατανομή, λόγω της συνδρομής της μεταβλητή αναφοράς.