Feature extraction for gas photoacoustic spectroscopy and content inverse based on overcomplete ICA bases

Optics & Laser Technology 48 (2013) 580–588

Contents lists available at SciVerse ScienceDirect

Optics & Laser Technology journal homepage: www.elsevier.com/locate/optlastec

Feature extraction for gas photoacoustic spectroscopy and content inverse based on overcomplete ICA bases Zhiying Wu a,b, Wang Zhang a, Jianwei Wang a, Qingxu Yu a,n a b

School of Physics and Optoelectronic Technology, Dalian University of Technology, No.2, Ling Gong Road, Gan Jing Zi District, Dalian 116023, PR China Information Engineering Department, Environmental Management College of China, Qinhuangdao 066004, PR China

a r t i c l e i n f o

a b s t r a c t

Article history: Received 14 November 2011 Received in revised form 26 June 2012 Accepted 27 July 2012 Available online 12 January 2013

In near-infrared band, macro representations of rovibrational-transition absorptions of some low molecule mass gases are ‘‘spiky’’ transient lines. The pressure broadening at ambient temperature and pressure make lineshape profiles at adjacent wavelength produce mutual projection or superposition, and it makes the ‘‘spiky’’ absorption features just be hidden in observable spectra. A blind source separation (BSS) model can extract features which hide within sampled spectra data, if a mixed gas phase photoacoustic signal is regarded as the weighted sum of feature absorptions by the spacing between adjacent feature wavelengths. A BSS model based on overcomplete ICA basis is proposed, and a weight truncation equation for transforming ICA bases based on the five-point sampling method is created so that a FastICA algorithm based on a negentropy approximation can be used to extract feature components. The resulting validity test using data from a real-time experiment of gas detection showed that the method detection limit can be improved from 16 ppb to 10 ppb with improved accuracy and signal-to-noise ratio. Crown Copyright & 2012 Published by Elsevier Ltd. All rights reserved.

Keywords: Photoacoustic spectroscopy Feature extraction Overcomplete ICA bases

1. Introduction Photoacoustic spectroscopy (PAS) is a tool for trace gas detection with minimal sample preparation. The photoacoustic (PA) effect, first discovered by Bell [1], has been used in many areas [2–9]. As light entering a PA cell is chopped, the pressure alternately increases and decreases, generating an acoustic signal. This signal is first detected by a microphone sensor, and then is amplified by a lock-in amplifier. When the measured gas sample is free of CO2, due to the vibration–translation relaxation time of the most gas molecules is fixed in around the order of magnitude of 10  6 s (at 1 atm) [10], the PA signal has the same phase with the probe laser, and thus the single gas PA signal can be expressed as sðlÞ ¼ FSmp P l ðlÞca

ð1Þ

where s(l) is the PA signal detected by microphone (mV/W), F is the constant of PA cell, Smp is the microphone sensitivity (mV Pa  1), Pl(l) is the incident-laser power (mW) at wavelength l (nm), c is the concentration of the detected gas, and a is the absorption coefficient at the gas absorption peak. The value of a can be described by a(n)¼Sg/[p(n  n0)2 þ g2], where n is the wavenumber in cm  1. n0 is the wavenumber of an absorption line at n Corresponding author at: School of Physics and Optoelectronic Technology, Dalian University of Technology, No.2, Ling Gong Road, Gan Jing Zi District, Dalian 116023, PR China. Tel.: þ 86 41184708379. E-mail address: [email protected] (Q. Yu).

its absorption peak. S is the line strength (cm  1/mol cm  2) at 296 K and 1 atm. g is the half width at half-maximum (HWHM) in cm  1. For real gas detection, Eq. (1) must be further modified as sðlc Þ ¼

n X

wj sj ðlc Þ

ð2Þ

j¼1

where lc is the center wavelength of an absorption peak in nm, sj(lc) indicates the amplitude of a certain single PA signal of the jth feature absorption at wavelength lc, s(lc) is the observable PA signal, which is the superimposition of all PA signals sj(lc) (j¼1,2,y,n) at lc, wj is superimposition coefficient, and n is the number of mixed gas components (n is not less than two). A solution to Eq. (2) has had several versions based on linear transformation, such as Levenberg–Marquardt fitting based on the least square method [10] or as linear equations [11]. However, the observed PA signal s(lc) on the left-hand side of Eq. (2) is not simply a linear combination of mixed components sj(lc) (j ¼1,2,y,n) on the right-hand side, but the blind combination of them by their cross-weighting contributions. Therefore, Eq. (2) cannot be solved by linear transformation unless some additional configurations are given [3,10,11]. Pichler et al. [12] used blind source separation (BSS) to analyze solid PAS, they measured the depth in the absorption sample with FTIR data acquired from two-layer PAS samples. Their BSS model was described by xðtÞ ¼ AsðtÞ

0030-3992/$ - see front matter Crown Copyright & 2012 Published by Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.optlastec.2012.07.028

ð3Þ

Z. Wu et al. / Optics & Laser Technology 48 (2013) 580–588

where x(t) ¼[x1(t), x2(t), y xk(t), y, xm(t)]T is the vector of observed sample, m is the number of observed variables xk(t) in the vector x(t) (k¼1, 2, y, m). s(t) ¼[s1(t), s2(t), y sj(t), y, sn(t)]T is the original vector with latent feature variables sj(t) (j¼1, 2, y, n), which are also called independent components (ICs) or source signals, n is the number of ICs, A is the unknown mixing matrix, and [.]T denotes matrix transposition. When m is not less than n, the solution to Eq.(3) is called the estimation of complete independent component analysis (ICA) basis [12], x(t)¼[x1(t), x2(t), y xk(t), y, xm(t)]T is called complete ICA basis. The methods for estimating Eq. (3) with these complete ICA bases to find original unknown variables sj(t) (j¼1, 2, y, n) have been determined by many researchers [12–14]. For a real low mass molecule gas detection at ambient temperature and pressure (usually at approximately 296 K, 1 atom), the feature absorptions are very ‘‘sharp’’ rovibrational–transition lines, but due to pressure broadening, they are hidden within the observable signal s(lc) [15]. So, the observed PA signal s(lc) in Eq. (2) can be regarded as the weighted sum of projections of the broadened absorption lineshape profiles by the spacing between adjacent feature wavelengths, the resultant projection coefficient or weight is wj (j¼ 1, 2, y, n). Then, another solution to Eq. (2) can also expressed as the BSS model similar to Eq. (3), and a group of equations can given as  T ð4Þ xn ðtÞ ¼ sðlc Þ Asn ðtÞ ¼

