Hongyan Gao

Hanping Xiaodong Mao

Zhang

Inspection of Lettuce Water Stress Based on Multisensor Information Fusion Technology

Keywords: Lettuce, Water stress, Information fusion 1.Introduction

Characteristics of reflection spectrum, multi-spectral images and temperature of lettuce canopy were gained to judge the lettuce's water stress condition which could lead to a precise, rapid & stable test of lettuce moisture and enlarged the models' universality. By the extraction of lettuce's multisensor characteristics in 4 different levels, quantitative analysis model of spectrum including 4 characteristic wavelengths, characteristic model of multispectral image and CWSI were established. These multi-sensor characteristics were fused by using the BP artificial neural network. Based on the fused multisensor characteristics, the lettuce moisture evaluation model was established. The results showed that the correlation coefficient of multi-spectral images model, spectral characteristics model and information fusion model were in turn increased, the correlation coefficients were respectively 0.8042、0.8547 and 0.9337. It was feasible to diagnose lettuce water content by using multi-sensor information fusion of reflectance spectroscopy, multi-spectral images and canopy temperature .The correct rate and robustness of the discriminating model from multi-sensor information fusion were better than those of the model from the single-sensor information.

proposes inspection of lettuce water stress based on the multi-sensor information fusion technology.

Multi-sensor information fusion might comprehensively process multi-source information which comes from some different sensors, so it can obtain more accurate and reliable conclusions [START_REF] He | Ning Peng:Multi-sensor data fusion and applications[END_REF]. Multi-sensor information fusion can greatly avoid the limitations of a single sensor and improves the performance of the system [START_REF] Wu Xiao | Near infrared spectroscopy and machine vision information fusion soil moisture detection[END_REF]. Because the sensors provide some uncertain information, multi-sensor information fusion technology is essentially a non-deterministic reasoning and decision-making process [START_REF] Ma | Based on Dempster Shafer evidence of multi-sensor information fusion technology and its application[END_REF].

Multi-sensor information fusion can be divided into three different layers which were the decision layer fusion, feature layer fusion and the raw data layer fusion [START_REF] Yan | Multi-sensor data fusion technology and its application[END_REF]. The research is the optimization and combination of spectra, images, and canopy temperature, because the spectrum analyzer, machine vision systems and other equipments are abstracted into different types of sensors, the features are different and collect different physical quantities, the information pattern and span are comparatively large, environment and objectives with time-varying features, crop characteristics are complex, so it does not suitable for the decision layer fusion and the raw data layer fusion. More practical option is feature layer fusion, the fusion not only retain a sufficient number of original information, but also achieve a level of data compression, contributes to real-time processing [START_REF] Chen | Based on near infrared spectroscopy and machine vision information fusion technology of multi-judge the quality of tea[END_REF].The concept of feature layer fusion is different features gather to form the new summary feature set, and then making decisions accordingly [START_REF]Runsheng Wang:Information fusion[END_REF].

The research obtains lettuce canopy spectrum, image information, canopy temperature and environmental temperature and humidity, etc. Then establish spectra model, mage model and the water stress index model. Ultimately, used BP neural network training samples and verification, water stress conditions on the lettuce is rapidly and non-destructively inspected.

Experimental Design and Sample Training

Instrument and Equipment

Spectrum measuring equipment is the United States ASD FieldSpec ® 3 handheld portable spectrum analyzer, the range 350~2500nm; at 350~1000nm, sampling interval is 1.4nm, resolution is 3nm; in 1000 ~ 2500nm, sampling interval is 2nm, resolution is 10nm. High-precision analytical balance weighs the quality of whole lettuce, accuracy is 0.1mg. Canopy multi-spectral image utilizes MS-3100 multispectral digital progressive scan camera, imaging spectral range is 400~1100nm, resolution is 1039×1392. Canopy temperature utilizes the TI50 infrared thermal imaging instrument, the range is -20~305℃, accuracy is 0.07℃.

Samples Training

Experiment started in Jiangsu university modern agricultural equipment and technology Venlo-type greenhouse. The variety is Italian anti-bolting lettuce. According to the Yamazaki Nutrient solution, the samples were divided into four levels, each level had 12 lettuce, so the predict set had 24 samples, calibration set had 20 samples. Four levels were:Group 1(W1) ensure adequate water were supplied throughout the growing season; Group 2 (W2), 3 (W3), 4 (W4) irrigated the standard formula of 75% 50%, 25% concentrate.

Experiment Design

Feature level fusion of based on BP neural network

Spectral model

Multi-spectral image model CWSI Lettuce samples (four levels)

Diagnostic assessment model of lettuce water

Results and Analysis

The Quantitative Analysis of Lettuce Moisture Content Based on Spectrum Technology

Fig. 2 showed the lettuce canopy reflectance spectra in the different water stress.

Combined with previous studies of discussion group, and referenced to the USDA researchers came to the main biochemistry components of the spectral absorption characteristics [START_REF] Pu | Peng Gong:Hyper spectral remote sensing and its application[END_REF], wavelength sensitive of water-related mostly concentrated in the near infrared band. As Fig. 2 shown, spectral reflectance of lettuce had significant difference in different water stress at the water sensitive bands. Wavelength/nm Spectral reflectance

Fig. 2. Lettuce canopy reflectance spectra under different water stress

In order to eliminate offset and drift which caused the spectral bias, highlight the hidden information and identification of samples. The first order derivative was carried out on the spectrum and conducived to extract characteristic wavelength. In first, all spectral points were divided into 4-sensitive bands: 1220-1300nm ,1410-1490nm ,1620-1700nm ,1900-1970nm, then removed the wavelength by stepwise regression [START_REF] Zheng | Based on stepwise regression of near-infrared spectral information extraction and model[END_REF] and got the sensitive wavelength of related to the lettuce canopy water stress:1267nm、1443nm、1661nm、1921nm. So as to eliminate the impact of multicollinearity, so the four wavelengths for partial least squares regression analysis [START_REF]Huiwen Wang:Partial least squares regression method and its application[END_REF], When two principal component score were extracted, cumulative contribution rate greater than 0.85. Ultimately obtained PLS regression model based on four sensitive wavelengths (X i ): Then 20 samples of spectral data tested the model, the correlation coefficient between dry water content of lettuce canopy measured and predictive value was 0.8547. The six channels lettuce canopy images were simultaneously acquired by the MS-3100 multi-spectral digital progressive scan camera (Fig. 4), they were R、G、B、 RGB、IR and CIR. This method not only contributed to extract image features of all the independent channels, but also easily achieved multi-spectral image pixel-level operation and integration(Image registration was not required).

Fig. 4. Lettuce multi-spectral images

Image processing based on MATLAB software. Median filtering method of 3 × 3 window would eliminated isolated noise points, reduced the image blurring, so it was used for image preprocessing. And two-dimensional maximum entropy segmentation was used for background segmentation; this method preserved more image information of crop canopy and contributed to the image feature extraction. Experiment used gray feature extraction method, finally AIR 810 and AIR 940 (AIR k is near infrared spectroscopy 810nm and 940nm canopy image mean gray value) significantly correlated to lettuce canopy water content. Establishing lettuce canopy water content forecast model by SPSS 13.0 for multiple linear regression analysis, including AIR 810 and AIR 940 image features variable:

Then 20 samples of spectral data tested the model, the correlation coefficient between dry water content of lettuce canopy measured and predictive value was 0.8042.

Canopy Water Stress Index (CWSI) Model Establishment

The TI50 infrared thermal imaging instrument obtained lettuce canopy temperature and real-timely monitored environmental temperature and humidity. According to the CWSI empirical model by Idso in the literature [START_REF] Cui | Based on the temperature of the summer maize canopy water stress index model study[END_REF], as follows:

( ) ( ) ( ) ( ) c a c a c a u l c a T T T T CWSI T T T T Π Π - -- = - -- (3)
( )

c a T T A B VPD Π - = + × (4) 1 7 .2 7 2 3 7 .3 0 .6 1 1 (1 ) 1 0 0 a a T T R H V P D e × + = × × - (5)
( )

c a u l T T A B VPG - = + × (6)
In the formula: Tc-the crop canopy temperature, °C; Ta-air temperature, °C; (Tc-Ta) Π -lower limit of the difference temperature between canopy and air, °C; (Tc-Ta) ullimit of the difference temperature between canopy and air, °C; VPD-air vapor pressure deficit, hPa. A, B-experience factor; VPG-the difference of air saturated vapor pressure between when the air temperature were Ta and Ta+ A, hPa. y = -1.4647x + 4.8287 R 2 = 0.8969 The model of difference temperature between canopy and air by data analysis:

VPD T T a c 4647 . 1 8287 . 4 - = - (7)
In the condition of full water supply, when VPD was 5.86, the difference temperature between canopy and air was minimum, it regarded as lower limit of CWSI. So (Tc-Ta) Π was -3.7544. When lettuce was severe water stress, canopy temperature reached the maximum during the experiment, it regarded as upper limit. So (Tc-Ta) ul was 1.3241(℃).

Lettuce CWSI model:

The Model Based on BP Neural Network Information Fusion

BP neural network has strong fault tolerance, distributed, storage, self-learning, adaptive, self-organization, nonlinear dynamic capabilities and handle complex environments [START_REF] Xiong | Based on neural network spectral recognition system design and analysis[END_REF].Therefore, using better self-learning and adaptive capacity, lettuce water stress condition was evaluated by BP neural network.

