Methods for determining a level of a medium in a container, computer readable media, and probing systems

Abstract

A computer implemented method for determining a level of a medium in a container comprises the following steps carried out by computer hardware components: sending a plurality of pulses through a probe in the container; sensing responses caused by reflections of the pulses; determining a statistical model related to the medium in the container based on the responses; and determining the level of the medium in the container based on the statistical model.

Description

The present disclosure relates to methods for determining a level of a medium in a container, computer readable media, and probing systems.

Determining the level of a medium (for example a fluid, in particular a liquid, such as oil or water, or a mixture of different fluid) in a tank or container with high accuracy may be an important task in various applications.

The fluid as referred to herein may be a liquid.

The level of the medium may be determined from Time Domain Reflectometry (TDR) sensor data, which may also be referred to as an echo curve.

FIG. 1 shows an illustration 100 of an echo curve 102 to illustrate the TDR concept. A horizontal axis 118 illustrates time (for example in ns), and a vertical axis 120 illustrates an amplitude, for example a voltage (for example in V). Determining the level 114 of the medium 108 in a container 106 (which may also be referred to as tank) may be achieved by sending a pulse down a TDR probe 104 (which may be straight and vertical), which is inserted into the container 106, and recording the resultant reflections, which collectively form the echo curve 102. The peak 122 at time τ₀ may correspond to the time of the sent pulse (corresponding to level 110 in the container 106). The peak 124 at time τ₁ may correspond to a reflection at the transition from the connector of the probe 104 to the container 106 (and may be a reflection due to impedance mismatch), and may correspond to level 112 in the container 106. The peak 126 at time τ₂ may correspond to a reflection at the surface of the medium 108 in the container 106, and may refer to the level 114 of the medium in the container 108. The peak 128 at time τ₃ may correspond to a reflection from the end of the probe 104, and may correspond to a level 116 in the container 108. The desired ToF 130 to determine the level 114 of the medium in the container 106 is τ₂−τ₀. By determining the ToF 130, and knowing the speed of the signal emitted by an emitter of the probe 104, the distance between the emitter of the probe 104 and the surface of the medium may be determined.

As explained with reference to FIG. 1, the medium level may be determined based on estimating the locations of positive reflection(s) in the echo curve, i.e. a Time of Flight (ToF) estimation task. A positive reflection may be a reflection related to a peak facing upwards (like peak 126). A negative reflection may be a reflection related to a peak facing downwards (like peaks 122, 124, 128). For example, the most positive reflection may occur at the medium surface. Time of flight (ToF) may then be mapped into a medium level, for example depending on the probe length. This however, may be challenging when the medium level is low (so that a reflection from the probe end may overlap with a reflection from the medium), and/or when there is foam above the medium surface, and/or the medium is layered (for example when the medium is a mixture of different fluids, in particular a liquids, with different densities lying layer on layer, which may lead to several positive reflections).

Presently used methods may use a thresholding procedure where a positive value is defined and a reflection in the echo curve exceeding this level is considered to be from the medium level. Its location (in other words: its time-instant) may be used to estimate the ToF. This may be combined by an interpolation method to improve the ToF estimation. Such a method, however, cannot reliably measure the medium level when the fluid level within the container is low (for example less than 7 cm) and/or the fluid has foam above it and/or the medium is a layered medium.

Presently used systems may address TDR performance issues for low fluid levels with TDR sensor probes having a vertical and horizontal part. Thus, mechanical changes to conventional TDR sensors with a straight probe may be made. This, however, may limit the applicability of the TDR sensors due to a complicated shape of the probe. For example, the horizontal part of the probe cannot be easily instrumented within typical containers. Moreover, it is a costly solution and difficult to manufacture.

Other presently used system may formulate the ToF estimation problem with sampling below Nyquist rate using finite rate of innovation approach. This, however, may lead to degradation in performance due to using sub-Nyquist rates which are no longer a limitation with other innovations, e.g. time-equivalent-sampling, enabling data acquisition at the required rates within longer time windows, and may require a substantial change to the hardware of the TDR sensors (for example for established legacy designs).

In presently used methods, the TDR sensor may send a narrow(er) pulse, i.e. wideband in the frequency domain, to address challenges due to overlapping pulses in the echo curve. This, may however complicate and increase the cost of the hardware involved, for example the processing chain, and may have limited impact on the achieved accuracy for low fluid levels.

It is an object of the present invention to provide methods, computer readable media and probe systems for determining a level of a medium in a container. In particular, methods, computer readable media and probe systems shall be provided which have the above mentioned properties of being able to reliably determine a level of a medium in a container even under difficult circumstances like low medium level, layered medium or foam on the medium.

The object is satisfied by a computer implemented method for determining a level of a medium in a container according to claim 1, a non-transitory computer readable medium according to claim 15, and a probing system according to claim 16. Embodiments are given in the subclaims, the description and the drawings.

In one aspect, the present disclosure is directed at a computer implemented method for determining a level of a medium in a container, comprising the following steps carried out by computer hardware components: sending a plurality of pulses through a probe in the container; sensing responses caused by reflections of the pulses; determining a statistical model related to the medium in the container based on the responses; and determining the level of the medium in the container based on the statistical model.

In other words, a statistical model may be used for determining the level of the medium. This may allow a robust detection of the level of the medium.

The statistical model may represent an echo curve. It has been found that instead of analysing a single echo curve, analysing a statistical model of the echo curve may increase robustness of the determination of the level of the medium.

The statistical model may comprise a probability distribution. The probability distribution may represent probabilities of peaks.

The statistical model may be determined based on fitting parameters of the probability distribution to the sensed responses. The parameters may comprise one or more of: mean, amplitude, and number of peaks.

The statistical model may be determined using a Markov Chain Monte Carlo interference technique.

Determining the statistical model may comprise adding peaks (which may be referred to as a “birth” of a peak), removing peaks (which may be referred to as a “death” of a peak), or maintaining a number of peaks. In any case (adding a peak, removing a peak, or maintaining the number of peaks), data such as location of each peak and amplitude of each peak may be updated based on measurement data.

