EP3101573B1

EP3101573B1 - Deconvolution method for emissions measurement

Info

Publication number: EP3101573B1
Application number: EP16177610.9A
Authority: EP
Inventors: Frank Berghof
Original assignee: AVL Test Systems Inc
Current assignee: AVL Test Systems Inc
Priority date: 2011-03-28
Filing date: 2012-03-14
Publication date: 2018-08-29
Anticipated expiration: 2032-03-14
Also published as: US20140019077A1; WO2012134815A3; EP2691901A2; EP2691901B1; CA2831593A1; US10520480B2; EP2691901A4; US20180321205A1; JP5932018B2; CN104303017B; EP3101573A1; CN104303017A; JP2014516404A; WO2012134815A2

Description

BACKGROUND

Emissions analyzers, or measurement instruments, measure certain gaseous constituents within a sample of exhaust, or aerosol, as a function of time, or are configured to measure particulate matter, such as soot, within an exhaust sample, as examples. The response of the instrument, however, may be uncorrected for the convolution of the measurement with some other signal representative of the transfer function, or the transient response, of the instrument. Deconvolution is a process used to reverse, or correct, the effects of convolution.
In one known method, the response of an instrument is recorded online in the time domain. Deconvolution of the recorded signal is performed offline in post-processing by (1) decomposing the recorded data, via a Fourier transform, into the frequency domain, (2) using a model to remove the effects convolution, and then (3) constructing a convolution corrected signal, via an inverse Fourier transform, back into the time domain.
WO 2011/004378 A1 discloses a method and system for reconstructing an input field where the latter is sensed by a measurement system. The method comprises: providing measured data corresponding to output field of said measurement system; providing data about sparsity of the input field, and data about effective response function of the measurement system; and processing the measured data based on said known data, the processing comprising: determining a sparse vector as a function of said measured data, said data about the sparsity of the input field, and said data about the effective response function; and using the sparse vector for reconstructing the input information.
US 2008/075347 A1 discloses a method for correcting x-ray scatter in projection radiography or in x-ray computer tomography, in which, for different observed objects, different scatter-correction convolution cores G are determined and in the examination of the respective objects object-specific scatter-correction convolution cores are applied to detector data which is created during the x-raying of these objects.

SUMMARY

Further disclosed is a method according to claim 1 for determining an inverse convolution function for use in correcting a response of an emissions measurement instrument. The method includes determining an idealized convolution function, the idealized convolution function being in the time domain. The idealized convolution function is transformed from the time domain to the frequency domain, and a regularizing filter function is divided by the transformed idealized convolution function. The result of the division is the inverse convolution function in the frequency domain. The inverse convolution function is then transformed from the frequency domain to the time domain.
These and other features of the present disclosure can be best understood from the following drawings and detailed description.

BRIEF DESCRIPTION OF DRAWINGS

The drawings can be briefly described as follows:

Figure 1 illustrates an example system including a measurement instrument configured to respond to exhaust.
Figure 2 illustrates an example method of correcting a response of a measurement instrument.
Figure 3 is representative of the details of the first step from Figure 2.
Figure 4 is representative of the details of the second step from Figure 2.
Figure 5 is representative of the details of the third step from Figure 2.

DETAILED DESCRIPTION