The research used a 3-layer structure of BP neural network for feature level data fusion. There were lettuce canopy reflectance spectra characteristic wavelength 1267nm, 1443nm, 1661nm, 1921nm and the image feature parameter AIR 810 , AIR 940 , CWSI as the input, dry basis moisture content measured as the output. Error index and the training step were respectively sited to 0.001 and 0.05, hidden nodes was 10. Using BP neural network to predict the same test set (Table 1). The average relative error of predicted and measured values was 9.24%, correlation coefficient R was 0.9337.

Conclusion

Object of study choose four different moisture content of lettuce, lettuce canopy water stress was evaluated by spectral characteristics, image feature information, canopy temperature and Environmental temperature and humidity. Spectral analysis model and the image model were established, then verifying the model, the correlation coefficient between water content of lettuce canopy measured and predictive value were 0.8547 and 0.8042.

The results showed that lettuce canopy water stress evaluation method based on the spectrum、 multi-spectral image and the CWSI of multi-sensor information fusion technology was feasible, and the correlation coefficient was 0.9337. Model of accuracy and stability were higher than a single information model. Results of the research for multi-sensor information fusion technology could regarded as reference to the rapid and accurate detection of lettuce water.

Fig. 1 .

 1 Fig. 1. The flow chart of lettuce water stress inspection based on multi-sensor information fusion technology

Fig. 3 .

 3 Fig. 3. Lettuce canopy reflectance spectra under different water stress

Fig. 5 .

 5 Fig. 5. the relationship of (Tc-Ta) and VDP

Table 1 . Samples predicted values and the measured values

 1

	Sample	Predictive	Measured	Relatively
	Number	Value (%)	Value (%)	Error (%)
	1	2984.269	2812.366	6.112
	2	2712.595	2508.462	8.138
	3	2440.624	2604.845	6.304
	4	2402.126	2598.000	7.539
	5	2290.337	2499.381	8.364
	6	1934.597	2161.101	10.481
	7	1852.465	2093.596	11.518
	8	1839.795	1721.290	6.885
	9	1791.997	1623.837	10.356
	10	1780.012	1968.000	9.552
	11	1710.706	1603.889	6.660
	12	1661.028	1798.400	7.639
	13	1690.231	1864.688	9.356
	14	1577.302	1704.605	7.468
	15	1603.096	1428.125	12.252
	16	1687.403	1490.943	13.177
	17	1575.115	1449.184	8.690
	18	1436.798	1625.636	11.616
	19	1206.933	1415.647	14.743
	20	1081.009	1174.524	7.962
	Average			9.240

Acknowledgements

This work was supported by a grant from the National High Technology Research and Development Program of China (863 Program)(No.2008AA10Z204 and 2008AA102208),"333 Talent Project" in Jiangsu Province


Min Zhang
email: zhangm@shou.edu.cn

Zhenhua Che

Jiahua Lu

Huizhong Zhao

Jianhua Chen

Zhiyou Zhong

Le Yang

Study on Thermal Conductivities Prediction for Apple Fruit Juice by Using Neural Network

Keywords: Thermal conductivity, Neural network, Apple fruit juice

Based on experimentally measured values by thermal probe method, the prediction model of thermal conductivities of apple fruit juice as a function of concentration and temperature was studied by neural network method. The optimal neural network was made of two hidden layers and every hidden layer had six neurons. The prediction result shows that the optimal model could predict thermal conductivity with a mean relative error of 0.11%, a mean absolute error of 0.00054W/mK, a mean standard error of 0.00039 W/mK, the linear relationship of 0.9993. The calculated precision was higher for BP neural network model than that for dual regression model. The presented results were proved that this model can be used with satisfactory accuracy for the prediction of thermal conductivity of apple fruit juice.

Introduction

Apple fruit juice as a traditional consumption is widely available on the Chinese market. It is a complex mixture of vitamins, pectin, sugar, mineral salts, organic acids andamino acid. Accurate values of thermophysical parameters over a wide temperature range are needed in the fruit juice industry in order to design and optimize the operation of handling and processing units [START_REF] Frandas | Thermal properties of fruit juices as a function of concentration and temperature determined using the photopyroelectric (PPE) method[END_REF]. Thermal conductivity of fruit juice is one of the important thermophysical parameters used to estimate the rate of heat transfer during processes, such as preservation, transportation and freezing. Generally, mathematical models which express the dependence of thermophysical properties on temperature and concentration are a very appealing alternative to experimentation, and a useful tool for the implementation of computer-aided routines for equipment design and process automation [START_REF] Telis-Romero | Thermophysical Properties of Brazilian Orange Juice as Affected by Temperature and Water Content[END_REF].

Over the years an extensive investigation on existing methods of measurement of thermal conductivity of fruits juice has been carried out [START_REF] Frandas | Thermal properties of fruit juices as a function of concentration and temperature determined using the photopyroelectric (PPE) method[END_REF][2] [START_REF] Shamsudin | Thermophysical properties of Thai seedless guava juice as affected by temperature and concentration[END_REF]. Some empirical equations have been proposed for the estimation of thermal conductivity of various fruits juice products as a function of composition [START_REF] Zainal | Effects of temperature on the physical properties of pink guava juice at two different concentrations[END_REF] [START_REF] Moreira Azoube | Effect of concentration on the physical properties of cashew juice[END_REF]. In these equations one can easily verify that thermal conductivity of fruits juice depends strongly on the water content, since water has the highest thermal conductivity of all food components [START_REF] Saravacos | Transport properties in processing of fruits and vegetables[END_REF].

In recent years, artificial neural networks as a kind of analysis method were obtained extensive attentions by food engineer. It has become increasingly popular during the last decade. The advantages of this method compared are its high speed, simplicity and large capacity which reduce engineering effort [START_REF] Kalogirou | Artificial neural networks in the renewable energy systems applications: a review[END_REF]. It can be used successfully in various fields of modeling and prediction in many engineering systems, mathematics, medicine, economics, meteorology and many others [START_REF] Kurt | Prediction of thermal conductivity of ethylene glycolwater solutions by using artificial neural networks[END_REF].

The objective of this study is to develop a BP artificial neural network model to predict the thermal conductivity of apple juice as a function of concentration and temperature and to compare the BP artificial neural network model established with dual regression model by the calculated precisions.

Methodology

The development of a neural network model involves: the generation of data required for training/testing, the training/testing of the BP artificial neural network model, the evaluation of the BP artificial neural network configuration leading to the selection of an optimal configuration, and validation of the optimal BP artificial neural network model with a data set not used in training before [START_REF] Shyam | Neural networks for predicting thermal conductivity of bakery products[END_REF].

Data generation of thermal conductivity of apple juice

The thermal conductivity of apple juice at different concentration and temperature were measured using a tiny thermal probe method. It has been described in detail in our previous studied [START_REF] Zhang | Experimental study on influence factors of thermal conductivities of apple juice[END_REF]. In this study, the thermal conductivity probe system was proposed, which needs short duration of the experiments, simplicity, high accuracy and relatively small sample requirement. The system worked satisfactorily for its application. The effects of concentration in the range of 0-70° Brix and temperature in the range of 3-50℃ on the thermal conductivity of apple juice were determined. The 192 data sets with different thermal conductivities affected by concentrations and temperatures were obtained through using this measurement system. The thermal conductivity of apple juice decreased with an increase in concentration while it increased with an increase in temperature. In addition, the influence of concentration on the thermal conductivity of apple juice was stronger than that of temperature. A double-factors regressive equation of thermal conductivity of apple juice as a function of concentration and temperature was obtained:

s t k 0035 . 0 0012 . 0 582 . 0 - + = (1)
Where k is the thermal conductivity of apple juice, W/m K; s is the concentration of the apple juice (° Brix) and t is the temperature of the apple juice, ℃.

2.2

Establishment of BP neural network

A multi-layer feed forward BP network structure with input, output and hidden layer(s) was used in this study as shown in Figure 1. It consists of the feed-forward pass and backward pass process. The input layer includes two neuron parameters of concentration and temperature. The output layer is the thermal conductivity of apple fruit juice, hidden layer suppose n neurons before training.

The network is processed through two stages: training stage and validation stage [START_REF] Yuan | Prediction for the Effect of Temperature and Water Content on the Soil Specific Heat by BP Neural Network[END_REF]. BP neural network achieves full connectivity between the upper and lower layers, and no connections between neurons in each layer. The neural cells of input and output layers were decided by the practical input and output parameters [START_REF] Zhang | Prediction model of thermal conductivities of fruits and vegetables based on BP neural networks[END_REF]. In the training stage, the network was trained to predict an output based on input data. In the validation stage, the network was validated to stop or continue training. The neurons of hidden layer were decided on the error decline in the training course. The learning process was composed of positive signal transmission and the error backpropagation. The actual output value was calculated by the process information flow from input layer, hidden layer to output layer by layer. If the actual output of the output layer does not match with the expected output, the error back-propagation stage was started [START_REF] Xu | Nonlinear fitting calculation of wood thermal conductivity using neural networks[END_REF]. The error signals of neural cells were based on this error. The each unit weight can be amended by the error signal values until the network output error was reduced to tolerant ranges.