n X

wj sj ðlc ÞT

ð5Þ

j¼1

xn ðtÞ ¼ Asn ðtÞ

ð6Þ

where A¼[w1, w2, y,wj, y, wn], which is the unknown matrix indicating superposition weight from every pressure broadening absorption lineshape; sn(t)¼[s1(lc), s2(lc), y,sj(lc), y, sn(lc)], which is the vector of original feature PA signals. The observable sample xn(t) is a mixture and an entire dataset with no given classification. That is, the number of the observable variables in the vector xn(t) is equal to 1 (m¼1) for a certain feature wavelength, and the number of sources in the vector

581

sn(t) is not less than two (nZ2). Therefore, xn(t) is an overcomplete ICA bases vector [16–19], and so the solution to Eq. (6) is an estimate of the BSS model based on overcomplete ICA basis. Consequently, the problem is first focused on a transformation from overcomplete ICA basis to complete ICA bases.

2. Estimation of the BSS model based on overcomplete ICA basis x*(t) 2.1. Transforming overcomplete ICA basis xn(t) into complete ICA bases xi  1(t) and xi(t) Herein, the analysis of two-component gas sample x*(t) is selected as an illustration. It completely can be promoted as multi-component gas sample with selecting small sample by adding window configuration based on the spacing of adjacent feature wavelength, see Appendix A. Fig. 1(a) is the selected synthesis sample of 6% H2O and 500 ppb NH3 feature absorptions (in 10  26 cm  1/mol cm  2) from HITRAN database, which are modulated by pressure broadening profiles, it is used for finding a transformation from overcomplete ICA basis xn(t) into the complete ICA bases x(t) by a numeric simulation. The sample xn(t) is empirically divided into two subsamples xi  1(t) and xi(t) by the spacing between adjacent feature wavelengths, the samples xi  1(t) and xi(t) have the same length N, i is an integer not less than 2. For convenience of explanation, xi  1(t) is designated as the feature component, also known as interesting sample, which is denoted by A, xi(t), which is to be removed, called as the projection sample, which is denoted by B, and A and B can be described as A

A ¼ ½xi1 ðtÞ ¼ ½sample1 , N

B

B ¼ ½xi ðtÞ ¼ ½sampleN , N

A

sample2 , B

sampleN1 ,

A

sampleNdm , B

sample1 

A

sampleNðdm 1Þ ,

A

sampleN 

ð7Þ

where the interesting sample A is corresponding to feature A component of 6% H2O, samplei denotes any an element in A, A i¼1, 2, y, N, sampleNdm ¼max([xi  1(t)]), which is the maximum of A and also is the maximum of the total sample xn(t), dm is the

Fig. 1. Selected sample from the HITRAN database at approximately temperature T¼296 K and pressure P¼ 1 atm. (a) Synthesis sample x*(t) in form of overcomplete bases and the processing chart for its transformation based on Dl-parameter. A complete-ICA-bases observation variable xi  1(t) is labeled as 2, and another variable xi(t) is labeled as 3. (b) The chart of the extracted feature components with Lorentz-pressure broadening line shape in sample s*(t) from the observable variable x*(t). The synthesis sample x*(t) of overcomplete ICA basis is labeled as 1, the original Lorentz-pressure broadening line shape of 500 ppb NH3 from the HITRN database is labeled as ^ The original Lorentz-pressure broadening line shape of 6% H2O from the HITRN database is 2, and the separated feature component of 500 ppb NH3 from 1 is labeled as 2. ^ The peak values of 2^ and 3^ are the expected source PA signals sj(lc) in s*(t), hidden labeled as 3, and the separated feature component of 6% H2O from 1 is labeled as 3. within x*(t) or x(t), where j¼ 1,2,y,n. In the test, the overcomplete ICA basis x*(t) is sampled at a full width (FW) of 0.2 nm. FW is far larger than FWHM (b).

582

Z. Wu et al. / Optics & Laser Technology 48 (2013) 580–588 A

A

distance from sampleNdm to sampleN , which is less than N but not equals to 0, otherwise, the separated feature PA signals have the same absorption peak without separated physically; the projection sample B is corresponding to feature component of 500 ppb B

NH3, samplei denotes any an element of the sample B, i¼1, 2, y,

corresponding recursive criterion rwi are respectively determined by    2 NG ðzÞ  k E GðzÞ EGðzGaussian Þ ð13Þ   GðzÞ ¼ exp z2 =2

A sampleN ¼max([xi(t)]),

which is the maximum of B. N, and In order to trim xn(t) into xi  1(t) and xi(t), a trimming criterion must be determined. Here, the criterion first is empirically determined by

Dl ¼ 9lsampleB lsampleA N

Ndm

9or Dn ¼ 9nsampleB nsampleA N

9

ð8Þ

Ndm

where Dl (unit: nm) in Eq. (8) is called Dl-parameter, lsampleA

Ndm

A

(unit: nm) is the wavelength of sampleNdm , lsampleB (unit: nm) is N

A

the wavelength of sampleN , lsampleA

Ndm

and lsampleB are also called N

quasi-feature wavelengths for making a distinction with the real feature wavelengths, here, and nsampleA and nsampleB are the wave

numbers

(in

A

cm  1)

Ndm

respectively