Figure 1 illustrates an example system 10 including an engine 12 and an exhaust pipe 14 downstream thereof. The engine 12 could be an engine of a vehicle, or could be a stand-alone engine in a lab, as examples. The engine 12 could further be any type of engine, including a diesel engine.
Exhaust 20 generated from the engine 12 flows downstream of the engine 12, is tapped at 22a, and a sample 24a of the exhaust 20 is directed to a sampling line 22b. A portion 24b of the sample 24a is directed toward a measurement instrument 26, whereas another portion 24c is directed toward a filter box 28 in parallel with the measurement instrument 26. The filter box 28 need not be present, however.
In this example, the measurement instrument 26 is a soot sensor, such as the AVL 483 Micro Soot Sensor (MSS), for example. The response (or, signal) from the measurement instrument 26 is indicative of a concentration of soot, as a function of time, within the portion 24b of the sample 24a.
A controller 30, which may be any type of known computer, is in communication with the measurement instrument 26 to record the response thereof. As those in the art would appreciate, the controller 30 could include a processor (or, CPU), screen, hard drive, mouse, keyboard, etc. The controller 30 is further configured to perform each of the calculations in the steps described below, and may be configured to communicate with other various components in the system 10.
A reference gas source 32 selectively in communication with the sampling line 22b by way of an adjustable valve 34. The controller 30, in one example, is configured to adjust the valve 34, however the valve 34 could be manually adjustable. In this example, the reference gas is a gas having a known soot concentration. The reference gas source 32 can include an appropriate reference gas, however, as will be appreciated from the below.
Notably, while a soot sensor is shown, this disclosure extends to other types of measurement instruments. For example, this disclosure extends to gas analyzers configured to measure a quantity (e.g., a concentration) of one or more gaseous constituents within a sample of exhaust, such as of CO₂, CO, NO, NO₂, NO_x, CH₄, HC, O₂, NH₃, and N₂O, as examples. The disclosed method can further be used to deconvolute data from any measurement instrument for which a convolution curve can be determined, such as temperature, pressure, flow rate, speed and torque measurements, as examples. The system 10 is likewise non-limiting, and this disclosure extends to other system set-ups, including those mounted for use on-road or in a lab.
Figure 2 shows a high-level overview of the steps in one example of the disclosed method. As shown, the response of the system 10 (specifically, the response of the measurement instrument 26) to a step input signal change is measured, at 100. Then, idealized and inverse convolution functions are then determined at 200 and 300, respectively. Steps 100, 200, 300 can be performed offline, before acquiring data during engine operation.
The results from steps 100-300 are then used in the fourth step, at 400, to deconvolute data acquired by the measurement instrument during engine operation. In one example, this data is acquired during an emissions test. The deconvoluted data can be further refined in an optional fifth step, at 500. Steps 100-500 are discussed in detail below.
As those in the art would immediately acknowledge, a function in the time domain is represented as n(t), for example, while the same function in the frequency domain would be represented as N(f). This notation is used throughout the application.
Figure 3 shows the detail of step 100. At 102-106 a sample of reference gas, which has a known quantity of a measurable exhaust component, is connected to the measurement instrument, via positioning of the valve 34, and an uncorrected response of the instrument x(t) is recorded. As noted, in the example where the measurement instrument 26 is a soot sensor, the reference gas would have a known soot concentration. Likewise, if the measurement instrument was configured to measure HC, a reference gas with a known HC concentration would be selected.
At 108, times T_A, T_B , and T_C are determined. As generally noted, these times
are times at which the amplitude of the recorded signal is at three different percentage values relative to the known signal. This is indicative of the attenuation caused by the measurement instrument and other measurement equipment. In this example, 10%, 50%, and 90% are used, for T_A , T_B , and T_C , respectively.
Figure 4 is representative of the details of step 200, the result of which is the determination of h(t), the idealized convolution function. This function generally represents an approximation of the real convolution function, using a model consisting of the Gauss function convoluted with the impulse response function: $h (t) = g (t) * i (t)$
where g(t) is the Gaussian function is defined as: $\frac{1}{σ \sqrt{2 π}} e^{- (\frac{{(t - μ)}^{2}}{2 σ^{2}})}$
and where i(t) is the impulse response function, defined as: $\frac{1}{τ} e^{- (\frac{t}{τ})}, for t \geq 0.$
At step 202 the ratio $\frac{τ}{σ}$
is determined, which is needed to calculate the normalized convolution function h_n(t) in step 204. In one example the ratio $\frac{τ}{σ}$
is determined from the following equation: $\frac{T_{B} - T_{A}}{T_{C} - T_{A}}$
In another example, a look-up table is used to determine the ratio. The inputs to an example look-up table are T_A, T_B, and T_C.
At step 204 the normalized convolution function h_n(t) is calculated. The normalized convolution function is: $h_{n} (t) = g_{n} (t) * i_{n} (t)$
where g_n(t) is the Gaussian function g(t) from above, with µ=0 and σ=1: $\frac{1}{\sqrt{2 π}} e^{- (\frac{t^{2}}{2})}$
and where i_n(t) is the impulse response function i(t) from above, with τ_n equal to the ratio $\frac{τ}{σ}$
determined in step 202: $\frac{1}{τ_{n}} e^{- (\frac{t}{τ_{n}})}, for t \geq 0$
A scaling factor k is determined at step 206, and is defined as: $k = \frac{T_{C} - T_{A}}{T_{C, n} - T_{A, n}}$
where T_A,n is the time at which ∫h_n (t) reaches A% of its maximum value (in this example, 10%), and where T_C,n is the time at which ∫h_n (t) reaches C% of its maximum value (in this example 90%).
At 208, scaling factor k can be used to determine the parameters σ, µ, and τ of the idealized convolution function h(t) based on the following equations: $σ = k$
$μ = k T_{B, n}$
$τ = k τ_{n}$
where T_B,n is the time at which ∫h_n (t) reaches B% of its maximum value (in this example, 50%). Having solved for these parameters, the idealized convolution function h(t) can then be determined by solving for g(t) and i(t), above.
As an alternative to step 200, the idealized convolution function h(t) could be approximated as the first derivative of the uncorrected instrument response x(t).
Figure 5 generally illustrates the steps for determining the inverse convolution function k(t). At 502 the idealized convolution function h(t) is transformed into the frequency domain by Fourier transformation, as follows: $H (f) = F (h (t))$
Next, as 304, a regularizing filter function R(f) is calculated from the following equation: $R (f) = \frac{{(H_{M A G} (f))}^{2}}{{(H_{M A G} (f))}^{2} + α}$
where H_MAG(f) is the magnitude, or absolute value, of H(f), and where α is a positive adjustable filter parameter. In one example, α is a constant, positive real value. In another example, α is a function of frequency, however a constant value is typically sufficient. As noted below, in the example where α is a constant, α can be tuned to adjust the convolution corrected instrument response y(t).
At 306 the inverse convolution function K(f) is calculated by: $K (f) = R (f) / H (f)$
Notably, R(f) and H(f) may include complex numbers, and thus, in one example, the above division follows the rules for division of two complex numbers and can be performed by dividing the magnitude of R(f) (e.g.,R_MAG(f)) by the magnitude of H(f) (e.g., H_MAG(f)) and subtracting the phase angle of H(f) (e.g., H_PHA(f)) from the phase angle of R(f) (e.g., R_PHA(f)).
At 308, the inverse convolution function K(f) is converted into the time domain by way of an inverse Fourier transformation to determine an initial inverse convolution function k_init(t): $k_{i n i t} (f) = F^{- 1} (K (f)) .$
The regularizing filter function R(f) depends from a positive adjustable filter parameter α, which may be a constant value, and need not be frequency dependent. The positive adjustable filter parameter α is generally representative of a signal to noise ratio.
Once k_init(t) is determined, the uncorrected instrument response x(t) recorded in step 100 is convolved with k_init(t) to construct an convolution corrected instrument response y(t), at step 310, as follows: $y (t) = x (t) * k_{i n i t} (t)$
The convolution corrected instrument response y(t) is then be evaluated relative to the known reference gas signal from step 100, at step 312. In one example this evaluation is performed by graphically comparing the two signals, however this could also be performed using a one-dimensional optimization algorithm to minimize the sum of squares of the deviations between the deconvoluted response and the signal representative of the known data.
As represented in steps 314-318, the positive adjustable filter parameter α can further be adjusted, or "tuned," to increase the accuracy of the inverse convolution function k_init(t), thus increasing the accuracy of the convolution corrected instrument response y(t) relative to the reference gas signal from step 100.
Tuning is dependent on varying the constant positive adjustable filter parameter α, from which k_init(t) depends. The dynamic response (or, slope) of y(t) is assessed at 314, while overshoots and undershoots (e.g., amplitude) of y(t) are accounted for at 316. As an example, increasing α would reduce the slope of y(t) (e.g., worse recovery of dynamic response) but also lower the over and undershoots. Once a desirable α is found (e.g., a value for α representing an acceptable compromise between error in slope and error due to over/undershoots is determined), the corresponding inverse convolution function is saved as k(t), at 320, for later use in step 400.
In step 400, the k(t) saved at 320 is used for deconvolution of the uncorrected instrument response m(t). In step 400, the system 10 would be arranged as shown in Figure 1, for example, such that valve 34 is adjusted so that the sample 24a sourced from the engine 12 is directed toward the instrument 26.
To construct the convolution corrected instrument response y(t), the uncorrected instrument response m(t) is convolved with k(t): $y (t) = m (t) * k (t)$
In one example, the controller 30 executes the convolution of m(t) with k(t) by way of the following Riemann sum: $y_{i} = \sum_{j = 1}^{n} k'_{j} m_{i - (j - 1)}$
where y_i is the i-th value of the convolution corrected instrument response vector, m _{i, (j-1)} is the i-(j-1)-th value of the uncorrected measured instrument response vector, k' is the flipped inverse convolution function in the time domain (as used herein, "flipped" means that the order of the values in the vector is reversed), n is the number of values in the inverse convolution function vector, j is the running index of the inverse convolution function vector, and i is the running index of the uncorrected instrument response vector.
Since the convolution corrected instrument response y(t) is calculated, at 400, entirely in the time domain using multiplication and summation, the calculation at step 400 can be done quickly and efficiently relative to other methods, such methods require transformations between the time and frequency domains. Post processing is thus not necessary with this disclosed method, and the convolution corrected instrument response y(t) can be determined online, during engine operation. Again, as noted above, the controller 30 can be used calculate the convolution corrected instrument response y(t).
At optional step 500, the convolution corrected instrument response y(t) can be further refined to eliminate deviations that may be present at step changes. In one example, this further refinement, called a derivative corrected instrument response p(t), can be calculated by solving for p(t) using the following equation: $p (t) + β (\frac{d p}{d t}) * k (t) = y (t)$
where β is a constant, k(t) is the inverse convolution function from 320, and y(t) is the convolution corrected instrument resulting from 400. In one example, p(t) is solved for iteratively, using y(t) as an initial estimate for p(t). Again, this fifth step is optional, and need not be included.
Although the different examples have the specific components shown in the illustrations, embodiments of this invention are not limited to those particular combinations. It is possible to use some of the components or features from one of the examples in combination with features or components from another one of the examples.
One of ordinary skill in this art would understand that the above-described embodiments are exemplary and non-limiting. That is, modifications of this disclosure would come within the scope of the claims. Accordingly, the following claims should be studied to determine their true scope and content.