A prior (in other words: a starting assumption on the probability distribution) and/or proposal distributions (in other words: (probability) distributions which may be used to update the estimate for the probability distribution) in the statistical model may be determined based on a peak proximity parameter (PPP, which may also be referred to as δ). It has been found that avoiding detecting peaks too close to each other may increase the robustness of the determination of the level of the medium.

The computer implemented method may further comprise determining point estimates for locations and amplitudes of (for example positive) pulses based on the statistical model. The level of the medium in the container may be determined based on the point estimates.

The computer implemented method may further comprise determining locations of the reflections based on the point estimates. The level of the medium in the container may be determined based on the locations of the reflections. The level of the medium in the container may be determined using at least one pre-determined rule based on the locations of the reflections.

The computer implemented method may further comprise tracking of the level of the medium in the container. This may increase robustness of the determination of the level of the medium by including a history of previous measurement data.

The computer implemented method may further comprise decomposing the statistical model into a plurality of reflections. The computer implemented method may further comprise identifying the medium based on the plurality of reflections.

The computer implemented method may further comprise identifying a medium situation based on the plurality of reflections. The medium situation may comprise at least one of a situation with one and only one layer of medium, a situation with a plurality of layers of medium, or a situation with presence of foam.

According to another aspect, a non-transitory computer readable medium may comprise instructions which, when executed by a computer system, make the computer system carry out the computer implemented method as described herein.

According to another aspect, a probing system may comprise: an emitter configured to send a plurality of pulses through a probe in a container; a sensor configured to sense responses caused by reflections of the pulses; and a determination circuit configured to determine a statistical model related to a medium in the container based on the responses, and to determine a level of the medium in the container based on the statistical model. The emitter and the sensor may be provided in a shared device, which may be referred to as “probe”, and may share parts; for example, some parts of the probe may be used while emitting, and the same parts may be used for detection.

The methods as described herein may also be applied to radar sensors (where also overlap of pulses may occur), for example freely emitting radar sensors, for example without any specific probe to transmit pulses. Similarly, a probing system using radar pulses may be used.

Exemplary embodiments and functions of the present disclosure are described herein in conjunction with the following drawings, showing schematically:

FIG. 1 an illustration of an echo curve to illustrate the TDR concept;

FIG. 2A an illustration of the fitted curves for the model using δ and the model not using δ over the data according to various embodiments;

FIG. 2B an illustration of the approximate posteriors, for the final data set according to various embodiments;

FIG. 3A an illustration of peaks from all samples of the model with δ according to various embodiments, overlaid with the data and the ground truth level;

FIG. 3B an illustration of peaks from all samples of the model without δ according to various embodiments, overlaid with the data and the ground truth level;

FIGS. 4A and 4B illustrations of the full posterior according to various embodiments for two peaks;

FIGS. 5A and 5B illustrations of the full posterior according to various embodiments for three peaks;

FIGS. 6A and 6B illustrations of the full posterior according to various embodiments for four peaks;

FIG. 7 a flow diagram illustrating a method for determining a level of a medium in a container according to various embodiments; and

FIG. 8 a computer system with a plurality of computer hardware components configured to carry out steps of a computer implemented method for determining a level of a medium in a container according to various embodiments.

According to various embodiments, an advanced Bayesian estimation method may be provided for the estimation of the locations of all reflections in the echo curve and their amplitudes. With these methods, the medium level may be estimated accurately even if it is very low and layered.

Accurate estimation of the medium level may be achieved by fitting the echo curve into a statistical model.

As will be described in more detail below, at any time instant t_k, the echo curve may be described by the sum of an unknown number N of Gaussian pulses ϕ(t_k|μ_n, σ²), wherein the means μ_n (in other words: location along time axis) and scaling values A_n are to be estimated, and the variance σ² (in other words, the width of the sent pulse) is known. Thus, the echo curve at a time instant t_k may be expressed as

$y_k=x\left(t_k\right)+{\epsilon }_k,x\left(t_k\right)=\sum _{n=1}^N{A}_n\varphi \left(t_k|{\mu }_n,{\sigma }^2\right)+{\epsilon }_k,x\left(t_k\right)$

where Gaussian noise ε(t_k) is added to represent measurements errors.

By defining statistical priors on the unknowns, i.e. A_n (for example Gaussian with a standard deviation of 1000), μ_n (for example uniform distribution with non-zero support from zero to the latest possible reflection from the probe end) and N (for example a discrete set of possibilities such as N∈{3,6,10,25,30}, etc.):

• • a) Markov Chain Monte Carlo (MCMC) inference technique may be applied to determine the full probability distribution, i.e. posterior p(a, μ, N|x) from the M available echo curve samples x=[x(t₁), . . . , x(t_M)]′ with a=[a₁, . . . , a_N]′ and μ=[μ₁, . . . , μ_N]′. According to various embodiments, various design choices (as described in more detail below), including selecting priors and sampling-pruning procedures, may be made to ensure the MCMC is robust and its implementation is computationally efficient.
  • b) Point estimates for the positive pulse locations and amplitudes may be obtained (e.g. means from the posterior); and
  • c) Reflection from the medium level may be obtained from the point estimates using fixed rules (for example assuming that the largest positive reflection earlier than the largest negative reflection, i.e. from probe end, pertains to the medium level).

According to various embodiments, a reversible jump Markov chain Monte Carlo method may be provided for efficient inference for the general problem of pulse fitting. In particular, the potential of an adopted parametric model overfitting to the (noisy) data via the inclusion of a peak proximity parameter may be reduced. This may facilitate learning a more representative underlying model and may reduce the computational cost. Synthetic and real data may be used to demonstrate the efficacy of the introduced Bayesian technique.

“Reverse Jump MCMC (RJ-MCMC) inference algorithm for TDR with a peak proximity parameter” may be a principled formulation for determining the number of pulses in the echo curve.