xðtÞ ¼ Asn ðtÞ

ð9Þ

where Eq.(9) is called as the BSS model based on complete ICA bases. It has been proven that when the estimated source variable is clearly clustered, or has a very ‘‘spiky’’ probability density function (PDF), a slightly modified version of the definition of differential entropy, which called negentropy in BSS [20], can be used for measuring this kind of ‘‘sharp’’ absorption peak by a principle stating that what is the most nongaussian is the most independent [14]. The use of the complete-ICA-bases vector x(t) as the input data to estimate the vector sn(t), some preprocessings for the x(t), such as normalizing the vector x(t), centering the normalized vector x(t) by zero mean, and whitening the centered vector x(t), will be made respectively ¼

xðtÞ ¼ ½xi1 ðtÞ, xi ðtÞT N

¼ xc ðt Þ ¼ xnorm ðtÞ zðtÞ N þ N whitening N þ N centering N þ N

ð10Þ

N

ð11Þ

where xðtÞ is the complete-ICA-bases vector with 2N samples, NþN

xnorm ðtÞ denotes the normalized complete-ICA-bases vector with NþN

2N samples, xc ðt Þ denotes the centered complete-ICA-bases vecNþN

tor with 2N samples, and zðtÞ denotes the whitened completeNþN

ICA-bases vector with 2N samples. Consequently, Eq.(9) can be described further by whiting vector z(t) zðtÞ ¼ Asn ðtÞ

where z is the short symbol for the vector z(t), NG(z) is negentropy of z, k is an constant, z Gaussian is a Gaussian variable with N(0,1) distribution, G(z) is a non-quadratic function, X is a symbol of grads, wi is called weight, g(.) is the derivative of G(z), E{.} is symbol of expectation, and b ¼E{G(z)} E{G(zGaussian)} is an estimation constant. It been proven that when the estimating progress for Eq. (12) goes into a steady state, then Eq. (15) equals to 0, i.e., rwi ¼0 [22]. In this work, recursive error e is fixed at 0.000000 in the simulation. 2.3. Determination of weight trunicate equation

After the trimming criterion based on Dl-parameter was fixed as shown in Eq. (8), a group of complete ICA bases xi  1(t) and xi(t) just is determined, let x(t)¼ [xi  1(t), xi(t)]T, and Eq.(6) can be described by

normalized N þ N

ð15Þ

to

A

2.2. Estimation of the BSS model based on complete ICA bases x(t)

NþN

  @N G ðzi Þ ¼ bE zgðwTi zÞ @wi

N

corresponding

sampleNdm and sampleN .

xnorm ðtÞ

rwi p

ð14Þ

ð12Þ

In order to estimate the BSS model as in Eq.(12) using FastICA recursive algorithm based on negentropy approximation [22], a recursive contrast function based on negentropy approximation NG(z) [26], the resultant non-quadratic function G(z), and the

Five-group samples are selected from Fig. 1(a), main test parameters are listed in Table 1, and superposition weights (projection coefficients) wi  1 and wi as shown in Fig. 1(b) from numeric simulation are listed in Table 2. Three separation-error curves as shown in Fig. 2 show that the tested data in Group 2 have a minimal separation error at its curve extremum point. Because a line width of typical vibration– rotation band is 0.05 cm  1 [21–23], then the contributions of source feature absorptions with pressure broadening lineshapes in the vector sn(t)0 to the observable xn(t) is not effective unless they distribute by the spacing less than 0.05 cm  1. In Table 1, the Dl parameter in Group 2 is 0.087 nm, i.e.,0.0376 cm  1, which less than 0.05 cm  1, and it is just an optimal at around 296 K and P¼1 atm. In order to further seek a limit relation to the transformation of ICA bases, a parameter LWT, which called weight truncation, is introduced into. Let LWT ¼ N  dm and dm a0(see the specification about dm in Section 2.1). It had been proven that excessive variation in the weights is harmful for estimation, and then weight truncation is frequently employed to control the weight variation and in turn control the variance inflation [24]. Weight truncation (more appropriately Winsorization) is a kind of weight modification strategy with often used for robust estimation in the presence of outliers [24–27]. With inspirited by the Golden section law, a series of mathematical inductions are made as follows: Group 1. N ¼14503, dm ¼ 1279, LWT ¼N  dm ¼13234 LWT/dm ¼10.4287; N/LWT ¼1.0959;(LWT/dm)  (N/LWT)¼9.5161; (LWT/dm)/(N/LWT)¼0.1833 [(LWT/dm)/(N/LWT)]/[(LWT/dm) (N/LWT)] ¼1.0196 Table 1 Five-group samples of x*(t) with the pressure broadening lineshape as shown in Fig. 1(a). The separation errors in Table 1 only correspond to the sample of interest A. Group

1 2 3 4 5

N (points)

dm (points)

Dl

Dn

(nm)

(cm  1)

NH3 separation error (10  26 cm  1/ mol cm  2)

14,503 14,513 14,523 14,543 14,573

1269 1279 1289 1309 1539

0.0086 0.0087 0.0088 0.0089 0.0091

0.0373 0.0376 0.0379 0.0385 0.0394

0.0258 0.0180 0.0621 0.1500 0.2797

Z. Wu et al. / Optics & Laser Technology 48 (2013) 580–588

583

Table 2 Weights (or projection coefficients): wi  1 and wi from the numerical simulation as illustrated in Fig. 1(b). Group 1

Group 2

Group 3

Group 4

Group 5

wi  1

wi

wi  1

wi

wi  1

wi

wi  1

wi

wi  1

wi

0.7951

0.6065

0.7950

0.6066

0.7948

0.6069

0.7943

0.6076

0.9974

0.0718

Fig. 2. Curves of feature separation error in each group x*(t), plotted respectively against the parameters. (a) N, (b) dm, and (c) Dl. Note that the separation error here is also the absolute deviation error.