Claims

A method of determining an inverse convolution function for use in correcting a response of an emissions measurement instrument (26), the method comprising:
determining an idealized convolution function, the idealized convolution function being in the time domain;

transforming the idealized convolution function from the time domain to the frequency domain;

dividing a regularizing filter function by the transformed idealized convolution function, the result being the inverse convolution function in the frequency domain; and

transforming the inverse convolution function from the frequency domain to the time domain.
The method as recited in claim 1, wherein the idealized convolution function is calculated by convolving a Gaussian function with an impulse response function.
The method as recited in claim 2, wherein the Gaussian and impulse response functions are based on a scaling factor.
The method as recited in claim 3, wherein the scaling factor is determined based on a normalized convolution function and on a response of the instrument (26) to a reference exhaust sample (24a).
The method as recited in claim 4, wherein the normalized convolution function is calculated by convolving a normalized Gaussian function with a normalized impulse response function.
The method as recited in claim 4, wherein the impulse response function is based on values from the response of the instrument (26) to a reference exhaust sample (24a).
The method as recited in claim 1, wherein the regularizing filter function is based on the transformed idealized convolution function and a positive adjustable filter parameter.
The method as recited in claim 7, wherein the positive adjustable filter parameter is a constant value independent of frequency.
The method as recited in claim 7, further comprising:
recording a response of the instrument (26) to a sample (24a) of reference gas;

