Roch Mollero

Xavier Pennec

Hervé Delingette

Nicholas Ayache

Maxime Sermesant

A Multiscale Cardiac Model for Fast Personalisation and Exploitation

Computer models of the heart are of increasing interest for clinical applications due to their discriminative and predictive abilities. However a single 3D simulation can be computationally expensive and long, which can make some practical applications such as the personalisation phase, or a sensitivity analysis of mechanical parameters over the simulated behaviour quite slow. In this manuscript we present a multiscale 0D/3D model which allows us to have a reliable (and extremely fast) approximation of the behaviour of the 3D model under a few simplifying assumptions. We first detail the two different models, then explain the coupling of the two models to get fast 0D approximation of 3D simulations. Finally we demonstrated how the multiscale model can speed-up an efficient optimization algorithm, which enables a fast personalisation of the 3D simulations by leveraging on the advantages of each scale.

Introduction

Three-D personalised cardiac models simulate the physical behaviour of a patient heart, in order to perform advanced analysis of the cardiac function. The movement of the myocardium is calculated under the influence of simulated electromechanics and haemodynamic condition, which requires complex and expensive calculations, particularly as the spatial and temporal resolution of the model increases (1 hour per heart beat, even several days for some models). This computational burden particularly impacts applications of the models where many simulations are required. The model personalisation phase (fitting the available clinical data) is particularly impacted, slowing the study of clinical applications such as disease caracterisation and prediction of case evolution

To tackle this computational burden, surrogate models of the cardiac function have been developed, for example for uncertainty quantification [START_REF] Neumann | Robust Image-Based Estimation of Cardiac Tissue Parameters and Their Uncertainty from Noisy Data[END_REF] or probability estimation [START_REF] Konukoglu | Efficient probabilistic model personalization integrating uncertainty on data and parameters: Application to eikonal-diffusion models in cardiac electrophysiology[END_REF] over the cardiac parameters. Those models usually rely on regression models which require many pre-computations to approximate nonlinear behaviours, thus scaling badly as the number of varying parameters increases.

Here we built a multiscale cardiac model by developping an novel coupling method between an original 3D cardiac model and a reduced lumped "0D" version of the same model. While the coupling of different scales has been used to improve boundary conditions in cardiac modelling and computational fluid dynamics [START_REF] Esmaily | A modular numerical method for implicit 0d/3d coupling in cardiovascular finite element simulations[END_REF], we propose here the simultaneous simulation of the same model at different scales, which enables extremely fast approximations of the behaviour of the 3D model with 0D simulations. It provides a natural approximation of all the desired cardiac values. Furthermore, we demonstrate the efficiency of the multiscale 0D/3D model when plugged into a personalisation process driven by the genetic algorithm CMA-ES. This creates a very fast and computationally efficient personalisation agent, successfully tested on the personalisation of 34 hearts with different geometries.

Multiscale Cardiac Model

In order to build a computer model of the myocardium, one needs different components to represent the passive stiffness, the active contraction and the boundary conditions (hemodynamics). We use here the Bestel-Clement-Sorine model (BCS) of the sarcomere contraction, in conjonction with the Mooney-Rivlin model for passive elasticity as described in [START_REF] Chapelle | Energy-preserving muscle tissue model: formulation and compatible discretizations[END_REF]. Hemodynamics are represented through global values of pressures and flows in the cardiac chambers.

Three-D Cardiac Model

Here we first use 3D biventricular mesh geometry extracted from MRI images, with synthetic myocardial fibres. The active sarcomere force and the passive force are computed at each point of the mesh with the Bestel-Clement-Sorine (BCS) model under the influence of the pre-computed electrophysiology. The boundary condition is a circulation model with the 4 phases of the cardiac cycle based using a 3-parameters Windkessel model of blood pressure as after-load. We also specify a resting mesh in addition of the initial mesh geometry to start the simulation with a geometry which already includes residual stress.

The complete details of the 3D electromechanical model as well as the C++ implementation of the finite element solver are discussed in [START_REF] Marchesseau | Fast parameter calibration of a cardiac electromechanical model from medical images based on the unscented transform[END_REF]. With meshes of 10 000 and 15 000 points and a time step of 5 ms, a single beat (around 0.9 s) takes between 25 minutes and 1 hour to compute. Finally, since we want to consider cardiac beats which are not affected by the initialization of the simulation, we compute 3 subsequent beats per simulation and only consider the last one, so the full 3D simulation takes between 75 minutes and 3 hours to be computed.

0D Cardiac Model

If we then make simplifying assumptions on the geometry and on the properties of the material, by approximating the ventricular shape as a sphere, we can derive the equations of a fast lumped model of the heart as detailed in [START_REF] Caruel | Dimensional reductions of a cardiac model for effective validation and calibration[END_REF].

The myocardial motion is then only described by the inner radius r of the ventricle (the "Resting Mesh" in 3D becomes simply a "Resting Radius" r0), the material is considered incompressible and has idealised tangential fibres. The electrophysiological activation is supposed synchronous and homogeneous over the sphere as well as the resulting active sarcomere force. Finally we use the same circulation model and Mooney Rivlin equations than in the 3D model.

We implemented the 0D model within the CellML framework, then converted into C code with OpenCOR [START_REF] Garny | Opencor: a modular and interoperable approach to computational biology[END_REF]. With a temporal discretisation of 0.1 milliseconds, we can simulate 10 to 20 heart beats per seconds.

Multiscale Coupling

Although with different mechanical assumptions, both models follow the same trends in behaviour when we change the global parameters. The idea of the coupling is to have those trends match within a specific domain of 3D parameters, in order to predict the behaviour of many 3D simulations by only computing 0D simulations. This is done by first matching several 3D simulations of the domain (called the sigma-simulations) with the 0D model, then learning a mapping between the 3D parameters and the 0D parameters.

Personalisation of the 0D Model to a 3D Simulation

The first step of the coupling is to make the two models match on a single simulated beat. Given one 3D simulated heartbeat, we want to accurately reproduce the pump behaviour of the heart with the 0D model (under the same boundary conditions). In particular we focus on reproducing 6 features of the volume curve by setting 5 main electromechanical parameters in the 0D model, both listed in Table 1.

We estimate the parameters by performing a derivative-free optimisation with the genetic algorithm CMA-ES [START_REF] Hansen | The cma evolution strategy: a comparing review[END_REF]. This powerful algorithm combines Bayesian principles of Maximum Likelihood with a local gradient descent on the members of each generation, by updating at each iteration both an internal covariance Ic and mean Im of the parameters to explore at the next iteration. We define the fitness (score) of each simulation by: f itness(x) = (V alues(x) -T argetV alues)/N ormalization which is the L 2 distance between the simulated and the objective values, normalized to be able to compare volumes (unit volume = 10 ml) with times (unit time = 20 ms), percentages (unit percent = 5 %) and flow (unit flow = 10 ml/s).

We use a population of 20 simulations per generation, and optimize over 100 generations. This takes less than 10 minutes on a multi-core computer with parallel computation of the simulations within each generation. Figure 3 shows the matching behaviours of a 3D simulation and its 0D counterpart.

Prediction of 3D behaviour from 0D Simulations

We want to be able to get reliable approximations of other 3D simulations within a domain Ω 3D of 3D parameters by doing 0D simulations, so we need to be able to convert a vector X ∈ Ω 3D into a relevant vector Y of 0D parameters which best approximates the 3D simulations. Here we focus in the mapping of N 1 =4 3D parameters toward N 2 =5 0D parameters listed in Table 2.

LV Volume (ml) over time t(s)

LV Flow (ml/s) over time t(s) The parameters do not behave necessarily the same for both models (and not in the same range of values). For example the 3D and 0D dampings are not calculated in the same way, and the "resting position" variable is a unitless value from 0 to 10 in the 3D case, whereas it is the "inner radius" (in cm) in the 0D mdoel. Also even for parameters coming from the same equations in both models (such as σ and c1) the values might be different for two matching 0D and 3D simulations. For example we can get the same volume curve with contractility and stiffness values which are 10 times smaller in the 0D implementation than in the 3D setting). The mapping is thus not trivial and is done through the following steps:

1. We first perform 2N 1 + 1 3D simulations at 2 points of each of the N 1 principal directions of the domain of interest Ω 3D and its center (called the sigma-simulations with parameters S i ). 2. We then fit the 0D model on each of those 2N 1 + 1 simulations with the method explained in the previous subsection, which gives the 2N 1 + 1 corresponding vectors of 0D parameters V i . 3. Then we compute a least square linear regression φ between the parameters of S i and the parameters V i . 4. Finally for each vector X of Ω 3D , we get the relevant vector Y = φ(X) of 0D parameters by using this linear regression.

Results of the Parameter Mapping

Here we evaluate the approximation of 80 other 3D simulations with parameters X j ∈ Ω 3D by their corresponding 0D simulations with parameters Y j = φ(X j ). We focus on the error made on 7 features, which we compare with the error made by a direct linear regression between the parameters of the sigma-simulations and the 7 features (which is the best regression model we can run with this small number of points). Results are shown in Table 3.