Group 2. N¼14503, dm ¼1279, LWT ¼N dm ¼13234 LWT/dm ¼10.4287; N/LWT ¼1.0959; (LWT/dm) (N/LWT) ¼9.4356; (LWT/dm)/(N/LWT) ¼0.1851 [(LWT/dm)/(N/LWT)]/[(LWT/dm)  (N/LWT)]¼1.0200 Group 3. N ¼14523, dm ¼1289, LWT ¼ N  dm ¼13234 LWT/dm ¼10.2669; N/LWT ¼1.0974;(LWT/dm) (N/LWT) ¼9.3557; (LWT/dm)/(N/LWT) ¼0.1862 [(LWT/dm)/(N/LWT)]/[(LWT/dm)  (N/LWT)]¼1.0203 Group 4. N ¼14543, dm ¼1309, LWT ¼N dm ¼13234 LWT/dm ¼10.1100; N/LWT ¼1.0989;(LWT/dm)  (N/LWT)¼9.0111; (LWT/dm)/(N/LWT) ¼0.2001 [(LWT/dm)/(N/LWT)]/[(LWT/dm)  (N/LWT)]¼1.0210 Group 5. N ¼14573, dm ¼1539, LWT ¼ N  dm ¼13034 LWT/dm ¼0.46914; N/LWT ¼1.11087;(LWT/dm) (N/LWT) ¼7.2511; (LWT/dm)/(N/LWT) ¼0.5748 [(LWT/dm)/(N/LWT)]/[(LWT/dm)  (N/LWT)]¼1.04464 where the expression [(LWT/dm)/(N/LWT)]/[(LWT/dm) (N/LWT)] is defined as truncation ratio, its values are listed in Table 3. From Table 3, the optimal truncation ratio of 1.0200 is obtained, and it remains in Group 2. So, the ratio is described as 1

ðLWT dm Þ=ðNL1 WT Þ 1 LWT dm NL1 WT

¼ 1:0200

ð16Þ A

A

where dm is the distance from sampleNdm to sampleN , LWT ¼N  dm is the defined weight truncation, which just expresses the truncating location of the global maximum A

sampleNdm in the interesting sample A. physically, LWT should be larger than dm (dm a0). Institute dm as in Eq. (16) with N  LWT, and a weight truncation equation can be created as L3WT 1:02NL2WT 1:02N2 LWT þ1:02N3 ¼ 0

ð17Þ

Table 3 1

1

1 N, LWT, dm, and truncation ratio ððLWT dm Þ=ðNL1 WT Þ=LWT dm NLWT Þ.

Group

N

LWT

dm

1 ððLWT dm Þ=ðNL1 WT Þ=LWT dm NLWT Þ

1 2 3 4 5

14,503 14,513 14,523 14,543 14,573

13,234 13,234 13,234 13,234 13,034

1269 1279 1289 1309 1539

1.01960 1.02000 1.02030 1.02100 1.04464

1

1

Table 4 Optimal parameters as computed by the created weight truncation equation of Eq. (17). Group

N (points)

LWT (points)

dm (points)

Dl (nm)

1 2 3 4 5

14,503 14,513 14,523 14,543 14,573

13,234 13,233 13,242 13,261 13,288

1269 1280 1281 1282 1285

0.0087 0.0087 0.0087 0.0087 0.0087

where N is the sample length of selected complete ICA bases and is unique known prerequisites. 2.4. Test of five samples using the created weight truncation equation Insert N¼14503, 14513, 14523, 14543, and 14573 into Eq. (15), a group of optimal parameters can be figured out respectively with listed in Table 4. A new feature separation procedure is described as (the recursive error e remain fixed at 0.0000001) N-Eq: ð17Þ-ðN,LWT Þ-xn ðtÞ

-

transformation

xðtÞ-Eq: ð9Þ

Estimation

-

FastICA

features ð18Þ

The superposition weights and separation errors (accuracy) for the interesting sample A are listed in listed in Table 5.

584

Z. Wu et al. / Optics & Laser Technology 48 (2013) 580–588

Table 5 Weights (wi  1 and wi) and separation errors with computed by the proposed BSS model of Eq. (6). Group 1

Group 2

Group 3

Group 4

Group 5

wi  1

wi

wi  1

wi

wi  1

wi

wi  1

wi

wi  1

wi

0.79499

0.60657

0.7950

0.6066

0.795001

0.60662

0.795001

0.6066

0.79489

0.60659

Separation error (10  26cm  1/mol cm  2) 0.0179 0.0180

0.01801

0.0180

0.0178

Fig. 3. Schematic diagram of the photoacoustic apparatus based on TEDFL (a near-infrared tunable erbium-doped fiber laser), where lock-in amplifier is used and microphone has sensitivity 22 (mV/Pa), the TEDFL emits from 1520 nm to 1610 nm.

2.5. Five-point sampling method for accurately truncating complete ICA bases x(t) Just as the discussed in subsections 2.1 to 2.3, overcomplete ICA basis xn(t) must accurately transformed into complete ICA bases x(t) (x(t) ¼[ xi  1(t), xi(t)]T ¼ [A,B]) if the proposed BSS model based on overcomplete ICA basis can be estimated successfully. Then, this transformation can be made by five-point sampling method based on weight truncation equation as shown in Eq. (17). After the weight truncation LWT was figured out with Eq. (17), the five-point sampling of vector x(t) can be described as   xðtÞ ¼ ½A,B ¼ Að1Þ,AðLWT Þ,AðNÞ,Bð1Þ,BðNÞ ð19Þ where A is the interesting sample, B denotes the projection A A A sample, A(1)¼sample1 , A(Lwt)¼ sampleNdm ¼sampleLWT , A(N)¼ A B B sampleN , B(1)¼ sample1 , B(N)¼ sampleN , see Fig. 1(a), and Procedure (18) can be regarded as the estimating process of features based on the five-point sampling method.