convolving the recorded response with the transformed inverse convolution function, the result being a convolution corrected instrument response; and

adjusting the positive adjustable filter parameter to adjust overshoots, undershoots, and a dynamic response of the convolution corrected instrument response relative to a reference gas signal.
The method as recited in claim 1, wherein the idealized convolution function is the first derivative of a response of the instrument (26) to a reference exhaust sample.
The method as recited in claim 1, wherein the instrument (26) is a gas analyzer configured to measure a concentration of a gaseous constituent of an exhaust sample (24a) as a function of time.
The method as recited in claim 1, further including calculating a derivative corrected instrument response to eliminate noise at step changes in the convolution corrected instrument response.
The method as recited in claim 1, further comprising:
recording a response of the instrument (26) to an exhaust sample as a function of time; and

convolving the recorded response with the inverse convolution function, the result being a convolution corrected instrument response.
The method as recited in claim 13, wherein a derivative corrected instrument response is calculated by solving for p(t) using the following equation: $p (t) + β \cdot (dp / dt) * k (t) = y (t)$
where p(t) is the derivative corrected instrument response, β is a constant, k(t) is the inverse convolution function, and y(t) is the convolution corrected instrument response.
The method as recited in claim 1, further comprising:
recording a response of the instrument (26) to a sample of reference gas, the reference gas having a known quantity of a measurable exhaust component;

convolving the recorded response with the transformed inverse convolution function, the result being a convolution corrected instrument response; and

comparing the convolution corrected instrument response to a reference gas signal that is representative of the known quantity.