Amin Khalaf
email: amin.khalaf@upmc.fr

Christian Camerlynck

Nicolas Florsch

Ana Schneider

Development of an adaptive multi-method algorithm for automatic picking of first arrival times: application to near surface seismic data

Keywords: automatic picking, adaptive algorithm, multi-nested window, kurtosis, Akaike Information Criterion, seismic tomography, look-like Monte Carlo

come

Development of an adaptive multi-method algorithm for

automatic picking of first arrival times: application to near surface seismic data

Information Criterion. They benefit from combining different properties of the signal in order to highlight first arrivals and finally provide an efficient and robust automatic picking. This strategy mimics the human first-break picking, where a global trend is first defined at the beginning of the picking procedure. The exact first-breaks are then sought in the vicinity of each point suggested by this trend. Three successive phases are combined in a multistage algorithm, each of them characterizing a specific signal property. Within each phase, the potential picks and their error range are automatically assessed and sequentially used as prior in the following phase picking. Since having realistic estimates of the error in picked traveltimes is crucial for seismic tomography, our adaptive algorithm automatically provides picked arrivaltimes with their associated uncertainties. We demonstrate the accuracy and robustness of the implemented algorithm with synthetic, pseudosynthetic, and real datasets presenting challenges for classical automatic pickers. The comparison of both manual and adaptive picking procedures demonstrates that our new scheme provides more reliable results even under different noisy conditions. All parameters of our multi-method algorithm are self-adaptive thanks to the sequential integration of each sub-algorithm results in the workflow.

INTRODUCTION

