import numpy as np import os import subprocess import time def Modelisation_VTEM(Wave,Gate,Hauteur,Res,Ep): """ Génère le fichier .cfl et exécute Airbeo pour modéliser la réponse 1D d'un modele synthétique du sous-sol Paramètres ---------- Wave : array [n x 2] Forme de l'onde du VTEM-plus. Gate : array [n x 2] Centre et largeur des Gate du VTEM-plus en ms. Hauteur : int Hauteur du vol en m. Res : array Resistances des couches du modèle du sous-sol. Ep : array Épaisseurs des couches du modèle du sous-sol. Returns ------- Data : ndarray [] Courbes dbz/dt en nV/m^2 (réponses du modèle) """ gg = Gate.shape vv = Wave.shape[0] Npar = len(Res) # dirs = os.path.dirname(os.path.realpath(__file__)) #%% Partie 1 : Création du fichier de commande .cfl AirbeoCfl = open('airbeo.cfl','w') # Entête du fichier Header = 'Modelisation VTEM-plus' AirbeoCfl.write("{}\n".format(Header)) AirbeoCfl.write("1 0 0 0 ") AirbeoCfl.write("!TDFD, DO1D, PRFL, ISTOP\n") AirbeoCfl.write('1 {} 0 1 {} 2 13.958000 '.format(vv, gg[0])) AirbeoCfl.write('!ISW, NSX, STEP, UNITS, NCHNL, KRXW, OFFTIME\n') #Forme du courant for i in range(0,vv-1): AirbeoCfl.write('{0:.6f} {1:.6f}\n'.format(Wave[i, 0], Wave[i, 1])) AirbeoCfl.write('{0:.6f} {1:.6f} '.format(Wave[-1, 0], Wave[-1, 1])) AirbeoCfl.write(' ! TXON, WAVEFORM(849)\n') #Fenêtre de temps for i in range(0,gg[0]-1): AirbeoCfl.write('{0:.6f} \n'.format(Wave[-1, 0]+Gate[i, 0])) AirbeoCfl.write('{0:.6f} '.format(Wave[-1, 0]+Gate[-1, 0])) AirbeoCfl.write(' ! TMS(44 centre)\n') for i in range(0,gg[0]-1): AirbeoCfl.write('{0:.6f}\n'.format(Gate[i, 1])) AirbeoCfl.write('{0:.6f} '.format(Gate[-1, 1])) AirbeoCfl.write(' ! WIDTH(44)\n') AirbeoCfl.write('0 13 0 ! TXCLN, CMP, KPPM\n') AirbeoCfl.write('531 4 ! TXAREA, NTRN\n') AirbeoCfl.write('0 0 0 ! ZRX, XRX, YRX(1)\n') AirbeoCfl.write( '1 1 0 1000 !NSTAT, SURVEY, BAROMTRC, LINE_TAG \n') AirbeoCfl.write( '0 0 {} 0 10 !EAST(1), NORTH(1), ALT(1), BEARING, DSTAT \n'.format(Hauteur)) AirbeoCfl.write( '{} 1 {} 0 ! NLAYER, QLYR, NLITH, GND_LEVEL\n'.format(Npar, Npar)) for i in range(0,Npar): AirbeoCfl.write('{} 1 1 1 0 0 0'.format(Res[i])) AirbeoCfl.write( ' ! RES, SIG_T, RMU, REPS, CHRG, CTAU, CFREQ(1) - Layer\n') for i in range(0,Npar-1): AirbeoCfl.write('{0:d} {1:.2f} '.format(i+1, Ep[i])) AirbeoCfl.write(' !LITH, THK\n') AirbeoCfl.write('{} !LITH, THK\n'.format(Npar)) del(i) AirbeoCfl.close() #%% Partie 2 : Modelisation et lecture des courbes dBz/dt modélisées # dirAirbeo = 'airbeo' subprocess.call('airbeo') time.sleep(0.1) Resultmf1 = open('Airbeo.mf1', 'r') count = 0 for line in Resultmf1: count +=1 if count == 16: lineMf1 = line return np.array([float(i) for i in lineMf1.split()[8:]]) def Modelisation_NanoTEM(Wave, Gate, Res, Ep): """ Génère le fichier de commande de Leroi et modélise la réponse 1D du NanoTEM en boucle centrale. Paramètres ---------- Wave : array [n x 2] Forme de l'onde du NanoTEM. Gate : array [n x 2] Centre et largeur des Gate du NanoTEM en ms. Res : array Resistances des couches du modèle du sous-sol. Ep : array Épaisseurs des couches du modèle du sous-sol. Returns ------- Data : ndarray [] Courbes dbz/dt en nV/m^2 (réponses du modèle) """ gg = Gate.shape Npar = len(Res) LeroiCfl = open('Leroi.cfl', 'w') # Entête du fichier Header = 'Modelisation NanoTEM' LeroiCfl.write("{}\n".format(Header)) LeroiCfl.write( '1 0 0 1 0 !TDFD, DO3D, ISYS, PRFL, ISTOP\n') LeroiCfl.write('0 4 31 2 {:.6f} {:.6f}'.format(3.9078, 3.9047)) LeroiCfl.write(' ! STEP, NSX, NCHNL, KRXW, REFTYM, OFFTIME\n') #Waveform LeroiCfl.write('{:.6f} {:.6f}\n'.format(Wave[0, 0], Wave[0, 1])) LeroiCfl.write('{:.6f} {:.6f}\n'.format(Wave[1, 0], Wave[1, 1])) LeroiCfl.write('{:.6f} {:.6f}\n'.format(Wave[2, 0], Wave[2, 1])) LeroiCfl.write('{:.6f} {:.6f}'.format(Wave[3, 0], Wave[3, 1])) LeroiCfl.write(' ! TXON, TXAMP(4)\n') #Time Gate for i in range(0,gg[0]-1): LeroiCfl.write('{:.6f} \n'.format(Gate[i,0])) LeroiCfl.write( '{:.6f} ! TMS(31, GATE center) \n'.format(Gate[-1, 0])) for i in range(0, gg[0]-1): LeroiCfl.write('{:.6f} \n'.format(Gate[i, 1])) LeroiCfl.write( '{:.6f} ! WIDTH(31)\n'.format(Gate[-1, 1])) #Survey Type LeroiCfl.write('1 ! SURVEY_TYPE\n') #Tx Position and Moment LeroiCfl.write( '1 1 1 1 4 1 ! NLINES,MRXL,NTX,SOURCE_TYPE,MXVRTX, TXMNT\n') LeroiCfl.write('4 0 ! N_VRTX(1) TxZ\n') LeroiCfl.write('0 20\n') LeroiCfl.write('20 20\n') LeroiCfl.write('20 0\n') LeroiCfl.write('0 0 !SXE(I,J), SXN(I,J)\n') #Rx Position and Moment LeroiCfl.write( '100 1 1 1 11 ! LINE, IDTX, RX_TYPE,NRX, UNITS \n') LeroiCfl.write( '3 0 0 3 0 250 ! CMP, SV_AZM, KNORM, IPLT, IDH, RXMNT\n') LeroiCfl.write('10 10 0 ! RXE, RXN RXZ\n') #Layer Earth Modele LeroiCfl.write('{} 0 {}'.format(Npar, Npar)) LeroiCfl.write(' !NPAR, NPLATE, NLITH\n') for i in range(0,Npar): LeroiCfl.write('{} -1 1 1 0 0 1'.format(Res[i])) LeroiCfl.write( ' ! RES,SIG_T(1),RMU, REPS, CHRG, CTAU, CFREQ(1) - Layer\n') for i in range(0,Npar-1): LeroiCfl.write('{} {} '.format(i+1, Ep[i])) LeroiCfl.write(' !LITH(i), THK(i)\n') LeroiCfl.write('{} !LITH(i), THK(i)\n'.format(Npar)) LeroiCfl.close() #%% Partie 2 : Modelisation et lecture des courbes dBz/dt modélisées subprocess.call('Leroi') time.sleep(0.1) Resultmf1 = open('Leroi.mf1', 'r') count = 0 for line in Resultmf1: count += 1 if count == 18: lineMf1 = line # lineMf1 = linecache.getline('Leroi.mf1', 18) return np.array([float(i) for i in lineMf1.split()[4:]]) def Inversion_VTEM(Data0, Wave, Gate, Hauteur, Res, Ep, flag = 0): """ Génère le fichier .cfl et .inv et exécute Airbeo pour inverser les données des courbes Data0. Paramètres ---------- Data0 : array_like Courbes dbz/dt en nV/m^2 à inverser. Wave : array_like Forme de l'onde du VTEM-plus. Gate : array_like Centre et largeur des Gate du VTEM-plus en ms. Hauteur : int Hauteur du vol en m. Res : array Resistances des couches du modèle du sous-sol. Ep : array Épaisseurs des couches du modèle du sous-sol. flag : int, optional flag=0 : inversion avec nombre minimum de chouches flag=1 : inversion lisse (multicouches) Default is 0. Returns ------- Data1 : ndarray Courbes dbz/dt en nV/m^2 qui ajuste les données (Data0) ResFinal : ndarray Résistivités rétournées par l'inversion EpFinal : ndarray Épaisseurs rétournées par l'inversion ProfFinal : ndarray Profondeurs rétournées par l'inversion """ vv = Wave.shape[0] gg = Gate.shape[0] Npar = Res.shape[0] nData = Data0.shape[0] ModeleFinal = np.zeros((Npar,4)) #%% Partie 1 : Création du fichier de commande .cfl AirbeoCfl = open('airbeo.cfl', 'w') # Entête du fichier Header = 'Inversion VTEM-plus' AirbeoCfl.write("{}\n".format(Header)) AirbeoCfl.write("1 1 0 0 ") AirbeoCfl.write("!TDFD, DO1D, PRFL, ISTOP\n") AirbeoCfl.write('1 {} 0 1 {} 2 13.958000 '.format(vv, gg)) AirbeoCfl.write('!ISW, NSX, STEP, UNITS, NCHNL, KRXW, OFFTIME\n') #Forme du courant for i in range(0, vv-1): AirbeoCfl.write('{0:.6f} {1:.6f}\n'.format(Wave[i, 0], Wave[i, 1])) AirbeoCfl.write('{0:.6f} {1:.6f} '.format(Wave[-1, 0], Wave[-1, 1])) AirbeoCfl.write(' ! TXON, WAVEFORM(849)\n') #Fenêtre de temps for i in range(0, gg-1): AirbeoCfl.write('{0:.6f} \n'.format(Wave[-1, 0]+Gate[i, 0])) AirbeoCfl.write('{0:.6f} '.format(Wave[-1, 0]+Gate[-1, 0])) AirbeoCfl.write(' ! TMS(44 centre)\n') for i in range(0, gg-1): AirbeoCfl.write('{0:.6f}\n'.format(Gate[i, 1])) AirbeoCfl.write('{0:.6f} '.format(Gate[-1, 1])) AirbeoCfl.write(' ! WIDTH(44)\n') AirbeoCfl.write('0 13 0 ! TXCLN, CMP, KPPM\n') AirbeoCfl.write('531 4 ! TXAREA, NTRN\n') #Rx Position and Moment AirbeoCfl.write('0 0 0 ! ZRX, XRX, YRX(1)\n') AirbeoCfl.write('1 ! NPASS\n') #Layer Earth Model AirbeoCfl.write( '{} 1 {} 0 ! NLAYER, QLYR, NLITH, GND_LEVEL\n'.format(Npar, Npar)) for i in range(0, Npar): AirbeoCfl.write('{} 1 1 1 0 0 0'.format(Res[i])) AirbeoCfl.write( ' ! RES, SIG_T, RMU, REPS, CHRG, CTAU, CFREQ(1) - Layer\n') for i in range(0, Npar-1): AirbeoCfl.write('{0:d} {1:.2f} '.format(i+1, Ep[i])) AirbeoCfl.write(' !LITH, THK\n') AirbeoCfl.write('{} !LITH, THK\n'.format(Npar)) #Inversion Control Parameter if flag == 0: AirbeoCfl.write( '90 1 0 1 2 ! MAXITS, CNVRG, NFIX, INVPRT,OUTPRT\n') elif flag == 1: AirbeoCfl.write( '90 1 {} 1 2 ! MAXITS, CNVRG, NFIX, INVPRT,OUTPRT\n'.format(nPar-1)) for i in range(0,nPar-1): AirbeoCfl.write( '1 {} 2 ! CTYPE, LYR_INDX, KPAR\n'.format(i+1)) AirbeoCfl.close() #Partie 2 : Création du fichier .inv et inversion AirbeoInv = open('airbeo.inv','w') AirbeoInv.write('/ Inversion de Data0\n') AirbeoInv.write( '1 2 0 3 1 ! NSTAT, SURVEY, BAROMTRC, KCMP, ORDER\n') AirbeoInv.write('0.1 ! DATA0 floor\n') AirbeoInv.write('0 0 0 ! N0STAT, N0CHNL, N0PTS\n') AirbeoInv.write('1000 0 0 {}'.format(Hauteur)) for i in range(nData-1): AirbeoInv.write(' {:.6E} '.format(Data0[i])) AirbeoInv.write(' {:.6E}\n'.format(Data0[-1])) AirbeoInv.close() #Partie 3 : Inversion et lecture du résultat subprocess.call('airbeo') ResultOut = open('Airbeo.out','r') for line in ResultOut: if 'Final Model after' in line: for i in range(4): ResultOut.readline() for j in range(Npar-1): ResultLine = ResultOut.readline() ModeleFinal[j,:] = np.array([float(i) for i in ResultLine.split()[1:]]) ResultLine = ResultOut.readline() ModeleFinal[Npar-1, 0] = np.array([float(i) for i in ResultLine.split()[1:]]) for i in range(6): ResultOut.readline() ResultLine = ResultOut.readline() Rms = float(ResultLine.split()[4]) for i in range(20): ResultOut.readline() ResultLine = ResultOut.readline() Data1 = np.array([float(i) for i in ResultLine.split()[:]]) ResFinal = ModeleFinal[:,0] ProfFinal = ModeleFinal[:-1,1] EpFinal = ModeleFinal[:-1,2] return(Data1, ResFinal, ProfFinal, EpFinal, Rms) def Inversion_NanoTEM(Data0, Wave, Gate, Res, Ep, flag=0): """ Génère le fichier .cfl et .inv et exécute Airbeo pour inverser les données des courbes Data0. Paramètres ---------- Data0 : array_like Courbes dbz/dt en nV/m^2 à inverser. Wave : array_like Forme de l'onde du VTEM-plus. Gate : array_like Centre et largeur des Gate du VTEM-plus en ms. Res : array Resistances des couches du modèle du sous-sol. Ep : array Épaisseurs des couches du modèle du sous-sol. flag : int, optional flag=0 : inversion avec nombre minimum de chouches flag=1 : inversion lisse (multicouches) Default is 0. Returns ------- Data1 : ndarray Courbes dbz/dt en nV/m^2 qui ajuste les données (Data0) ResFinal : ndarray Résistivités rétournées par l'inversion EpFinal : ndarray Épaisseurs rétournées par l'inversion ProfFinal : ndarray Profondeurs rétournées par l'inversion """ gg = Gate.shape[0] Npar = Res.shape[0] nData = Data0.shape[0] ModeleFinal = np.zeros((Npar, 3)) #%% Partie 1 : Création du fichier de commande .cfl BeowulfCfl = open('beowulf.cfl', 'w') BeowulfCfl.write('Inversion NanoTEM\n') BeowulfCfl.write('1 0 0 ! TDFD, ISYS, ISTOP \n') BeowulfCfl.write( '0 4 31 2 3.907800 3.904700 ! STEP, NSX, NCHNL, KRXW, REFTYM, OFFTIME\n') BeowulfCfl.write('') #Forme du courant BeowulfCfl.write('{:.6f} {:.6f}\n'.format(Wave[0, 0], Wave[0, 1])) BeowulfCfl.write('{:.6f} {:.6f}\n'.format(Wave[1, 0], Wave[1, 1])) BeowulfCfl.write('{:.6f} {:.6f}\n'.format(Wave[2, 0], Wave[2, 1])) BeowulfCfl.write('{:.6f} {:.6f}'.format(Wave[3, 0], Wave[3, 1])) BeowulfCfl.write(' ! TXON, TXAMP(4)\n') #Time Gate for i in range(0, gg-1): BeowulfCfl.write('{:.6f} \n'.format(Gate[i, 0])) BeowulfCfl.write( '{:.6f} ! TMS(31, GATE center) \n'.format(Gate[-1, 0])) for i in range(0, gg-1): BeowulfCfl.write('{:.6f} \n'.format(Gate[i, 1])) BeowulfCfl.write( '{:.6f} ! WIDTH(31)\n'.format(Gate[-1, 1])) #Survey Type BeowulfCfl.write('1 ! SURVEY_TYPE\n') #Tx Position and Moment BeowulfCfl.write( '1 1 1 1 4 1 ! NLINES,MRXL,NTX,SOURCE_TYPE,MXVRTX, TXMNT\n') BeowulfCfl.write('4 0 ! N_VRTX(1) TxZ\n') BeowulfCfl.write('0 20\n') BeowulfCfl.write('20 20\n') BeowulfCfl.write('20 0\n') BeowulfCfl.write('0 0 !SXE(I,J), SXN(I,J)\n') #Rx Position and Moment BeowulfCfl.write( '100 1 1 1 11 ! LINE, IDTX, RX_TYPE,NRX, UNITS \n') BeowulfCfl.write( '3 0 0 3 0 250 ! CMP, SV_AZM, KNORM, IPLT, IDH, RXMNT\n') BeowulfCfl.write('10 10 0 ! RXE, RXN RXZ\n') #Layer Earth Modele BeowulfCfl.write('{} {}'.format(Npar, Npar)) BeowulfCfl.write(' !NPAR, NLITH\n') for i in range(0, Npar): BeowulfCfl.write('{} 1 1 0 0 1'.format(Res[i])) BeowulfCfl.write( ' ! RES,RMU, REPS, CHRG, CTAU, CFREQ(1) - Layer\n') for i in range(0, Npar-1): BeowulfCfl.write('{} {} '.format(i+1, Ep[i])) BeowulfCfl.write(' !LITH(i), THK(i)\n') BeowulfCfl.write('{} !LITH(i), THK(i)\n'.format(Npar)) #Inversion Control Parameter if flag == 0: BeowulfCfl.write( '0 90 1 2 ! NFIX, MAXITS, CNVRG, INVPRT\n') elif flag == 1: BeowulfCfl.write( '{} 90 1 2 ! NFIX, MAXITS, CNVRG, INVPRT\n'.format(Npar-1)) for i in range(Npar-1): BeowulfCfl.write( '1 {} 2 ! CTYPE, LYR_INDX, KPAR\n'.format(i+1)) BeowulfCfl.close() #Partie 2 : Création du fichier .inv et inversion BeowulfInv = open('beowulf.inv','w') BeowulfInv.write('0 ! FD_ORDER\n') BeowulfInv.write('100 1 3 ! LINE_CHK, NSTAT, KMP(J)\n') BeowulfInv.write('{:.6E} ! Data floor\n'.format(Data0[-1])) BeowulfInv.write('1 ') for i in range(nData-1): BeowulfInv.write(' {:.6E} '.format(Data0[i])) BeowulfInv.write('{:.6E}\n'.format(Data0[-1])) BeowulfInv.close() #Partie 3 : Inversion et lecture du résultat subprocess.call('Beowulf') ResultOut = open('Beowulf.out', 'r') for line in ResultOut: if 'Final Model After' in line: Rms = float(line.split()[8]) for i in range(3): ResultOut.readline() for j in range(Npar-1): ResultLine = ResultOut.readline() ModeleFinal[j, :] = np.array( [float(i) for i in ResultLine.split()[1:4]]) ResultLine = ResultOut.readline() ModeleFinal[Npar-1, 0] = np.array([float(ResultLine.split()[1])]) for i in range(9): ResultOut.readline() ResultLine = ResultOut.readline() Data1 = np.array([float(i) for i in ResultLine.split()[2:]]) ResFinal = ModeleFinal[:, 0] ProfFinal = ModeleFinal[:-1, 1] EpFinal = ModeleFinal[:-1, 2] return(Data1, ResFinal, ProfFinal, EpFinal, Rms)