Pulse fitting may refer to a problem occurring in many applications, such as Nuclear Magnetic Resonance spectroscopy and Time Domain Reflectometry (TDR). The goal of pulse fitting may be to estimate the parameters of underlying pulses from which the observed data are generated. For example, amplitudes may be used to inspect the faultiness of wires. In TDR, which follows the same principles as guided radar, an electromagnetic pulse may be sent down a probe and the combination of attained reflections yields the echo curve. Another TDR instrumentation, with numerous uses in the food, beverage, and pharmaceutical industries, may be employed to measure the medium/fluid levels in sealed industrial containers. Here, the location of the pulse reflected by the fluid surface may correspond to its level in the container (so can be seen as a Time of Flight (ToF) inference problem). Estimating this pulse location may be especially challenging when the reflection pertaining to the fluid surface is severely obfuscated by other reflections, for instance, from the probe's end when the fluid level is low (i.e., the peak at τ₂ can be too low to detect or be concealed by that at τ₃, as illustrated in FIG. 1) and/or for mediums with a low dielectric constant, which are less reflective and/or in case the reflections of the fluid overlap with reflections from the mechanical container construction.

For simplicity, and without loss of generality, details may be provided for a one-dimensional case. It will be understood that a higher dimensional case may be provided accordingly.

Each datum, y_k, may be assumed to be a discretised observation of a continuous curve of pulses, x(t), with additive Gaussian noise. To further simplify, pulses may be illustrated to be restricted to (unnormalised) Gaussian shapes of equal and known variance, σ². Other pulse shapes, such as log-normal or Lorentzian may also be used. Thus, for all k=1, . . . , K,

$\begin{array}{cc}{y}_k=x\left(t_k\right)+{\epsilon }_k,x\left(t_k\right)=\sum _{i=1}^N{A}_i\varphi \left(t|{\mu }_i,{\sigma }^2\right),& \left(1\right)\end{array}$ ${ϵ}_k\stackrel{\mathrm{idd}}{~}\left(0,{\sigma }_y^2\right),\varphi \left(t|\mu ,{\sigma }^2\right)=e^{-\frac{{\left(t-\mu \right)}^2}{2{\sigma }^2}},$

where t_k is the k^th sampling time of t=(t₁, . . . , t_K)^T and σ_y² is the measurement noise variance. The quantities of interest are the peak locations, μ=(μ₁, . . . , μ_N)^T, and the amplitudes, A=(A₁, . . . , A_N)^T, with N the number of peaks, which is unknown. The methods according to various embodiments may determine μ, and/or A, and/or N.

As used herein, vectors are printed in bold, and scalars and matrices are printed non-bold.

According to various embodiments, a peak proximity parameter (PPP), which may be referred to as δ, may be introduced into both the prior and the proposal distributions of a Markov Chain Monte Carlo (MCMC) method. Markov Chain Monte Carlo are known as such and are widely used.

According to various embodiments, a Reversible Jump MCMC (RJ-MCMC) inference scheme for pulse fitting may be provided. As described herein, methods both with and without use of the PPP are compared to demonstrate the value of this addition on both synthetic and real data.

In the frequency domain (and/or in the frequentist domain), pulse fitting may be handled by a non-linear least squares approach, or with a variation thereof, such as a weighted least squares method, which focuses on fitting the most meaningful peaks. Within a Bayesian framework, a RJ-MCMC scheme may be utilised, where the reversible jump may be used to generate samples with different numbers of peaks. To achieve a variation in peak numbers, a birth-and-death MCMC scheme may be used, where additional peaks may be introduced (corresponding to “birth”) or where peaks may be removed (corresponding to “death”). Due to the data-driven proposal distributions used, dimension-changing steps (for example “birth” step or “death” step) being rejected too often (which would make the sampling inefficient) may be largely circumvented, as proposals may have high likelihood contributions.

RJ-MCMC may be an extension of MCMC which allows so-called ‘reversible jump’ moves. These may yield a parameter's sample having different dimensionality to the previous. As described herein, (single) peak birth and death move types may be used, as well as a fixed dimension Metropolis-within-Gibbs move type, derived from the standard Metropolis-Hastings (MH) and Gibbs sampler schemes, as will be described in more detail below.

Vectorising equation (1) yields y=FA+ε where F:=F(t, μ) is a matrix satisfying F_ki=ϕ(t_k|μ_i, σ²). Linearity in A is beneficial, but is not available for μ. A Metropolis step may be used to sample μ, which may require a proposal distribution to be chosen. The prior, p(μ|N), may typically be uninformative, and the likelihood complex, so the proposal is normally a random walk MH. Slice sampling may be used to propose new samples for each μ_i. This method may provide a replacement for the random walk MH that does not require tuning. It may depend on the current μ_i sample, and so can produce several consecutive samples which yield a low acceptance probability. According to various embodiments, a proposal distribution may be provided which avoids this dependence, which may allow for better mixing of the RJ-MCMC sampling.

One way to avoid the complexity of allowing dimension-changing steps is to assume excessive dimensionality, in this case, assuming a fixed number of peaks which is more than the curve could comprise. The hope is that the amplitude samples of excess peaks are close to zero. As will be shown further below, peaks may cancel or combine to fit the data better as it can use the extra complexity to better explain the noise. However, this results in the samples poorly capturing the true curve. This remains a concern even when using RJ-MCMC schemes. However, according to various embodiments, the PPP may be provided to avoid this overfitting. Another method to overcome it is to propose peak merging and peak splitting steps. The strategy proposed here is simpler, by just preventing peaks from being near enough that a split or merge step would be necessary.

The RJ-MCMC scheme according to various embodiments will be described in the following. This involves specifying the reversible jump steps, which may be of birth and death varieties, as well as the fixed dimension step type. Within this last, which is responsible for the majority of the exploration of the parameter space of each dimension, a Metropolis-within-Gibbs scheme may be provided.

Let j index the S samples, {μ^j, A^j}, the first B of which are discarded. Then, with bracketed tildes being optional, define

$\begin{array}{cc}{\stackrel{(~)}{\mu }}_i^j={\left({\mu }_{1:i-1}^{j^T},{\stackrel{(~)}{\mu }}_i^j,{\mu }_{i+1:N^j}^{j-1^T}\right)}^T,{\tilde{\mu }}_b^j={\left({\mu }^{j-1^T},{\tilde{\mu }}_b^j\right)}^T,& \left(2\right)\end{array}$