Seismic arrival-time picking is crucial for many seismological studies. Using picked arrivaltimes, we can locate earthquakes, understand the interior structure of the Earth, study the geotectonic context and evaluate seismic risk for a region [START_REF] Billings | Errors in hypocenter location: Picking, model, and magnitude dependence[END_REF][START_REF] Hafez | Earthquake onset detection using spectro-ratio on multithreshold time-frequency sub-band[END_REF][START_REF] Ross | Automatic picking of direct P, S seismic phases and fault zone head waves[END_REF][START_REF] Stewart | Real-time detection and location of local seismic events in central california[END_REF]. More specifically in reflection and refraction seismology, first-arrival times are picked to guess the near-surface structure [START_REF] Blias | Optimization approach to automatic first arrival picking for three-component threedimensional vertical seismic profiling data[END_REF][START_REF] Lanz | Refraction tomography over a buried waste disposal site[END_REF][START_REF] Lawton | Computation of refraction static corrections using first-break traveltime differences[END_REF]. The quality of inverted models (near-surface layers' velocity and thickness) therefore depends on the accuracy of the picked first-arrival times [START_REF] Senkaya | A semi-automatic approach to identify first arrival time: the Cross-Correlation Technique[END_REF] since a small error in times may result in a non-negligible error in velocity and depth calculations (Ramananantoandro and Bernitsas, 1987).

The first events can be classically defined as the first remarkable onset of seismic energy, i.e.

"the first noticeable departure of the seismic pulse from a background signal" [START_REF] Hatherly | A computer method for determining seismic first arrival times[END_REF]. However, they do manifest themselves as a noticeable variation in amplitude, frequency content, statistical properties, and/or polarization properties in time and/or frequency domain. The quality of these events depends on the source type and near-surface conditions [START_REF] Yilmaz | Seismic data analysis: Processing, Inversion, and Interpretation of Seismic Data[END_REF]. The ambient noise may over-shadow the early part of the seismic phase and delay its apparent occurrence.

Traditionally, the first arrival times are manually picked by human operators, but the huge volume of data (especially in active-source seismic applications) makes this visual handpicking tedious, time-consuming and thus inefficient. Furthermore, the subjectivity of human-eye picking [START_REF] Hatherly | A computer method for determining seismic first arrival times[END_REF][START_REF] Jiao | Detection of seismic refraction signals using a variance fractal dimension technique[END_REF] puts alongside with factors relative to the operator's physiological and/or psychological state make the results of manual picking results often biased or inconsistent [START_REF] Saragiotis | PAI-S/K: A robust automatic seismic P phase arrival identification scheme[END_REF][START_REF] Sabbione | Automatic first-breaks picking: New strategies and algorithms[END_REF].

On the other hand, automatic picking is a more objective tool since it allows rapid analysis of large data sets while providing homogeneous and consistent results. Computer-based picking methods can be divided into semi-automatic (e.g. [START_REF] Lou | AIMBAT: A Python/Matplotlib Tool for Measuring Teleseismic Arrival Times[END_REF][START_REF] Vandecar | Determination of teleseismic relative phase arrival times using multichannel cross-correlation and least squares[END_REF]) and full-automatic (e.g. [START_REF] Coppens | First arrival picking on common-offset trace collections for automatic estimation of static corrections[END_REF][START_REF] Küperkoch | Automated determination of P -phase arrival times at regional and local distances using higher order statistics[END_REF][START_REF] Lomax | Automatic Picker Developments and Optimization: FilterPicker--a Robust, Broadband Picker for Real-Time Seismic Monitoring and Earthquake Early Warning[END_REF] algorithms. These algorithms can be further classified into two main approaches: single-trace and multi-trace. Singletrace algorithms can be either used as a "detector" in order to recognize a seismic phase arrival, or as a "picker" in order to determine a precise onset time of a detected phase [START_REF] Allen | Automatic phase pickers: Their present use and future prospects[END_REF]. Multitrace algorithms can be seen as an extension of single-trace algorithm, working on several traces at once. They are usually used for picking passive and active-source seismic data, benefiting from the similarity of waveforms from nearby events. Recently, [START_REF] Tselentis | Strategy for automated analysis of passive microseismic data based on S-transform, Otsu's thresholding, and higher order statistics[END_REF] provided an excellent classification of the common event-detection and phase-picking methodologies. Another exhaustive review with a new classification (i.e. single-level, hybrid and multi-level algorithms) was recently provided by Akram and Eaton (2016a) with discussion about the key parameters of each algorithm on downhole microseismic data. Among these algorithms, the most reported for single-trace (P-, Sphase, or both) picking are energy ratio criteria (STA/LTA) [START_REF] Allen | Automatic phase pickers: Their present use and future prospects[END_REF][START_REF] Allen | Automatic earthquake recognition and timing from single traces[END_REF][START_REF] Baer | An automatic phase picker for local and teleseismic events[END_REF][START_REF] Earle | Characterization of Global Seismograms Using an Automatic-Picking Algorithm[END_REF], autoregressive-Akaike information criteria (AIC) [START_REF] Maeda | A method for reading and checking phase times in auto-processing system of seismic wave data[END_REF][START_REF] Sleeman | Robust automatic P-phase picking: an on-line implementation in the analysis of broadband seismogram recordings[END_REF][START_REF] Takanami | A new efficient procedure for the estimation of onset times of seismic waves[END_REF], higher order statistics (HOS) [START_REF] Saragiotis | Automatic P Phase Picking Using Maximum Kurtosis and k-Statistics Criteria[END_REF][START_REF] Saragiotis | PAI-S/K: A robust automatic seismic P phase arrival identification scheme[END_REF], polarization information [START_REF] Cichowicz | An automatic S-phase picker[END_REF], and wavelet transform (WT) [START_REF] Anant | Wavelet transform methods for phase identification in threecomponent seismograms[END_REF][START_REF] Bogiatzis | Continuous Wavelet Decomposition Algorithms for Automatic Detection of Compressional-and Shear-Wave Arrival Times[END_REF][START_REF] Zhang | Automatic P -Wave Arrival Detection and Picking with Multiscale Wavelet Analysis for Single-Component Recordings[END_REF].

Although a multitude of algorithms exists, the aforementioned automatic methods are less commonly used in the literature for first-arrival picking on near-surface refraction records. [START_REF] Peraldi | Digital processing of refraction data study of first arrivals*[END_REF] proposed an algorithm based on cross correlation of the seismic trace with a reference "model" trace. However due to incorrect assumption about waveform stationarity, this algorithm does not pick the first arrivals accurately [START_REF] Gelchinsky | Automatic picking of first arrivals and parameterization of traveltime curves*[END_REF]. [START_REF] Hatherly | A computer method for determining seismic first arrival times[END_REF] proposed an algorithm combining the correlation method with a linear least-squares prediction method. A drawback of this algorithm is that it is non-robust against random noise.

Gelchinsky and Shtivelman (1983) used a method based on spatial correlation properties of the refracted signals. To prevent an apparent cycle skipping problem with the cross correlation approach, [START_REF] Zelt | Picking noisy refraction data using semblance supplemented by a Monte Carlo procedure and spectral balancing[END_REF] implemented a Monte Carlo technique in which the cross correlation function is transformed into a probability distribution using semblance optimization in order to pick low amplitude arrivals in crustal refraction data. Like other automatic picking methods, the accuracy of cross-correlation based methods decreases with the signal-to-noise ratio (SNR). However, the accuracy of these algorithms may be improved when combined with (i) semi-automatic processes in which human operators can interactively guide the algorithm to the correct pick (e.g. [START_REF] Senkaya | A semi-automatic approach to identify first arrival time: the Cross-Correlation Technique[END_REF], (ii) iterative stacking/filtering procedures that leads to improve SNR signals (e.g. Akram and Eaton, 2016b;[START_REF] Iqbal | Iterative interferometry-based method for picking microseismic events[END_REF], (iii) or both (e.g. [START_REF] Lou | AIMBAT: A Python/Matplotlib Tool for Measuring Teleseismic Arrival Times[END_REF][START_REF] Pavlis | Array processing of teleseismic body waves with the USArray[END_REF].

Data pre-processing can also improve the detectability of first arrivals and increase the accuracy of automatic picking [START_REF] Mousa | A new technique for first-arrival picking of refracted seismic data based on digital image segmentation[END_REF][START_REF] Mousa | Enhancement of first arrivals using the τ-p transform on energyratio seismic shot records[END_REF]. In strong heterogeneous velocity fields, the first arrivals do not necessarily correspond to the most energetic wavefield [START_REF] Geoltrain | Can we image complex structures with first-arrival traveltime?[END_REF], especially in the case of refracted arrivals at far offsets. Therefore, seismic noise can considerably impede automatic picking of weak [START_REF] Tselentis | Strategy for automated analysis of passive microseismic data based on S-transform, Otsu's thresholding, and higher order statistics[END_REF]. One way to address this problem is to filter the data before automatic picking. But when the first arrivals and noise share the same frequency band, data filtering becomes ineffective for automatic picking. [START_REF] Allam | Seismic Imaging of a Bimaterial Interface Along the Hayward Fault, CA, with Fault Zone Head Waves and Direct P Arrivals[END_REF] showed the potential effects of causal and non-causal filters on head wave picking quality. Non-causal filters can introduce small artifacts ringing to the front of high-amplitude arrivals, whereas causal filters can produce a phase shift. Additionally, [START_REF] Aldersons | Toward a Three-Dimensional Crustal Structure of the Dead Sea region from Local Earthquake Tomography[END_REF] and [START_REF] Senkaya | A semi-automatic approach to identify first arrival time: the Cross-Correlation Technique[END_REF] discussed the troublesome effects inherent to the choice of corner frequencies and slopes of a filter .

The combination of multi-methods was also proposed to increase the accuracy of automatic techniques when picking more complicated phases (e.g. S-phase) [START_REF] Ross | Automatic picking of direct P, S seismic phases and fault zone head waves[END_REF][START_REF] Diehl | Automatic S-wave picker for local earthquake tomography[END_REF][START_REF] Küperkoch | Automated determination of Sphase arrival times using autoregressive prediction: application to local and regional distances[END_REF]. Since each method is based on identifying a specific property of the part of the seismic trace where the first arrival occurs, the hybrid algorithms combine the advantageous properties of each individual method.

In near-surface seismic tomography, picking far offset traces remains a challenge. When picked, the overall bad quality of far offsets traveltimes usually degrades the inversion reliability though they are required to increase the investigation depth. To mitigate such effects, first arrivaltimes can be used with their associated uncertainty estimate in order to balance the trade-off between investigation depth and inversion quality. According to [START_REF] Aldersons | Toward a Three-Dimensional Crustal Structure of the Dead Sea region from Local Earthquake Tomography[END_REF], estimating picking uncertainties is the second most important goal in an automatic picking algorithm. Di [START_REF] Di Stefano | Automatic seismic phase picking and consistent observation error assessment: application to the Italian seismicity[END_REF] automatically estimated the uncertainties of first arrivals by using quality weighting algorithm which takes into account the waveform sampling rate, spectral density analysis and SNR. [START_REF] Diehl | Automatic S-wave picker for local earthquake tomography[END_REF] obtained a robust estimate of error intervals for automatically picked S-arrivals by combining picking information from different methods (STA/LTA detector, polarization detector, and autoregressive picker). [START_REF] Küperkoch | Automated determination of P -phase arrival times at regional and local distances using higher order statistics[END_REF] introduced an automatic quality estimation of picked P-onset from local slope and SNR of the characteristic function (CF) used in their picking procedure.

Although using combined algorithms has become a standard approach for automatic picking in deep seismology, it is still not common for near surface seismic data.

Here we propose a new adaptive algorithm based on the combination of three picking methods, respectively Multi-Nested Windows (MNW), Higher Order Statistics (HOS), and Akaike Information Criterion (AIC). This work focuses on overcoming most of the shortcomings inherent to the use of a single picking method. It also aims at providing accurate P-wave arrivaltimes along with their uncertainties, which are required information for reliable near surface tomography processing.

Exploiting all sub-algorithms in a single scheme proves to be beneficial in order to avoid filtering. The proposed new scheme is thus auto-adaptive and nearly parameter-free algorithm. Furthermore, all sub-algorithms are applied on raw single-component records in time domain in order to improve computational efficiency.

METHODOLOGY Principle and methods

Most of automatic picking techniques rely on different noise and signal properties in time or/and frequency domain. The transition point (noise to signal) can be identified with an algorithm that enhances and detects a specific attribute. Our proposed strategy is based upon an adaptive algorithm combining three picking methods (MNW, HOS, and AIC). The implemented algorithm mimics the human first-break picking, where at the beginning a global trend is defined, and then the precise first arrival times are searched in the vicinity of this trend. The picking is performed using a three stage algorithm that combines sequentially the three sub-algorithms. In this section, we describe the mathematical background of the three methods and their implementation in adaptive workflow either for a trace-by-trace or a shot-gather picking. To better illustrate and facilitate the comparison, each trace is normalized to its maximum amplitude before any process.

Multi-Nested Windows (MNW) method

The particles motion increases gradually when the disturbance of the seismic waves arrives, so usually no sudden take-off in the trace at the first break can be noticed [START_REF] Ricker | The form and laws of propagation of seismic wavelets[END_REF]. Therefore, the early advent of the first arrival (here, P-wave) is usual low, and the later one would be more energetic. However, in near surface prospection, first arrivals could be obscured by ambient noise. This case is commonly met at far offsets due to the spherical divergence of the front wave, the attenuation and the interference with background noise [START_REF] Mallinson | Enhanced refractor imaging by supervirtual interferometry[END_REF].

Let

represent the amplitudes of trace at a time t ( [ ] of samples). The average energy in the before, after and delayed windows (respectively BEA, AEA, and DEA) are defined by Eq. (1):

∫ ∫ ∫ (1)
where , and are the lengths (here in samples) of the corresponding windows in terms of dominant period length ( ) of the first arrival wave; d is the time delay for shifting the DEA window. [START_REF] Trnkoczy | Understanding and parameter setting of STA / LTA trigger algorithm 1 Introduction[END_REF] quantified the tradeoffs between the windows size and trigger sensitivity of the short-term average/long-term average (STA/LTA; [START_REF] Allen | Automatic earthquake recognition and timing from single traces[END_REF] algorithm. It is claimed that the preceding window should be a few times longer than the later window for microseismic data picking [START_REF] Chen | A multi-window algorithm for automatic picking of microseismic events on 3-C data[END_REF][START_REF] Han | Time picking and random noise reduction on microseismic data[END_REF][START_REF] Wong | Automatic time-picking of first arrivals on noisy microseismic data[END_REF]. In this work, the lengths of the windows are . The energy ratios and are defined by Eq.

(2).

(

) 2
where β is a constant used to avoid numerical instability in the energy ratio computation. In this study, β was adopted as 0.005 being suitable for the datasets used. It is noticed that is slightly equivalent to energy ratio (ER=STA/LTA) methods (e.g. [START_REF] Han | Time picking and random noise reduction on microseismic data[END_REF]. However, when ER exceeds a pre-defined threshold, the first arrival is considered as detected. But the onset picking may need further processing. [START_REF] Wong | Automatic time-picking of first arrivals on noisy microseismic data[END_REF] showed that the onset of the first arrival is very close to the maximum derivative of ER. In the case of noisy traces, the first-break times derived from ER technique are somewhat later [START_REF] Gaci | The use of wavelet-based denoising techniques to enhance the first-arrival picking on seismic traces[END_REF]. To overcome the problems due to first arrival weakness, we propose a new characteristic function (CF mnw ) by summing the two energy ratios ( ), calculated in the nested windows shown in Fig. 1a, as follows.

(3)

Fig. 1b demonstrates the ability of incorporated energy ratios to enhance the presence of the first arrival onset, where it is shown that the new CF is greater in value and sharper than the traditional STA/LTA ratio at the first arrival time. Thereby, the accuracy in picking should be enhanced by comparing to traditional energy ratio methods.

Beside the drawbacks concerning windows size, threshold choice is another crucial problem of energy ratios methods (Akram and Eaton, 2012). We assume that the noise average does not depend on time, i.e. the seismic noise satisfies local stationary in short observation-scale [START_REF] Ait Laasri | Automatic detection and picking of Pwave arrival in locally stationary noise using cross-correlation[END_REF]. Consequently, when all windows are merely comprised of noise part, and

(given by Eq. (2) are circa one, resulting in CF mnw around two (Eq. ( 3)). Based on these considerations, the theoretical threshold is given by Eq. (4):

(4)

where is the standard deviation computed within the window preceding the instant t.

By adding term in Eq. (4), the impact of most noise fluctuations should be mitigated, making this threshold slightly time dependent. However, when the later and delayed windows encompass a part of the seismic phases, the CF mnw becomes greater than 2. The beginning of the potential zone (BPZ; P-wavetrain trend) is determined once the CF mnw exceeds Thr(t). In order to minimize spike fluctuations on the picking procedure, a smoothing process of CF mnw is performed via local regression tools (using for example loess MatLab function) within a window of half dominant period ( ).

Since CF mnw characterizes the waveform change of the P-wavetrain, the potential first arrival times can be allocated to the two first local maxima of the smoothed CF mnw (SCF mnw ) after BPZ within a predefined searching window. Therefore, one of the two initial picks that correspond to the maximum SNR represents the first break (tP 1 ) (Fig. 1). The approximated interval error (tE 1 ) can be assessed by where is the absolute difference between BPZ and the first local maximum, and is the absolute difference between the two maxima if they exist. The results of MNW method can be summarized by three parameters: the early beginning of the potential zone (BPZ), the potential first arrival time , and the approximated interval error .

Higher Order Statistics (HOS) method

The seismic noise follows closely a Gaussian probability distribution [START_REF] Ait Laasri | Automatic detection and picking of Pwave arrival in locally stationary noise using cross-correlation[END_REF], while at the opposite, the seismic signal is non-Gaussian [START_REF] Persson | Statistical tests for regional seismic phase characterizations[END_REF][START_REF] Saragiotis | PAI-S/K: A robust automatic seismic P phase arrival identification scheme[END_REF].

Higher statistics (as skewness and kurtosis which are normalized central moments of order three and four, respectively) are well suited to check the gaussianity of time series [START_REF] Giannakis | Time-domain tests for Gaussianity and time-reversibility[END_REF]. Moving skewness and kurtosis reflect in different degrees the changes in distribution shape around a chosen position of the time series [START_REF] Press | A computer algorithm for automatic picking of refraction first-arrival time[END_REF]. However, the ability of HOS to pick first arrivals was successfully tested by [START_REF] Saragiotis | PAI-S/K: A robust automatic seismic P phase arrival identification scheme[END_REF] for real seismic applications.

Furthermore, the superior efficiency of kurtosis over other metrics (such as skewness) to detect and pick seismic phases has been shown by many authors (e.g. [START_REF] Küperkoch | Automated determination of P -phase arrival times at regional and local distances using higher order statistics[END_REF][START_REF] Nippress | Optimized automatic pickers: application to the ANCORP data set[END_REF].

In this study, kurtosis-based attribute (CF k ) was used to characterize the amplitude distribution, where CF k is calculated using Eq. ( 5) with a n k -samples moving window over the seismic trace.

∑ ( ̅ ) (5)
where a i represents a vector of seismic amplitudes inside a window of size n k , ̅ and σ t are respectively the mean and the standard deviation over n k samples. The kurtosis-value for each window is ascribed to the end sample of the window (Fig. 2a). CF k increases when the moving window integrates the first part of the seismic phase, due to the transition from Gaussian to non-Gaussian behavior (Fig. 2b). In case of complicated or emergent onsets, the kurtosis increases progressively before reaching its maximum value. The very beginning of this increase matches well the human-operator picking [START_REF] Baillard | An Automatic Kurtosis-Based Pand S-Phase Picker Designed for Local Seismic Networks[END_REF]. The onset time of the considered phase can be determined using the maximum CF k derivative [START_REF] Langet | Continuous Kurtosis-Based Migration for Seismic Event Detection and Location, with Application to Piton de la Fournaise Volcano, La Reunion[END_REF][START_REF] Lois | A new automatic S-onset detection technique: Application in local earthquake data[END_REF]. To address this issue, we applied the successive transformations proposed by [START_REF] Baillard | An Automatic Kurtosis-Based Pand S-Phase Picker Designed for Local Seismic Networks[END_REF] to CF k followed by smoothing; in a similar fashion as in MNW method, we can also apply low pass filtering to CF k by using the same function and the smoothing parameters. Hence, the first arrival time (tP 2 ) is defined by the minimum of the smoothed Baillard's function (SCF Bk ) (Fig. 2b). We notice that when the nonimpulsive arrivals are submerged by the noise, the kurtosis local maxima identifying this arrivals shift up in time. The approximated error (tE 2 ) can therefore be estimated by the time interval between the maximum of the CF k function and the picked onset. For the clear and impulsive phases this time interval is very small, while in case of noisy and non-impulsive phases this interval becomes significantly larger. The results of HOS method can be summarized in two parameters: the potential first arrival time ( ) and the time sensibility range ( ).

Akaike Information Criterion (AIC) method AIC [START_REF] Akaike | Markovian representation of stochastic processes and its application to the analysis of autoregressive moving average processes[END_REF] is widely used in autoregressive (AR) techniques (e.g. [START_REF] Leonard | Multi-component autoregressive techniques for the analysis of seismograms[END_REF][START_REF] Leonard | Comparison of Manual and Automatic Onset Time Picking[END_REF][START_REF] Takanami | A new efficient procedure for the estimation of onset times of seismic waves[END_REF]. This approach is based on the assumption that the seismic trace can be split into locally stationary segments, where each of them is satisfying a different autoregressive model. Hence, the point where the AIC is minimized indicates the optimal separation of two time series (noise and signal plus noise) where the first break arrival [START_REF] Sleeman | Robust automatic P-phase picking: an on-line implementation in the analysis of broadband seismogram recordings[END_REF]. In this study, instead of using AR coefficients for calculating AIC function, the CF based AIC (CF aic ) is calculated via Maeda's version (1985). Let a(t) a seismic trace of length N, the CF aic is computed by Eq. ( 6) directly from amplitudes sequences splitting in two segments (Fig. 3a).

(

) 6
where and are the variance of a in the backward and forward windows.

Under ideal conditions, the minimum of the CF aic is searched for, and its corresponding time is interpreted as first arrival time. In some circumstances when noise and other later wavefield interfere, or when the seismic trace contains many significant features, AIC picker may fail to define the arrival as a global minimum. It is so necessary to apply AIC picker within a time limited window (searching window) containing only the useful phase [START_REF] Zhang | Automatic P -Wave Arrival Detection and Picking with Multiscale Wavelet Analysis for Single-Component Recordings[END_REF]. Taking into account these points for improving the accuracy of AIC picking, we implemented the multi-model inference approach [START_REF] Li | An improved automatic time-of-flight picker for medical ultrasound tomography[END_REF]. This approach incorporates the entire set of models in the searching window to compute a weighted average model, based on normalized Akaike weights as following:

-Define the minimum of the CF aic in the searching window (CF_aic min ).

-Calculate the differences (Δ i ) between CF aic values and CF_aic min :

(7)

-Compute the normalized Akaike weights :

⁄ ∑ ⁄ (8)
-Set all values of out the searching window to zero.

-The potential first break time is obtained by using the weighted average:

∑ (9)
where t i is the time corresponding to each sample in the time series.

For traces showing high SNR, the difference between the improved AIC picker based on weighted model averaging scheme and ordinary AIC picker based on the best model (corresponding to the minimum AIC value) is negligible (Fig. 3b). On the other hand, when the traces show low SNR (where the global minimum of AIC function is not very sharp), the new AIC picker is more accurate since it integrates the other candidate models based on its relative importance estimated by Akaike weights [START_REF] Li | An improved automatic time-of-flight picker for medical ultrasound tomography[END_REF].

The sharpness of the AIC minimum often reflects the onset clarity [START_REF] Diehl | Automatic S-wave picker for local earthquake tomography[END_REF][START_REF] Toomey | Tomographic imaging of the shallow crustal structure of the East Pacific Rise at 9°30′N[END_REF]. A prominent AIC minimum is typically produced by an impulsive event, whereas the broader local minima are related to non-impulsive or noisy events. We estimate the time uncertainty of AIC picker from normalized Akaike weights function within searching window. The lower and upper limit of error interval correspond with, respectively, to the first and last sample where is greater than an appropriate threshold. We found that a threshold of 10% of the maximum function gives satisfying estimate of the uncertainty interval. The results of the AIC method can be summarized also by two parameters: the potential first-break of P-wave ( ), and the error range ( ).

However, in free-noisy traces one method is enough to pick correctly the first arrivals.

Whereas in low SNR traces, each method picking has its own CF characterizing a specific property of the trace, their results may therefore be significantly different. Hence, using multi-method algorithm helps to fill gaps in each method maintaining a good performance, and further, to guide picking algorithm toward the interesting zone.

Algorithms implementation

The main idea of the implemented strategy is that multi-algorithm integration should enhance the picking capability of weak events and reduce picking errors whatever the noise level. In this three stage algorithm the information provided by each sub-algorithm is used either to orient the following phase picking (defining a target window), or to setup the essential parameters for the next phase, or both. The sequential integration of each sub-algorithm results in the implemented workflow provides auto adaptive parameter estimation. Therefore, one input parameter is needed only for the first picking algorithm (in MNW). Subsequently, all potential picks derived from subalgorithms are used to determine the first arrival times and their corresponding relative errors. In the following subsections we describe the workflow of the complete picking algorithm and its parameters for a signal trace (trace-by-trace algorithm) and for a set of traces (entire-shot gather algorithm).

Single-trace algorithm picking

We extracted a vertical component seismic trace from a simulated record by finite differences modeling in order to describe the internal stage operations of the picking algorithm (inputoutput parameters and its respective CFs). A significant noise level (random and bandlimited noise) is added ( Fig. 4a). This trace presents a double serious challenge to be picked automatically, since the first arrival has a weak onset and is completely embedded within added noise. The first arrival was picked manually on clean trace before adding the noise and is marked by a black vertical solid line ( Fig. 4a).

Phase 1 (MNW sub-algorithm): the main goal of this phase is to enhance the presence of weak event and to determine the very beginning of the considering seismic phase. The CF characterizes the energy variations as a function of time. The first period length of the P-wave (T d ) is the only input parameter that needs fine-tune by the user. However, it can be obtained directly from the input data trace. For this test trace, we estimated T d to 75 samples. The picking procedures, as described previously, are applied on the smoothed CF (SCF mnw ). The searching window (SW1) is located just after BPZ. We find a length equal to one and a half period (1.5T d ) is suited to pick the two local minima for which one of them is potentially the first arrival (tP 1 ).

Fig. 4b shows the SCF mnw and the related results.

Phase 2 (HOS sub-algorithm): the size of the moving window (N k ) used in calculating the kurtosis-based (CF K ) can be estimated automatically from the approximated sensitivity interval of the previous sub-algorithm (tE 1 ). For most tested traces, we found that N k twice the (tE 1 ) provides satisfactory results. In order to avoid the under-or over-estimating of CF K , we added a security condition in which when estimated N k less than 0.5T d or greater than 2T d , the algorithm setup the N k to one period (T d ). The searching window (SW2 = tE 1 +T d ) is

automatically set up around the potential first arrival derived from MNW sub-algorithm (tP 1 ); tE 1 before tP 1 and T d after. The potential first arrival (tP 2 ) and associated time uncertainty (tE 2 ) are also determined according to the descriptions in the related section (HOS method).

The related CFs and results of this sub-algorithm are illustrated on the Phase 3 (AIC sub-algorithm): Since the performance of AIC methods depends strongly upon the searching window configuration, the results of previous sub-algorithms (MNW and HOS based kurtosis) were jointly used to lead the AIC picker, and further to setup its searching window size. We define the size of this window and its location (respectively, SW 3 and C aic ) according to Eq. (10).

(10)

Finally, the potential first arrival ( ) and its time uncertainty ( ) are also obtained according to the related explications in AIC method section. The related attribute functions and the results are illustrated on the where is a quality factor of the picked phase calculated into a predefined time windows and is given in dB [START_REF] Küperkoch | Automated determination of P -phase arrival times at regional and local distances using higher order statistics[END_REF]:

( ), (12)
where is the RMS-amplitude of the noise calculated within a time window of length (nois_wind = 3Td) before the determined . The is the RMS-amplitude of the P-wave evaluated on the signal window (sig_wind = 1Td) starting immediately after .

In general, Q-values can yield crude information on the quality pick. If the arrival is correctly picked, Q has high and positive value. On the other hand, Q-value close to zero or negative indicates, respectively, very early picking (noise signal or dead trace) or late picking. All picks corresponding with negative Q values were excluded out the final picking process.

Error and quality assessment

While an ultimate goal of any automatic picker is to assure the highest accuracy possible, an associated uncertainty estimation is of great importance in a posteriori processing of first arrival time (e.g. in tomographic inversion). For these purposes, we suggest an approximated estimation of final pick error by using a simple statistical tool. From a practical and physical point view, the standard deviation can be considered as a good error estimator of repeated measurements of a physical parameter (here, the arrival time). Therefore, the final time uncertainty ( ) is produced by the standard deviation of three potential picks resulting from the aforementioned subalgorithms.

(

) 13
Since the first arrival onset is completely embedded by the random noise, the final pick is slightly shifted ( Fig. 4e). Therefore, that demonstrates the high sensibility of the three picking methods to the higher amplitudes. Although the considerable added noise, the correct first arrival time is still in the range error of the implement algorithm ( Fig. 4e).

Final picking reliability of single-trace algorithm:

To analyze the reliability of the final picking in the single-trace algorithm, we use a pseudosynthetic dataset. These dataset was generated by adding 1000 noise scenarios to a real trace characterized by high SNR. Since the first-order Gauss-Markov process is the best candidate that can model the seismic noise [START_REF] Baziw | Microseismic Event Detection Kalman Filter: Derivation of the Noise Covariance Matrix and Automated First Break Determination for Accurate Source Location Estimation[END_REF], we use this process to create the seed of our synthetic noise realizations as following:

generation of a long time series (the seed series) of random noise via Gauss-Markov process.

apply S-transform [START_REF] Stockwell | Localization of the complex spectrum: the S transform[END_REF] on this one-dimensional noise. As expected, the algorithm performance deteriorates with decreasing SNR. The quality of final picking and/or the normalized final error could be good criteria for assessing picked results. For this dataset, we can easily define 5 dB as a threshold for primary filtering out the bad picks with no good picks below this threshold. However, that is not evident for bad picks corresponds to earlier or later high amplitude produced by addition noise or posterior arrivals, respectively (Fig 6d).

To clarify the MNW effect on the final picking we perform a cross-analysis of MNW picking results. We define the acceptable error for MNW as ±T d . This error-value seems a reasonable choice because it corresponds to a part of the searching window of HOS that follows the initial pick estimates from MNW. A drawback of the three-stage algorithm occurs when the MNW stage fails to successfully define the BPZ, leading to sequential failure as sub-algorithms are inter-connected by auto adaptive parameters transmission, and to false picking. In such case, high-fidelity filtering (e.g. [START_REF] Aldersons | Toward a Three-Dimensional Crustal Structure of the Dead Sea region from Local Earthquake Tomography[END_REF] or denoising (e.g. [START_REF] Parolai | Denoising of Seismograms Using the S Transform[END_REF][START_REF] Tselentis | Strategy for automated analysis of passive microseismic data based on S-transform, Otsu's thresholding, and higher order statistics[END_REF][START_REF] Gaci | The use of wavelet-based denoising techniques to enhance the first-arrival picking on seismic traces[END_REF] are still possible during preprocessing. In near surface seismic surveys, exploitation of multi-trace redundant data and events coherency is another candidate solution.

Entire-shot record algorithm picking

The basic idea of this algorithm is to integrate all picks provided from preceding subalgorithms by exploitation the properties of coherency of the first arrivals and its consistency over all offsets. This algorithm is an extension of the single-trace algorithm, and it is adapted for full automatically picking of common shot-gathers. In order to better determine the trend of Pwavetrain and to work under different noise conditions we optimized MNW sub-algorithm by adding/involving two tools; variable-threshold engine and look-like Monte Carlo approach.

Using a fixed MNW-threshold might lead to false picking due to the high variability of the noise and the P-waveforms. In order to elevate this drawback, we proposed a variable-threshold engine combined with an iterative process. The latter is a look-like Monte Carlo approach; an optimal solution is built from the database picks that arise from MNW sub-algorithm via variable thresholding. Here, the solution is the set of picks (BPZ picks) forming the first-breaks trend over the whole shot gather. Using two subroutines: In the first one, a matrix of all possible solutions is generated by a multi threshold engine picking (the threshold varies between 0.5 and 10). The second one is an iterative system consisting of: (1) randomly building a solution, (2) calculating the corresponding cost function, (3) repeating the process for a predetermined number of iterations, (4) choosing the best solution which corresponds to maximum of the cost function. Our cost function combines three main properties in order to highlight the presence of the first arrivals:

Energy cost function ( ): we defined (Eq. ( 14)) as a function of the afterenergy ratio (AER) modulated by given picks quality (Q) and normalized by an error ( ). AER is CF mnw average within a limited window just after the built solution. Q is estimated by Eq. ( 12).

is given by the standard deviation of the solutions matrix over threshold dimension ( ). This cost function has to maximize over the whole shot record.

∑ ( )

Smooth cost function ( ): by definition, a curve is smooth when all the derivatives exist and are continuous. The simplicity and efficiency of centered second-order derivative ( ) make it an attractive candidate for analyzing smoothness curve. The second order derivative of the travel time shape has to minimize (i.e. the solution is hold as close as possible to a smooth curve).

The shape of the trend is normalized by an attributed error ( ), which can be defined as the standard deviation of the selected solution ( ). The is then given through the Eq. (15):

∑ ( | | ) (15)

SNR cost function (

): the is defined by Eq. ( 16), and has to be maximized.

∑ ( )

where is calculated for raw traces in predetermined windows around the built solution.

The total cost function is the sum of the normalized costs functions as follows:

(17)
where , β and γ are weights in [0, 1], which allow controlling the importance of individual objective functions. In this study, we give an equal contribution of the three cost functions ( ). The trend of the first arrivals is then to find the optimal solution ( ) which maximize the total cost function ( ). The problem can be simply formulated as,

The general flow of our iterative system is somewhat similar to Monte Carlo approach (randomizing selection and then minimizing a likelihood function). In our approach, the random selection is limited or somehow guided, and the likelihood function is absent (i.e. full directly searching).

Taking in consideration that the first arrivals are consistent and coherent over all offsets, their trend should be as smoothed as possible. In order to eliminate outliers and high first order derivative points, the final solution should be passed through a smoothing process. We used rloess matlab function, which is based on local regression using linear least squares; it assigns lower weights to outliers during the regression. The smoothing window length is chosen depending on the complexity of subsurface and therefore on the expected form of travel time curve (in percent of total traces number).

In very noisy case, the correct first arrivals could be found in the shadow of strong ambient noise. Hence, they might be invisible for the variable-threshold engine. In such case, it is preferred to rebuild the solutions matrix, but within a small tolerance window centered on the smoothed solution. The tolerance window length of two time dominant period (2T d ) is suitable for this purpose.

The new candidate matrix is processed again though the iterative system. The final smoothed solution can then be interpreted as the trend of the potential zone of the P-wavetrain arrivals.

As a whole, the rest of the algorithm has the same workflow as the single-trace algorithm. In particular, the input parameters of each sub-algorithm (especially the searching windows size) were adjusted by using statistical tools. The abundance of each sub-algorithm outputs permits also to set forwardly the SW location for the next sub-algorithm (i.e. guiding mode).

Since the trend might be delayed due to use of multi-threshold and/or smoothing process in MWN sub-algorithm, it is appropriate to add a confidence interval (0.5 T d ) before the estimating smoothed trend. Therefore, the P-wavetrain zone is located around the smoothed trend rather than just after it (as in single-trace algorithm).

Table 2. Since FD modelling requires to grid the model into a finite number of cells [START_REF] Carcione | Seismic modeling[END_REF], the model (500 x 150 m) was discretized into 1000 x 300 grid mesh (i.e. the space is equally sampled in x and z at 0.5 m). The Perfected Matched Layer (PML) absorbing boundary with 120 grid points was implemented along at three sides, with free surface. We generate many inline shot gathers; in each one, a total of 0.6 s with time increment of 0.5 ms was modelled for the vertical point force located at the free surface. Fig. 8b shows an example end-on spread shot-gather

(the source at 198 m) of the vertical component of the wavefield simulated, in which the theoretical P-wave arrival times are denoted in small green circles.

We added the noise adding to the simulated wavefields, which represents another difficulty inherent to weak events, along complex changes in the waveforms, and nonlinear travel-time curve in presence of a fault. Two types of added noise (random zero-mean Gaussian and band-limited noises) were computed using the function suaddnoise of the Seismic Unix (SU) package [START_REF] Cohen | CWP/SU: Seismic Un*x Release No. 43R1: an open source software package for seismic research and processing[END_REF]. The band-limited noise has approximately the double frequency range of the simulated source.

We applied the entire-shot picking algorithm to pick the first arrival times in the shot gather shown in picking performed satisfactorily, an over-estimating of the arrival times is noticed after adding the noise. This shift-down of the picks increases with the offset and it is remarkable at the events corresponding to the fault.

To assess our algorithm robustness, we compare the picks with those determined theoretically for two cases: noise-free and noisy records. In the case of noise-free shot gather, 60%

of the automatic picks are within [0,-2] ms of the actual picks, and 90% are within [0, -10] ms (Fig.

9a, red histogram

). This analysis can be considered as a calibration step of the algorithm parameter setting; the median value of the automatic-manual differences is about -1 ms. A possible explanation of this deviation is that the P-wave dominant period length varies with the offset. For noisy data (random and band-limited noises), the blue histogram in Fig. 9a shows that about 35% of the differences are found at 10 ms, and 75% are within [0, 20] ms. Fig. 9b shows a statistical analysis per quality class, which are determined in analogy to the empirical uncertainty categories (A-D) of [START_REF] Lois | A new automatic S-onset detection technique: Application in local earthquake data[END_REF]. Clearly, all distributions of the four classes are found on the positive side of the difference between automatic and actual picks. The visual check of individual sub-algorithm results and corresponding waveforms shows that delay time picking is mainly due to their high sensibility to the high amplitudes in the presence of high noise levels. The dispersion symmetry around the median and the absence of outliers show again good consistence of the two picking methods.

With noise addition, the first arrivals quality was degraded by about 20% of their level at the noise-free gather (Fig. 10; see red dotted curves evaluation comparing to the red solid ones). The accuracy of our adaptive algorithm versus the first arrivals quality is presented in the tables included in the Fig. 10. The Root-Mean-Square (RMS) of the differences in case of zero noise is 3.6 ms, and the maximum RMS in noisy cases is 16 ms, which normally corresponds with C and D quality classes.

Generally, increased noise levels result in anomalously late picks, because parts or even complete tends to become larger at far-offsets. Following the initial inspections on individual shot gathers, six periods ranging from 11 to 35 samples in length were examined. This key parameter was calibrated through a trial-and-error procedure in order to unify its value over all offsets. We can see that a value of 25 samples is a reasonable choice because it provides the best compromise between minimum RMS and maximum traces picked within the predefined margin error of manual picking (±5 ms) (Fig. 11). This error-value corresponds to approximately a quarter period of the 40-50 Hz first arrival signals [START_REF] Zelt | 3D seismic refraction traveltime tomography at a groundwater contamination site[END_REF]. The other required parameters are set up automatically during the process (i.e. auto adaptive parameters) as described in the previous sections.

In order to reduce calculation time, only the first 250 ms of the time record were taken in consideration. Fig. 12 shows an example of the picking results by our adaptive algorithm (blue line)

with related error, and those picked manually (green line). Unlike the far-offset traces, an excellent agreement is observed between manual and automatic picks of near-offset traces.

Applying the proposed algorithm to the whole profile data (more than 12500 traces), more than 88% of traces were picked in the implicit presumed error interval of the manual picks (±5 ms), with only 7 ms as RMS. A more detailed view of the statistical analysis of the automatic picking results is shown in Fig. 13. The residuals between our automatic picks and manual picks (automaticmanual) are shown as a function of offset (Fig. 13a) and first arrivals quality (SNR) (Fig. 13b), as well as their distribution in a histogram (Fig. 13c). The residuals decrease quickly and concentrate (close to zero) for the high quality traces, which largely reflect the high performance of the proposed algorithm. Unfortunately, the residuals depart from zero after offset 200 m and are dispersed in the positive side, which suggests that our algorithm picks systematically late arrivals at far-offsets.

Another possible explanation of this deviation is that the increasing dominant period length with offset might have a significant effect on algorithm accuracy. The signal-to-noise ratio (SNR) is however the main source of error for even a manual picking [START_REF] Zeiler | Seismogram Picking Error from Analyst Review (SPEAR): Single-Analyst and Institution Analysis[END_REF]. Fig. 13d shows the normalized cumulative histogram of the related picking errors, with 90% of the available traces picked with standard errors less than ±3 ms.

However, the tendency of late picking at those offsets reinforces the need to enhance these arrivals before automatic picking, or at least to omit automatically or manually their picks before any posterior processing. Additionally, there might be some dead traces where their picks are naturally wrong and must be also ignored [START_REF] Sabbione | Automatic first-breaks picking: New strategies and algorithms[END_REF].

A further advantage of the proposed algorithm is that the individual error derived automatically can be employed to define the pick quality, thus weighting the data in travel time inversion (e.g. [START_REF] Bauer | Tomographic P wave velocity and vertical velocity gradient structure across the geothermal site Groß Schönebeck (NE German Basin): Relationship to lithology, salt tectonics, and thermal regime[END_REF]. Furthermore, the performance of the developed algorithm versus other automatic pickers and with many datasets recorded by various sources-types was stated by [START_REF] Khalaf | Développement d'une nouvelle technique de pointé automatique pour les données de sismique réfraction[END_REF].

From our experience, the new algorithm processes 72 traces shot gather in about 10 seconds along with thorough quality check. Since this study is primarily focused on accurate picking, no attention is paid to time calculation improvement that may be a subject of further research.

CONCLUSIONS

We have introduced an adaptive multi-method algorithm to pick automatically the first arrival times of near surface seismic data by combining three methods (MNW, HOS, and AIC). The MNW method was combined with iterative look-like Monte Carlo approach supplying a good localization of the first arrivals trend, even if the arrival onset embedded by a high noise level.

Integrating the three methods results as a combination of the first arrival properties (energy, gaussianity, and stationarity) yielded fruitfully as a precise, nearly parameter-free, and full automatic picking algorithm.

Testing our algorithm on synthetic and real data demonstrated good results in comparison to those calculated theoretically and picked manually, respectively. The synthetic example presented in this paper is realistic for near-surface studies in terms of the velocities distribution, depth investigation, source-receiver geometries, and noise data. The noise robustness tests (with no added noise and with additive random Gaussian and limited-band noises) demonstrated the overall reliability of proposed algorithm, resulting in maximum RMS of 16 ms in the case of exaggerated noise. For real data, the effectiveness of the proposed picking method is highly dependent on the quality of the first arrivals and their dominant period, which is the key parameter in the overall methodology.

The proposed algorithm provides automatically not only the first arrival times, but also their uncertainties. These uncertainties can be used in a posteriori processing such as in tomographic inversion or data phase-editing prior inversion.

The number of user-tuning parameters was kept as low as possible as in order to reduce parameters errors. In fact, the main period of the first arrivals is the only required parameter. All the others parameters are auto-estimated and adapted across the processing workflow.

Finally, in cases of very weak events with high noise presence, local enhancement of the first arrival signals may help to improve our algorithm accuracy. For instance, like this enhancement may be achieved by local stacking or very high fidelity denoising. Our algorithm is appropriate for common shot gather, but it is expected to provide a basis for more effective picking on other seismic data collections, such as common-offset gather (COG). However, their detailed advantages and drawbacks analysis is beyond the scope of the present paper.

Fig. 4c .

 4c Fig. 4c. It is worth mentioning that it is sufficient to compute the kurtosis attribute only for

Fig. 4d .

 4d Fig. 4d.

-

  for each scenario, random selection (the bandlimited frequency and the start time) on the S-Transform time-frequency map , zeroing out of the selected zone and inverse S-Transform, normalization the output time series by its maximum and multiplication it by a random value (0.1 -0.9). The resultant noise scenarios can be considered more realistic in terms of random and bandlimited noises and S/N values. Fig 5 illustrates the pseudo-synthetic dataset generation via one noisy trace realization. The only input parameter (T d ) of the single-trace algorithm is about 30 ms. By assuming that the acceptable error is a quarter of the dominant period (T d /4 = 7.5 ms), 70% of final picks falls in the interval error [-7.5, 7.5] with RMS value equal to 3.66 ms. A detailed view of the statistical analysis of the automatic picking results is shown in Fig 6. The time differences (automatic -manual) and final error normalized by T d are shown as a function of first arrival quality (SNR) (Fig 6a), as well as their histogram distribution (Fig 6b).

  Fig 6c shows time difference pie chart for the four possible crosscombinations of MNW and Final Picking (FP) results. Among 70% of successful FP in its interval error, only 2% of picks were not initialized by MNW. Otherwise, despite than 9% of picks were well initialized by the MNW algorithm, the FP could not pick them correctly due to the relative highquality weighting of the other pickers (HOS and AIC). Most of these picks are often associated with later arrivals. The MNW could not determine the potential zone of the last 21%. Since energy ratio algorithms are sensitive to high noise fluctuations (Akram and Eaton, 2016a), MNW was probably trapped by high-amplitude noise fluctuations far from the first arrival. The same behavior is observed for FP. The corresponding picks are characterized by low final quality picking (Fig 6d). Normalizing the pick quality by the corresponding errors issued from final picking can be considered as a good criterion to evaluate the three stage algorithm consistency. High values indicate a good arrival quality, and the sub-algorithms will therefore pick in a narrow time interval (small error), whereas lower values (close to 0) correspond to erroneous first arrival (noise) or/and dispersive picks algorithms (large error). Fig 7 shows the consistence of three algorithms on the categories in the Fig 6.

  //fr.mathworks.com/matlabcentral/fileexchange/24531-accurate-fast-marching, last accessed March 1st, 2018). The automatic process was repeated for the same shot-gather in Fig.8abut with additive random Gaussian noise of 20%, with limited-band noise by [40-180] Hz, followed by the two previous noises together (Fig. 8b, Fig. 8c, and Fig. 8d), respectively. The picks resulting from the automatic algorithm have a good overall consistence with the theoretical picks. Although automatic

FiguresFig. 1 .Fig. 2 .Fig. 3 .Fig. 4 .Fig 5 .Fig 6 .Fig 7 .Fig. 8 .Fig. 9 .Fig. 10 .Fig. 11 .Fig. 12 .Fig. 13 .Fig. 14 .Table 1 :Fig. 1 .

 123456789101112131411 Figures Fig. 1. Illustration of MNW method: (a). normalized free-noisy trace and the three windows (red = before instant t; black = just after instant t; green = after a delay instant t) used on the calculation of CF mnw . (b). Energy ratios (ER1 = energy ratio calculated between the before window and after window; ER2 = calculated energy ratio between delayed window and before window; CF mnw = Characteristic Function of the incorporated energy ratio; CF mnw = Smoothed Characteristic Function; Thr = Threshold). The black line indicates the first break...................................................... 46

Fig. 2 .Fig. 3 .

 23 Fig. 2. Illustration of HOS based kurtosis method: (a). normalized free-noisy trace and windows (red rectangle) used for calculating CF k . (b). the resulting kurtosis CF (blew solid line), and transformed kurtosis based on Baillard procedure (red solid line). The black dashed line indicates the manual pick of first arrival of P-wave.

Fig. 4 .

 4 Fig. 4. Illustration of the Single-trace algorithm applied to a synthetic trace: (a). Noise-free synthetic trace (black), after adding noise (blue), manual first arrival (vertical black line). The first period of P-wave is framed by the red rectangle (b). MNW sub-algorithm topics: threshold (dashed red line), BPZ (vertical red line), picked potential first arrival and the associated error (respectively, vertical pointed and horizontal green line). (c). HOS sub-algorithm topics: the guide picker (vertical pointed green line), maximum of kurtosis (small green circle), picked potential first arrival and associated error (respectively, vertical pointed and horizontal red line). (d). AIC sub-algorithm topics: the picker guide (vertical pointed blue line), the red part of CF aic is considered for calculating CF wa , the picked potential first arrival and associated error (respectively, vertical and horizontal purple lines). The gray zones in b, c, and d represent the searching windows SW1, SW2, and SW3, respectively. (e). Zoom in of the trace in (a) around the first arrival (vertical solid line). The final pick and its error are marked by vertical pointed and horizontal solid lines, respectively. The potential first arrivals and their associated errors picked by the three sub-algorithms are redelineated on the noisy trace.

Fig. 8 .

 8 Fig. 8. A synthetic vertical components record with the first arrivals times obtained theoretically (in green) and automatically with their interval errors using the proposed adaptive algorithm (in blue and the error barre in red): (a).Noise-free shot-gather. (b). With 20% random white noise. (c). With random band-limited noise [40-180] Hz. (c). With the noise in (a) and in (b) (random and band-limited noise).

Fig. 11 .

 11 Fig. 11. Trial-and-error results for determining the best first arrival period-length: The curves show values expressed as RMS between automatic and manual picks (left axis) and percentages of hit-rate picking in the manual range-error (right axis).

Fig. 12 .

 12 Fig. 12. Shot-gather shown as an example of automatic picking result.

ACKNOWLEDGEMENTS AND DATA

The synthetic data will be available for download upon request. The authors are grateful to Dr. Antoine Chevalier for a number of helpful comments and suggestions concerning to Monte Carlo approach, and to Dr. Sylvain Pasquet, Dr. Mustafa Senkaya and another anonymous reviewer for their very constructive comments, which helped to significantly improve the present manuscript.

The parameters setting of all sub-algorithms are recapitulated in the Table 1 that presents a comparative summary of the two algorithms' parameters.

VALIDATION THE PROPOSED ALGORITHM Test on synthetic data

Synthetic datasets were produced by finite-difference (FD) modelling in the time-space domain using FDELMODC program developed by [START_REF] Thorbecke | Finite-difference modeling experiments for seismic interferometry[END_REF]. The geometry of the acquisition system consists of an array of 72 receivers with a spacing of 4 m in flat topography.

The seismic sources are simulated by an explosive Ricker wavelet with maximum-frequency of 90 Hz, spreading every 4 m between the receivers. The geological model used to create wavefield data is assumed to be an elastic media incorporating realistic subsurface structures. This model is composed of many layers (C1-C3) and a substratum (S) fractured by a normal fault (F); three layers at the left side of the fault and two layers on the right one (Fig. 8a). All velocities and densities are assumed to be constant in each layer, except the fault zone with a linear velocity gradient about 12.5 m/s per m, with v 0 = 625 m/s. Densities ( ) are set using the Grardner's rule for P-wave velocities

) for each layer. The details about the model geometrical parameters are shown in Fig. 8a and phases of the recorded signal can be masked by the noise. The implemented algorithm kept a good performance despite the decreasing of quality arrival due to additive noise, where the maximum RMS in the very noisy case did not exceed the acceptable error.

Test on real data

Exploration-scale active-source seismic dataset was chosen to validate the implemented algorithm on real data. The dataset is characterized by a variable quality level of first arrivals when increasing source-receiver distances. The performance assessment is achieved by comparing the algorithm result with those manually picked. The Douglas's concept (1997) "the less filtering, the better picking algorithm" has been applied, where no data preprocessing was performed before the automatic picking process, but all traces are scaled with respect to their maximum value.

Field data

As a test, we have applied our algorithm to a 800-m-long seismic profile conducted in the Chaux Forest in the Jura region (France). An accelerated 40 kg weight (hitch-mounted at the front of a tractor) dropped onto a steel plate, is used as a seismic source for generating the energy, with 4 m source interval along the profile. Four shots at each source location were averaged to increase the signal-to-noise ratio. The recorder system consists of 72 14-Hz vertical component geophones, with 4 m receiver spacing along 318 m of maximum offset. The recording time was 1 s with a delay of -20 ms. The data were sampled at 1 ms. A roll-along technique acquisition was used to perform the total profile survey. In this study more than 12500 recorded traces were considered from 174 shots.

The dominant period (T d ) of the P-waves can be directly estimated by visual inspection of the seismic section. Since the earth acts as a minimum phase low pass filter [START_REF] Aki | Quantitative seismology[END_REF], T d

List of captions

Tables Table 1 :


Yaoyao Ning

Xiaoqing Wang

Kangjia Sheng

Lili Yang

Wei Han

Chaoni Xiao

Jianli Li

Yongmin Zhang

Shaoping Wu

A novel colorimetric and fluorescent turn-on pH sensor with a notably large Stokes shift for its application

A novel naked-eye colorimetric and fluorescent turn-on pH sensor based on naphthalenone scaffold was rationally designed and facilely synthesized. The probe exhibited strong alkaline pH-dependent (pKa=9.92) behavior, high selectivity, sensitivity, high quantum efficiency (Ф= 0.64) and rapid respond to pH fluctuations (within 1 min) in the range of 9.0 -14.0. In addition, probe DDTM displayed a notably large Stokes shift of 198 nm. The response mechanism of the fluorescent probes was proposed to be caused by OH --induced structure changes from the nonfluorescent DDTM form to the highly emissive addition product (more than 150-fold fluorescence enhancement) along with ICT changes. Live cell imaging data revealed that probe DDTM could selectively monitor pH changes with low cytotoxicity and cell membrane permeability. Therefore, this probe could be acted as an effective intracellular pH imaging sensor under extreme alkaline condition in the biomedical and biological fields.

Introduction

It is well known that the pH values are uniquely of high importance among the few chemical parameters. Specially, intracellular pH plays an important role in cell function and regulation such as cell volume regulation [ 1 ], vesicle trafficking, cellular metabolism [2], cell membrane polarity [3], cellular signalling, cell activation, growth, proliferation and apoptosis [4]. As an important indicator of cellular health, abnormal pH values of intracellular could lead to cellular functional disorders and serious diseases such as cancer [5,6], stroke [ 7 ], Alzheimer's disease [ 8 ], neurodegenerative disorders [9] and cardiovascular disease [START_REF] Vaughan-Jones | Intracellular pH regulation in heart[END_REF]. Therefore, it is extremely crucial to sense and monitor pH changes in living cells for investigating cellular functions and further providing insight into physiological and pathological processes.

A variety of analytical techniques have been reported to measure pH value including acid-base indicator titration [11], potentiometric titration [12], nuclear magnetic resonance [13], Raman sensor [14] and so on. Hence, it is essential for the health professional to have simple and effective methodologies for the rapid, selective, and sensitive detection of pH values in biological specimens. Fluorescent probes have been found widespread application for sensing and imaging pH values in living cells and tissues due to its high sensitivity, easy operation, real-time and noninvasive imaging properties [15]. However, a lot of reported pH probes still have some drawbacks such as low fluorescence quantum yields, poor photostability and narrow Stokes shifts. For the sake of overcoming these drawbacks, it is important to develop new longwave fluorescent probes with large Stokes shifts, and high turn-on fluorescence brightness.

Naphthalenone derivatives have been found to exhibit a variety of biological activities [16,17] and have also been surveyed for their suitability as dyes scaffold [ 18,19]. In addition, aromatic heterocyclic compounds containing a five-membered furan or thiophene unit have been extensively studied in recent years due to their impressive optical properties and excellent charge-transport properties [20,21]. Therefore, construction of naphthalenone derivatives containing heterocyclic unit will be vital in synthesis of novel organic materials with appealing photonic and electronic functionality.

As part of our ongoing program in the study of novel fluorescent dyes [22,23], in the current work, a novel pHresponse fluorescent probe DDTM was rationally designed and facilely synthesized by connecting a DTC fluorophore and a dicyanomethylidene group with double bond. In DDTM, the dicyanomethylidene group was used as an acceptor moiety while the amino group acts as an electron donor to form pushpoll system. This probe displayed that ICT was conjugated organic π-systems with acceptor and donor subunits. To our surprise, the dicyanomethylidene group of probe DDTM, one of cyanide receptors, has no reaction with cyanide ion except OH - [24]. Compared with other reported pH probes [25][26][27][28][29][30][31][32], probe DDTM exhibited strong pH-dependent behaviour, high selectivity, sensitivity and responded linearly and rapidly to respond pH fluctuations within the range of 9.0 -14.0. Specially, probe DDTM displayed a notably large Stokes shift of 198 nm and high fluorescence quantum yields. The response mechanism of the fluorescent probes could be caused by the OH --induced structure changes from the nonfluorescent DDTM form to the highly emissive addition product along with ICT changes. Furthermore, live cell imaging data revealed that probe DDTM could selectively monitor pH changes with low cytotoxicity. Therefore, this probe could be used to visualize extreme alkaline within the biomedical and biological fields (Scheme 1). Scheme 1. Strategies to response of pH probe.

Materials and methods

Materials

All chemical reagents were purchased from commercial suppliers and used without further purification. All organic solvents used for this work were of analytical grade also without further purification. All reaction processes were monitored through thin layer chromatography (TLC) which was performed by using Merck F254 gel-60 plates. Silica gel (200 mesh) was used as the solid phase for column chromatography by Qingdao Ocean Chemicals (China). All metal ions were prepared from corresponding salts. All samples were prepared at room temperature before UV-vis and fluorescence determination. The phosphate-buffered saline (PBS) was prepared from sodium dihydrogen phosphate and disodium hydrogen phosphate. Probe standard solution dissolved in DMSO and diluted to 50 mL volumetric flask. Distilled water was prepared from a Millipore water purification system. All glassware was cleaned with ultrapure water for three times and dried before use.

Instrumentation and methods

A pH meter (Sartorius PB-10) was used to determine the pH value. UV-vis absorption spectra were measured at room temperature with a Shimadzu UV-2550 spectrophotometer in 1 cm quartz. The fluorescence spectra were taken on a Hitachi F-7000 Fluorescence spectrophotometer (Hitachi, Japan) using a 5.0 nm slit width in a 1×1 cm quartz cell. The cells fluorescence imaging was obtained by laser scanning confocal fluorescence microscopy Olympus FV1000 (Olympus Corporation, Tokyo, Japan). High resolution mass spectra (HRMS) were recorded on a micro TOF-QII mass spectrometer. 1 H NMR and 13 C NMR spectra were obtained by Varian Gemini 2000 DMX600 MHz FT NMR spectrometers using CDCl3/DMSO-d6 as solvent and TMS as internal standard. The chemical shifts were represented by ppm. The quantum yield was measured at room temperature referenced to Quinine sulphate dehydrates in aqueous solution of 0.1 M sulfuric acid. The structures of synthesized compounds were characterized by HRMS, 1 H NMR, 13 C NMR and the related spectra are shown in the Supplementary Information.

2.3.

Synthesis of compound 2-((7-(dimethylamino) -4,5dihydronaphtho [1,2-b] thiophen -2-yl) methylene) malononitrile (DDTM) A solution of DTC [33] (0.0500 g, 0.19 mmol, 1.0 equiv.) and malononitrile (0.0130 g, 0.19 mmol, 1.0 equiv.) in anhydrous CH2Cl2 (5 mL) was treated with triethylamine (0.2 mL) at room temperature under argon. The reaction mixture was stirred at room temperature for 3 h. The solvent was evaporated under reduced pressure, and the residue was purified by silica gel column chromatography. (Yield: 40.5%, Rf = 0.5, PE: EtOAc= 3:1). 1 H NMR (600 MHz, CDCl3) δ 7.58 (s, 1H), 7.45 (s, 1H), 7.36 (d, J = 8.5 Hz, 1H), 6.62-6.48 (m, 2H), 3.04 (s, 6H), 2.96-2.87 (m, 2H), 2.81 (dd, J = 23.5, 16.5 Hz, 2H). 13

Titration experiments of probe DDTM

For UV-vis and fluorescence titrations, stock solution of probe DDTM (1.0 mM) was prepared in DMSO solution. Different portions of the stock solution were then diluted to different concentrations for further usage. Generally, 50 µL of probe DDTM solution was added into 5 mL of buffer solution (10 mM PBS buffer, 20% DMSO) in a colorimetric tube. The other ions were dissolved in deionized water and added to the probe solution under the same condition. The pH variations of the solution were adjusted by adding the different volumes of NaOH (0.1 M or 1.0 M). The resulting solution was shaken well and incubated for 5 min at room temperature before measurements. The spectra of these solutions were recorded by means of fluorescence method. Fluorescence measurements were carried out with a slit width of 5 nm (λex = 380 nm) in 10 mm quartz cuvettes at room temperature and the scan rate was 1200 nm/min.

Cytotoxicity experiments

HeLa cells were purchased from the ATCC Cell Bank. Cell counting kit-8 (CCK-8) was used to detect the cytotoxicity of probe DDTM. Cells were seeded in 96-well plates at a concentration of 2 × 10 4 cells per well and cultured for 24 h in DMEM (supplemented with 10% fetal bovine serum (FBS)) in an incubator (37 ℃, 5% CO2). After the cells were incubated with probe DDTM at different concentrations (0, 5, 10, 15, 20 and 25 µM) for 24 h, CCK-8 (10 µL) was added to each well of the 96-well plate for 3 h at 37 ℃, the absorbance was measured at 450 nm using a microplate reader. Five replicates were done for each treatment group.

Cell imaging of probe DDTM

HeLa cells were placed on a 20 mm diameter glass bottomed culture dish and allowed to adhere for 24 h for intracellular imaging of OH -. Before confocal scanning imaging, these cells were incubated with 10 µM of probe DDTM at 37 ℃ for 30 min followed by washing three times with phosphate-buffered saline (PBS, pH = 7.4). Subsequently, the cells were washed with PBS for three times to remove excess DDTM, then in situ treated with different pH value solution (pH 10.0 and 12.0, respectively), cells were then visualized with a FV1000 laser scanning confocal microscope. The excitation wavelength was set to 405 nm and 543 nm.

Results and discussion

Synthesis of probe DDTM

The general synthetic route of probe DDTM was shown in Scheme 2. DTC was conveniently synthesized according to our previous work. Then DTC reacted with malononitrile to obtain the probe DDTM by Knoevenagel condensation at moderate yield in one step. The structure of probe DDTM was characterized and confirmed by 1 H NMR, 13 C NMR and HRMS data.

Photophysical properties of probe DDTM

Probe DDTM good solubility in common organic solvents such as CH2Cl2, CHCl3, THF and DMSO. On one hand, probe DDTM was slightly soluble in water but DMSO could help increase its water solubility. On the other hand, the respond of probe DDTM toward OH -has the maximum fluorescence emission wavelength in DMSO solvent (Fig. S1). So DMSO was preferred to be co-solvent in the photophysical properties experiment. In addition, the maximum emission wavelength of probe DDTM toward OH -was obtained in PBS buffer (Fig. S2). Hence, probe DDTM could dissolve completely in a mixture of DMSO/PBS buffer (10 mM PBS buffer, 20% DMSO) and the corresponding solution exhibited an efficient fluorescent turn on mode for pH value change within 1 minute. The colour of probe DDTM in a mixture of DMSO/PBS buffer pH 7.4 was light pink with a distinct absorption band at 430 nm. On the contrary, the colour of the mixture was changed from light pink to yellow and the absorbance of probe DDTM drastically increases at 430 nm when pH value increased to 11.0 (Fig. S3). The change of colour was readily observed by naked-eye. In addition, no emission wavelength of probe DDTM was appeared from 450 nm to 650 nm with the fluorescence quantum efficiency of 0.01 (Quinine sulphate dehydrate in 0.1 M H2SO4 was used as the main standard, ФS = 0.54, λex=365 nm). The low quantum efficiency was occurred by ICT process from the amino group of naphthalenone scaffold as an electron acceptor group to the dicyanovinyl moiety as an electron acceptor group by its π-conjugation. To our delight, the remarkable fluorescence enhancement was observed more than 150-fold and the quantum efficiency increased to 0.64 due to the disruption of ICT phenomena of probe DDTM when pH value increased to 11.0, the fluorescence colour was changed from colourless to orange to yellow with the maximum emission wavelength at 578 nm. The dicyanovinyl group could be modulated by reacting with OH -, and this nucleophilic addition reaction could finally interrupt the πconjugation leading to the turn on mode of this detection.

UV-vis absorbance and fluorescence response of probe DDTM to pH

The UV-Vis and fluorescence properties of probe DDTM were measured in a mixture of DMSO/PBS buffer system (10 mM PBS buffer, 20% DMSO) for detecting pH changes. As shown in Fig. 1(a), the absorbance at 430 nm was gradually increased and had a good linear relationship with pH value in the range from 10.5 to 13.5 (R 2 = 0.9946) along with the pH value increase. At the same time, the solution colour was drastically changed from light pink to yellow accompanied with these absorption changes. So probe DDTM could be used as nakedeye-visible pH indicators in aqueous media due to their sensitive colour response to solution pH change. Further for fluorescent responses, as depicted in Fig. 2, probe DDTM itself had no fluorescence intensity at 578 nm, however, the fluorescence intensity increased with a concomitant enhancement in pH changes. With the variation of pH from 8.0 to 14.0, the probe showed turn-on fluorescence. Notably, a good linearity (R 2 = 0.9958) between the florescence intensities with pH in the range of 9.5-10.7 was obtained, the regression equation was Y=2324.5X-21250. Additionally, on the basis of the Henderson-Hasselbalch equation (log [(Imax-I)/ (I-Imin)]=pH -pKa), the pKa value of probe DDTM was determined to be 9.92 by analysis of the dependency of fluorescence intensities at 578 nm (Fig. 3). These results demonstrated that probe DDTM could be used to detect pH value quantitatively in extreme alkaline solution.

Selectivity experiment of probe DDTM to pH

The intracellular environment is quite complex and many biological molecules are present, such as amino acids, cysteine (Cys), homocysteine (Hcy) and glutathione (GSH), metal ions and reactive oxygen species (ROS) and so on. In order to verify whether these biomolecules could disturb the pH measurement efficiently in living cells, the selectivity of DDTM to alkaline over these species was investigated at pH 11.0, the results were shown in Fig. 4(a), metal ions and reactive oxygen species showed very weak fluorescence intensity at 578 nm. In contrast, the strong fluorescence signal of probe DDTM was observed at 578 nm in alkaline solutions. Meanwhile, all biomolecules of fluorescence response in neutral and alkaline were evaluated to survey competition experiments. As shown in Fig. 4(b), addition of biomolecules has no interference on the fluorescent spectroscopic changes to pH under the testing conditions. All these results indicated that probe DDTM could serve as a highly selective chemosensor for pH change in extreme alkaline solution.

Time-dependent effect and reversibility study of probe DDTM

The time-dependent course analysis showed that probe DDTM could quickly respond to pH value change for less than one minute at pH 11.0. As shown in Fig. 5(a), after addition of DDTM to PBS buffer solutions at pH 7.4 and 11.0 at room temperature, respectively, the fluorescence intensity of the probe instantly reached equilibrium and remained unchanged for at least 30 min. The experimental results proved that probe DDTM has rapid respond to pH changes in extreme alkaline solution.

As for pH probes, the reversibility was one of the vital requirements for practical applications. Therefore, the reversibility of probe DDTM was explored that the pH value of the solution was switched back and forth between pH 4.0 and 11.0 using concentrated hydrochloric acid and aqueous sodium hydroxide. As shown in Fig. 5(b), probe DDTM exhibited a highly reversible response to pH and the response time in different pH values were less than 10 s. Switching between the fluorescence on/off states could be repeated along with the colour change repeatedly between colourless (pH 4.0) and yellow (pH 11.0). These results demonstrated that probe DDTM could show quick and reversible responses to pH changes, and also be potentially useful to real-time monitor pH value changes in extreme alkaline solution.

Proposed detection mechanism and theoretical computation

The mechanism for the fluorescence enhancement of probe DDTM was proposed as shown in Fig. 6. In alkaline solution, hydroxyl ion reacted as nucleophilic reagent on the double bond of dicyanovinyl group in probe DDTM to give the addition intermediate, which could finally interrupt the πconjugation leading to produce strong fluorescence. Then hydroxyl of the addition intermediate could be occurred in protonation process in acid condition, finally probe DDTM is regenerated and the water molecule is released. To further confirm the proposed detection mechanism of probe DDTM towards hydroxyl, probe DDTM was reacted with tetrabutyl ammonium hydroxide at room temperature for 10 min in MeOH solution, then the reaction solution was diluted appropriately and tested by mass spectroscopy, an expected signal at m/z 322.9321 corresponding to addition intermediate [probe DDTM + OH] -was observed (Fig. S4), and these data were consistent with the proposed mechanism. To better understand the optical responses of probe DDTM upon adding with OH -, density functional theory (DFT) calculations were performed using the CAM-B3LYP/6-311+G(d, p) level of Gaussian 09 program [34,35]. The solvent effects were taken into account in water through the polarizable continuum model (PCM) [36]. The optimized structures of DDTM and DDTM-OH -were shown in Fig. 7. The molecular backbone of the DDTM was nearly coplanar while the dihedral angle of thiophene ring and dicyanovinyl group in the DDTM-OH -was 48.9°, indicating the dicyanovinyl group was twisted in response to OH -. Molecular orbital surfaces of the HOMO and LUMO for compound DDTM and DDTM-OH -were displayed in Fig. 8. For DDTM, the electron density in the HOMO was distributed over the whole molecular while that in the LUMO was mainly concentrated on the thiophene ring and dicyanovinyl group, respectively. As for DDTM-OH -, the electron in HOMO and LUMO orbitals was mostly spread over the moiety of DTC. Obviously, there was no ICT process occurred in the DDTM-OH -, in that a sp 3 hybridized carbon atom existed in vinyl moiety of dicyanomethylidene and dicyanomethylidene twisted. The interruption of the ICT process also resulted in the increased HOMO-LUMO energy gap of DDTM-OH -(6.39 eV) in comparison with that of the DDTM (4.63 eV) and is also explained for the appearance of a new absorption band at 430 nm in the absorption spectrum of DDTM-OH -. Therefore, these DFT calculations verified that OH-was bound to the vinyl group of probe DDTM and are in good agreement with the experimental results.

Cytotoxicity and fluorescence imaging for living cells

Firstly, the cytotoxic effects of probe DDTM were studied in order to evaluate its imaging performance in vivo. MTT, a standard cell viability protocol, was employed to evaluate the cytotoxicity of probe DDTM on HeLa cells by treating HeLa cells with the probe for 48 h. The probe showed low cytotoxicity and excellent compatibility with more than 80% viability at concentration of 25 μM under experimental conditions (Fig. S3).

Next, the cellular imaging assays were performed to verify whether probe DDTM could be applied for determining pH changes in the biological application. HeLa cells were incubated with probe DDTM at pH 7.4, 10.0 and 12.0 for 5 min, after that intracellular fluorescence was monitored by confocal microscopy (Fig. 9). Under neutral condition, no fluorescence was observed. On the contrary, moderate fluorescence was observed when cells were incubated with probe DDTM at pH 10.0. Furthermore, the strong blue fluorescence was observed at pH 12.0. At the same time, the similar results were obtained in red channel when excitation wavelength was 543 nm. These results indicated that probe DDTM displays sensitive fluorescence enhances in a pHdependent manner in living cells system.

Conclusions

In summary, a novel pH-sensitive colorimetric sensors based on naphthalenone scaffold was designed and synthesized.

Probe DDTM showed a naked eye visible colour changes from light pink to yellow, and the fluorescence colour was changed from colourless to orange to yellow in the range pH 9.0-14.0. In addition, the probe exhibited strong pH-dependent behaviour, high selectivity and rapidly respond to pH fluctuations within the range of 9.0 -14.0. Furthermore, probe DDTM displayed a notably large Stokes shift of 198 nm. Live cell imaging data revealed that probe DDTM could selectively monitor pH changes with low cytotoxicity and cell membrane permeability. Therefore, the probe could be acted as an effective intracellular pH imaging agent under extreme alkaline conditions in the biomedical and biological fields.

Scheme 2 .

 2 Scheme 2. Synthesis of probe DDTM.

Fig. 1 .

 1 Fig. 1. (a) UV-Vis absorption spectra of probe DDTM (10 μM in 10 mM PBS buffer, 20% DMSO) in accurate pH values ranging from 5.0 to 14.0 at 25℃. (b) Linear relationship of UV-Vis absorption at 430 nm and pH values from 10.5 to 13.5. Inset: the colour of probe DDTM in the different pH values under daylight lamp.

Fig. 2 .

 2 Fig. 2. (a) Emission spectra of probe DDTM (10 μM in 10 mM PBS buffer, 20% DMSO) in accurate pH values ranging from 5.0 to 14.0 at 25℃ (λex=380 nm, slit: 5.0 nm/5.0 nm). Inset: Fluorescence intensity at 578 nm by pH values according to the fluorescence titration (pH 5.0~14.0). (b) The linear relationship of fluorescence intensity at 578 nm and pH values from 9.5 to 10.7. Inset: the color of probe DDTM in the different pH values under UV 365 nm.

Fig. 3 .

 3 Fig. 3. Linear regression relationship between the pH value and log [(Fmax-Fx)/(Fx-Fmin)].

Fig. 4 .

 4 Fig. 4. (a) Emission spectra of probe DDTM in the presence of different metal ions and biologically relevant analytes (10 μM) in the pH 11.0 solution. (b) Fluorescence respond of probe DDTM (10 μM) in the presence of different metal ions and biologically relevant analytes 1.0 mM in the two pH values. 1: Blank; 2: AcO -; 3: ClO4 -; 4: CO3 2-; 5: ClO -; 6: NO3 -; 7: ONOO -; 8: SO3 2-; 9: Cys; 10: GSH; 11: H2O2 (0.5 mM); 12: Ba 2+ ; 13: Ca 2+ ; 14: CN -; 15: Cu 2+ ; 16: Fe 3+ ; 17: Hg 2+ ; 18: K + ; 19: 2+ ; 20: Na + ; 21: DTT; 22: H2S. (λex=380 nm, slit: 5.0 nm/5.0 nm).

Fig. 5 .

 5 Fig. 5. (a) Time course of fluorescence intensity of probe DDTM (10 μM in 10 mM PBS buffer, 20% DMSO) in the two pH values. (b) pH reversibility of probe DDTM (10 μM in 10 mM PBS buffer, 20% DMSO) between pH 4.0 and pH 11.0 at 25 ℃ (λex=380 nm, slit: 5.0 nm/5.0 nm).

Fig. 6 .

 6 Fig. 6. Proposed mechanism of probe DDTM towards pH changes.

Fig. 7 .

 7 Fig. 7. Optimized structures of DDTM and DDTM-OH -at the CAM-B3LYP/6-311+G(d, p)/PCM theoretical level.

Fig. 8 .

 8 Fig. 8. HOMO-LUMO energy levels and frontier molecular orbitals of probe DDTM and the adduct DDTM-OH -.

Fig. 9 .

 9 Fig. 9. Confocal fluorescent imaging of HeLa cells incubated with probe DDTM (10 μM) at pH 7.4 (a, b, c), 10.0 (d, e, f) and 12.0 (g, h,

Fig. 10 . 3 . 8 .

 1038 Fig. 10. Dynamic fluorescence images of living HeLa cells on incubation with 10 μM DDTM for different times 10 min (a, b, c) and 30 min (d, e, f) in pH 12.0 medium. Left: Excitation wavelength was 405 nm. Right: Excitation wavelength was 543 nm. Scale bar: 40 μm.3.8. Practical application of probe DDTM in pH test stripsIn order to further explore the practical application of probe DDTM, pH test strips were prepared and utilized to sense different pH according to the excellent properties of probe DDTM. The test strip was prepared by immersing filter papers in probe DDTM solution (10 μM in 20% DMSO). After drying in the air, the strip displayed different colour change when the test strips immersed in different pH solution. As shown in Fig.11, the changes of colour in different pH solution were observed in naked eyes and fluorescent light, respectively. Therefore, the distinct colour changes promoted probe DDTM to be conveniently used for precise pH paper preparation in real time applications.

Fig. 11 .

 11 Fig. 11. Photographs of colorimetric responses of probe DDTM test strips with different pH solutions. (a) Naked eye in visible light. (b) Fluorescence colour changes under 365 nm UV lamp.

  C NMR (151 MHz, CDCl3) δ 152.53, 151.97, 149.76, 139.29, 138.62, 136.61, 131.25, 130.63, 126.94, 118.68, 114.93, 111.54, 71.69, 40.62, 32.28, 30.39, 30.06, 23.90. HRMS (C18H15N3S): calcd. for [M+H] + 306.1059; found: [M+H] + 306.0987.

Acknowledgements

We sincerely thank the Program for Changjiang Scholars and Innovative Research Team in University (No. IRT_15R55), the Science Foundation of Northwest University (No. 15NW17), the International Science & Technology Cooperation Program of Shaanxi Province (No. 2016KW-003) and Opening Foundation of Key Laboratory of Resource Biology and Biotechnology in Western China (Northwest University), Ministry of Education.

Notes and references

Conflicts of interest

There are no conflicts of interest to declare.