3. Validity test for feature extraction and concentration inverse using the samples from real-time breath ammonia detection 3.1. Sample selection of overcomplete ICA basis vector xn(t) Four-group measured data from the experiment on a real-time detection of breath ammonia with high concentration of CO2 and H2O are selected [11]. It is a second harmonic gas PAS detection at 37 1C and 1 atm, the gas system is made up of a certified mixture containing 1 ppm of NH3 buffered in nitrogen (N2) and four different content of certified mixture of about 6% H2O. The experimental setup is schematically depicted as shown in Fig. 3, where a near-infrared tunable erbium-doped fiber laser (TEDFL) is used as probe laser. The TEDFL output is divided into two branches with a 1  2 fiber splitter. About 1% of the optical power enters a hydrogen cyanide (Wavelength References, Standard Gas Cells) gas reference cell, which is used for calibrating laser wavelength to the selected multi-target gases absorption lines, and about 99% of the laser power is amplified by a computer-controlled erbium-doped fiber

amplifier (EDFA, Amonics, AEDA-27-BFA). The output of available EDFA has up to around 500 mW with the same linewidth and wavelength as the seed laser. 3.2. From xn(t) to the extracted feature PA signals in the vector sn(t) Wavelength-modulation spectroscopy is a well established technique to increase the signal-to-noise ratio (SNR) [28–31]. Wang et al. used second harmonic gas PAS detection with probed by tunable wavelength laser to further improve the SNR of the detected PA signal and to increase detection limit up to 16 ppb [11]. However, for Wang et al.’s second harmonic gas PAS detection, PA signal profile at least has two modulations, one is the pressure broadening modulation at about 300 K 1 atm, and another is second harmonic modulation. A tunable laser (TEDFL) is used for probing gas PAS, this causes the broadened absorption peak to absorb the energy from the wavelength modulation and to convert it into the increased-amplitude modulation [19], i.e., a higher SNR can be obtained in the second harmonic gas PAS detection. For a consideration on utility, the sample xn(t) should be more close to the feature absorption peak, and thus it can be selected within the HWHM of second harmonic profile as shown in Fig. 4(a), which is a small sample space. Fig. 4(b) illustrates the extracting process of feature PA signals using procedure (18), including the transformation from overcomplete ICA basis sample xn(t) into the complete ICA-basis sample x(t) based on five-point sampling method as shown in Eq. (19), the extraction of feature PA signals, the removal of background noise, and two demodulations of the modulations. 3.3. Second harmonic and pressure broadening modulations According to by the publication document [15], a single and normalized feature PA signal in second harmonic gas PAS detection can be further deduced as sj ðlc Þ ¼

g L ðnc ÞA2 ðnc Þ 2 nc þ ðk1 nc þ k2 ÞLa Ia ðlc Þcac Þ a2

ð20Þ

where sj ðlc Þdenotes the separated PA signal with modulated both by pressure broadening and by second harmonic. A2(nc) is the

Z. Wu et al. / Optics & Laser Technology 48 (2013) 580–588

585

Fig. 4. (a) The small sample space ‘‘a ’’ of overcomplete ICA basis x*(t), where ‘‘a ’’is far less than the HWHM and complete-ICA bases x(t) ¼[xi  1(t), xi(t)]T. (b) The five-point sampling for determining complete-ICA bases x(t) (x(t) ¼[xi  1(t), xi(t)]T ¼ [B,A]¼ {B(1),yy, B(N),A(1), yy, A(LWT), yy, A(N)}) within the sample space ‘‘a’’ and the resultant processing steps; tag 1 and tag 2 are two separated second harmonic PA signals with strong background noise; tag 3 and tag 4 are two PA signals with removed the background noises, but the fluctuating etalon fringe noises of about 0.25 mV/W from the optical instrument remain superposed on their peaks; tag 5 is the separated feature PA signal by weight w2 and it remain a mixture of NH3 and H2O PA signals at wavelength lc, tag 6 is the separated feature PA signal by weight w1 and it also remain a mixture of NH3 and H2O PA signals at wavelength l0 c, and they were removed the second harmonic modulation.

amplitude of separated PA signal with modulated only by second harmonic at feature frequency nc, see tags 3 and 4 as shown in Fig.4(b), a denotes modulation depth, Ia is the absorption intensity of the detected gas, La is the absorption path length, i.e. the length of PA cell, c is the concentration of detected gas, ac is the absorption coefficient at center frequency nc , which is also known as feature frequency of the detected gas, k1 and k2 are constants depending on a given experiment, and g L ðnc Þ is the maximum of normalized lineshape function. The normalized Lorentzian function is the singly peaked function with described by [32]. According to Eq. (20), it can be seen that the PA signal sj ðlc Þ contains, apart from a base frequency component nc , a second harmonic modulation n2c , a pressure broadening modulation g L ðnc Þ, and the pressure broadening modulation can be expressed as g L ðnÞ ¼ Z

1

p

gððnnc Þ2 þ g2 Þ1

ð21Þ

þ1

g L ðnÞdn ¼ 1

ð22Þ

1

where g is the HWHM of Lorentzian lineshape function, and the maximum gL(nc) of the function gL(n) is 1/pg.

Sample xn(t)

After the extracted feature PA signal sj ðlc Þ with the two modulations as in Eq. (20) was removed the second harmonic modulation accordingly and a molar absorption coefficient was introduced into, and then the signal can be normalized by a single laser line in unit volume as

NH3 1 ppm (%) N

LWT

dm

NH3 þH2O (w2) H2Oþ NH3 (w1)

H2O 6.15 6.29 6.43 6.56

807 807 807 807

79 79 79 79

74.9998 76.7157 76.9096 76.9799

886 886 886 886

67.1301 68.6511 70.1819 71.6042

without removal. In this detection, T¼310 K, lc ¼1522.44 nm,