∑ = - = N c p m k k N MAE 1 1 (2) ∑ = - = N c m p m k k k N MRE 1 1 (3) 1 ) ( 1 2 - - = ∑ = N k k SE N c p m (4) ∑ = - - = N c m p m k k k R 1 2 2 2
) ( 1 [START_REF] Moreira Azoube | Effect of concentration on the physical properties of cashew juice[END_REF] Where N is the number of training samples; k m and k p are the actual measurement values and the neural network fitting calculated values of thermal conductivity.

Results and discussion

The thermal conductivity data set of the 144 cases in our study were divided into two groups; training and testing sets. The BP neural network model was trained using 116 randomly selected data (accounting for 80% of the total data) while the remaining 28 data (accounting for 20%) was utilized for testing of the network performance. In order to improve training characteristics, the thermal conductivity data set experimentally obtained were normalized using a simple normalization method. The Maximum and minimum values of the inputs and outputs and normalization values ere given in Table 1. The modeling program was implemented under MATLAB. The BP learning algorithm has been used in feed-forward, which included an input layer, a hidden layer and an output layer network. There are two input parameters in the input layer, namely the temperature difference and concentration of the apple juice. The output is the thermal conductivity of apple juice. The number of neurons in the input layer is equal to the number of input parameters and the number of neurons in the output layer is the number of output parameters. Optimal number of the neurons in the hidden layer was determined by trying different networks [START_REF] Kurt | Prediction of thermal conductivity of ethylene glycolwater solutions by using artificial neural networks[END_REF]. The number of neurons was increased from 1 to 8 based on the trial and error method in the hidden layer.

Input layer and hidden layer use the sigmoid transfer function for the transfer function, the output layer use the purelin function, in order to improve the prediction accuracy, the test results would be normalized, so that concentration

C i ∈[0,1], temperature t i ∈[0,1], thermal conductivity k mi ∈[0,1].
In the process of model, the network is initialized by means of Newff function in Matlab neural network toolbox, the trained model parameters set as follows: the maximum epoch of 10000, network performance goal error of 10 -6 , learning rate of 0.02, momentum constant of 0.9. The errors associated with desired output response are adjusted in the way that reduces these errors in each neuron from the output to the input layer. When the network output error was reduced to the required tolerance, the training automatically stopped. According to equation ( 2) -( 5) the neural network established was training and the error of the fitting calculation results of each network was analyzed. It is found that best structure of the model has six neurons in the first and second hidden layer presented in Table 2. The MRE, MAE, SE and R 2 for this optimal configuration were 0.92%, 0.0047 W/mK, 0.000525 W/mK, and 0.9957, respectively. So the optimized BP neural network model architecture has a configuration of 2-4-4-1 neurons. The prediction thermal conductivity of apple juice for the optimized BP neural network model with a configuration of 2-4-4-1 neurons were compared with the experimental ones as shown in Figure 2 for the 116 training data set. The results demonstrated good agreement between the predicted and the experimental values of thermal conductivity (R 2 =0.9957).

MAE (W/m•K) SE (W/m•K) R 2
The performance of the optimal neural network (two hidden layers and six neurons in each hidden layer) was validated using a smaller data set consisting 28 cases. Correlation between predicted and measured values of thermal conductivity was shown in Figure 3.This network predicted thermal conductivity values with an MRE of 0.11%, MAE of 0.00054 W/mK and SE of 0.00039W/mK, respectively. It still kept a good agreement between the predicted and the measured values of thermal conductivity (R 2 =0.9993). The results presented that this model can be used with satisfactory accuracy for the prediction of thermal conductivity of apple fruit juice. The mean relative error of 28 prediction values of thermal conductivity cases in this model is smaller than that in double-factors regressive equation whose mean relative error is 1.73% [START_REF] Zhang | Experimental study on influence factors of thermal conductivities of apple juice[END_REF].

Conclusion

In this work, the prediction model of thermal conductivities of apple fruit juice as a function of concentration and temperature was studied by neural network method. The optimal neural network was made of two hidden layers and every hidden layer had six neurons. The optimal neural model could predict thermal conductivity with MRE of 0.11%, MAE of 0.00054W/mK, SE of 0.00039 W/mK, R 2 of 0.9993 while the mean relative error in double-factors regressive equation was 1.73%. The prediction precision of BP neural network model was higher than that of dual regression model. The presented results were proved that this neural network model can be used with satisfactory accuracy for the prediction of thermal conductivity of apple fruits juice and further applications.

Fig. 1 .

 1 Fig. 1. Model of BP neural network. The indicators of thermal conductivity of apple fruit juice including the mean absolute error MAE, mean relative error MRE, standard deviation SE, coefficient of determination R 2 can judge the right type of performance of the network and optimize neural network model. The four indicators used to compare the performance of various BP artificial neural network configurations were calculated by the following expressions[9]:

Fig. 2 .

 2 Fig. 2. Comparison of measured and predicted values of thermal conductivity for the training sets.

Fig. 3 .

 3 Fig. 3. Comparison of measured and predicted values of thermal conductivity for the validation sets.

Table 1 .

 1 Input and output parameters ranges and normalization values.

	Parameters	Minimum value	Maximum value	Normalization value
	Concentration,%	0	70	70
	Temperature,℃	3	50	40
	Thermal conductivity, W/mK	0.341	0.647	0.646

Table 2 .

 2 Error parameters with different BP neural network structure.

	Hidden layers	Neurons	MRE (%)

Acknowledgments. Funding for this research was provided by the Natural Science Foundation of China, Grant No 30771245 and Leading Academic Discipline Project of Shanghai Municipal Education Commission, Project No J50704.


Xiaoli Wang

Guomin Zhou

Study on Pretreatment Algorithm of Near Infrared Spectroscopy

Pretreatment of near-infrared spectral data is the basis of feature extraction, quantitative and qualitative analysis and establishment of models, it plays a significant role in obtaining the data and get reliable results. The purpose of the paper is compared the advantages and disadvantages of the S-G, derivative and multiple algorithm methods of spectral preprocessing through the example of apple leaves. S-G algorithm can smooth the data relatively better, but we must according to the specific circumstances of the case while chose the width of the window and the order of polynomial; Kernel smoothing is better than S-G in two-ends data processing, but its processing speed is slower than S-G.

Derivative algorithm can get more stable reflectance, but it is sensitive to noise, so it need to be used with the smoothing algorithm. Multiple scatter correction can be used effectively to eliminate the translation and offset of baseline. All of above algorithms have been applied in the system of near infrared spectroscopy processing system of leaves and satisfactory result was obtained.

Keywords:Near Infrared Spectroscopy;Spectral pretreatment;Smoothing algorithm;Derivative algorithm;Multiple scatter correction algorithms

Introduction

Near infrared spectroscopy is electromagnetic waves whose wavelength is between 780 ~ 2500nm, and it is intermediate between visible light and infrared, NIR refers primarily to low-energy electronic transitions and the clusters containing hydrogen atoms. Various objects' components information can be obtained through the near-infrared spectroscopy.

However, when we survey and evaluate the components information of the objects, the measurements may be quite different for the difference of the weather, environment, humidity, temperature and human factors, as well as the existence of dark current. Near infrared spectroscopy is mainly the spectral information of the measured object, but also includes all sorts of interference and noise from outside. Therefore, the spectral data must be preprocessed before establish the model to reduce noise and interference information, lay the roots for simplifying the operational process of modeling and improve the accuracy of the analysis. The main idea of Savitzky-Golay [START_REF] William | Numerical Recipes in C++: The Art of Scientific Computing[END_REF][START_REF] 法の大ざっぱな説明 | [END_REF][START_REF] Gorry | General Least-Squares Smoothing and Differentiation by the Convolution (Savitzky-Golay) Method [J][END_REF] (later referred to as S-G) filter algorithm is that use higher order polynomial to approximate fit the fixed number of data (also called the active window). For each point i, use the least square method, fit all the data of this activity window to a polynomial expression. Then replace the original value with the value of polynomial (g i ) at point i. When the window moves to the next point i+1, we use least-squares fitting algorithm to fit the original data rather than the data of the polynomial in the new window. Using the S-G algorithm requires that horizontal axis x i have a uniform spacing (x i+1 -x i ≡ ∆ x).

Materials and Methods

experimental materials

We

The core formula of the S-G algorithm:

∑ - = + = r n c g l n n n i n i f , T T T M n A A A e c 1 1 ) ( - + = ( r l n n n ≤ ≤ - ) ;
Here f i is the original data, g i is the data after smoothing, is M+1

1 + M e dimensional unit
vector, n r is the number used to the right, i.e., later than i, n l is the number of points used to the left of a data point i, i.e., earlier. [START_REF] Wand | Kernel Smoothing [M[END_REF] . Same as S-G algorithm, suppose the data after smoothing is g(i), window width is h, weighted function is S(t), then the expression of fitting value in time t (such as the spectrum reservation time) is:

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ - - - - - - - - - - - - = - - - - 1 1 ) 1 ( ) 1 ( ) 1 ( 1 0 0 0 1 ) 1 ( ) 1 ( ) 1 ( 1 ) ( ) ( r 1 r 1 1 1 n n n n n n n n n n n n A m m r r m r m r l m l m l l m l m l L L M M L M M L M M L M M L L ) 1 ( ) 1 ( R + × + + ∈ M n n r l

Kernel smoothing algorithm

∑ = = n j j j y t S t g 1 ) ( ) ( ) (
The weighted function proposed Gasser and Muller is:

* ( ) * ( 1) 1 ( ) ker ( ) j j t j t u t S t n h h - - = ∫
In the two formulas above, kern (x) function is the core of smooth; it has three different weighted functions:

Uniform function: 0.5 1 ker ( ) 0 x n x otherwise ⎧ ≤ ⎪ = ⎨ ⎪ ⎩ Quadratic function: 0.75(1 ) 1 ker ( ) 0 x x n x otherwise ⎧ - ≤ ⎪ = ⎨ ⎪ ⎩ Gaussian function: 1/ 2 2 ker ( ) (2 ) exp( / 2) n x u π - = -

Derivative algorithm

Usually, derivative [START_REF] Ni | Progress and effect of spectral data pretreatment in NIR analytical technique[END_REF] spectrum use different methods to approximate calculate:

First derivative spectrum:

λ λ ρ λ ρ λ ρ ∆ - = - + 2 )] ( ) ( [ ) ( 1 1 ' i i i
Second derivative spectrum:

2 2 2 ' ' 4 )] ( ) ( 2 ) ( [ ) ( λ λ ρ λ ρ λ ρ λ ρ ∆ + - = - + i i i i

Multiple scatter correction

The concrete account step of multiple scatter correction [START_REF] Yong Jun | Research on the Method of Choosing Optimum Wavelengths Combination by Using Multiple Scattering Correction Technique[END_REF] is: First of all, calculate the average spectrum of all samples as a standard spectrum; secondly, fit the spectral data for each sample with standard spectra through one-dimensional linear regression, and then get the offset (regression coefficient) and translational (regression constant) of each spectrum relative to standard spectrum; finally, the original spectra data of each sample minus regression constant then divided by regression coefficient and you get spectrum after correction. Specific process is as follows:

(1)calculate the average spectrum:

n A A n i j i j i ∑ = = 1 , , (2)
the linear regression:

i j i i i b A a A + = (3)multiplicative Scatter Correction: i i i msc i a b A A ) ( ) ( - =

Experiment and Analysis

We use MATLAB software to process data.

After programming and debugging,spectral data of apple leaves were processed and compared and we obtained the following results: We find that the original spectrum, in the 1000nm and 2000-2500nm, has relatively large noise. The spectrum also has large noise after S-G (9, 9, 1) processing. After S-G (99, 99, 1) processing, the smoothing effect is obvious, but many spectral information is lost. With the increase of the fitting polynomial order, the noise is also increased. (55, 55, 1) is the relatively ideal parameter. In practice, parameter selection is very important and we have to analyze specific issues. By the comparison of Table 1, we can find, the discreteness of the data after S-G smooth is smaller than the original data; by the same smoothing window, kernel smoothing effect is better than the S-G smoothing; the effect of smooth is better when the kernel weighting function is not uniform but quadratic; and after repeated treatment, we find that variance and range change a little after 55 smoothing window.

Thus, we selected kernel smoothing algorithm (the window width is 55, the polynomial weighted function is quadratic function) as the smoothing method of apple spectral data.

Derivatives

Discussion and outlook

After the application of the above algorithm to the processing system of Near Infrared Spectroscopy in leaves, we found that:

1、S-G smoothing can be better to smooth the spectral data, while the window width and degree of polynomial selection is very important, the more data points within the window, the more seriously the spectral resolution drops, also the more severely the spectrum distorted; But if too few data points, the smoothing result is not satisfied; the higher polynomial fitting degree, the larger data noise and smoothing effect is not very good, but the low-order polynomial fitting can not fit the data very good; so, the selection of the window width and the degree of polynomial need to be based on the actual situation. There, we choose the (55, 55, 1)as the smooth parameters for apple spectral data.

2、As the same parameters, after the kernel smoothing algorithm processing and then modeling analysis, the predicted result is more accuracy than S-G, but the time complexity is higher, therefore, we should select the smoothing algorithm according to the size of the data in practice. In summary, we can get the following conclusions:

The suitable pretreatment for apple leaf spectral data was found, and provided an important reference for the establishment of models.

But for the difference of the apple leaves' physical and chemical properties, spectroscopic may show some differences, so when we establish the models, we have to change the data preprocessing method at any time according to the different data.

Figure1. 2 . 2 . Data smoothing algorithm 2 . 2 . 1 .

 22221 Figure1. Near infrared spectroscopy of apple leaves

Figure 2 .

 2 Figure 2. S-G smooth Figure 2 is the original spectral and the spectrum which was treated by different S-G smoothing of parameters (except the original spectrum, in order not to overlap spectrum, ordinate value was added by 0.1 in proper order).

Figure 3 .Table 1 .

 31 Figure 3. Kernel smooth Figure 3 is the spectrum which was processed by kernel smoothing (except the original spectrum, in order not to overlap spectrum, ordinate value was added by 0.1 in proper order). We find, compared with the S-G algorithm, kernel algorithm has a better smooth result at the both ends of the data (Note the smoothing effect in about 400nm in Figure 2 and Figure 3). The disadvantage of kernel algorithms is the high time complexity. When the data is big, the processing speed is slow. Standard deviation (Std Dev) is a measure of the degree of data dispersion, indicated the deviation of the data from the mean, it is also a target of variable stability. Range is the difference between the maximum and minimum data which described the dispersion of data. The range is bigger, the dispersion of data is bigger.Table 1. The comparison of the raw data and the data which was processed by smoothing algorithm the raw data

Figure 4 .Figure 5 .Figure 6 .Figure 7 .Figure 4 and 3 . 3 .Figure 8 .Figure 9 .Figure 10 .

 45674338910 Figure 4. The first derivative spectra of the original data

3 、

 3 Derivative algorithm processing results are not very satisfactory. A cause of this phenomenon is that the absorption peaks are clear, but its relatively narrow peaks are not conducive to model diagnosis. 4、Multiple Scatter Correction can improve the modeling accuracy better, and the accuracy of the results can greatly improve.


Lianqing Zhou
email: lianqing@zju.edu.cn

Guiying Pan

Zhou Shi

Methodology Comparison for Effective LAI Retrieving based on Digital Hemispherical Photograph in Rice Canopy

Keywords: Rice Canopy, Digital Hemispherical Photograph, Effective LAI, Methodology Comparison

The paper investigates methodology comparisons for retrieving effective leaf area index (LAI) using digital hemispherical photograph(DHP) in rice canopy. A set of self-making DHP instrument equipped with a fish-eye Len is utilized to acquire DHP in rice canopy, and some self-developing DHP processing procedures are utilized to pre-process DHP and extract effective LAI from DHP (LAI DHP ) rapidly. Based on Beer-Lambert's law and gap fraction that computed from DHP, four methods of single zenith angle (SZA), Lang, Mill formula(Mf) and iterative formula(IF) are used to derive LAI DHP (effective LAI from them are called LAI SZA , LAI Lang LAI Mf and LAI IF ,respectively). LAI SZA , LAI Lang , LAI Mf and LAI IF are inter-compared and are also compared with from AccuPAR LP-80(LAI APAR ) and direct manual method (LAI direct ). It is found that, in general, LAI SZA , LAI Lang LAI Mf and LAI IF are similar to LAI APAR , but slightly lower than LAIdirect. During their intercomparsion, LAI Lang is more similar to LAIAPAR than other three and LAI Mf is more similar to LAI direct , while LAI SZA and LAI Lang are almost the same. It is implied that Lang method outperformed the other three when compared with AccuPAR LP-80 and Mill formula method outperformed the other three when compared with direct manual measure.

Introduction

Rice is one of the uppermost crops in China, and are increasingly important worldwide-both to China and to South Asian countries, and there is an obvious need to obtain accurate estimates of its leaf area index (LAI). LAI is an essential input into many models of rice growth and yield estimation as well as being an essential component of comparative studies of many canopy-level attributes such as carbon cycle, transpiration and water use efficiency and understory synthetic photon flux density (PPFD) capture [START_REF] Macfarlane | Adams: Estimation of leaf area index in eucalypt forest using digital photography[END_REF].So LAI is a dimensionless quantity characterizing the canopy of an ecosystem and a key component of biogeochemical cycles in ecosystems. Also managers (farmers and foresters), ecologists, site and global modelers, request information about canopy leaf area index. Unfortunately, this interface between ecosystem and atmosphere is very difficult to quantify, due to its spatial (horizontal and vertical) and temporal variability: annual cycles and interannual variability interact with the crop structure, stratification and heterogeneity [START_REF] Wu | Advance in ground based LAI measurement methods[END_REF].

Leaf area index can be measured directly by destructive harvest or allometric approaches [START_REF] Norman | Canopy structure, in Plant physiological ecology: field methods and instrumentation[END_REF][START_REF]Daughtry: Direct measurements of canopy structure[END_REF].Direct measurement may be practical with short canopies such as those of many arable crops, but is usually laborious and timely consuming. Indirect methods for determining LAI are now commonly used to overcome this problem [START_REF] Welles | Some indirect methods of estimating canopy structure, in Instrumentation for Studying Vegetation Canopies from Remote Sensing in Optical and Infrared Regions[END_REF][START_REF] Welles | Canopy structure measurement by gap fraction analysis using commercial instrumentation[END_REF] and have an additional advantage of being non-destructive. These methods include direct measurement of intercepted radiation using line quantum sensors [START_REF] Pierce | Running: Rapid estimation of coniferous forest leaf area index using a portable integrating radiometer[END_REF] or radiometers, inclined point quadrant [START_REF] Wilson | Inclined point quadrats[END_REF], gap fraction techniques [START_REF] Norman | Canopy structure, in Plant physiological ecology: field methods and instrumentation[END_REF][START_REF] Welles | Some indirect methods of estimating canopy structure, in Instrumentation for Studying Vegetation Canopies from Remote Sensing in Optical and Infrared Regions[END_REF] and capacitance sensors [START_REF] Vickery | An improved electronic capacitance meter for estimating herbage mass[END_REF]. Of these techniques, it is the group based on gap fraction data that is now the most widely applied.

Digital hemispherical photographs (DHPs) have been widely used to measure canopy structure with the recent rapidly advancement in digital cameras and may have brought us a needed tool at a reasonable cost to estimate LAI. Digital hemispherical photography system(DHPS) made with a fisheye lens allow the acquisition of DHPs without the need of scanners to digitize images and can be quickly inspected on the camera's viewer, or on laptop screen and in a timely fashion in the field. Moreover, for a given camera, the hemispherical images are of consistent size and position on the digital array. This retrieval can be done in a consistent manner at many zenith angles [START_REF] Chen | Evaluation of hemispherical photography for determining plant area index and geometry of a forest stand[END_REF][START_REF] Frazer | Canopy openness and leaf area in chronosequences of coastal temperate rainforests[END_REF][START_REF] Fournier | Modeling light obstruction in three conifer forests using hemispherical photography and fine tree architecture[END_REF] .

The model commonly used with indirect methods (including DHPS) to determine the LAI is the Poisson law. It assumes that leaves are uniformly and randomly distributed, which may be valid for homogeneous canopies [START_REF] Levy | Direct and indirect measurements of LAI in millet and fallow vegetation in HAPEX-Sahel[END_REF],but does not hold for canopies with aggregative patterns [START_REF] Lemeur | A critical review of light models for estimating the shortwave radiation regime of plant canopies[END_REF][START_REF] Nilson | Inversion of the frequency of gaps in plant stands[END_REF].To allow the use of the Poisson law, the concept of effective LAI is proposed [START_REF] Chen | Black: Measuring leaf area index of plant canopies with branch architecture[END_REF][START_REF] Chen | Cihlar: Plant canopy gap-size analysis theory for improving optical measurements of leaf area index[END_REF][START_REF] Chen | Optically-based methods for measuring seasonal variation of leaf area index in boreal conifer stands[END_REF] as a result of the contribution of woody elements to the total plant cover, which results in overestimation of LAI, and clumping of foliage, which results in underestimation of LAI. So, in this paper, a set of self-made DHPS is applied to acquire DHPs of rice canopy, and some self-developing DHP processing procedures are utilized to process DHP and to compute gap fraction rapidly from them. Then four methods, single zenith angle (SZA), Lang, Mill formula(Mf) and iterative formula(IF),are utilized to derive effective LAI based on those gap fractions(effective LAI from those four are called LAI SZA , LAI Lang , LAI Mf and LAI IF , respectively). At last methodology comparison will be investigated in detail.

Materials and methods

DHPS

In this paper, a DHPS based a set of self-made digital photography sensor with a fisheye lens is utilized to capture fisheye photographs from rice canopy, which contents of fisheye lens, filters changer, CMOS camera and laptop computer, as shown in Figure1. In Figure 1, the optical component of DHPS is the fisheye lens with a view zenith angle of 180° mounted on the CMOS camera, which can acquire visible band DHP of rice canopy. These DHPs are stored in format of JPG on the laptop computer hard disk.

Site description and experiment design

The study area is located in China National Rice Research Institute (CNRRI), Fuyang city of Zhejiang province. DHPs are acquired began in Sept., 2009, 30 days after early rice seedling being transplanted, and once per week in cloudy day till rice tassels. Four sampling plots are set up in the selected paddy field, and AccuPAR LP-80, DHPS and directly manual measure are carried out sequentially in each sampling plot. X parameter of AccuPAR LP-80 is initialized to 1.0 and detector is set up along two directions: obeying to rice seedling row and 45° clockwise. Mean value of two LAI-readings from two directions respectively is taken as the effect LAI by AccuPAR LP-80,which is called LAI APAR .

The DHPS is set levelly in the center of each sampling plot, and inside rice canopy with its fisheye lens oriented upward to 0.5m beneath rice canopy top, or above rice canopy with fisheye Len downwards to 0.5m away from top of rice canopy. Gap fraction is computed from DHPs as the input of calculating LAI. During direct manual measure, five rice seedlings or sixty pieces of rice leaves are picked out. Length (L)and width(W) of each leaf are measured by ruler, and leaf area(LA) is calculated by LA=L×W×0.83 [8] 。Effect LAI from direct manual measure (LAI direct ) is equal to all rice LAI in a sampling plot divided by its area.

Methodologies

Working flow

The working flow for extracting effect LAI from DHPs of rice canopy is shown as Figure 2.

Fig. 2. The flow sheet

Where SZA is single zenith angle, Mf is Mill formula, IF is Iterative formula.

Digital image pretreatment

In DHPS method, effect LAI is derived from gap fraction which computed from binary image of rice canopy. So a series of image pretreatment steps must be carried out to transfer original color photo to gray and to binary consisting of two types of pixels: leaves and background (non-leaves).

In color DHP from fisheye lens upward, rice leaves and background (sky) can be distinguished more clearly in the blue band than in other two. As we know, in visible band, rice leaves absorb more blue than red and green ray, transmission and reflection of blue is the least in rice canopy. So the blue band of DHP is selected as original data for binary classification.

On the other hand, in color DHP from fisheye lens downwards, rice leaves and background can be learn more clearly in the green band than in other two. In this situation, fisheye lens of DHPS is setup downwards above rice leaves, the green gray is reflected mostly by rice leaves, which is acquired by fisheye lens, and absorbed by non-leaves. So the green band of DHP is selected as original data for binary classification.

After color DHP being grayed, a typical two-peak image can be obtained: pixels of foliage making for one and the background make for another in the histogram, and there is little pixel which gray value is between the two peak values, which leads to the valley between two peaks. So the valley value can be applied as the threshold to classify foliage and background pixels clearly and reasonably. In this paper, in gray image from fisheye lens upward, pixels with gray value being lower than the threshold are classified as foliage and higher as background(sky or non-foliage),and in fisheye lens downward, the higher as foliage and lower as background(sky or non-foliage) [START_REF] Zhu | Trans): Digital Image Processing[END_REF].

Detection of gap fraction

Gap fraction was estimated using an overlay defining n annuli covering the hemispherical image with their positions determined by a range of zenith angles, such as n is 9 and zenith angles ranging from 5° to 85° with a step of 10. The location of the midpoint of each annulus can be fixed so that annuli are positioned with equal zenith angles. Gap fractions are determined for each annulus using image analysis and computed by Eq. ( 1):

0 01 () () ( ) ( ) P P PP      (1)
Where P(θ) is gap fraction at zenith angle θ,P 1 (θ) the fraction of foliage and P 0 (θ) background (none-foliage).

Methodology

Theory Basis.

In DHP method effect LAI (LAI eff ) is computed from the gap fraction P(θ) following the Poisson law [START_REF] Monsi | Uber den Lichtfaktor in den Pflanzengesellschaften und seine Bedeutung für die Stoffproduktion[END_REF][START_REF] Welles | Instrument for indirect measurement of canopy architecture[END_REF]:

ln ( ).cos () eff P LAI G      (2)
Where LAI eff is effect LAI, P(θ) is gap fraction at zenith angle θ, G(θ) is the mean projection of a leaf area unit in a plane perpendicular to direction θ which is directly dependent of the leaf angle distribution.

There are two distinct characters about G function [START_REF] Weiss | Coppin: Review of methods for in situ leaf area index (LAI) determination Part II. Estimation of LAI, errors and sampling[END_REF][START_REF] Lang | Leaf area and average leaf angle from transmission of direct sunlight[END_REF]: ① at a view angle of 57.5°, G(θ) can be considered as almost independent of leaf inclination, e.g. θ≈57.5 , G(θ) ≈0.5,as shown in Fig. 3(a); ② for a uniform leaf azimuth distribution and a constant leaf normal angle, G(θ) can be approximated as a linear function of θ in the 25-65° range. The slope of the regression (∂G(θ)/∂θ) was then related to the average leaf inclination angle(ALIA) by polynomial fitting, as shown in Fig. 3 Where T(57.5°) is gap fraction computed within 55-60° zenith angles. For this particular direction, G(θ) is almost independent of leaf inclination simplifying the LAI retrieval process [START_REF] Wilson | Estimation of foliage denseness and foliage angle by inclined point quadrats[END_REF].

Lang model.

On the assumption that a uniform leaf azimuth distribution and a constant leaf normal angle, G(θ) can be approximated as a linear function of zenith angle θ in the 25-65• range. The slope of the regression (∂G(θv)/∂θv) is then related to the ALIA by polynomial fitting. Using an initial estimate of LAI based on gap fraction measurements at a 55° (close to 57.5°) zenith angle, the slope ∂G(θ)/∂θ can be estimated and then ALIA can be derived following Eq.( 4):

S S S SS               (4)
x=-3+(ā÷9.65) -0.6061 (5)

2 2 1/ 2 2 1/ 2 1 1 1 1 2 2 1/ 2 2 1/ 2 2 2 2 2 ( tan ) cos ( ) , ( 1
) , x 1 (sin ) / ( tan ) cos ( ) , ( 1
) , x 1 1 ln[(1 ) / (1 )] 2 x G x x x G x x x                                        > (6)
Where ā is ALIA,S is the slope ∂G(θ)/∂θ [START_REF] Welles | Instrument for indirect measurement of canopy architecture[END_REF][START_REF] Lang | Leaf area and average leaf angle from transmission of direct sunlight[END_REF]. So the route for Lang method is followed as: detecting gap fraction P(25°) -P(75°) from DHP with zenith angle ranging from 25° to 75°; ②retrieving effect LAI by SZA; ③computing G(25°) -G(75°) by introducing gap fractions obtained in ① and LAIs in ② into Eq.( 2);④the ALIA ā is calculate by Eq.( 4) and then filled in Eq.( 5) in order to obtain x; ⑤x then is utilized in Eq.( 6) and assuming θ is 57.5°,a new G(57.5°) is worked out and then filled in Eq.(3) to calculate a new LAI by SZA.

Therefore, the LAI estimates can be refined, and the process is iterated several times until convergence.

Miller formula.

Chen and Black [START_REF] Chen | Black: Measuring leaf area index of plant canopies with branch architecture[END_REF] derived LAI from the gap fraction measured in all directions using the formula of Miller [START_REF]Miller: A formula for average foliage density[END_REF], which assumes that gap fraction depends only on the view zenith angle θ:

   2 / 0 sin cos )) ( ln( 2      d T LAI (7)
A practical method was proposed to compute the integral of Eq. ( 7) from gap fraction measurements in several directions [START_REF] Welles | Instrument for indirect measurement of canopy architecture[END_REF]:

        i i n i i T LAI sin cos )) ( ln( 2 (8)
Where n is denotation of using an overlay defining n annuli covering the hemispherical image with their positions determined by a range of zenith angles, T(θ i ) is gap fraction in the ith zenith angle ring, Δθ is the zenith resolution of annulus expressed by radian.

One of the main limitations with this technique is the necessity to sample the entire directional range of gap fraction variation, which might prove difficult for larger zenith angles.

Iterative formula.

In this method, a gap fraction model contents of LAI and average leaf inclination angle(ALIA) simultaneity is built by inputting Eq. ( 5) to Eq.( 6) then to Eq.( 2). Then the gap fraction model is inversed by using an iterative optimization technique [START_REF] Perry | Indirect sensing of plant canopy structure with simple radiation measurements[END_REF]. The Poisson model is used and the leaf angle inclination is assumed to be azimuthally isotropic with an ellipsoidal zenith angle distribution. Starting with a series of initially given LAI and ALIA, the gap fraction model is run in the forward direction to simulate the gap fraction. The variables LAI and ALIA are then iteratively changed, using the simplex algorithm [START_REF] Nelder | Mead: A simplex method for function minimization[END_REF], until a good agreement is met between the simulated and measured gap fraction values.

Results and discussion

Variety of gap fraction with zenith angle

The trend of gap fraction varies with zenith angle θ is illustrated in Figure 4,where there are five curves created by gap fractions from fisheye lens upward(labeled as "LenUp") and one for downward (labeled as "LenDown") ,n is the number of concentric annulus. As shown in Fig. 4,gap fractions decrease while zenith angle θ increasing, which might be caused by that rice leaves captured by fisheye lens of DHPS increase while θ increasing.

Impact of fisheye lens orientation on gap fraction.

① Gap fractions from fisheye lens downward is averagely higher than upward. It might be that the reflection of leaves underlaying rice canopy cannot reach to fisheye lens of DHPS because of being sheltered from the topper when lens of DHPS is downward, which leads to the shadowed leaves being misclassified as background in DHP to increase gap fraction. ②Gap fractions from fisheye lens downward is lower than upward when zenith angle is up to the zenith(zenith angle is close to 0°).It could be the sky inverted reflections in water, captured when fish eye lens downward, are misclassified as "leaves" because of almost the same gray value as rice leaves, which result decrease of gap fractions.

③Gap fraction from fisheye lens upward decrease more quickly than upward while zenith angle increasing. It is may be because of the shape of rice leaf: rice leaf is gladiate and its tip is captured in DHP in a greater probability when fisheye lens downward, which leads to gap fraction decreasing. And it is on the contrary when fisheye lens of DHPS is oriented upward.

Impact of zenith angle resolution on gap fraction.

There involve all gap fraction information from full view of zenith angle from 0° to 90° in DHP. In this paper DHP is divided into n annuli, and n gap fractions are detected respectively. When n grows view zenith resolution reaches more higher which leads to more detail gap fraction being detected and curves undulates more sharply ,especially in zenith angle of 0<θ<10°,as shown on Fig. 4. As we know, there are few pixels contented by a single zenith annulus in a relatively small zenith angle such as 0<θ<10°,so a little changing of background pixels number might bring a biggish changing of gap fraction .

In this paper, gap fraction is estimated using an overlay defining n (n = 9,18,27,36,45,54,63,72,81 and 90) annuli covering the hemispherical image, and 10 zenith resolutions are reached. There are only 5 gap fraction curves of fisheye lens upward where n = 9,18,27,36, 72, 90 and 1 downward of n=18 are picked out as representative zenith resolutions. For these curves of fisheye lens upward, there is little overlapping among the curve of n=9 and the others which are closed to each other. There exists almost the same trend while fisheye lens is oriented upward.

LAI retrieving

Comparison of four methods.

In this paper LAI SZA is calculated based on gap fraction detected from the annulus of 55°60°. LAI SZA from four sampling plots are 2.52, 2.53, 2.49 and 2.50 while fisheye lens is oriented upward and 2.33, 2.35, 2.30 and 2.33 while downward. They are close to LAI Lang which is 2.5 when lens downward and 2.31 when upward. But actually G(θ) is derived by LAI SZA in Lang way, which means that Lang way is built upon SZA.LAI from other three approaches are illustrated by Fig. 5,in which, LAI Mf is the minimum and LAI Lang is as much as LAI IF .

Impact of fisheye lens orientation on effect LAI.

LAIs from fisheye lens upward are always higher than downward. As mentioned in 3.1.1 that shadow of rice leaf is misclassified as background to increase gap fraction in DHP from fisheye lens downward. And a higher gap fraction may leads to a relatively lower LAI according to Eq.2.

Impact of zenith angle resolution on effect LAI.

Miller formula is integral, so increasing zenith resolution may do a great deal of good to improve its retrieving precision of LAI. As demonstrated by Fig. 5, LAIsMf from two Len orientations keep on increasing while n increases until convergence. LAI Lang vary little with variety of zenith resolution except for n=9, which indicates that Lang method is a relatively stable one. LAI IF will fluctuate around a value when zenith resolution reaches to a certain level, which means that Iterative formula method is a sensitive one.

So, in the following text, zenith resolution of n=90 is applied in Miller formula and n=18 in the others.

Comparisons among DHP and AccuPAR and direct manual measure

LAI DHP from two lens orientations and four methods varying from 2.14 to 2.50 are lower than 2.75 of LAI direct and close to LAI APAR from AccuPAR LP-80,as shown in Table 1. LAIs from DHPS and AccuPAR PL-80 are both lower than from direct measure. It may be that DHP and AccuPAR PL-80 are both based on the Poisson law which assumes that leaves are uniformly and randomly distributed, which may be valid for homogeneous canopies [START_REF] Levy | Direct and indirect measurements of LAI in millet and fallow vegetation in HAPEX-Sahel[END_REF] but does not hold for crops with aggregative patterns such as rice [START_REF] Lemeur | A critical review of light models for estimating the shortwave radiation regime of plant canopies[END_REF][START_REF] Nilson | Inversion of the frequency of gaps in plant stands[END_REF].The heterogeneousness and aggregation of rice canopy must be taken into account if measuring a LAI with more higher precision than now.

Conclusions

In this paper it is investigated that effect LAIs of rice canopy are calculated based on gap fractions detected from DHPs taken by a set of self-made DHPS with fisheye lens and comparisons are carried out among LAIs from DHP and AccuPAR PL-80 and direct manual measure. It is found that LAIs from DHPs are lower than from direct manual measure and close to from AccuPAR.

In this paper the precision of LAI from DHP is not as well as from direct manual measure, but DHP is a quick, nondestructive and simple way to obtain LAI in situ. In addition, more parameters of rice canopy structure will be detected simultaneously from DHP if more models are applied in DHP pretreatment. So it is believed that DHP will be a widely utilized way to measure rice canopy structure parameters

In the next study for detecting rice canopy structure parameters it is argued that: gap fractions will vary with rice leaves shape during calculating because of relatively lower rice canopy;② The heterogeneousness and aggregation of rice foliage will impact on the precision of extracting LAI.

Fig. 1 .

 1 Fig. 1. Digital hemispherical photography system (DHPS) for rice canopy

Fig. 3 . 2 . 4 . 2

 3242 Fig.3. The projection function G[18]

Fig. 4 .

 4 Fig. 4. Gap fractions vary with the view angles

Fig. 5 .

 5 Fig. 5. LAIs from different lens direction, zenith angle resolution and method

Table 1 .

 1 Comparisons of LAIs from DHP and AccuPAR and direct method

				L
		Methods		AI
		Single	zenith	2
		angle(SZA)		.49
	DHP Fisheye Len Upward	Miller formula Lang		2 .32 2 .50
		Iterative formula	2 .42
		Single	zenith	2
		angle(SZA)		.31
	DHP Fisheye Len Downward	Miller formula Lang		2 .14 2 .31
		Iterative formula	2 .37
	Direct manual measure		2 .75
	AccuPAR LP-80			2 .32

1 Project supported by the National High Technology Research and Development Program (863) ( 2008AA10Z206).


Liu Zhao
email: zhaoliu@nercita.org.cn

Zhihong Ma

Yunxia Luan

Anxiang Lu

Jihua Wang

Ligang Pan

{liu Zhao

Molecular Methods of Studying Microbial Diversity in Soil Environments

Keywords: Microbial diversity, SIP, T-RFLP, RAPD, DGGE, Soil

Microorganisms live in all parts of the biosphere and are critical to nutrient recycling in ecosystems. In recent years, the development of methodologies for the analysis of microorganisms and microbial ecology, at the molecular level, has progressed phenomenally. This review introduces and compares the various molecular methods for studying microbial diversity in soil environments, and the advantages and disadvantages of current methods are proposed as well.

Introduction

Microorganisms are ubiquitous in the environment and fulfill a range of important ecological functions, particularly those associated with nutrients cycling processes and maintenance of ecosystem health in soil [START_REF] Nannipieri | Role of stabilised enzymes in microbial ecology and enzyme extraction from soil with potential applications in soil proteomics[END_REF] . Soil contain an estimated 109 prokaryotes and more than 2000 genome types per gram of soil, with an average genome type representing less than 0.05% of the soil community [START_REF] Torsvik | Microbial diversity and function in soil: from genes to ecosystems[END_REF][START_REF] Gans | Computational improvements reveal great bacterial diversity and high metal toxicity in soil[END_REF] .

Until a few decades ago, soil microorganisms could only be studied by microscopic observation or culture-dependent methods. In less than a decade, using nucleic acid probe technique for the detection of microorganisms had exploded. Nowadays, molecular ecology techniques based on sequence comparisons of nucleic acids can be used to provide molecular characterization while providing a classification scheme that predicts natural evolutionary relationships [START_REF] Handelsman | Conversations with the silent majority[END_REF][START_REF] Schutte | Advances in the use of terminal restriction fragment length polymorphism (T-RFLP) analysis of 16S rRNA genes to characterize microbial communities[END_REF] .

These laboratory-based works have been spectacularly successful in revealing details of soil microbial diversity. So we attempt to draw together some of these studies, with emphasis on the advantages and disadvantages of current molecular methods, for studying microbial diversity in soil environments.

Stable Isotope Probing

Stable isotope probing (SIP) has been coupled with nucleic acid methods to provide a culture-independent method of linking the identity of bacteria with their function in the environment [START_REF] Mcdonald | Stable isotope probing of nucleic acids in methanotrophs and methylotrophs: a review[END_REF][START_REF] Dumont | Stable isotope probing-linking microbial identity to function[END_REF] . Soil is either incubated after adding a 13 C-labeled substrate or a plant is labeled with [START_REF] Dunfield | Seasonal changes in the rhizosphere microbial communities associated with field-grown genetically modified canola (Brassicanapus)[END_REF] C-CO 2 . Soil DNA or RNA is then extracted and centrifuged in a density gradient to separate [START_REF] Dunfield | Seasonal changes in the rhizosphere microbial communities associated with field-grown genetically modified canola (Brassicanapus)[END_REF] C-labeled nucleic acids from those containing. Then separated, labeled DNA can be amplified using PCR. Analysis of the PCR products, through cloning and sequencing for example, allows the microbes that have assimilated the labeled substrate to be identified [START_REF] Griffiths | CO 2 )-13 C pulse labelling of plants in tandem with stable isotope probing: methodological considerations for examining microbial function in the rhizosphere[END_REF][START_REF] Leake | Carbon fluxes from plants through soil organisms determined by field (CO 2 )-13 C pulse-labelling in an upland grassland[END_REF] . SIP-based approaches hold great potential for linking microbial identity with function, but at present a high degree of labeling is necessary to be able to separate labeled from unlabeled marker molecules.

Terminal-Restriction Fragment Length Polymorphism

Terminal-restriction fragment length polymorphism (T-RFLP) has multiple advantages for its rapid gain in popularity: data are quantitative and comparable between laboratories, the final electrophoresis step can be performed with automated sequencing equipment at core sequencing facilities [START_REF] Bent | The tragedy of the uncommon: understanding limitations in the analysis of microbial diversity[END_REF] . This method provides a picture of the community including incorporates diversity and phylogenetic details. The profiles can be generated by using the procedure from physical-capture [START_REF] Blackwood | Evaluating the physical-capture method of terminal restriction fragment length polymorphism[END_REF] , fluorescence scanner [START_REF] Ikeda | Microbial community analyses using a simple, rapid detection method for DNA fingerprints with a fluorescence scanner[END_REF] or 32 P-labelled primers [START_REF] Dunfield | Seasonal changes in the rhizosphere microbial communities associated with field-grown genetically modified canola (Brassicanapus)[END_REF] . Because of its relative simplicity, T-RFLP has been applied to the analysis of soil microbial diversity, for instance, fungal ribosomal genes [START_REF] Genney | Fine-scale distribution of pine ectomycorrhizas and their extramatrical mycelium[END_REF][START_REF] Kennedy | Soil bacterial and fungal community structure across a range of unimproved and semi-improved upland grasslands[END_REF] , bacterial 16S rRNA genes [START_REF] Hullar | Recurring seasonal dynamics of microbial communities in stream habits[END_REF][START_REF] Thies | Soil microbial community analysis using terminal restriction fragment length polymorphisms[END_REF] , and archaeal 16S rRNA genes [START_REF] Leybo | Effect of hydrogen concentration on the community structure of hydrogenotrophic methanogens studied by T-RFLP analysis of 16S rRNA gene amplicons[END_REF] . In addition, T-RFLP has been used for the analysis of functional genes such as those encoding for nitrogen fixation and methane oxidation [START_REF] Mohanty | Differential effects of nitrogenous fertilizers onmethane-consuming microbes in rice field and forest soils[END_REF][START_REF] Zak | A molecular dawn for biogeochemistry[END_REF] .

Fluorescence in Situ Hybridization

Fluorescence in situ hybridization (FISH) is a relatively new technique utilizing fluorescently labeled DNA probe to detect genes of microorganisms in soil samples. Apart from allowing direct visualization of bacteria in the environment, FISH also has the added advantage of being able to detect active cells by targeting rRNA. This method has been reliably used for identification and quantification of ammoniaoxidizing bacteria [START_REF] Briones | Influence of different cultivars on populations of ammonia-oxidizing bacteria in the root environment of rice[END_REF][START_REF] Meyer | Challenges for simultaneous nitrification, denitrification, and phosphorus removal in microbial aggregates: mass transfer limitation and nitrousoxide production[END_REF] . Several studies have used FISH coupled with microautoradiography (FISH-MAR), this combined approach allows in situ identification and provides information on substrate utilization in complex microbial communities [START_REF] Ouverney | Combined microautoradiography-16S rRNA probe technique for determination of radioisotope uptake by specific microbial cell types in situ[END_REF][START_REF] Zwirglmaier | Fluorescence in situ hybridisation (FISH)-the next generation[END_REF] .

Microarray

Microarray for microbial community analysis has been classified into three main categories depending on the combination of probe types and target molecules exploited: community genome arrays, rRNA-based oligonucleotide microarrays and functional gene arrays [START_REF] Zhou | Challenges in applying microarrays to environmental studies[END_REF][START_REF] Zhou | Microarrays for bacterial detection and microbial community analysis[END_REF] . Tiquia constructed a 50-mer oligonucleotide microarray using 763 genes involved in nitrogen cycling [START_REF] Tiquia | Evaluation of 50-mer oligonucleotide arrays for detecting microbial populations in environmental samples[END_REF] . The increased use of cultivationindependent metagenomic approaches employing large-insert cloning could lead to an important extension of the CGA approach for large genomic fragments from uncultivated microorganisms [START_REF] Handelsman | Metagenomics: application of genomics to uncultured microorganisms[END_REF][START_REF] Riesenfeld | Metagenomics: genomic analysis of microbial communities[END_REF] . In the meantime, adapting enzymatic signal amplification approaches for microarray analysis [START_REF] Small | Direct detection of 16S rRNA in soil extracts by using oligonucleotide microarrays[END_REF][START_REF] Denef | Validation of a more sensitive method for using spotted oligonucleotide DNA microarrays for functional genomics studies on bacterial communities[END_REF] , the use of high sensitivity microarray hybridization detection devices hold much promise for enhancing the sensitivity of direct detection of extracted rRNA [START_REF] Mckendry | Multiple label-free biodetection and quantitative DNA-binding assays on a nanomechanical cantilever array[END_REF][START_REF] Epstein | High-density, microsphere-based fiber optic DNA microarrays[END_REF][START_REF] Thompson | Label-free detection of nucleic acid and protein microarrays by scanning Kelvin nanoprobe[END_REF] .

Random Amplified Polymorphic DNA

Random amplified polymorphic DNA (RAPD) use DNA products by PCR which are based on random priming of the target DNA, these primers are usually 10 bases in length and are designed to a number of random target sites [START_REF] Podila | Functional genomic approaches for studies of mycorrhizal symbiosis[END_REF] . Huai studied the genetic variation and spatial distribution of the ectomycorrhizal fungus, the 33 sporophores studied belonged to distinct genotypes [START_REF] Huai | Genetic diversity of an ectomycorrhizal fungus Tricholoma terreum in a Larix principis-rupprechtii stand assessed using random amplified polymorphic DNA[END_REF] . The advantage of this approach is that it requires very little sample material and obtains results rapidly. However, this method is less susceptible to base changes in the target DNA. Thus usually the PCR is performed at low stringency for the first few cycles. This technique allows the generation of product when mismatches between template and primers occur. Hence, similar patterns generated using this technique better reflect phylogeny of the phages. One disadvantage of this technology is the results may be difficult to repeat by other users [START_REF] Micheli | Reproducible DNA fingerprinting with the random amplified polymorphic DNA (RAPD) method[END_REF][START_REF] Henry | Quantification of denitrifying bacteria in soils by nirK gene targeted real-time PCR[END_REF] .

DNA Fingerprinting

DNA fingerprinting is used to distinguish differences in the genetic makeup of microbial populations from different samples. The advantage of this technique is that it enable high sample throughput and can be used to target sequences that are phylogenetically or functionally significant. Nowadays, the most common used techniques are denaturing gradient gel electrophoresis (DGGE) or temperature gradient gel electrophoresis (TGGE). They have been used to profile fungal microbial communities from many diverse environments [START_REF] Möhlenhof | Molecular approach to the characterisation of fungal communities: methods for DNA extraction, PCR amplification and DGGE analysis of painted art objects[END_REF][START_REF] Nikolcheva | Determining diversity of fungi on decaying leaves: comparison of traditional and molecular approaches[END_REF] . DGGE fingerprints of total DNA from rhizosphere revealed that there was a relationship between fungal community composition and rhizosphere development [START_REF] Gomes | Dynamics of fungal communities in bulk and maize rhizosphere soil in the tropics[END_REF] . Fungal community diversity was studied by PCR-DGGE followed by sequence analyses of ITS fragments in soil samples [START_REF] Anderson | Potential bias of fungal 18S rDNA and internal transcribed spacer polymerase chain reaction primers for estimating fungal biodiversity in soil[END_REF] . Although DGGE is a promising technique, it can still underestimate fungal diversity. The number of bands depends on the resolution of the gels, this takes time to optimize and is difficult to reproduce [START_REF] Fromin | Statistical analysis of denaturing gel electrophoresis (DGGE) fingerprinting patterns[END_REF] . This has already been demonstrated in previous studies that phylogenetically distant taxa can have comigrating bands and that one band does not mean one unique phylotype [START_REF] Sekiguchi | A single band does not always represent single bacterial strains in denaturing gradient gel electrophoresis[END_REF][START_REF] Kisand | Limited resolution of 16S rDNA DGGE caused by melting properties and closely related DNA sequences[END_REF] .

A new DNA-based fingerprinting approach, two-dimensional polyacrylamide gel electrophoresis (2D-PAGE), can be used to separate PCR amplicon of the ITS regions [START_REF] Jones | Soil microbial community analysis using two-dimensional polyacrylamide gel electrophoresis of the bacterial ribosomal internal transcribed spacer regions[END_REF] . But this technique does not lend itself to high sample throughput. However, its improved ability to discriminate between soil communities and retrieve sequence information make it a powerful technique for elucidating key differences in community structure.

Conclusions and Future Directions

Soil has been dubbed "The Final Frontier". Current knowledge pertaining to the diversity and distribution of soil ecosystem is still rudimentary. In recent years, obviously improvement in traditional approaches combined with various molecular techniques, such as molecular and phylogenetic inventorying via clone libraries [START_REF] Loy | Microarray and functional gene analyses of sulfate-reducing prokaryotes in low sulfate, acidic fens reveal co-occurrence of recognized genera and novel lineages[END_REF] , quantitative real-time PCR [START_REF] Kolb | Quantitative detection of methanotrophs in soil by novel pmoA-targeted real-time PCR assays[END_REF] , pyrosequencing [START_REF] Huber | Microbial population structures in the deep marine biosphere[END_REF] , quantitative membrane hybridization [START_REF] Hristova | Desulfotomaculum genus-and subgenus-specific 16S rRNA hybridization probes for environmental studies[END_REF] has provided new data on these aspects. Nonetheless, we feel that these sophisticated methods are highly relevant to the existing knowledge of the role of microorganism in ecosystem processes, but the application of the techniques to this end is far from complete.

From recent molecular ecological studies, we have seen that most carbon-and nitrogen-cycling gene sequences are divergent from those of the model organisms on which most of our existing appreciation of biogeochemical cycles are based. However, our understanding of ecologically relevant microorganism involved in these cycles is limited. Hence we are aware that further understanding of how they fit into the complexities of ecosystems will require both bottom-up and top-down approaches.

Acknowledgments

. This work was supported by the public foundation from Beijing Municipal Rural Affair Committee (20080901), National High Technology Research and Development Program 863 (2010AA10Z403, 2007AA10Z202) and Beijing Municipal Science and Technology Commission Program (Z09090501040901).


Renwang He

Dandan Xie

Yuling Zhao

Yibo Yang

Simulation Study of a Novel Algorithm for Digital Relaying Based on FPGA

Keywords: radix-4 FFT, Qua-input/output, digital protection, MATLAB

The paper proposes a new fast Fourier transform (FFT) protection algorithm for digital relaying based on field-programmable gate array (FPGA), constructs detailed models and performs corresponding simulations. Simulation results show that the harmonic components of voltages and currents can be conveniently calculated, and the digital protection functions can be effectively fulfilled with enhanced accuracy and reduced storage requirements.

Introduction

The power system fault signals contain a large number of harmonics and a DC component [START_REF] Babu | Accelerated Fourier Algorithm based Digital Multifunction Relay for Transmission Line Protection[END_REF], it needs to extract necessary information to estimate whether the digital protection should operate. Most of the digital protections need to extract the fundamental signal, and filter the other harmonics to avoid malfunction [START_REF] Zhang | Application and research of new algorithm for microcomputer-based transformer relaying protection[END_REF]. Fast Fourier Transform (FFT) algorithm can extract every harmonic without attenuation, guaranteeing the other integer harmonics and the DC component attenuating to zero, but can not eliminate the non-integer harmonics and the capacity for suppressing nonperiodic components in low-frequency is poor. The radix-4 decomposition, which divides the input sequence recursively to form four-point sequences, has the advantage that it requires only trivial multiplications in the four-point DFT. This results in the highest throughput decomposition, while requiring non-trivial complex multiplications in the post-butterfly twiddle-factor rotations only. Complex multiplication decrease and the four inputs and outputs parallel processing improves processing speed greatly.

Fourier digital protection algorithm

When ignore the non-integer periodic components, the voltage and current signals in power system can be expressed in the form of Fourier series,

1 0 ( ) cos( ) . k k x t X k t k ω φ ∞ = = ∑ + (1)
where expresses the DC component; Expand expression (1) and accord the orthogonality of the trigonometric functions in [0, ], it can be derived:

1 T 1 1 0 1 2 cos ( ) cos( ) . T R k k k X X x t k T φ = = ∫ t d t ω (2) 1 1 0 1 2 sin ( )sin( ) . T I k k k X X xt k td T φ = - = -∫ t ω (3) Rk X
, are the real and imaginary parts of the harmonic component respectively. , , are analog while , , are digital.

Ik X th k k X Rk X Ik X ) (k X ) (k X R ) (k X I
Each cycle samples N points and the integral is approximate to the sum [START_REF] Liu | Application of Fourier algorithm based on FFT in microcomputer-based relay protection[END_REF].

1 0 2 2 ( ) ( ) cos . N R k R n X X k x n nk N N π - = ⎛ ≈ = ⎜ ⎝ ⎠ ∑ ⎞ ⎟ (4) 1 0 2 2 ( ) ( ) sin .
N Ik I n X X k x n nk N N π - = ⎛ ≈ =- ⎜ ⎝ ⎠ ∑ ⎞ ⎟ (5)
where is the sampled value of .

) (n x th n ) (t x

Radix-4 FFT algorithm

The sampled power system voltage or current signals are discrete periodic sequences. Their Discrete Fourier Transform (DFT) expressions are:

2 1 1 0 0 ( ) ( ) ( ) . N N j nk nk N N n n X k x n e x n W π - - - = = = = ∑ ∑ (6)
where 5) and ( 6), it can be obtained

N j N e W π 2 - = From (
1 0 2 2 ( ) ( ) ( ) ( ) . N nk Rk Ik R I N n X jX X k jX k x n W X k N N - = + = + = = ∑ (7)
Each harmonic needed in digital protection can be calculated by using DFT algorithm.

FFT algorithm is an improvement of the DFT algorithm. This paper mainly introduces the radix-4 FFT algorithm according to the DIF (decimation in frequency) [START_REF] Cheng | Digital Signal Processing Tutorial[END_REF].

The sampled points during every cycle are , and expression (7) can be written as:

L N 4 = 0 1 1 3 3 3 1 2 1 0 0 0 0 ( ) ( , , , ) . L nk L L N n n n X k x n n n n W - - - = = = = ∑∑ ∑ (8)
where

∑ - = = 1 0 4 L i i i n n , , ∑ - = = 1 0 4 L i i i k k 3 ,
, , ,,) n

L L L L L L L L W W W W X n n x n n X n n x n n W W W W X n n x n n W W W W X n n x n W W W W - - - - - - - - ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦
⎡ ⎤ ⎡ ⎡ ⎤ ⎤ ⎢ ⎥ ⎢ ⎢ ⎥ -- ⎥ ⎢ ⎥ ⎢ ⎢ ⎥ = ⎥ ⎢ ⎥ ⎢ ⎢ ⎥ - - ⎥ ⎢ ⎥ ⎢ ⎢ ⎥ -- ⎣ ⎦ ⎥ ⎣ ⎦ ⎣ ⎦

ω

  is the amplitude and phase of the harmonic component respectively.

  th k

Acknowledgments

The authors would like to thank the Research Fund of East China Jiaotong University which supports the project (09DQ08).

The multiplication of -j and j just needs only to exchange the real part with the imaginary and adds the necessary plus or minus sign. Expression (9) shows that the computer storage requirements are greatly reduced and the processing speed is enormously enhanced.

MATLAB simulations

Sample 64 points per cycle. The 64 points will be divided into four groups, paralleled qua-input and qua-output natural order in the same address, the most critical factor in the calculation is to determine twiddle factors and data reordering. The data flow figure of the basic processor, A0, A1, A2, A3 is the four sets of data. T0, T1, T2 are the corresponding twiddle factors in multiplications each four points group. The processing flow is shown in Figure 1. MATLAB simulation results verify the correctness and effectiveness of the proposed algorithm.

Conclusions

In this paper, the 64-point data in four groups have been reordered twice, and both of the input and the output of the qua-input/output radix-4 FFT algorithm are in natural order. All the necessary power system harmonics are employed as information and criterion of the digital protection. MATLAB simulation results show that the proposed algorithm can reduce the computer storage requirements and increase the processing speed of the digital protection system.