where {tilde over (μ)}_i^j is the j^th proposed sample for μ_i as is {tilde over (μ)}_b^j, but for a new peak (in other words: a new-born peak, indicated by index “b”). Define these similarly for A (where relevant). Also, let ω_−I=ω\ω_I, for any vector ω. Finally, the number of peaks, N, may have a Poisson prior, whose parameter is λ: p(N)˜Po(λ).

The choice of prior for the amplitudes is chosen to yield tractability, and so is chosen to be Gaussian,

$\begin{array}{cc}A❘N~\left(m_A,{\sum }_A\right).& \left(3\right)\end{array}$

By standard Gaussian theory, the posterior becomes

$\begin{array}{cc}p\left(A❘\mu ,y,N\right)=\left(A❘{\tilde{m}}_A,{\sum ^{~}}_A\right),& \left(4\right)\end{array}$ ${\tilde{m}}_A=m_A+{\sum }_A{F}^T{\left(F{\sum }_A{F}^T+{\sigma }_y^2{I}_K\right)}^{-1}\left(y-{\mathrm{Fm}}_A\right),$ ${\sum ^{~}}_A={\sum }_A-{\sum }_A{F}^T{\left(F{\sum }_A{F}^T+{\sigma }_y^2{I}_K\right)}^{-1}F{\sum }_A,$

where I_∘ is the ∘×∘ identity matrix. Then, p(A_i|A_−i, μ, y, N), which can be derived by more Gaussian properties, becomes

$\begin{array}{cc}p\left(A_i❘A_{-i},\mu ,y,N\right)=𝒩\left(A_i❘m_i,{\sigma }_i^2\right),m_i={\tilde{m}}_{A_i}+{\tilde{\Sigma }}_{{𝒜}_{-\mathrm{ii}}}^T{\tilde{\Sigma }}_{{𝒜}_{-i-i}}^{-1}\left(A_{-i}-{\tilde{m}}_{A_i}\right),{\sigma }_i^2={\tilde{\Sigma }}_{{𝒜}_{\mathrm{ii}}}-{\tilde{\Sigma }}_{{𝒜}_{-\mathrm{ii}}}^T{\tilde{\Sigma }}_{{𝒜}_{-i-i}}^{-1}{\tilde{\Sigma }}_{A_{-\mathrm{ii}}}.& \left(5\right)\end{array}$

It is from this that most amplitude samples will be generated.

The prior on the locations may be chosen to be independently uniform. However, in another embodiment, a modification to this may be provided, for example by using of the peak proximity parameter, δ. For the interval, L=[t₁, t_K],

$\begin{array}{cc}p\left(\mu ❘N\right)\propto 1\left(\mu \in S_N\right),S_N=\left\{{\mu }^{\prime }\in L:❘{\mu }_i^{\prime }-{\mu }_{i^{\prime }}^{\prime }❘>\delta ,\forall i,i^{\prime }\le N\right\}& \left(6\right)\end{array}$

with 1(∘) being the indicator function.

Using this (with or without δ), the posterior is intractable, jointly or conditionally on the other locations. Thus, it may be desirable to define a proposal distribution for peak i, q_μ, from which to get samples. For computational efficiency whilst retaining sample quality, use of the data within the proposal may be essential. As such, q_μ may be based on a discrete convolution between a discretised Gaussian peak and the ‘unexplained’ data, placing high probability of peak placement near to areas of insufficient data explanation. At all t_k within δ of any μ_i′, i′≠i, the proposal, to match the prior, is set to zero. Then,

$\begin{array}{cc}{q}_{\mu }\left({\mu }_i=t_k❘{\mu }_{-i},A_{-i},y,N\right)\propto ❘\sum _{k^{\prime }=1}^K{r}_{k^{\prime }}^i\varphi \left(t_{k^{\prime }}❘t_k,{\sigma }^2\right)❘\times 1\left({\left({\mu }_{-i}^T,t_k\right)}^T\in S_N\right),& \left(7\right)\end{array}$

where r_k′ⁱ=y_k′−F(t_k′, μ_−i)A_−i is the k′^th ‘unexplained’ datum without using peak i. As possible locations are continuous, q_μ is defined continuously by linear interpolation. With a slight abuse of set notation, let {tilde over (t)} be, in ascending order, t and all points exactly δ from a current peak: {tilde over (t)}=t∪{τ∈L|∃i′≠i:τ=μ_i′+δ}. Hence, q_μ({tilde over (t)}_k| . . . )=q_μ(t_k′| . . . ) if {tilde over (t)}_k=t_k′, and zero if {tilde over (t)}_k∉t, is well defined. Thus, for t∈[{tilde over (t)}_k, {tilde over (t)}_k+1), we have

$\begin{array}{cc}{q}_{\mu }\left({\mu }_i=t❘{\mu }_{-i},A_{-i},y,N\right)=\frac{{\tilde{t}}_{k+1}-t}{{\tilde{t}}_{k+1}-{\tilde{t}}_k}{q}_{\mu }\left({\tilde{t}}_k❘\dots \right)+\frac{t-{\tilde{t}}_k}{{\tilde{t}}_{k+1}-{\tilde{t}}_k}{q}_{\mu }\left({\tilde{t}}_{k+1}❘\dots \right).& \left(8\right)\end{array}$

The normalising constant for (7), Q*, may be computed as

$\begin{array}{cc}{Q}^{*}=\frac{1}{2}\sum _{k=2}^K\left({\tilde{t}}_k-{\tilde{t}}_{k-1}\right)\left(q_{\mu }\left({\tilde{t}}_{k-1}❘\dots \right)+q_{\mu }\left({\tilde{t}}_k❘\dots \right)\right).& \left(9\right)\end{array}$

Now the RJ-MCMC step types may be described: birth, death, and fixed dimension. For the fixed dimension type at step j (thereby accepting N^j=N^j-1), each peak i may be sampled,