g ¼0.0376 cm  1 as in Table 1. The every separated PA signal sj(lc) as in Eq. (23) is listed in Table 6, j¼ 1(in column 7 in Table 6), 2 (in column 6 in Table 6). According to tag 5 and tag 6 as shown in Fig. 4(b), it can be given as 8 0 0 > sNH3 ðlc Þ þ sH2 O ðlc Þ ¼ w1 > > > > cNH IaNH ðl0c Þ 0 0 > > < sNH3 ðlc Þ  sH2 O ðlc Þ ¼ cH O3IaH O3 ðl0 Þ c 2 2 3 ð24Þ > sNH3 ðlc Þ þ sH2 O ðlc Þ ¼ w2 > > > > cNH IaNH ðlc Þ > > : sNH3 ðlc Þ  sH2 O ðlc Þ ¼ cH O3IaH O3 ðlc Þ 2 3

8 0 cNH IaNH ðlc Þ 0 2 0 > < sðlc Þ w1 sðlc Þ þ cH O3IaH O3 ðl0c Þ ¼ 0 2

2 3

c I ðl Þ > : sðlc Þ2 w2 sðlc Þ þ c NH3I aNH3 ðlc Þ ¼ 0 c H O aH O 2

ð25Þ

2 3

Solving the two quadratic equations as Eq. (25), it can be obtained

ð23Þ

where sj(lc) is the PA signal magnitude of the jth feature absorption (unit: mV/W), F is the constant of PA cell, cj is the concentration of the jth detected gas, Ia(lc) is absorption intensity of the detected gas (cm  1/mol cm  2) at lc, k is Boltzmann constant (1.3806503  10  23 m2 kg s  2 K  1), T is Kelvin temperature(unit: K), and gL(lc)¼gL(nc)¼1/pg, which the pressure broadening modulation

Quasi-feature PA signal (mV/W)

Optimal parameter (Points)

2

3.4. Concentration inverse of the detected gas

cj Ia ðlc Þ sj ðlc Þ ¼ F g L ðlc Þ kT

Table 6 Optimal parameters and the separated quasi-feature PA signals ‘‘w1’’ and ‘‘w2’’ as shown in Fig. 4(b), which both remain as a mixture like tag 5 and tag 6 as shown in Fig. 2(b).

0

sH2 O ðlc Þ ¼

0

sNH3 ðlc Þ ¼

ffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 cNH IaNH3 ðlc Þ w1  w21 4 c 3I ðl0 Þ H2 O aH2 O3

c

2 w1 þ

ð26Þ

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 cNH IaNH3 ðlc Þ w21 4 c 3I ðl0 Þ H2 O aH2 O3

2

c

ð27Þ

586

Z. Wu et al. / Optics & Laser Technology 48 (2013) 580–588

Table 7 Resulting inverse concentrations of pure-feature PA signal sj(lc) which are hidden in x*(t) as shown in Fig. 4(a), where w01 ¼ w1(l0 c) þ w1(lc), and w02 ¼w2(l0 c) þw2(lc). Sample xn(t)

Feature PA signal(mV/W)

Gas concentration and inverse error

NH3

1 ppm (%)

NH3 (w02 :mV/W)

H2O (w01 :mV/W)

NH3 (ppm)

Error (%)

H2O (%)

Error (%)

H2O

6.15 6.29 6.43 6.56

23.6707 23.9505 23.9087 23.55731

74.9277 75.7153 75.8096 75.9371

0.9899 1.0017 1.0009 0.9847

1.01 0.17 0.09 1.53

6.1504 6.3067 6.4200 6.5227

0.065 0.266 0.156 0.570

sNH3 ðlc Þ ¼

sH2 O ðlc Þ ¼

MDL ¼ st ðn1, a ¼ 0:01Þ

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cNH IaNH ðlc Þ w2  w22 4 cH O3IaH O3 ðlc Þ 2

2 3

2 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi cNH IaNH ðlc Þ w2 þ w22 4 cH O3IaH O3 ðlc Þ 2

2

2 3

ð28Þ

ð29Þ

4. Discussion 4.1. Signal to noise ratio (SNR) and method detection limit (MDL) Signal processing aims at extracting information from the raw signal. The difficulty in reaching this goal depends both on the characteristics of the noise-free signal and the noise. As an estimation of analytical accuracy and confidence level, SNR is introduced into signal processing and detection. When the signal is a transient line (e.g., a feature absorption peak of small molecule), the ratio of the mean around the maximum (xmax ) and the standard deviation (s) of the measured signal is defined as SNR [33] ð30Þ

In the numeric simulation for separating the feature absorption with hidden in pressure-broadening profile in Section 2, the analytical H2O feature absorption standard deviation s¼0.010462437, the analytical H2O feature absorption SNR¼675.5498, the analytical NH3 feature absorption s¼0.023179119, the analytical NH3 feature absorption SNR¼90.1886, and he lowest analytical accuracy is 0.0178% as in Table 5; while in the validity test for extracting feature PA signal with real-time data based on second harmonic detection in Section 3, the extracted NH3 feature PA signal standard deviation

ð31Þ

where s is the standard deviation of the replicate analyses, and n is the number of replicate analyses, and t(n  1, a ¼ 0.01) is the Student’s t-test value for the 1  a (99%) confidence level for n replicate analyses [34] t¼

where sH2 O ðlc Þ and sNH3 ðlc Þ are respectively the amplitude of H2Ofeature-PA signal with demodulated the second-harmonic modulation and the amplitude of NH3-feature-PA signal with demodulated the second-harmonic modulation at wavelength of 1522.44 nm, IaNH3 ðlc Þ is the absorption intensity of NH3 at 1522.44, IaH2 O3 ðlc Þ is IaNH3 ðlc Þ is the absorption intensity of 0 0 0 H2O at 1522.44; lc ¼1522.40 nm, sH2 O ðlc Þ and sNH3 ðlc Þ are respectively the amplitude of H2O-feature-PA signal with demodulated the second harmonic modulation and the amplitude of NH3feature-PA signal with demodulated the second harmonic modulation at the weighted-mean wavelength of 1522.40 nm, 0 IaNH3 ðlc Þis the absorption intensity of NH3 at 1522.40, 0 IaH2 O3 ðlc Þ is the absorption intensity of H2O at 1522.40, and cH2 O and cNH3 are respectively the volume–ratio concentration of H2O and NH3 in this detection. By querying HITRAN database, it can been found that CO2 has very weak contribution to xn(t) at wavelength of 1522.44 nm, so, it was kept out of the consideration in feature separation, if any, it can be fixed by adding windows configuration as Appendix A in the behind. Table 7 gives the inversed concentration of the resultant feature PA signals sj(lc)(j ¼1,2) with removed the pressure broadening modulation.