Both methods provide good approximations of the features. However if we exclude tVmin (which is harder to approximate because of the temporal discretisation of the models), the simulations with the 0D model provide overall better approximations of the values, particularly of those with the most nonlinear behaviour such as the percentage of filling (FVat) and the flow (Qat) when the atrium opens (resp by 43% and 40%). This is because the 0D model naturally captures the nonlinear behaviour of those features which would require many simulations to be well approximated by a nonlinear regression model.

Efficient Personalisation of the Multiscale Model

Personalisation is the first step for the application of 3D cardiac modelling to patient-specific data. After extracting a patient-specific mesh geometry from the images, we want the 3D simulation to match specific clinical data. This can be done as in Section 3.1, by optimizing some parameters of the 3D cardiac model with CMA-ES to minimize a score related to the distance between the simulation and the target clinical data. Since the convergence speed of CMA-ES increases with the size of the population evaluated at each iteration, we want to set this size as high as possible, but the computational cost of each 3D simulation is the limiting factor. However with the coupled model, we show that we can evaluate the scores of the members of a large population with a limited number of 3D simulations, which allows to benefit from the increased speed of convergence while keeping a reasonable computational burden.

Full Personalisation Pipeline

In the personalisation process of the 3D model, we want to reproduce 5 features of the volume curve (listed in Table 4) by optimizing the 4 parameters described in Section 3.2. We initialize a CMA-ES optimization process with a high population size m=80, and derive the Coupled-CMA optimization method by adapting the original CMA-ES as follows:

At each iteration, the CMA-ES algorithm asks for the scores of m 3D simulations X j whose parameters are drawn from a multivariate distribution, ruled by its own internal mean Im and covariance Ic.

Instead of performing these m 3D simulations, we calculate a mapping φ and m surrogate 0D simulations with parameters Y j = φ(X j ), as explained in 3.2. Only the 2N 1 + 1 = 9 3D sigma-simulations are required, which are drawn at the center and within the principal axis of the multivariate distribution (two simulations per axis at +/-standard deviation of the axis from the mean).

Then we calculate the scores of these m surrogate simulations Y j = φ(X j ) and provide the scores to the CMA-ES algorithm "as if" they were the scores of the m 3D simulations.

Results

In order to compare the efficiency of the Coupled-CMA optimization, we compare it with the full CMA optimization with a population size of 80 (CMA-80), as well as with a direct CMA optimization with a population size of 9 (CMA-9) which has the same computational burden than the Coupled-CMA. This is done by comparing the scores of the 3D simulations made with the algorithm's internal mean Im at each iteration, which is the algorithm's current best estimate of the minimal parameters (In the case of Coupled-CMA, it is the score of the central sigma-simulation). Results are shown in Figure 5 (left).

The CMA-80 expectedly converges a lot faster than the CMA-9, but the main result is that the Coupled-CMA optimization converges at the same speed than the CMA-80. From a computing viewpoint, the Coupled-CMA lowers the number of 3D simulations from 80 to 9 at each iteration while keeping the accuracy, thus saving between 50 and 150 hours of CPU time at each iteration.

This similarity of convergence between CMA-80 and Coupled-CMA was confirmed on other geometries and values to fit, and overall we applied this method to the personalisation of 34 patient hearts, sucessfully matching the features to their target clinical values, within 6% of the mean value of each feature in the population, for at least 90% of the cases, as show on Figure 5 (right).

Conclusion

In this manuscript we detailed how a fast reduced "0D" model relying on the same equations but simplifying assumptions can be coupled to the original 3D model to approximate computationally-costly 3D simulations with fast 0D simulations within a domain of 3D parameters. For this we developed a novel method to learn a mapping between the parameters of 3D sigma-simulations within the domain, and the parameters of corresponding 0D simulations with the same mechanical behaviour. In learning how to go from 3D parameters to 0D parameters, we learned the difference between the 0D (lumped) hypothesis and the 3D hypothesis. We then used the coupled model to build a fast personalisation agent, by substituting many 3D simulations in the optimization process by their 0D counterparts. This enabled us to benefit from the higher speed of convergence while considerably reducing the computational burden.

Several aspects of this method can be directly extended to more (nonlinear) features and parameters, such as the pressure, heart rate and the boundary conditions. We could also improve the computational gain by reusing the mapping from one iteration of the CMA-ES to the next one when the domains are similar enough instead of recomputing it. Finally, we expect the personalisation method to scale well to more complex optimization problems with a larger parameter space, since CMA-ES is well suited for those problems.

Fig. 1 :

 1 Fig. 1: Multiscale model and personalisation pipeline

Fig. 3 :

 3 Fig. 3: Comparison between a volume curve simulated with the 3D model (blue) and the matched 0D simulation (red). Orange lines and points outline the 6 fitted values

Fig. 5 :

 5 Fig. 5: Left: Score per iteration for the 3D methods. Right: Fitting error of each feature for the 34 patients (in percentage of the mean observed value of each feature)

Table 1 :

 1 Features and parameters in the 0D model personalisation over a 3D simulation

	Features	0D Parameters
	Maximal Volume Vmax	Maximal Contraction σ
	Minimal Volume Vmin	Stiffness c1
	Percentage of filling when the atrium contracts FVat Resting radius r0
	Flow when the atrium contracts Qat	damping eta
	Time where the volume is minimum tVmin	Atrio-Ventricular Delay AV.
	Time when ejection starts teject

Table 2 :

 2 Three-D and 0D Parameters of the Mapping

	Three-D Parameters	0D Parameters
	Maximal Contraction σ Maximal Contraction σ
	Stiffness c1	Stiffness c1
	Resting Mesh RM	Resting radius r0
	3D damping µ	0D damping eta
		Atrio-Ventricular Delay AV

Table 3 :

 3 Feature approximation by the Mapping: Range of values of each feature, mean absolute error with the mapping (M AEm) and mean absolute error with the direct regression (M AEr)

	Feature tVmin (ms) Vmax (ml) Vmin (ml) EF (%) SV (ml) FVat (%) Qat (ml/s)
	Range 30.3	9.97	34.6	27.7	29.7	11.1	93.9
	M AEm 6.35	2.03	1.06	0.17	0.98	0.78	3.62
	M AEr 3.86	2.09	0.89	0.15	1.21	1.38	5.99

Table 4 :

 4 Features and parameters in the 3D model personalisation to a real volume curve

	Features	3D Parameters
	Maximal Volume Vmax	Maximal Contraction σ
	Stroke Volume SV	Stiffness c1
	Percentage of filling when the atrium contracts FVat Resting Mesh RM
	Flow when the atrium contracts Qat	3D damping µ
	Time where the volume is minimum tVmin

Ackowledgements This work has been partially funded by the EU FP7-funded project MD-Paedigree (Grant Agreement 600932) and contributes to the objectives of the ERC advanced grant MedYMA (2011-291080).


Dong Wang
email: wangd@nercita.org.cnwangjh@nercita.org.cn

Donghai Han

Zhihong Ma

Ligang Pan

Liu Zhao

Ping Han

Researcher, Ph.D Jihua Wang
email: wangd@nercita.org.cn

The Application of Classical Least Square Algorithm in the Quantitative Analysis of Lime in Wheat Flour by ATR-MIR Spectroscopy

Keywords: Classical Least Square, Attenuate Total Reflection, Mid-Infrared Spectroscopy, Wheat Flour, Lime

In this thesis, classical least square (CLS) regression was applied in the attenuate total reflection mid-infrared (ATR-MIR) spectra data processing of lime in wheat flour to develop the corrected peak height model of the raw spectra, the corrected peak area model of the raw spectra, the corrected peak

Introduction

Wheat flour is one of the important food sources for human beings, of which the quality safety is highly focused by all quarters of society. However, the quality safety incidents about wheat flour have been reported frequently in recent years, such as lime added in it, which results in not only negative influence but also being harmful to the customers' health. Therefore, the development of rapid, simple, practical and low-cost determination method for detecting lime in wheat flour becomes particularly urgent and very important.

Among the common instrument analytical methods, mid-infrared (MIR) spectroscopy analysis is widely applied in the fields of the identification and characterization of microbes [START_REF] Santos | Fourier transform infrared as a powerful technique for the identification and characterization of filamentous fungi and yeasts[END_REF], detection and identification of bacteria [START_REF] Hamzah | Rasco: Fourier transform infrared spectroscopy, detection and identification of Escherichia coli O157:H7 and Alicyclobacillus strains in apple juice[J][END_REF], quality analysis of edible oil [START_REF] Gurdeniz | Detection of adulteration of extra-virgin olive oil by chemometric analysis of mid-infrared spectral data[END_REF][START_REF] Rubé N | Monitoring of fatty acid composition in virgin olive oil by Fourier transformed infrared spectroscopy coupled with partial least squares[END_REF], quantitative determination of anthocyanins in cherry [START_REF] Pappas | Quantitative determination of anthocyanins in three sweet cherry varieties using diffuse reflectance infrared Fourier transform spectroscopy[END_REF] with the advantages of distinguishable peaks, strong characteristics, high sensitivity, rapid analysis and non-invasive analysis.