$\begin{array}{cc}{\tilde{\mu }}_i^j~q_{\mu }\left({\mu }_i❘{\mu }_{-i}^j,A_{-i}^j,y,N^j\right),{\tilde{A}}_i^j~p\left(A_i❘A_{-i}^j,{\tilde{\mu }}_i^j,y,N^j\right).& \left(10\right)\end{array}$

Both μ_i^j={tilde over (μ)}_i^j and A_i^j=Ã_i^j may be accepted with probability α_i^j=min(1, {tilde over (α)}_i^j), where, after cancelling all distributions of the unchanged peak parameters, and the location priors, one finds

$\begin{array}{cc}{\tilde{\alpha }}_i^j=\frac{p\left(y❘{\tilde{\mu }}_i^j,{\tilde{A}}_i^j,N^j\right)}{p\left(y❘{\mu }_{i-1}^j,A_{i-1}^j,N^j\right)}\frac{p({\tilde{A}}_i^j❘A_{-i}^j,{\tilde{\mu }}_i^j,N^j}{p\left(A_i^{j-1}❘A_{-i}^j,{\mu }_{i-1}^j,N^j\right)}\times \frac{p\left(A_i^{j-1}❘A_{-i}^j,{\mu }_{i-1}^j,y,N^j\right)}{p\left({\tilde{A}}_i^j❘A_{-i}^j,{\tilde{\mu }}_i^j,y,N^j\right)}\frac{q_{\mu }\left({\mu }_i^{j-1}❘{\mu }_{-i}^j,A_{-i}^j,y,N^j\right)}{q_{\mu }\left({\tilde{\mu }}_i^j❘{\mu }_{-i}^j,A_{-i}^j,y,N^j\right)}.& \left(11\right)\end{array}$

To clarify, a new peak i is proposed, and this proposal is accepted (or rejected) before proposing the following peak.

For a birth type at step j, Ñ^j=N^j-1+1 is proposed, as is the new peak, together with all amplitudes jointly

$\begin{array}{cc}{\tilde{\mu }}_b^j~q_{\mu }\left({\mu }_b❘{\mu }^{j-1},A^{j-1},y,{\tilde{N}}^j\right),{\tilde{A}}^j~p\left(A❘{\tilde{\mu }}_b^j,y,{\tilde{N}}^j\right).& \left(12\right)\end{array}$

The apparent redundancy of sampling all amplitudes, rather than just that of the new peak, A_b^j, may be justified by the desire for all A_i^j to be conditioned on this new peak's location. This choice may increase both the birth and death acceptance rates.

For a death step, Ñ^j=N^j-1−1 may be proposed, along with the peak to be deleted, I, which may be done with probability p(i|A^j-1, N^j-1)∝|A_i^j-1|⁻¹ (placing higher probability on less significant peaks), and all remaining amplitudes. That is,

$\begin{array}{cc}{I}^j~p\left(i❘A^{j-1},N^{j-1}\right),{\tilde{A}}^j~p\left(A❘{\mu }_{-I^j}^{j-1},y,{\tilde{N}}^j\right).& \left(13\right)\end{array}$

For the (generalised) acceptance probabilities of these step types, we define a useful function ρ=ρ(μ, A, a, i),

$\begin{array}{cc}\rho =\frac{p\left(y❘\mu ,A,N\right)}{p\left(y❘{\mu }_{-i},a,N-1\right)}\frac{p\left(A❘\mu ,N\right)}{p\left(a❘{\mu }_{-i},N-1\right)}\frac{q_N\left(N-1❘N\right)}{q_N\left(N❘N-1\right)}\times & \left(14\right)\end{array}$ $\frac{p\left(a❘{\mu }_{-i},y,N-1\right)}{p\left(A❘\mu ,y,N\right)}\frac{p\left(i❘A,N\right)}{q_{\mu }\left({\mu }_i❘{\mu }_{-i},a,y,N\right)}\frac{1}{R}\frac{\lambda }{N},$

where N=dim μ, q_N=(N|N′) is the proposal of the step type for N and N′, and R=∫_L1((μ_−i^T, t)^T∈S_N)dt is the length of all points >δ from all peaks. The acceptance probabilities (which include accepting N^j=Ñ^j) for birth and death, β^j=min(1, {tilde over (β)}^j) and γ^j=min(1, {tilde over (γ)}^j), respectively, are then

$\begin{array}{cc}{\tilde{\beta }}^j=\rho \left({\tilde{\mu }}_b^j,{\tilde{A}}^j,A^{j-1},{\tilde{N}}^j\right),{\tilde{\gamma }}^j={\rho \left({\mu }^{j-1},A^{j-1},{\tilde{A}}^j,I^j\right)}^{-1}.& \left(15\right)\end{array}$

For q_N, we choose q_N(N+1|N)=q_N(N−1|N)=c, for some c, when N>0, and q_N(1|0)=1. Finally, initialization may be done with N⁰=0 peaks.

The RJ-MCMC method according to various embodiments may be summarized as follows, wherein S may be the number of sampling steps to be performed:

Input: σ2, λ, S, B, mA, ΣA, σy2, c, δ
for 1 ≤ j ≤ S do
sample step type by sampling Ñj ~ qN(N | Nj−1)
If Ñj = Nj = Nj−1 (fixed dimension) then
for 1 ≤ i ≤ Nj do
sample {tilde over (μ)}ij and Ãij, as in (10)
accept μij = {tilde over (μ)}ij and Aij = Ãij with probability αij,
using (11)
else If Ñj = Nj−1 + 1 (birth) then
sample {tilde over (μ)}bj and Ãj, as in (12)
accept Nj = Ñj, μj = {tilde over (μ)}bj and Aj = Ãj with
probability βj, using (14)
else if Ñj = Nj−1 − 1 (death) then
sample Ij and Ãj, as in (13)
accept Nj = Ñj, μj = μ−Ij and Aj = Ãj with
probability γj, again using (14)
Output: {μj, Aj, Nj}j=B+1S

In the following, experimentation results will be provided.

The method according to various embodiments, using the PPP, is compared against the same model, only without the PPP; this may illustrate the effect of δ. The model without the PPP may only differ from its counterpart in the prior on locations, which is now simply uniform, and the proposal distribution, q_μ, which simplifies by avoiding the need to define {tilde over (t)}, instead being just the linear interpolation of q_μ(t| . . . ).

Comparisons are made on both synthetic 1D (one-dimensional) data and real measurements from a TDR sensor.

Each of the 100 synthetic data sets is generated from the priors of the model without the PPP. From the samples generated, a ‘fitted curve’ can be found:

$\begin{array}{cc}\hat{x}\left(t\right)=\frac{1}{S-B}\sum _{j=B+1}^{S}F\left(t,{\mu }^j\right)A^j.& \left(16\right)\end{array}$

Five metrics are considered. Firstly is the Root Mean Square Error (RMSE) between the true and fitted curves. The next two, ‘overlap’ and ‘cancellation’ averaged over samples, aim to quantify overfitting, by comparison to the same on the true curve. The ‘overlap’ for i, i′ same-sign peaks (i.e., sign(A_i^j)=sign(A_i′^j) is

$\begin{array}{cc}\frac{{\left(\sqrt{2\pi }\sigma \right)}^{-1}}{{\sum }_{i\le N^j}❘A_i^j❘}{\int }_R\underset{x}{\min }\left\{A_i^j\varphi \left(x❘{\mu }_i^j,{\sigma }^2\right),A_{i^{\prime }}^j\varphi \left(x❘{\mu }_{i^{\prime }}^j,{\sigma }^2\right)\right\}\mathrm{dx},& \left(17\right)\end{array}$

including standardisation via division of the total peak area. The ‘overlap’ metric is (17) summed over all same-sign peaks, and averaged over samples. Cancellation is simply the area under the linear interpolation of |F^jA^j|, standardised as for overlap, and averaged over samples. The final metrics are the average number of peaks per sample, and the average run time, in seconds (machine: Apple M1, CPU@3.2 GHz, 16 GB RAM). If a method has much higher overlap and cancellation than the true curve, but fits the data well, then this indicates probable overfitting.

Both methods are given the same hyperparameters, which, where relevant, are the same as those used to generate the data. Namely, (S, B, σ², A, c, σ_y²)=(300, 100, 20, 5, 0.1, 1). All peaks are treated symmetrically and independently, so (m_A, Σ_A)=(0_N, σ_A²I_N), ∀N, with σ_A²=100. Additionally, for the model using it, the PPP is

$\delta =\frac{\sigma }{2}=\sqrt{5}.$

Its relatively small size compared to the peak ‘width’, ≈6σ′, is noteworthy.

TABLE 1
Averaged results across all 100 synthetic data sets.
Method	RMSE	Overlap	Cancel	N	Time (s)
δ used	3.88	0.04	0.07	4.39	2.25
No δ	4.31	0.15	0.39	12.4	6.40
Truth	0	0.04	0.11	4.86	—

The averaged results for all 100 data sets are shown in Table 1.

FIG. 2A shows an illustration 200 of the fitted curves for both models (the model using δ illustrated by line 208; and the model not using δ illustrated by line 210) over the data (illustrated by dashed line 206) according to various embodiments. A horizontal axis 202 illustrates a time index, and a vertical axis 204 illustrates a return (for example corresponding to the echo curve). Both methods fit the data similarly well, with use of δ giving only marginal improvement. The four true peaks have μ≈(94, 158, 168, 174)^T and A≈(−2, 5, 11, 3)^T. Three are very close, but the PPP model (i.e. the model using δ) is still able to fit the data well. However, its use allows for seemingly more accurate representation of the true curve; the overlap and cancellation metrics are far nearer to those for the true curve than without the PPP (i.e. the model not using δ).

FIG. 2B shows an illustration 250 of the approximate posteriors p(N|y), for the final data set according to various embodiments. A horizontal axis 252 illustrates the number of peaks, and a vertical axis 254 illustrates the frequency. The average estimated number of peaks with δ (illustrated by blocks 258) is much closer to the ground truth of 4 (illustrated by dashed line 256) than without (illustrated by blocks 260), slightly underestimating rather than largely overestimating. This underestimation may be reasonable, since true peaks with A≈0 are unlikely to be modelled. Despite the additional complexity of using the PPP, the average run time is much faster, due to the reduction in peaks being modelled. It is to be noted that the code has not been optimised, however should be comparable between the two methods, due to sharing of functions and structure. The most obvious drawback of using the PPP is if true peaks are closer than δ together. This does not appear to be an issue, at least in this experiment (the RMSE may be inflated were it so). If this is systematically the case in an application, then use of δ may be carefully considered there, and the method may be used without the use of δ.

For further evaluation, real unprocessed TDR measurements may be used. This data may be for an oil medium with a dielectric constant of approximately 2, at a low level (24.7 mm) in a sealed container; TDR vertical probe is of length ≈900 mm. The resultant echo curve may be provided with a sampling period of T≈60.9 ps.

TABLE 2
Estimating the fluid level of a TDR set-up, both with
and without use of δ. The known fluid level generates
a peak at 113.5, which is considered the ground truth.
	Method	{circumflex over (μ)}	Error	N	Time (s)
	δ used	111.3	2.2	2.4	1.0
	No δ	109.9	3.6	4.5	2.8

Inference on this data set yields FIG. 3A and FIG. 3B, showing individual peaks which model the data, and performance metrics in Table 2.

FIG. 3A shows an illustration 300 of peaks from all samples of the model with δ according to various embodiments, overlaid with the data and the ground truth level. A horizontal axis 302 illustrates a time index, and a vertical axis 304 illustrates a return (for example corresponding to the echo curve). The peaks are illustrated by solid lines 310, the data is illustrated by dotted line 308, and the time index corresponding to the fluid level is illustrated by dashed line 306.

FIG. 3B shows an illustration 350 of peaks from all samples of the model without δ according to various embodiments, overlaid with the data and the ground truth level. A horizontal axis 352 illustrates a time index, and a vertical axis 354 illustrates a return (for example corresponding to the echo curve). The peaks are illustrated by solid lines 360, the data is illustrated by dotted line 358, and the time index corresponding to the fluid level is illustrated by dashed line 356.

Whilst the ‘true curve’ is unknown, the fluid level's peak is known in the experiment. An estimate of this, {circumflex over (μ)}, is found for both methods and, as the synthetic experiment suggested may be the case, the use of the PPP results in a better representation of the curve, shown by the more accurate fluid level estimate. It is noted that the pulse reflected from the medium level (with positive amplitude) is not clearly visible in FIG. 3A and FIG. 3B, since it heavily overlaps with the stronger negative reflection from the probe's end. As before, the average number of peaks is lesser with the PPP than without, which increases run time. A nonlinear least squares method may struggle to get a sensible estimate for this data set, and does not learn the number of peaks as standard.

When making a crude (but reasonable) approximation to the full posterior

$p\left(\mu ,A,N❘y\right)\approx p\left(\mu ,A❘\tilde{N},y\right)1\left(N=\tilde{N}\right)$

by fixing the number of peaks to the best estimate, Ñ, only one result (for the fixed number of peaks) may be provided.

The methods according to various embodiments may use Reversible Jump MCMC, which may make a better approximation of the full posterior (albeit at higher cost) by varying N in a more complex way.

FIGS. 4A, 4B, 5A, 5B, 6A, 6B show illustrations 400, 450, 500, 550, 600, and 650 of the full posterior according to various embodiments for two peaks (N=2, in FIG. 4A and FIG. 4B), for three peaks (N=3, in FIG. 5A and FIG. 5B), and for four peaks (N=4, in FIG. 6A and FIG. 6B). In FIGS. 4A, 4B, 5A, 5B, 6A, 6B, the plots have means (illustrated along axis 402) and amplitudes (illustrated along axis 404), and a density (illustrated by axis 406, and by grey scale 414). The synthetic data 408 based on which the inference was made is shown for reference.

Each surface is shown with two plots (for two peaks FIG. 4A and FIG. 4B, for three peaks FIG. 5A and FIG. 5B, and for four peaks in FIG. 6A and FIG. 6B), where in FIGS. 4A, 5A, and 6A the height of the peaks can be seen, and in FIGS. 4B, 5B, and 6B a view from above is shown. However, FIG. 4A and FIG. 4B show the same plot (in other words, the same data). Likewise, FIG. 5A and FIG. 5B show the same plot (in other words, the same data). Likewise, FIG. 6A and FIG. 6B show the same plot (in other words, the same data).

FIGS. 4A and 4B show two peaks 410 and 412.

FIGS. 5A and 5B show three peaks 502, 504, and 506.

The four peaks in FIGS. 6A and 6B are so low that they cannot be easily seen.

FIG. 7 shows a flow diagram 700 illustrating a method for determining a level of a medium in a container according to various embodiments. At 702, a plurality of pulses may be sent through a probe in the container. At 704, responses caused by reflections of the pulses may be sensed. At 706, a statistical model related to the medium in the container may be determined based on the responses. At 708, the level of the medium in the container may be determined based on the statistical model.

According to various embodiments, the statistical model may represent an echo curve.

According to various embodiments, the statistical model may include or may be a probability distribution.

According to various embodiments, the statistical model may be determined based on fitting parameters of the probability distribution to the sensed responses.

According to various embodiments, the parameters may include or may be one or more of: mean, amplitude, and number of peaks.

According to various embodiments, the statistical model may be determined using a Markov Chain Monte Carlo interference technique.

According to various embodiments, determining the statistical model may include adding peaks, removing peaks, or maintaining a number of peaks.

According to various embodiments, a prior and/or proposal distributions in the statistical model may be determined based on a peak proximity parameter.

According to various embodiments, the method may further include determining point estimates for locations and amplitudes of (for example positive) pulses based on the statistical model. According to various embodiments, the level of the medium in the container may be determined based on the point estimates.

According to various embodiments, the method may further include determining locations of the reflections based on the point estimates. According to various embodiments, the level of the medium in the container may be determined based on the locations of the reflections.

According to various embodiments, the method may further include tracking of the level of the medium in the container.

According to various embodiments, the method may further include decomposing the statistical model into a plurality of reflections.

According to various embodiments, the method may further include identifying the medium based on the plurality of reflections.

According to various embodiments, the method may further include identifying a medium situation based on the plurality of reflections.

Various embodiments provide a method for pulse fitting with RJ-MCMC via the use of a peak proximity parameter, δ. This prevents pulses from encroaching on one another, reducing overfitting by not considering an excessive number of underlying peaks. The benefits of this approach are demonstrated on both synthetic and real TDR data. Results show that the formulation according to various embodiments not only significantly improves the learning of the latent data model, but also substantially reduces the computational complexity, as captured by the run times. These derive from the reduction of peaks being modelled, thus circumventing the aforementioned overfitting and requiring fewer computations. Various embodiments may be applied to high dimensional signals.

The methods as described herein, including MCMC, may effectively find the best model that fits the observed data (in other words: best fits the TDR echo curve), including the number of pulses it contains. Thus, the method according to various embodiments may be a pulse fitting (optimisation) method, which may also be referred to as pulse separation method.

The methods according to various embodiments provide a solution with the statistical nature, i.e. involving statistical approaches. The possibility of both negative and positive reflections may be addressed by the Bayesian pulse fitting methods according to various embodiments by enabling to learn N (i.e. the number of pulses). The methods according to various embodiments provide an efficient MCMC scheme due to both the marginalised sample acceptance approach (i.e. accept/reject each proposed {tilde over (μ)}_n individually, rather than accepting/rejecting their collective, {tilde over (μ)}), and a data-driven proposal distribution used to sample said pulse means, improving efficiency.

The methods according to various embodiments may accurately estimate medium levels under different scenarios, including for low level fluids and layered mediums. The methods may not require a mechanical change to the TDR sensor probe. The methods may not entail any changes to the processing chain or hardware. The methods may permit incorporating, in a principled way, known prior information (for example based on the probe length and/or type of fluid) to deliver accurate TDR sensing.

According to various embodiments, the sensed medium may be fully characterised by decomposing the echo curve into individual pulses. This may permit additional TDR sensing capabilities such as the automatic identification of the medium (for example identification of the type of the fluid) and the measurement scenario (for example whether a medium with one layer is present or whether a multiple layers medium (in other words: a medium with more than one layer) is present, and/or whether foam is present).

According to various embodiments, the medium level may be tracked over time via extending the MCMC formulation to include a dynamical model on the medium level. This may improve the robustness of the results (for example to mitigate outliers).

FIG. 8 shows a computer system 800 with a plurality of computer hardware components configured to carry out steps of a computer implemented method for determining a level of a medium in a container according to various embodiments.

The computer system 800 may include a processor 802, a memory 804, and a non-transitory data storage 806. A probe 808 (which may include an emitter and a sensor) may be provided as part of the computer system 800 (like illustrated in FIG. 8), or may be provided external to the computer system 800.

The processor 802 may carry out instructions provided in the memory 804. The non-transitory data storage 806 may store a computer program, including the instructions that may be transferred to the memory 804 and then executed by the processor 802. The probe 808 may include (or may act as) an emitter configured to send a plurality of pulses through a probe in a container. The probe 808 may further include (or may further act as) a sensor configured to sense responses caused by reflections of the pulses.

The processor 802, the memory 804, and the non-transitory data storage 806 may be coupled with each other, e.g. via an electrical connection 810, such as e.g. a cable or a computer bus or via any other suitable electrical connection to exchange electrical signals. The probe 808 may be coupled to the computer system 800, for example via an external interface, or may be provided as parts of the computer system (in other words: internal to the computer system, for example coupled via the electrical connection 810).

The terms “coupling” or “connection” are intended to include a direct “coupling” (for example via a physical link) or direct “connection” as well as an indirect “coupling” or indirect “connection” (for example via a logical link), respectively.

It will be understood that what has been described for one of the methods above may analogously hold true for the computer system 800.

REFERENCE NUMERAL LIST

100     illustration of an echo curve
102     echo curve
104     probe
106     container
108     medium
110     level in the container
112     level in the container
114     level of the surface of the medium in the container
118     horizontal axis
120     vertical axis
122     peak at time τ0
124     peak at time τ1
126     peak at time τ2
128     peak at time τ3
130     desired ToF
200     illustration of fitted curves
202     horizontal axis
204     vertical axis
206     dashed line
208     line
210     line
250     illustration of approximate posteriors
252     horizontal axis
254     vertical axis
256     dashed line
258     blocks
260     blocks
300     illustration of peaks
302     horizontal axis
304     vertical axis
306     dashed line
308     dotted line
310     solid lines
350     illustration of peaks
352     horizontal axis
354     vertical axis
356     dashed line
358     dotted line
360     solid lines
400     illustrations of full posterior with two peaks
402     axis
404     axis
406     axis
408     synthetic data
410     peak
412     peak
414     grey scale
450     illustrations of full posterior with two peaks
500     illustrations of full posterior with three peaks
502     peak
504     peak
506     peak
550     illustrations of full posterior with three peaks
600     illustrations of full posterior with four peaks
650     illustrations of full posterior with four peaks
700     flow diagram illustrating a method for determining a level of a
        medium in a container according to various embodiments
702     step of sending a plurality of pulses through a probe in the container
704     step of sensing responses caused by reflections of the pulses
706     step of determining a statistical model related to the medium in the
        container based on the responses
708     step of determining the level of the medium in the container based on
        the statistical model
800     computer system according to various embodiments
802     processor
804     memory
806     non-transitory data storage
808     probe
810     connection
SICK AG S29654PUS - Br/Mx

Claims

1. Computer implemented method for determining a level of a medium in a container, the method comprising the following steps carried out by computer hardware components:sending a plurality of pulses through a probe in the container;sensing responses caused by reflections of the pulses;determining a statistical model related to the medium in the container based on the responses; anddetermining the level of the medium in the container based on the statistical model.

2. The computer implemented method of claim 1, wherein the statistical model represents an echo curve.

3. The computer implemented method of claim 1, wherein the statistical model comprises a probability distribution.

4. The computer implemented method of claim 3, wherein the statistical model is determined based on fitting parameters of the probability distribution to the sensed responses.

5. The computer implemented method of claim 4, wherein the parameters comprise one or more of: mean, amplitude, and number of peaks.

6. The computer implemented method of claim 1, wherein the statistical model is determined using a Markov Chain Monte Carlo interference technique.

7. The computer implemented method of claim 1, wherein determining the statistical model comprises adding peaks, removing peaks, or maintaining a number of peaks.

8. The computer implemented method of claim 1, wherein a prior and/or proposal distributions in the statistical model are determined based on a peak proximity parameter.

9. The computer implemented method of claim 1, further comprising: determining point estimates for locations and amplitudes of pulses based on the statistical model;wherein the level of the medium in the container is determined based on the point estimates.

10. The computer implemented method of claim 9, further comprising: determining locations of the reflections based on the point estimates;wherein the level of the medium in the container is determined based on the locations of the reflections.

11. The computer implemented method of claim 1, further comprising: tracking of the level of the medium in the container.

12. The computer implemented method of claim 1, further comprising: decomposing the statistical model into a plurality of reflections.

13. The computer implemented method of claim 12, further comprising: identifying the medium based on the plurality of reflections.

14. The computer implemented method of claim 12, further comprising: identifying a medium situation based on the plurality of reflections.

15. Non-transitory computer readable medium comprising instructions which, when executed by a computer system, make the computer system carry out the computer implemented method of claim 1.

16. Probing system, comprising: an emitter configured to send a plurality of pulses through a probe in a container;a sensor configured to sense responses caused by reflections of the pulses;a determination circuit configured to determine a statistical model related to a medium in the container based on the responses, and to determine a level of the medium in the container based on the statistical model.