SNR ¼ xmax =s

s¼0.008363412, the extracted NH3 feature PA signal SNR¼118.8869, and the lowest analytical accuracy is 1.53% as in Table 7. The method detection limit MDL, as defined by the USEPA, is the minimum concentration that can be measured and reported as greater than zero at the 99-percent confidence level (US Environmental Protection Agency, 1997). The MDL was calculated using the formula [34]

xm0 pffiffiffi s= n

ð32Þ

where x is sample mean,m0 is the specified value in certain detection. In this study, for the concentration inverse of gas NH3 based on feature PA signal extraction, the MDL was calculated as 10 ppb (m0 ¼1 ppm, n¼ 4,t ¼1.18046, s(NH3) ¼ 0.008363412). By contrast of Wang et al. [11], the MDL is further improved from MDL¼ 16 ppb to 10 ppb. SNR and MDL depend on each other, SNR describes the effect of random error on a particular measurement, and estimates the expected accuracy of a series of measurements. Samples spiked in the appropriate range for a MDL determination typically have a SNR of not less than 2.5, and thus a highly analytical accuracy can be reached in for the spiked signal. Generally, SNR may always be far higher than 10 for a spiked signal, and a SNR less than 2.5 indicates that the random error in a series of measurements is too high, and the determined MDL is probably high [33,34]. The calculated results of SNR in this section showed that a highly improved SNR was obtained when the measured samples were processed with the proposed BSS model and the five-point sampling method based on the created weight truncation equation, and consequently, the MDL also was improved accordingly. 4.2. Noise immunity Shown as tags 1–4 in Fig. 4(b), not only are there many random background noise but some outliers. By digital filter, many more noises can be removed. However, about 0.25 mV/W fluctuating etalon fringe noises (like tag 3 and tag 4 as shown in Fig. 4(b)) from the optical instrument in real detection, although very weaker than the extracted feature PA signal, are superposed on the extracted feature PA signal by further modifying its maximal peak value. The test results show that these etalon fringe noises can not affect on the feature extraction. So, the FastICA algorithm based on negentropy has excellent noise immunity than other linear equations or the similar analytical method based on linear transformation. 4.3. Dl-parameter and subsample length N as well as weight truncation LWT The replicate test results in this study show: at about 296 K 1 atm, the optimal Dl-parameter is 0.087 nm (equivalent to 0.0376 cm  1), which just is the representation of typical pressure broadening[21–23]. According to the Shannon sampling theorem, for the separation of feature absorption components as in Section 2, if amplitude difference between the separated sample A and B is 0.1, a downsampling need to be made, and then the subsample length N can be 13 and the weight truncation LWT at least is 12; if the separated magnitude difference is 0.5040, an upsampling need to be made, and then N can be 70 and LWT can be 64, at least; For the feature extraction as in Section 3 whose samples against

Z. Wu et al. / Optics & Laser Technology 48 (2013) 580–588

587

weight truncation equation show a strong noise immunity, and thus the method detection limit of detected gas has been largely improved under a larger SNR.

Acknowledgments This research was supported by a grant from the National Natural Science Foundation of China (No. 60677010).

Appendix A

Fig. 5. Schematic drawing of the feature extraction for multi-component gas mixture (three or more than three gases), by adding windows based on the adjacent-feature-wavelength, the determined sample x*(t) in window 1 can be processed using the proposed processing procedures as procedure (18), the determined sample x*(t) in window 2 can similarly be processed by procedure (18), and the features in three-component gases sample just are extracted. This method based on adding window configuration is available for more than three gases’ sample.

different concentrations background (H2O), N can be 100–1000. However, when N¼886, the inverse concentration accuracy is almost approximate to saturation, and thus, as a compromise between analytical accuracy and processing speed, N is fixed 886 and the resulting weight truncation LWT is 79. 4.4. Feature extraction for intricate gas PAS sample Strictly speaking, when the measured sample of gas PA signal is free of CO2, due to the vibration–translation relaxation time of the most gas molecules is fixed in around the order of magnitude of 10  6 s (at 1 atm) [9], the produced PA signals have the same phase with the probe laser, and they can be expressed as the linear sum described as in Eq. (1) or Eq. (2). But when the detected gas sample contains a certain concentration of CO2, the so-called kinetic-cooling effect makes the phase of produced PA signal lag behind that of the probe laser, or even phase flip [9], the proposed BSS model need to be estimated with FastICA algorithm involving convolution processing, and the created weight truncation equation, the proposed BSS model base on overcomplete ICA basis and the resulting processing procedures remain available after the convolution calculation was added in the FastICA algorithm.

5. Conclusions In this work, a BSS model and a truncation equation for processing overcomplete ICA basis are built up. The resulting procedure is designed for applicability to make a data analysis to the low mass molecule gas samples with the ‘‘sharp’’ rovibrational–transition-absorption peak [21], by the scaleinvariant theory [35], this proposed BSS model and the created weight truncation equation can effectively extract the features remained in the statistic sample’s ‘‘entropy’’ or ‘‘negentropy’’ though these features are hidden due to pressure broadening or harmonic modulation. Especially, for the second harmonic gas PA signal detection with the large background noise, the proposed model and the resultant five-point sampling method based on the