However, the traditional MIR analysis involves the sample preparation such as tabletting with KBr, which is difficult for the operators. Attenuated total reflection (ATR) technology can acquire the information of samples according to the reflective signals of the sample surface, which is free of sample preparation. ATR technology has the advantages of high sensitivity, distinguishable peaks, user-friendly control and non-interference by humidity. Based upon the advantages mentioned above, ATR technology is now being widely applied in the fields of analysis of juice sweetness [START_REF] Jha | Authentication of sweetness of mango juice using Fourier transform infrared-attenuated total reflection spectroscopy[END_REF], biology structure studies [START_REF] Vassiliki | Hamodrakas: Dogfish egg case structural studies by ATR FT-IR and FT-Raman spectroscopy[END_REF], material study [START_REF] Sune | In situ ATR-IR study of nitrite hydrogenation over Pd/Al2O3[J][END_REF] and so forth. Therefore, it is convenient to collect the MIR spectra of powder samples, such as wheat flour, by ATR technology.

It is known by the experiment data that lime, such as calcium oxide, calcium hydroxide, has the characteristic absorption peaks in the MIR range, which is the theoretical foundation of the qualitative and quantitative analysis to it. The characteristic absorption peak in MIR range of wheat flour is far away from that of lime, i.e. the interference from wheat flour is not very obvious. Based on the MIR spectral characteristics of wheat flour and lime, CLS algorithm has a certain potential application in the quantitative analysis of lime in wheat flour by MIR spectra data.

Classical least square (CLS) regression is a data processing method applied in the data regression of analytical chemistry commonly [START_REF] Li | Spectrophotometric determination of cefotaxime sodium with potassium ferricyanide[END_REF][START_REF] Yu | Determination of tungsten in rare earth-magnesium master alloy by spectrophotometry with thiocyanate and acid digestion[END_REF]. The basic principle of CLS is that to calculate the evaluated values of the regression coefficients between the independents and dependents according to a group of obvious values in order to decrease the errors between the specified values and the predicted values furthest. The obvious errors of CLS obey three hypotheses, i.e. the obvious values are independent each other with the same precision, the obvious errors are non-systematical, the obvious errors are as the law of normal distribution [START_REF] Lu | Modern Near-infrared Spectroscopy Analysis Technology[END_REF]. Compared to partial least square regression and principal component regression, CLS regression is easy to understand relatively. Meanwhile, it is easier to achieve CLS algorithm on the special MIR detection instruments, which is in favor of the promotion of this technology.

In this research, ATR-MIR spectra of the wheat samples mixed with lime were collected. The quantitative determination of lime in wheat flour by ATR-MIR combined with CLS algorithm was studied. The corrected peak height / area models of the raw spectra and the corrected peak height / area models of the 2 nd derivative spectra were developed respectively. The external validation set was employed to evaluate the prediction performance of the models. This method is of great importance to the rapid quantitative screening of lime in wheat flour. This research, to some extent, will provide some references to not only the design and manufacturing the special MIR spectrometer for the quality safety screening of wheat flour, but also the quantitative determination of other banned additives in wheat flour.

Experiment

Instrument and Materials

Instrument: FT-MIR spectrometer, type: Spectrum 400, Perkin Elmer, U.S..

The ATR-MIR spectra are collected by FT-MIR spectrometer combined with ATR accessory, of which, the reference spectrum is collected by air, the wavenumber range is 4000 cm -1 -540 cm -1 , the wavenumber resolution is 4 cm -1 , the time of accumulation is 240. The reference spectrum is collected once every 30 minutes.

Materials: wheat flour milled by the winter wheat from Changping experimental base of Beijing, without any additives. Calcium oxide and calcium hydroxide: A.R..

Method

The quantitative calibration models of lime are developed by the corrected peak height and corrected peak area of lime in ATR-MIR spectra combined with the algorithm of classical least square. The external validation set, of which the average and standard deviation are similar to the calibration set, is employed to evaluate the prediction performances of the models. The contents of lime in both the calibration set and external validation set are distributed in the range of 0.6 % -16.0 % uniformly.

The mass of each sample is about 50 grams. The mass data of wheat flour and lime in each sample are recorded according to analytical balance accurately. The information of the content of lime in the samples is showed in table 1.

Tab. In order to decrease the error of peak height and peak area measurement, the corrected peak height / area models of the 2 nd derivative spectra are discussed despite the corrected peak height / area models of the raw spectra. The 2 nd derivative [START_REF] Shiroma | Application of NIR and MIR spectroscopy in quality control of potato chips[J][END_REF] is one of the data preprocessing method that is applied in the data preprocessing in MIR data commonly, which can eliminate the baseline drifting and baseline translation.

The corrected height is the distance between the peak top of the target peak and the Corrected peak height of raw spectra y = 1.0020x + 0.0061 0.9649 1.17 378

Corrected peak area of raw spectra y = 0.9981x + 0.0198 0.9699 1.07 444

Corrected peak height of 2 nd derivative spectra y = 0.9988x + 0.0106 0.9647 1.17 376

Corrected peak area of 2 nd derivative spectra y = 0.9947x + 0.0497 0.9606 1.23 334

It can be seen from figure 3 that for all of the regression equations by self-prediction, the slops are all close to 1, the intercepts are close to 0, the correlation coefficients are no less than 0.96, the F-test values are much larger than that looked up in F-test value table (F (1,28) = 7.64 when the significance level α = 0.01), which indicates that, for each model, a very remarkable correlation exists between the specified values and estimated values by the self-prediction of the calibration model, i.e. CLS regression has a certain potential application in the quantitative analysis of lime in wheat flour by ATR-MIR spectra data.

External Validation

The external validation set is employed to evaluate the prediction performance of the Comparing the four models, it can be seen that the differences among the four models are not significant in regression equations, correlation coefficients and the prediction results of the external validation set, the reason of which is that the characteristic peak of lime in ATR-MIR spectrum (3642 cm -1 ) is far away from the absorption peak of wheat flour in ATR-MIR spectrum (3300 cm -1 ), therefore, the interference from the spectral background becomes weak. According to that the detection limit is three times the RMSEC [START_REF] Lu | Modern Near-infrared Spectroscopy Analysis Technology[END_REF], the detection limits of the corrected peak height / area models of the raw spectra and the corrected peak height / area models of the 2 nd derivative spectra are 3.51 %, 3.21 %, 3.51 % and 3.69 % respectively, which are much less than that added in wheat flour for the purpose of weight increment. This method can fulfill the demand of rapid quality screening of wheat flour in the market.

Conclusion

Classical least square algorithm was applied in the development of the quantitative calibration models of lime in the wheat flour samples. The result indicated that the CLS regression has a certain potential application in the quantitative determination of lime in wheat flour by ATR-MIR spectra data. The detection limits of the corrected peak height / area models of the raw spectra and the corrected peak height / area models of the 2 nd derivative spectra are 3.51 %, 3.21 %, 3.51 % and 3.69 % respectively, which are much less than that added in wheat four for the purpose of weight increment. This method can fulfill the demand of rapid quality screening of wheat flour in the market, which is of great and positive meaning to the monitoring of wheat flour in the market and the guarantee of the customers' health. This research, to some extent, will also provide some references to the quantitative determination of other banned additives and the design and manufacturing of the special MIR spectrometer for the quality screening of wheat flour in the market. The method to decrease the detection limit is still being researched.

Fig. 1 .Fig. 2 .Fig. 3 .

 123 Fig. 1. the sketch maps of the measurement of the corrected peak height (A1) and the corrected peak area (A2) of the raw spectra and the corrected peak height (B1) and the corrected peak area (B2) of the 2 nd derivative spectra

= 7 .Fig. 4 .Tab. 3 .

 743 Fig. 4. The correlations between the estimated and specified values of the calibration set

Fig. 5 .

 5 Fig. 5. The correlations between the estimated and specified values of the external validation set A1: the corrected peak height model of the raw spectra, A2: the corrected peak area model of the raw spectra, B1: the corrected peak height model of the 2 nd derivative spectra, B2: the corrected peak area model of the 2 nd derivative spectra

table 4 .

 4 Tab. 4. The regression equations, correlation coefficients (r p ), SEP, RPD and F-test values (F p )

	of the external validation set
	Models	Regression equation	r p	SEP RPD	F p
	Corrected peak height of raw spectra	y = 1.4658x -1.1709	0.9954 2.59 1.73	862
	Corrected peak area of raw spectra	y = 1.4471x -1.1580 0.9954 2.47 1.81	862
	Corrected peak height of 2 nd derivative spectra	y = 1.4425x -1.0632 0.9964 2.48 1.81 1119
	Corrected peak area of	y = 1.4224x -1.0096 0.9935 2.42 1.85	607

Acknowledgment. This thesis is supported by Beijing Postdoctoral Research Foundation; China Postdoctoral Science Foundation (49 th ) (No. 20110490318); Postdoctoral Science Foundation of Beijing Academy of Agriculture and Forestry Science; Special Research of Public Welfare Industry, the Ministry of Agriculture of People's Republic of China (No. 201203046) and National High Technology Research and Development Program 863 (2010AA10Z201).

Supported by Beijing Postdoctoral Research Foundation; China Postdoctoral Science Foundation (49 th ) (No. 20110490318); Postdoctoral Science Foundation of Beijing Academy of Agriculture and Forestry Science; Special Research of Public Welfare Industry, the Ministry of Agriculture of People's Republic of China (No. 201203046) and National High Technology Research and Development Program 863 (2010AA10Z201).


Hao Ma

Haiyan Ji

Xue Liang

Zhenhong Rao

Study of Quantitative Analysis for Moisture Content in Winter Wheat Leaves Using MSC-ANN Algorithm

Keywords: MSC, ANN, Moisture Content, Reflectance spectrum, Quantitative analysis 1

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Instruction

The moisture content is one of the most important factors to effect photosynthesis, respiration and biomass in plant [1][START_REF] Elvidge | Visible and near infared reflectance characteristics of dry plantmaterials[END_REF][START_REF] Mooney | Exchange of materials between terrestrial ecosystems and the atmosphere[END_REF]. Shortage of moisture would effect the physiological process, biochemical process and morphostucture of plant, and further effect the growing, output and quality of crops [4]. For the past few years there are many flood and drought disasters arising, which make the crop cutting into large area.

Therefore, the scientific management of moisture content in crops increasingly becomes one of the most important measures in crop production.

Presently, the moisture content of large area crops has been researched by means of remote sensing, which indicates the deficit of moisture by measuring the temperature of plant canopy [5]. Another common method is to establish an inverse model about reflectance of leaves and moisture content, which bases on the characteristic of moisture absorbing infrared light. From the 80's, many famous scholars have proved that: the spectral reflectance of plant be used to detect status of moisture content is feasible [START_REF] Penuelas | Estimation of Plant Water Concentration by the Reflectance Water index WI (R900/R970)[END_REF][START_REF] Penuelas | Reflectance Indices Indicative of Changes in Water and Pigment Contents of Peanut and Wheat Leaves[END_REF][8][9][10][11]. In this paper we would investigate the mehtod of processing spectrum to promote the result of modeling.

Multiplicative scatter correction(MSC) is a very effective spectrum preprocessing method of dealing with living leaves of plant. It is a multi-variable scatter correction technique [12]. The original spectrum processed with the method can effectively eliminate the scatting caused by the physical factors, enhance the spectral information related to the tested ingredients, reduce the impact of baseline translation and drift between samples caused by the scatting, and can greatly improve the signal-to-noise ratio of correlative spectra.

In this paper we gathered the reflectance spectra of a new collection of winter wheat leaves, then analyzed the spectral data with the method of principal component analysis after the spectra preprocessed with the method of MSC, and finally made study of the model with the training system of artificial neural network. From the result of modeling, we can get the better result by the model of MSC-ANN, especially dealing with the living leaves in the field.

The algorithm 2.1 The algorithm of MSC

The use of MSC firstly requires to create a "ideal spectrum" of the tested samples, which is used to correct the others as a standard spectrum [13]. Actually, we usually calculate the average spectrum of all the samples as the standard spectrum, operate each sample spectrum with the method of unary linear regression to get their regression constant and regression coefficients, and then correct the original spectra of each sample to improve the signal-to-noise ratio. The concrete process of arithmetic is as follow:

(1) The average spectrum:

n n i j i j i A A    1 , , (2)
Unary linear regression:

b m A A i i i i   (3) MSC: m b A A i i i MSC i ) ( ) (  
In the formula, i is the number of samples, j is the number of wavelength.

The algorithm of ANN

The neural network developed in recent years is a very active interdiscipline, which refers to many disciplines such as biology, electronics, computers, mathematics, physics and so on, and it has a extensive future [14]. Artificial neural network is an effective means of researching spectroscopy, which is widely applied in the study of spectrum and NMR. Currently the most active model is the back propagation(BP) algorithm model, proposed by Rumelhart [START_REF] Tang | Academic Periodical of Farm Products Processing[END_REF], which involves parallel distributed processing theory and multi-layer network. BP neural network model is a powerful learning system to achieve the nonlinear mapping between input and output in high degree [START_REF] Bao | [END_REF], and it has proved the model can approximate to any continuous nonlinear curves [17]. The preprocessor spectra is operated with the method of principal component analysis to get appropriate principal components, which are used as the input of multi-layer feedforward neural network, and finally establishing the training model of BP neural network.

Experiment

Test method and instrument

we collected respectively 40 leaves as the samples in booting and milk stages of winter wheat, and then measured their reflectance spectra by using the FieldSpec@Pro FR portable spectrometer, which produced by the U.S. ASD company, and 1800 integrating sphere that produced by LI-COR company. The true values of moisture content were gathered by drying method [18]. We first surveyed the fresh weight of the leaves with AB104N precision electronic balance, and put the leaves in 105℃ to dry until their weight were constant, the moisture content is as follow:

Moisture content(%) = (fresh weight -dry weight) / fresh weight × 100%

The information of sample spectra

Scanned range is from 400nm to 2100nm, the data are collected in each nanometer of the spectral region, so we can get 1700 data points from each spectrum.

Actually, each sample was scanned 5 times, the average spectrum is calculated to be used as the true spectrum of sample. That in figure 1 is three spectrograms choosed.

From the figure we can see there are three characteristic peaks in each spectrum, 1160nm, 1440nm and 1930nm, especially the peak point of the 1440nm is the most notable.

Error Analysis and Model Evaluation

There are several factors in source of the error. The first reason is that the smoothness, cleanliness of leaves and the external environment changing make a great influence on the reflectance spectra. And in addition, water lost in the process of measurement and hard to achieve constant weight in drying attribute to another reasons. These measurement errors accumulate eventually to cause differences in the results of the quantitative analysis. However, the artificial errors can be minimized by MSC-ANN algorithm. From the figure of prediction model, we can see that the prediction values are very close to the chemical values measured by standard method.

So it is feasible to quickly determine the moisture content of winter wheat leaves by NIR reflectance spectroscopy and the model of MSC-ANN algorithm.

Conclusion

In this paper, a new modeling method for quantitative analysis moisture content in winter wheat leaves is studied: Multiplicative scatter correction -artificial neural network(MSC-ANN). In the research, we establish the quantitative analysis model of moisture content in different growing seasons of winter wheat leaves, and then make an analysis of the model. By analyzing the experimental data, MSC-ANN algorithm in dealing with chemical and biological living components is better than the traditional algorithm, according to the source of errors, the measurement errors mainly due to human factors, and another factors are influenced by the noise of instrument and background noise caused by environment changing. These errors can be effectively decreased by the re-training of the artificial neural network, which is better than other algorithm, so the MSC-ANN algorithm proved that it is reliable as a method of modeling.

Fig. 1 4 . 1

 141 Fig.1 Reflectance spectra of winter wheat leaves

Fig. 2

 2 Fig.2 Relationships between predicted and ture values of calibration set in booting stage

Fig. 3

 3 Fig.3 Relationships between predicted and ture values of predicted set in booting stage

Fig. 4

 4 Fig.4 Relationships between predicted and ture values of calibration set in milk stage

Fig. 5

 5 Fig.5 Relationships between predicted and ture values of calibration set in milk stage

1 ACKNOWLEDGMENT. This work was supported by 863 National High-Tech Research and Development Plan (No. 2007AA10Z211).


Jun Luo

Tao Liu

Yande Liu

FT-NIR and Confocal Microscope Raman spectroscopic studies of sesame oil adulteration

Keywords: FT-NIR, Confocal Micro-Raman, Sesame Oil, Adulteration, PCA, PLSDA

Abstract:Adulteration of edible vegetable oils is a very serious problem affecting its commercial value and customers' health. FT-NIR and Confocal Micro-Raman have been explored for discriminating adulteration of sesame oils. The soybean, corn, and peanut oil were mixed into sesame oils in the range of 5~20%(v/v),respectively,and their Confocal Micro-Raman and FT-NIR spectra were collected. Principal component analysis (PCA) and Partial least squares discriminant analysis(PLSDA) were used for the discrimination and classification of adulteration in sesame oil based on spectral data. In addition, the spectra were subjected to the pretreatements (eg, Savitzky-Golay smoothing, Multiplicatnive Scatter Correction, first derivatives and second derivatives) before developing PCA and PLS-DA models. FT-NIR was found to give better efficient in classification of adulteration oils by using PLS-DA with ca.100% classification accuracy, and Raman also gave 100% classification accuracy. The results demonstrated that FT-NIR and Confocal Micro-Raman both can be used to determine the authenticity of edible oils rapidlly.

Introduction

As is well known, authenticity is a very important quality criterion for edible oils, and adulteration is a major problem for the specialist food markets. Many edible oils maybe become polluted through careless handling, or through unintentional or deliberate intervention. Many expensive oils are adulterated with cheaper alternatives to increase profit margins, or expensive ingredients are replaced with a cheaper alternatives. In some cases, this action simply means reducing the quality of the product, but such practices may even have moral and ethical implications. Numerous methods for assessing the authenticity of edible oils have been proposed, such as pyrolysis-mass spectrometry and gas chromatography [START_REF] Yang | Oil fingerprinting analysis using commercial solid phase extraction (SPE) cartridge and gas chromatography-mass spectrometry (GC-MS)[END_REF], focused on the existence or absence of certain chemical compounds in the adulterated sample. However, the  Corresponding author:Jun Luo,Male, Master,Lecturer,Born in 1983,major in Agricultural Mechanization Engineering traditional methods for determination of authenticity based on titration are expensive, complex, laborious, time-consuming, require skilled personnel and sophisticated instrument capability. So, there is a need to develop a simple, cheap, and rapid alternative method to determine adulterants.

NIR spectra assesses chemical structures through the analysis of the molecular bands in the NIR spectrum (eg.C-H) [START_REF] Cozzolino | Verification of silage type using near-infrared spectroscopy combined with multivariate analysis[END_REF], but the peaks in the near-infrared (NIR) region (1,100~2,500nm or 9,091~4,000 cm -1 ) are broad and weak due to the overtone and combination vibrations originated from the chemical composition of samples. NIR spectroscopy was the first applied to determine the presence of adulterants (corn oil, sunflower oil, and raw olive residue oil) in olive oils [START_REF] Wesley | Measurement of Adulteration of Olive Oils by Near-Infrared Spectroscopy[END_REF].Kasemsumran, Sumaporn [START_REF] Kasemsumran | Partial Least Squares Processing of Near-Infrared Spectra for Discrimination and Quantification of Adulterated Olive Oils[END_REF] made a new processing based on partial least squares (PLS) algorithm for the discrimination and determination of adulterants in pure olive oil using nearinfrared (NIR) spectroscopy. Balabin,RomanM. [START_REF] Balabin | Near-infrared (NIR) spectroscopy for motor oil classification: From discriminant analysis to support vector machines[END_REF] tried to evaluate the efficiency of different methods for motor oils classification by using near-infrared (NIR) spectroscopy and SVM and PNN chemometric techniques were found to be the most effective ones for classification of motor oil.

The potential of Raman spectroscopy for the analysis of oils has been recognized for some decades. Raman spectroscopy has already been successfully used to determine important composition or physical structure parameters in a wide range of lipid sample types. Ming-Qiang Zou [START_REF] Zou | Rapid Authentication of Olive Oil Adulteration by Raman Spectrometry[END_REF]established a method based on Raman spectroscopy vibration bands to authenticate genuine/fake olive oil and this method can reliably distinguish the genuine olive oils from the olive oils containing 5% (volume percentage) or more of other edible oils, such as soybean oil. Ru-Hui Weng [START_REF] Weng | Authentication of Camellia oleifera Abel Oil by Near Infrared FourierTransform Raman Spectroscopy[END_REF] determined the authenticity of Camellia oleifera Abel Oil by Near Infrared Fourier Transform Raman Spectroscopy and a linear relationship with a high correlation coefficient (R 2 = 0.9938) was obtained. Rasha M. El-Abassy [START_REF] Rasha | Rapid Determination of Free Fatty Acid in Extra Virgin Olive Oil by Raman Spectroscopy and Multivariate Analysis[END_REF] used a visible Raman spectroscopic method for determining the free fatty acid (FFA) content of extra virgin olive oil with the aid of multivariate analysis.

This work shows the results of using Confocal Micro-Raman and FT-NIR spectroscopy combined with principal component analysis (PCA) and partial least squares discriminant analysis (PLS-DA) to detect the classification of adulteration in sesame oil and compares the performance of two spectroscopic techniques.

2

Materials and methods

Sample description

The samples used in this investigation were sesame (Shenzhen, Jinlongyu brand), soybean (Shenzhen,jinlongyu brand), corn (Shanxi, Xiwang brand) and peanut oil (Shenzhen,jinlongyu brand) that obtained from a local supermarkets. The concentration of adulteration was in the range of 5~20% by weight respectively, and the same was repeated with each concentration. All weighing was done by use of a analytic balance AUX120 (SHIMADZU, Germany) with the accuracy of 0.0001g. Totally, 66 adulterations of sesame oil mixed with soybean oil(n=22), corn oil(n=22) or peanut oil(n=22) were prepared, and their spectra were collected.

FT-NIR analysis

The NIR measurements were performed using transmittance mode in the region of 12,000~4,000 cm -1 by means of TENSOR37 FT-NIR spectrometer(Bruker Optics Inc Germany), equipped with the RT-InGaAs detector. All the spectral data were collected with a 4 cm -1 spectral resolution and a total of 32 scans were co-added to ensure a sufficient signal-to-noise ratio. The oil samples were scanned in a rectangular quartz cuvettes. The quartz cuvette was cleaned with pure ethyl alcohol after successive measurements and dried to ensure the best possible sample spectra. To avoid any spectral variation due to some factors, air background spectrum was collected before each sample spectrum.

Raman analysis

The Confocal Micro-Raman Spectrometer ( BRUKER, SENTERRA ) combined with a spectrometer module and a microscope module was used in this experiment .

The spectrometer module was equipped with a CCD detector (thermo-electrically cooled, 1025×256 pixel) and a semiconductor laser operating at of 785 nm excitation (50mW at sample). The total spectral range (90~3,500 cm -1 ) was recorded within 15 s of total integration time at 9 cm -1 resolution and every measurement consisted of five averaged signal accumulations each with an actual exposure time of 3 s. The oil samples analyzed were poured into a silver spoon and placed on the sample stage (motorized stage 100mm×80mm) .

Discriminant analysis

The OPUS software (verion. 6.5; Bruker Optik GmbH, Germany) was employed for the spectral data collection. PCA and PLS-DA were performed by using the Unscrambler software (verion. 10.1; CAMO AS, Trondheim, Norway). In order to minimize the effect of light scattering, The spectra were subjected to Savitzky-Golay smoothing (7-point, 2orders),Multiplicative Scatter Correction, first derivatives(7point Savitsky -Golay filter) and second derivatives (7-point Savitsky-Golay filter) before developing PCA and PLS-DA models for the discrimination. For the PLS-DA analysis, a set of values were given to the spectra of examined adulterant arbitrarily for the dummy variable "species", such as corn=1.0, soybean=2.0, peanut=3.0. All samples with an error around predicted value more than or equal to 0.5 were classified as not being members of the examined adulterant type. On the other hand, those with an error value lower than 0.5 were classified as being members of the adulterant type.

The results are presented as percentage , the expression is given as correct discrimination rate (CDR) of models.

N N-n CDR 

(1)

where N represents total number of calibration (validation)samples for each model, n is the number of erroneous judgement of samples in every model.

Principal component analysis (PCA) is a nondirected, nonsupervised method of data reduction, which determines the main sources of variation within the spectra, allowing transformation of several related data points into a single number [START_REF] Beattie • Steven | Classification of Adipose Tissue Species using Raman Spectroscopy[END_REF]. Discriminant analysis is a supervised classification technique where the number of groups and the samples that belong to each group are previously defined [START_REF] Cozzolino | Verification of silage type using near-infrared spectroscopy combined with multivariate analysis[END_REF]. Partial least squares discriminant analysis (PLS-DA) differs from PCA in that each source of variation is not ranked in relation to the overall variation, but according to how it accounts for variation between data subsets(groups).

For chemometric analysis, five replications of each sample were used for calibration and validation models. A total of 110 spectra of NIR or Raman were acquired for each adulterant. The same sets of samples were used for enabling comparison of the spectral techniques. For each sample, four of 5 measurements were used for the development of calibration model and rest for validation. Hence a total of 264 individual samples were used for calibration and 66 for validation.

Results and discussion

As we know, every oil differs in composition, length and unsaturated degree of the fatty acids as well as their positions in the chain. Raman spectra represent a combined fingerprint pattern unique to each oil and are used for discriminant analysis. Many analysis were conducted using regions specific to the C=C bond vibration. The double bond considered as an unsaturated bond and the intensities, areas or heights of its peak in this region might indicate the degree of unsaturation in fatty acids, sterols, and vitamins.

FT-NIR spectra

Fig. 1 shows some FT-NIR spectra of sesame oil adulterated with corn oil、soybean oil and peanut oil in the 12,000~4,500 cm -1 region. The spectra shows three major band groups at around 8,263, 7,186, and 5,788 cm -1 . The region between 5,500 and 6,000cm -1 is the first overtone of the C-H stretching from-CH 2 , -CH 3 and -CH=CH-functional groups. The region between 7,000~7,400 cm -1 contains bands assigned to the combination modes of the CH stretching and CH deformation vibration modes. The region between 8,000-8,500 cm -1 consists of bands due to the second overtone of CH stretching modes of the CH 3 and CH 2 groups stretching modes. It is noted in Fig. 1 that those mean NIR spectra are very similar to each other, and it is not easy to classify the spectra into each adulterant type. Therefore, the PCA and PLS-DA method are applied in the study to identify the unknown adulterations in sesame oil samples. An overall PCA was carried out using all samples of each oil. Fig 2 shows the classification results of the PCA on data sets containing three kinds of the oil spectra in score plots. The PCA model was achieved using 1PCs and explaining 97% of variance in the spectral data. The plots shows that the peanut oil is clearly separated from the others. But corn oil has so many portion in common with soybean oil, which gives no boundaries. As with all non-supervised methods for classification, the plot built by PCA reveals and explains the patterns of the samples studied, but the assignation of new samples to a class may be a difficult task, as the method does not calculate a rule to generate boundaries or regions for the groups of samples obtained. Hence, supervised method (PLS-DA) was introduced in order to improve the results. The statistical results of PLS-DA based on FT-NIR spectra are summarized in Table 1. The discrimination models are developed by using two selected spectral ranges and five spectra types: original and four pretreatment spectra. The rate of correct classification obtained from discriminant analysis of the spectra between 6,410.4 and 8,917.4 cm -1 is found to be better than the data from the spectral region between 4,000 and 12,000 cm -1 . In the region 6,410.4~8,917.4 cm -1 , the models using any pretreatment yield 100% correct classification of the prediction set (%CP), especially the second derivative spectra with only the lowest number of PLS factor 2. As is well-known, a small number of PLS factor required for developing a model that provides acceptable results indicates the stability and accuracy of such model [START_REF] Martens | Multivariate Calibration[END_REF].However, when the wide spectral range(4,000~12,000 cm -1 ) is used, the model developed by using the second derivative spectra shows the worst classification performance. This is because the wide spectral range includes a region that contains more noise than relevant information to the model development for the concentration of an analyte. Therefore, the optimal models could be selected from the models giving the best predictive accuracy with the smallest number of PLS factor shows that many major bands of the oleic acid were detected in the Raman spectra of the oil samples, and are attributed to the main components in the oil, which are fatty acids. It should be noted that the region between 1,600 and 1,700 cm -1 contains the stretching vibration of the C=C bond which is very distinctive around 1,655 cm -1 in Raman spectra. It demonstrates that spectra can display non-polar groups of samples which can not be observed in the FT-NIR spectra. According to Fig 4, we can see that the discriminant models built with Raman spectra by PCA are the same weak to decide for samples as the models built with FT-NIR spectra. Although corn oil can be separated from others successfully,giving a value of good prediction 100%, Soybean and peanut oil are both distributed in the same circle. Therefore, we need better methods to make models perfect.

Fig.4: PCA plots for classification of oils using Raman

The spectra regions (90~1,795, 2,818.5~3,159 and 90~3,500cm -1 ) were used for discriminant analysis. The results obtained from PLS-DA for the adulteration discrimination in sesame oil are summarized in Table 2. The rate for correct classification with PLS-DA methods was found to be more than 95.5% for both the calibrationand validation samples due to the degree of unsaturation. It is noted in Table 4 that the excellent prediction results are obtained mostly from the models built with the regions of 90~1,795 and 2,818.5~3,159cm -1 at the same time. It also can be seen that when the spectral data are pretreated by MSC, the model gives the highest rate(100%) for correct classification. However, the results from Confocal Micro-Raman measurements were not as good as these from FT-NIR. Vibrations from polar groups, such as C=O and O-H, are very weak in Raman spectra, which result in the missing of some important vibrations and these vibrations maybe play very important roles in discriminant analysis. 4

Conclusions

This study compared Confocal Micro-Raman and FT-NIR spectroscopy for detecting the classification of adulteration in sesame oil. The results revealed that either of the spectral methods is easier to perform and the FT-NIR method was found to be more efficient. However, Raman spectroscopy techniques could provide exquisite structural insights into functional groups of oils for discriminant analysis. In addition, it is clear from the results above section that the PLS-DA is a very powerful method to identify adulterant types in the adulteration of sesame oil samples

1 Fig. 1 :Fig. 2 :

 112 Fig.1:FT-NIR spectra in the 12,000~4,500cm -1 region of sesame oil with adulteration

  -Raman spectra of pure sesame oil are shown in Fig 3. The Fig 3

Fig. 3 :

 3 Fig.3: Raman spectra of pure sesame oil in the 90~3,500cm -1 region

Table 1

 1 Discrimination analysis of FT-NIR spectra of the adulterant types in sesame oil

		pretreatment	Factors	% CDRC a	% CDRV b
	6,410.4-	None	6	100	100
	8,917.4cm -1	smoothing	6	100	100
		MSC	4	100	100
		1st	3	100	100
		2nd	2	100	100
	4,000-	None	10	99.6	90.9
	12,000cm -1	smoothing	10	99.6	89.4
		MSC	10	99.6	95.5
		1st	4	86.4	63.6
		2nd	5	41.7	45.5
	CDRC a =correct discrimination rate of calibration, CDRV b = correct discrimination rate of
	validation

Table 2

 2 Discrimination analysis of Raman spectra of the adulterant types in sesame oil CDRC a =correct discrimination rate of calibration, CDRV b =correct discrimination rate of validation

		pretreatment	Factors	%CDRC a	%CDRV b
	90~1,795	None	10	100	96.9
	2,818.5~3,159cm -1	smoothing	9	99.2	97.0
		MSC	10	100	100
		1st	5	96.7	96.7
		2nd	7	100	97.0
	90~3,500cm -1	None	11	100	98.5
		smoothing	10	99.2	98.5
		MSC	9	99.6	100
		1st	6	96.2	95.5
		2nd	5	99.2	95.5

ACKNOWLEDGEMENTS

The authors gratefully acknowledge the financial support provided by National Science and Technology Support Program (31160250,61178036), Natural Science Foundation of Jiangxi Province (2008GQN0029, 2007GZN0266), Special Science and Technology-Support Program for Foreign Science and Technology Cooperation Plan (2009BHB15200), Technological expertise and academic leaders training program of Jiangxi Province (2009DD00700).


Guangrong Shen

Xiumei Huang

Pei Zhou

Lumei Wang

Yuee Zhi

Non-Point Source Pollution Characteristics of Agriculture-derived Nitrogen in Groundwater in Suburban Area of Shanghai Based on Models

Keywords: DNDC model, L-THIA model, nitrogen content, NPS pollution, peach orchard

This paper analyzes the characteristics of the NPS nitrogen pollution of the groundwater in Guoyuan Village, Pudong District, Shanghai. And the associated effects on the surface and groundwater around the study area are discussed in detail based on the successive observed data, DNDC and L-THIA model. The results show that both of the surface and groundwater are polluted so seriously that they are not suitable to drink. The average content of total nitrogen in surface water is 6.3 mg/L and 16.85 mg/L in the groundwater, and both of them attribute to Grade Ⅴ surface water standard (≤2.0 mg/L) according to the national standard (GB 3838-2002). It is concluded that the nitrogen pollution comes mainly from the fertilizer of the peach orchard based on the further modeling analysis. Therefore, reasonable adjustment of fertilization measures and project may be an effective and practical approach to control the nitrogen pollution around the peach orchard.

Introduction

Non-point source (NPS) pollution is the main factor of surface and groundwater pollution as point source pollution is becoming under control. Among all the non-point pollution, the agriculture-derived nitrate NPS pollution of groundwater has become an environmental issue which may cause algal bloom and eutrophication in aquifers, and even produce potential hazards to human health though the link is still disputable [START_REF] Liu | Regional differentiation of non-point source pollution of agriculture-derived nitrate nitrogen in groundwater in northern China[END_REF]. The heavy use of nitrogen (N) fertilization for intensive farming and cropping system with low N use efficiency is often responsible for nitrate overloading into groundwater. The role of agricultural NPS in water quality degradation is the most prominent [START_REF] Boers | Nutrient emissions from a culture in Netherlands: causes and remedies[END_REF][START_REF] Tim | Evaluating agriculture NPS pollution using integrated geographic information systems and hydrologic/water quality mode1[END_REF], especially in the suburbs of Shanghai where ecological agriculture is predominant. In 2001, the application of N in fields was up to 612.75 kg/hm 2 , which is 59.6% higher than the average level in China, and with 1.3×10 6 kg N•a -1 flowing to the aquifers [START_REF] Zong-Hao | Soil pollution countermeasure in Nanhui district[END_REF]. Characteristics of agriculture-derived NPS nitrogen pollution and the load estimation are the important part for regional water quality control. Concerning studies on the pollution load include many aspects, such as estimation of the relationship between total precipitation and pollution load, estimating NPS pollution load of the long-term average output by using the model (L-THIA) [START_REF] Ding-Qiang | Law of NPS pollutants losses in a typical small watershed of Dongjiang Drainage Basin of Guangdong Province[END_REF], farmland runoff pollution load model [START_REF] Jian-Yu | Research on NPS nitrogen pollution model in farmland of Shanghai suburbs[END_REF], the use of NPS Watershed SWAT model system for different scenarios of watershed management [START_REF] Bao-Lin | Identification of NPS pollution in the Miyun Reservoir watershed[END_REF] and agricultural NPS pollution load. There are some studies on the characteristics of agricultural NPS pollution focusing on analyzing the factors such as the land use types, topography and its associated impacts on the NPS pollution. The spatial and temporal variability of agriculture-derived NPS pollution is the subject of intensive research recently.

In this study, the characteristics of the agriculture-derived NPS nitrogen pollution of groundwater is analyzed based on the biogeochemical process model DNDC, the L-THIA model and the simultaneous observation data in soil, surface and groundwater. Guoyuan Village, Xinchang Town of Pudong District, located in the southern Shanghai, China, traditionally a high-yield area for growing peaches, is selected as the study area. We aim at exploring the optimum efficiencies of water and fertilizer usage in peach orchard by analyzing the relationship between agricultural management measures and agriculture-derived NPS nitrogen pollution. Results are discussed in terms of the multi-objective management for peach orchard in order to reduce or avoid the nitrogen pollution in groundwater.

Study area

The research is conducted in the peach orchard of Guoyuan Village (31°03'N, 121°39'E), located in Xinchang Town of Pudong District on the east of Huangpu River in Shanghai, China (Fig. 1(a)). The site situates in a subtropical monsoon zone with 4 m altitude. The average annual temperature is 15.7 ℃ and 1100 mm of the average annual rainfall, ranging from 300 to 1300 mm with most falling during the period from July to September [START_REF] Hong-Yan | Characteristics of phosphorus NPS pollution on groundwater in urban agricultural region of Shanghai[END_REF]. There has been a long history of peach planting in this area since 1621. The planting area reached 5336 hm 2 in 2005 and peach became the symbol of this district [9]. The study area, which is in roughly rectangular shape with the area of about 33.35 hm 2 , is adjacent to living area of Guoyuan Village (Fig. 1(b)). Peach trees are planted at a density of an amount of 900 per hectare. The peach orchard is surrounded by river and some ditches distribute regularly inside. A main ditch from east to west divides the peach orchard into northern and southern parts. The north of the main ditch is segmented vertically by 11 ditches with equal space while several ditches distribute irregularly on the other side (Fig. 1(b)). The sources of the water pollution are mainly from the peach orchard and residential area because there are no other land use types in the district. The compound planting mode that green vegetables are planted under the peach trees when the peach is ripe (from June to August) is employed in the peach orchard. 70% of the peach orchard is interplant one-season vegetables, others of two-season vegetables [START_REF] Jian-Fu | Benefit analysis of fruit--crops compound management model--take peach-brassica chinensis compound model in Shanghai city Nanhui area as example[END_REF]. This compound planting and management mode in study area increases soil nitrogen content and associates with the surface and groundwater pollution. The location of Guoyuan Village & ditches in peach orchard 3 Material and method

Data Collection

Measurements of N in surface, groundwater and soil have been carried out since June, 2009. The monitoring locations in soil, surface and groundwater were chosen respectively according to their different features. GPS was used to fix the monitoring location. It is worthy of mentioning that the monitoring locations of surface and groundwater were collected precisely according to the topography and the flowing direction of the ditches to make sure that monitoring locations can present the spread of nitrogen in the soil runoff as much as possible (Fig. 2).

 Sixteen soil samples were collected from different depth (0 ~ 15 cm, 15 ~ 30 cm, 30 ~ 45 cm) in June, 2009. Total nitrogen in soil was analyzed by conventional method and nitrate (in terms of nitrogen) by Bremner method.

Model Parameters

Field management including fertilizer application and tillage system can influence soil fertility directly through fertilizer inputs and nitrogen circulation in agro-ecological process. Models can simulate the processes responsible for production, consumption and transport of N 2 O in both long and short terms, and also allow spatial simulation [START_REF] Willams | NO x and N 2 O emissions from soil[END_REF]. Simulation mechanistic models consider all factors affecting N 2 O circulation in the agriculture ecosystem. In this study, these factors include the planting mode, fertilizer application, irrigation strategy and the tillage measures of peach orchard in recent ten years (Tab.1), which were acquired through field survey, visiting local farmers and references [START_REF] Hong-Yan | Characteristics of phosphorus NPS pollution on groundwater in urban agricultural region of Shanghai[END_REF]. Fertilizer (compound fertilizer 545 kg/hm 2 ) (The N concentration of the compound fertilizer is 15%, the organic material content in the organic fertilizer is 25%, and the N concentration is 0.45%.)

In addition, the spatial data such as the land use maps of different years, soil type maps of Guoyuan village are necessary, which are used to calculate the pollution load in agricultural runoff based on hydrologic impact assessment model. These data were provided by local village committee or obtained through remote sensing image interpretation. Meteorological data such as the temperature, rainfall data of recent thirty years were provided by the weather station of Pudong district.

Study approach

The study area is a relatively independent agro-ecosystem. The compound management mode and agricultural management practices directly impact on the soil-water NPS nitrogen pollution in the area around peach orchard. Modeling allows the complex links among soil physical, chemical and microbial processes that underpin nitrification, de-nitrification and decomposition to be examined. Therefore, this paper quantitatively calculates and analyzes the mechanism of nitrogen circulation in soil based on DNDC model, and estimates the nitrogen content in soil runoff by using NPS pollution load model L-THIA. The measure values of monitoring locations are compared with the simulating values in order to analyze in detail the mechanism of nitrogen dynamic change in agro-ecosystem and the relationship between the compound management model and the pollution characteristics of the surface and groundwater.

Biogeochemical model

The DNDC model is a mechanism model, which is used to simulate and assess the emission of N 2 O, NO, N 2 and CO 2 from agro-ecosystem [START_REF] Jian-Fu | Benefit analysis of fruit--crops compound management model--take peach-brassica chinensis compound model in Shanghai city Nanhui area as example[END_REF][START_REF] Willams | NO x and N 2 O emissions from soil[END_REF][START_REF] Krobel | Modeling water dynamics with DNDC and DAISY in a soil of the North China Plain: A comparative study[END_REF][START_REF] Smith | A comparison of performance of nine soil organic matter models using datasets from seven long term experiments[END_REF]. Although this model was initially created according to the growing conditions of America, it has been verified in many agro-ecosystems such as grassland, farmland and forest ecosystem in many countries such as America, Canada, India, European, China and New Zealand and so on. DNDC model, which is considered as one of the best biogeochemical models [START_REF] Li-Gang | Apply DNDC model to analysis long-term effect of soil organic carbon content under different fertilization and plough mode[END_REF], is able to simulate the emission of a fixed point and also can assess the emission of a region. The model has reasonable data requirement and is suitable for simulation at appropriate temporal and spatial scales. Tab.2 shows the parameters and associated values required by DNDC model in this study. DNDC model can simulate the biogeochemical process under a certain management of production and learn about soil nitrogen pollution and the nitrogen flux during the process. It also can analyze the effects of different management measures such as fertilization and tillage on soil organic carbon content [START_REF] Rong-Ta | Potential relationship between soil nitrogen & phosphorus content and water pollution in Dianchi Lake[END_REF]. This provides the basis for further quantitative analysis of nitrogen load in agricultural runoff.

NPS pollution load model

Agricultural runoff is the key factor to quantify and simulate agricultural NPS pollution. The concentration of output pollutants in runoff and the pollution load per unit area are related with the regional climatic conditions [START_REF] Tao | Spatial distribution characteristics of non-point pollution of Miyun reservoir areas based on L-TH IA model[END_REF]. Long-term hydrologic impact assessment model (L-THIA), which is integrated with GIS, can estimate the average annual runoff and NPS pollution load based on land use type, soil type and daily rainfall for at least thirty consecutive years [START_REF] Kyoung | Development of the Long Term Hydrologic Impact Assessment (LTHIA) WWW Systems[END_REF].

The adjacent western and southern areas of the study area are also peach orchards, and the northern and the eastern areas are residential areas. Survey and experimental data show that the soil nitrogen is carried over by rainfall and irrigation water to enter the river via surface runoff, resulting in different levels of water eutrophication. The nitrogen content in the runoff is analyzed via simulating by L-THIA model. The migration of nitrogen is estimated comprehensively and quantitatively by using DNDC model. The results are tested and compared with the measured values to explore the mechanism of the nitrogen circulation in agro-ecosystem, and learn about the nitrogen pollution load in surface and groundwater under current field management practices and local climate conditions.

Results and Discussion

The analysis of Water samples

The total nitrogen and nitrate nitrogen contents measured in water of the peach orchard are showed in Fig. 3. Obviously, the nitrogen content is widely different between surface and groundwater. The nitrogen content in groundwater is about 1 to 2 times higher than that in surface water. The standard error of nitrogen content in surface water is smaller than that in groundwater, which represents that the nitrogen content in surface water changed little among all the samples. In particular, as shown in Fig. 3, total nitrogen in surface water measured in May 2010 (14.57 mg/L) is significantly higher than that in other times, which is probably due to large amount of fertilizer application in May. Besides, the rainfall of the day when we sampled reached 45.5 mm (22 nd May). However, the variance of nitrogen content in groundwater is very distinct: the total nitrogen content ranges from 0.24 to 49.23 mg/L; the nitrate content from 0.18 to 45.41 mg/L. Furthermore, we learn about the spatial variation features of nitrogen content in water body reference to Fig. 2 and Fig. 3. The nitrogen content in surface water samples around the orchard is relatively low, while those near the main ditch are higher. This distribution pattern of nitrogen content indicates that fertilization probably is the main cause resulting in superfluous nitrogen in water. As for groundwater, the nitrogen content is lower in residential area and higher in the orchard. These characteristics are in accordance with that agricultural production activities contribute to NPS nitrogen pollution in water body.

The measured results state that monthly average total nitrogen content of surface water reaches 6.3 mg/L, higher than the value of surface water Grade Ⅴ (≤ 2.0 mg/L) according to Chinese surface water quality standard (GB 3838-2002) (in terms of total nitrogen), which means the surface water in study area is polluted seriously and belongs to surface water Grade Ⅴ based on the current water quality assessment method in China [START_REF] Ke-Wang | Groundwater quality assessment of urban drinking water sources in China[END_REF]. For groundwater, 16.85 mg/L of monthly average total nitrogen is also much higher than that of the surface water Grade Ⅴ(≤ 2.0 mg/L). The monthly average nitrate nitrogen in groundwater is 15.24 mg/L which falls in to groundwater Grade Ⅲ (≤ 20 mg/L) according to Chinese groundwater quality standard (GB/T14848-93) (in terms of nitrate). However, 20% of groundwater samples belong to water Grade Ⅴ (> 30 mg/L) in terms of nitrate nitrogen content. It is summarized that groundwater in the district belongs to groundwater GradeⅤ, far lower than surface water GradeⅤ. The water pollution in this area is serious, which is consistent with the conclusion by Huang Hongyan [START_REF] Hong-Yan | Characteristics of phosphorus NPS pollution on groundwater in urban agricultural region of Shanghai[END_REF].

Characteristics of NPS soil nitrogen pollution

To explore the characteristics of NPS nitrogen pollution in the peach orchard, DNDC and L-THIA model are employed to simulate C, N circulation in this agro-ecosystem. Climate, management measures such as fertilization application and soil input variables for DNDC are listed in Tab.2. Field data used for all of the variables are listed except for clay fraction and depth of the soil water retention layer. Here default values are used. In this study, the nitrogen input that come mainly from the organic fertilizer, fertilizer, nitrogen sedimentation in rainfall, totalizes to 548.53 kg N/(hm 2 •a) according to DNDC. The output nitrogen is 339.48 kg N/(hm 2 •a) in this agro-ecosystem including crop absorption, the volatilization of NH 3 , N 2 O, NO, N 2 . So, it is calculated that the yearly nitrogen flux change or the net increased nitrogen in soil is about 209.05 kg N/ (hm 2 •a) based on simulation.

In addition, the simulated results of DNDC model show that average total nitrogen content from different soil depths is 1887.2 kg/hm 2 , while the field measured value in June 2010 is 1728.7 kg/hm 2 . The fitness of the two parts is up to 91.6%, which suggests that the simulated results of the C, N circulation in peach orchard by DNDC is credible and valuable for quantificational C, N monitoring.

L-THIA model is used to simulate the runoff depth, total nitrogen & nitrate nitrogen load based mainly on the local climate and land use type data. Fig. 4 displays the runoff depth map, total nitrogen & nitrate nitrogen load maps of the peach orchard derived from L-THIA. Simulations show that the average annual runoff depth of the peach orchard is 190 mm. As it shows in Fig. 4, the maximum nitrogen pollution load is in orchard when compared with other land use types, where the output of total nitrogen is 0.0085 kg/cell (each cell is 3 m * 3 m), or 9.44 kg/hm 2 , and the output of nitrate nitrogen is 0.0031 kg/cell, or 3.44 kg/hm 2 .

The estimated values from L-THIA model represent the average nitrogen output in a region for a long time, reflecting the average nitrogen pollution level. While the measured values of nitrogen in surface water depend on the crop growing stages, farm management practices and climatic conditions. Comparing the measured values with the simulation results from L-THIA model, the results indicate that simulated total nitrogen content in runoff accounts for 78.08% of measured values in surface water, and nitrate nitrogen covers 52.28% of measured values. Hence, it may be reasonable that there are other nitrogen sources affecting the surface water pollution besides the soil we have known. In addition, the compound management mode of interplanting vegetables in peach orchard promotes the soil nitrogen circulation and increases nitrogen content to some extent. However, the mode is not considered in L-THIA model, which may result in the simulated values littler than it is. Shanghai, situated in the Yangtze River Delta region, owns a plentiful water source and its soil is with higher nitrogen content, all of which cause serious NPS nitrogen pollution [26]. The average annual rainfall is 1100 mm with most falling from July to September in Shanghai and nitrate could descend 1 cm in soil with about per 2 ~ 3 mm precipitation [START_REF] Yan-Lin | China N NPS pollution of farmland and empirical estimation: Ⅳ various types of pollution levels and trends[END_REF]. According to this rule, 209.05 kg N/(hm 2 •a) of the net increased nitrogen content in the peach orchard soil would go down 330 ~ 550 cm at differently speed. In this study, the samples in groundwater are from the wells whose depth is 4 m. Therefore, the descending nitrogen in soil from agricultural activities might be one of the main reasons resulting in the higher nitrogen concentration in groundwater. Furthermore, taking 11% as the soil seepage and runoff nitrogen loss occupying total nitrogen input [START_REF] Jian-Fu | Benefit analysis of fruit--crops compound management model--take peach-brassica chinensis compound model in Shanghai city Nanhui area as example[END_REF], the detailed analysis based on L-THIA model and DNDC simulation shows that 1.7% of annual net increased nitrogen goes to the surface water via the soil runoff, 5.8% to the groundwater via soil leaching loss. Thus, 3.5% is cumulative in the soil leachate. In this case, here about 19.2 kg N/(hm 2 •a) would enrich nitrogen in groundwater so much as resulting in pollution.

Conclusion and Discussion

The NPS nitrogen pollution in surface and groundwater is quantitatively analyzed in detail in this paper based on the monitoring location measured data, DNDC and L-THIA.The surface water falls into Grade Ⅴ in terms of nitrogen content and the groundwater quality belongs to GradeⅤ surface water standard, while the nitrate content in groundwater belongs to GradeⅤ groundwater standard according to national associated standard,

The fertilizer application in agricultural production is the main pollution source of surface and groundwater in the study area. The mechanism of migration and output of nitrogen in the process is directly related to NPS pollution of the groundwater where nitrogen content variation is wide. It cannot but mention that, there is the same agro-ecological environment around the study area with the same peach orchards and residential areas. The study area is a relatively independent region for the river separates it from the surrounding environment. This study analyzes the nitrogen characteristics of NPS pollution by taking the agro-ecosystem of peach orchard as the case. It is difficult to absolutely divide the boundary between the study area and the conjoint area because of the complexity of the environment of agricultural soil and water, which may affect the accuracy of results. Meanwhile, the same agro-ecosystems surrounding the peach orchard suggests that the results of this study can reflect the status of NPS nitrogen pollution for such agro-ecosystem. In particular, the simulated values, the measured values and the speculated or predicted values of soil nitrogen leakage should include or cover the NPS pollution from some possible fuzzy boundaries. In addition, soil drainage changing with rainfall intensity and duration is one of the main reasons to cause groundwater contamination although the critical cause of nitrogen pollution is heavy fertilizer application. Therefore, to improve fertilization strategy is an effective means to reduce nitrogen pollution. Some measures such as adjusting the fertilization structure, the ratio between compound fertilizer and organic fertilizer, the proportion of nitrogen, phosphorus and potassium in fertilizer and so on should be done to control and stable the amount of nitrogen fertilizer.

Fig. 1 .

 1 Fig. 1.

Fig. 2 .

 2 Fig. 2. The distribution map of monitoring locations

Fig. 3 .

 3 Fig. 3. Nitrogen content in water of the peach orchard

Fig. 4

 4 Fig. 4 The distribution maps of runoff depth, TN & NO 3 -N

Table 1 .

 1 The management measures of the peach orchard

	Month	Management
	1	Fertilizer (compound fertilizer 750 kg/hm 2 ,
		organic fertilizer 1.5×10 4 kg/hm 2 ), deep tillage
	3	Fertilizer (compound fertilizer 600 kg/hm 2 )
	5	Fertilizer (compound fertilizer 600 kg/hm 2 ), shallow tillage
	9~11	Interplanting green vegetable
	10	Irrigation, deep tillage
	11

Table 2 .

 2 The parameters and partial values in DNDC model

		Parameter		Value
	Climate	Temperature	&	Daily data (omission)
	information	precipitation
		N content in rainfall (mg	2.53
		N/L)
		NH3 content in the air (ug	0.06
		N/m 3 )
	Soil	CO 2 content (ppm) Bulk density (g/cm 3 )		350 1.33
	background	Texture		clay

Acknowledgment This study is founded by Shanghai Leading Academic Discipline (B209), Foundation of science and engineering cross of Shanghai Jiaotong University(2010) and Projects (08dz1900404) supported by Major project of Shanghai and Technology Commission.