The created weight truncation equation as Eq. (17) is derived from the superposition projection of two-component feature components. Multi-component gas sample (three or more than three gases) can be processed by adding windows configuration as shown in Fig. 5 before using the created truncation equation and then the xn(t) in each window just can be further processed with the equation and resultant processing procedures.

References [1] Alexander Graham Bell. Upon the production of sound by radiant energy. Philosophical Magazine 1881;11:510–28. [2] Sun Mingjian, Feng Naizhang, Shen Yi, Li Jiangang, Ma Liyong, Wu Zhenghua. Photoacoustic image reconstruction based on Bayesian compressive sensing algorithm. Chinese Optics Letters 2011;9:061002–5. [3] Zhang Wang, Wu Zhiyingu, Qingxu Y. Photoacoustic spectroscopy for fast and sensitive ammonia detection. Chinese Optics Letters 2007;5:677–9. [4] Rosencwaig A, Gersho A. Photoacoustics and photoacoustic spectroscopy. Journal of Applied Physics 1976;47:64–9. [5] McDonald FA, Wetsel Jr. GC. Generalized theory of the photoacoustic effect. Journal of Applied Physics 1978;49(4):2313–22. [6] McDonald FA. Photoacoustic effect and the physics of waves. American Journal of Physics 1980;48:41–7. [7] McClelland JF. Analytical Chemistry 1983;55:89A–105A. [8] McClelland JF, Luo S, Jones RW, Seaverson LM, Bicanic D, editors. 1992;69:113–24. [9] McClell JF, Jones RW, Luo S, Seaverson LM. In: Coleman PB, editor. Practical sampling techniques for infrared analysis. Boca Raton: CRC Press; 1993. p. 107–44. [10] Yang Xiao long, Xu Yu Qing, Li Shao cheng. Photoacoustic measurement of multi-component gas mixtures applying nonlinear fitting algorithm. Optoelectronics  Laser 2003;14(2):164–7. [11] Wang J, Zhang W, Li L, Yu Q. Breath ammonia detection based on tunable fiber laser photoacoustic spectroscopy. Applied Physics B 2011;103:263–9. [12] Pichler Arthur, Sowaa Michael G. Independent component analysis of photoacoustic depth profiles. Applied Spectroscopy 2005;59(2):164–72. [13] Amari S-I, Cichocki A, Yang HH. A new learning algorithm for blind source separation. Advances in neural information processing systems 8 (NIPSn95). Cambridge, MA. 1996: pp. 757–763. [14] Bell AJ, Sejnowski TJ. An information-maximization approach to blind separation and blind deconvolution. Neural Computation 1995;7:1129–59. [15] Daylight Solutions, Inc., all rights reserved. A brief comparison of 2f and broadly swept detection techniques. /www.daylightsolutions.netS; 2007. pp. 1–5. [16] Bell AJ, Sejnowski TJ. The ‘independent components’ of natural scenes are edge filters. Vision Research 1997;37:3327–38. [17] Pajunen P. Blind separation of binary sources with less sensors than sources. In: Proceedings of the International Conference on Neural Networks. Houston (Texas); 1997. [18] Olshausen BA, Field DJ. Sparse coding with an overcomplete basis set: a strategy employed by V1? Vision Research 1997;37:3311–25. [19] Lewicki M, Sejnowski TJ. Learning overcomplete representations. Neural Computation 2000;12:337–65. [20] Hyv€arinen A, Oja E. Independent component analysis: algorithms and applications. Neural Networks 2000;13(4):411–30. [21] Goody RM, Yung YL. Atmospheric radiation. Theoretical basis. New York: Oxford University Press; 1989 [chapter 2]. [22] Liu LH, Jiang J. Inverse radiation problem for reconstruction of temperature profile in axisymmetric free flames. Journal of Quantitative Spectroscopy and Radiative Transfer 2001;70:207–15. [23] Garcı a-Cuesta E, Galva n IM, de Castro AJ. Neural networks and spectral feature selection for retrieval of hot gases temperature profiles. In: Proceedings of the International Conference on Computational Intelligence for Modelling, Control and Automation; 2005. pp. 81–86.

588

Z. Wu et al. / Optics & Laser Technology 48 (2013) 580–588

[24] Lee H. Outliers in business surveys. In: Cox B, Binder D, Chinnappa B, Christian-son A, Colledge M, Kott P, editors. Business survey methods. New York: John Wiley; 1995. p. 503–26. [25] Elliott RE, Little RJA. Model-based alternatives to trimming survey weights. Journal of Official Statistics 2000;16:191–209. [26] Chinnappa N, Christianson A, Colledge M, Kott P, editors. New York: Wiley; 1996. [27] Potter F. A study of procedures to identify and trim extreme sample weights. In: Proceedings of the Survey Methods Section, American Statistical Association; 1990. pp. 225–230. [28] Arndt Rolf. Analytical line shapes for Lorentzian signals broadened by modulation. Journal of Applied Physics 1965;36(8):2522–4. [29] Reid J, Labrie D. Second-harmonic detection with tunable diode lasercomparison experiment and theory. Applied Physics B 1981;26(3):203–10.

[30] Dharamsi AN, Lu Y. High-resolution spectroscopy using high order derivative techniques. Applied Physics Letters 1994;65(18):2257–9. [31] Dharamsi AN, Lu Y. Sensitive density-fluctuation measurements using wavelength-modulation spectroscopy with high-order-harmonic detection. Applied Physics B 1996;62(3):272–8. [32] /http://mathworld.wolfram.com/LorentzianFunction.htmlS. [33] /http://www.statistics4u.info/fundstat_eng/cc_signal_noise.htmlS. [34] Analytical detection limit guidance and laboratory guide for determining method detection limits. Wisconsin Department of Natural Resources Laboratory Certification Program; April 1996. [35] Brown M, Lowe DG. Invariant features from interest point groups. In: Proceedings of the British machine vision conference. Cardiff (Wales); 2002. pp. 656–665.