ON A MODEL OF A RADIOMETRIC ISOTOPE FLAW DETECTOR

V. I. Gorbunov, F. M. Zavyalkin, M. S. Kvasnitsa

(Presented by the scientific seminar of the Research Institute of Electronic Introscopy)

One of the possible directions for improving the quality of operation of radiometric instruments is the use, in their development, of methods of the theory of statistical decisions, in particular, the optimal synthesis of isotope flaw-detection systems. In the present case, along with detection of a signal (defect), it is also necessary to solve the problem of estimating its principal parameters (for example, volume). A specific feature of this circle of questions is the indeterminacy of the form of the useful signal. Consequently, when solving the problem of estimation and detection, it is necessary to single out from the set of all possible signals such a class of signals in which their common features can be established. The quality of operation of the system and its structure depend on the form of the signal. This means that, when choosing the form of the useful signal, it is necessary to take into account the stability of the system with respect to small changes in the signal form. On the other hand, it is necessary to choose such a signal form that would be characterized by the smallest possible number, preferably integral, of parameters.

Taking the above into account, the following model representations and assumptions are proposed.

1. The form of the signal \(s(t)\) is approximated by a Poisson form. This assumption is based on the comparability of the length of the collimating aperture with the extent of the defect in the direction of inspection, and also on the comparability of the time during which the defect passes through the field of view of the radiation detector with the value reciprocal to the passband of the smoothing filter.

Normalizing the signal \(s(t)\) on the data-sampling interval \((0, T)\) and assuming that the most essential parts of the signal are contained inside the interval \((0, T)\), we write the functional form of the signal \(s(t)\) as

\[ s(t)=\frac{4A}{a}(t-\tau)\exp\left(-\frac{t-a}{a}\right)1(t-\tau), \tag{1} \]

where \(\tau\) is the moment of appearance of the signal.

The amplitude \(A\) of the signal \(s(t)\) may be regarded as a random variable uniformly distributed on the interval \((0, A_m)\). The parameter \(A\) may be uniquely related to the volume of the defect.

The quantity \(a\) in equation (1), which determines the duration of the signal, is also a random variable and may be related by a one-to-one dependence to the extent of the defect in the direction of inspection. However, the above-indicated form of the signal \(s(t)\) is most justifi-

for local defects. Consequently, the parameter \(a\) in this model is undesirable. We take its probability distribution to be uniform on the interval \([a_0, ka_0]\), \(k \ge 1\).

1. The time \(\tau\) at which the signal \(s(t)\) appears is also a random variable, uniformly distributed over the observation interval \((0,T)\). This assumption may be justified as follows: the frequency of occurrence of defects in the item being inspected is small, i.e., the mean time interval between two signals is significantly greater than the interval \((0,T)\). The flow of signals in this case may be represented as Poisson. Then, neglecting the probability of the appearance of two signals on the interval \((0,T)\), \(\tau\) has a uniform probability-density distribution.

2. The noise \(n(t)\) at the output of an inertial system described by a first-order linear differential equation may be represented by an asymptotically normal process with covariance function

\[ B(\tau)=\frac{\lambda_0\gamma}{2}e^{-\gamma|\tau|}, \tag{2} \]

where \(\lambda_0\) is the intensity of pulses at the output of the radiation detector;
\(\gamma\) is the passband of the inertial system.

Thus, under the indicated model assumptions, our problem consists in determining the structure of a device estimating the parameters \(A,\tau\), on the basis of which a decision is made as to the presence or absence of a defect, and also in evaluating the quality of operation of the given system.

The received oscillation \(r(t)\) under the two hypotheses is written as

\[ r(t)= \begin{cases} s(t)+n(t), & 0\le t\le T:H_1,\\ n(t), & 0\le t\le T:H_0. \end{cases} \]

To find the optimal structure of the radiometric-information processing system, we use its direct derivation on the basis of a sufficient statistic \(r_1\) [1]

\[ r(t)=r_1s(t)+y(t), \tag{3} \]

where \(r_1\) is a random variable obtained by means of certain operations on \(r(t)\);

\(r_1\) and \(y(t)\) are statistically independent;

\(y(t)\) does not depend on which hypothesis is true.

The statistic \(r_1\) is defined by

\[ r_1=\int_0^T r(u)g(u)\,du. \tag{4} \]

The further problem consists in finding a function \(g(t)\) possessing the required properties. From expressions (3) and (4) we obtain

\[ r(t)= \begin{cases} (s_1+n_1)s(t)+y(t):H_1,\\ n_1s(t)+y(t):H_0, \end{cases} \]

where

\[ s_1=\int_0^T As(u)g(u)\,du, \tag{5} \]

\[ n_1=\int_0^T n(u)g(u)\,du. \tag{6} \]

Let us introduce the constraint

\[ 1=\int_0^T s(u)g(u)\,du. \tag{7} \]

Then from expressions (5) and (7), \(s_1=A\). If \(n_1\) is a random variable with zero mean, then

\[ n(t)=n_1s(t)+y(t);\quad 0\leq t\leq T . \]

Let us determine the condition that the function \(g(t)\) must satisfy in order that \(\overline{n_1y(t)}=0\), or \(\overline{n_1n(t)}=n_1^2s(t)\). Using expression (6), this condition can be written as

\[ \int_0^T B(t-u)g(u)\,du = s(t)\int_0^T\int_0^T B(\theta-u)g(\theta)g(u)\,du\,d\theta . \tag{8} \]

Equality (8) is satisfied if

\[ \int_0^T B(t-u)g(u)\,du=As(t). \tag{9} \]

Thus, the sufficient statistic \(r_1\) is obtained by correlating \(r(t)\) with \(g(t)\). The structure of the information-processing system is completely determined by the function \(g(t)\), which is the solution of the integral equation (9) and is nothing other than the impulse response function of the system being determined. Methods for solving equations of this type (Fredholm equations of the first kind) are widely known [2, 3]. The general solution of equation (9) for \(B(\tau)\), determined by expression (2), has the form

\[ g(t)=\frac{A}{2\gamma^3} \left[ -\frac{\partial^2 s(t)}{\partial t^2}+\gamma^2 s(t) \right] +C_1\delta(t)+C_2\delta(t-T), \tag{10} \]

where

\[ C_1=\frac{\gamma s(0)-s'(0)}{\gamma^2}; \quad C_2=\frac{\gamma s(T)+s'(T)}{\gamma^2}. \]

It follows from expression (1) that \(C_1=C_2=0\). Neglecting terms of second order of smallness in expression (10), we obtain

\[ g(t)= \frac{2A[2a-(t-\tau)]}{\gamma^2 a^5} \exp\left\{-\frac{t-\tau}{a}\right\}1(t-\tau). \tag{11} \]

Let us estimate the parameters \(A\) and \(\tau\) of the signal \(s(t)\). The logarithm of the likelihood function is written in general form as

\[ \ln \Phi[r(t)] = \int_0^T r(u)g(u)\,du - \frac{1}{2}\int_0^T s(u)g(u)\,du . \tag{12} \]

The averaged likelihood ratio for our conditions is

\[ \overline{\Phi}[r(t)] = \frac{1}{a_0(k-1)} \int_{a_0}^{ka_0} \exp\left\{ \frac{2A}{\gamma^2} \int_0^T \left[ \frac{2a-(t-\tau)}{a^5} \right] r(t)\times \right. \]

\[ \left. {}\times \exp\left(-\frac{t-\tau}{a}\right) -\frac{A^2}{\gamma^3a^5} \right\}\,da\,dt, \tag{13} \]

For weak signals, expression (13) can be simplified, also assuming that the observation \(r(t)\) is carried out over the interval \((0,T)\), comparable with the maximum possible signal duration \(\bigl(T\sim(3\div5)a_0k\bigr)\).

\[ \Phi_1=A\sum_{n=0}^{m}\gamma_n L_n+\frac{\gamma_1}{3}A^2, \tag{14} \]

where

\[ L_n=\int_{\tau}^{T}(t-\tau)^n r(t)\,dt; \]

\[ \gamma_1=\frac{3(k^3+k^2+k+1)}{4a_0^5k^4};\quad \gamma_n=(-1)^n\frac{(n+2)(k^{n+2}+k^{n+1}+\ldots+1)}{n!a^{n+4}k^{n+3}}; \]
\[ \Phi_1=(\Phi[r(t)]-1)\frac{\gamma\sigma^2}{2}. \tag{15} \]

The estimates \(\hat A\) and \(\hat\tau\) of the parameters of the signal \(s(t)\) are found from the system of equations obtained from expression (14) by differentiating it with respect to the given parameters. We determine the estimate \(\hat A\), and then maximize expression (14) with respect to \(\tau\). The existence condition has the form

\[ \left[\sum_{n=0}^{m}\gamma_nL_n-\frac{2}{3}\gamma_1A\right]_{A=\hat A}=0, \]

and the estimate itself is

\[ \hat A=\frac{3}{2\gamma_1}\sum_{n=0}^{m}\gamma_nL_n. \tag{16} \]

Since the solution is unique and the second derivative is negative, the absolute maximum of the likelihood function exists for any value of \(\tau\). We note that, in nonlinear estimation, in practice there almost always exist many local maxima of the likelihood function. From expressions (14) and (16) we obtain

\[ \left\{\left[\sum_{n=0}^{m}\gamma_nL_n\right]\left[\sum_{n=0}^{m-1}(n+1)L_n\gamma_{n+1}\right]\right\}_{\tau=\hat\tau'}=0. \tag{17} \]

Putting \(m=1\), we obtain the first crude approximation

\[ \hat\tau'= \frac{\displaystyle\int_{0}^{T}tr(t)\,dt-\displaystyle\int_{0}^{\hat\tau'}tr(t)\,dt} {\displaystyle\int_{0}^{T}r(t)\,dt-\displaystyle\int_{0}^{\hat\tau'}r(t)\,dt} -\frac{2}{3}a_0\frac{k(k^2+k+1)}{k^3+k^2+k+1}. \tag{18} \]

The second and subsequent approximations are found as the roots of polynomials of the third and higher degree in \(\tau\) from equation (17).

Let us determine expressions for calculating the lower bound of the mean-square error of the estimates \(\hat A\) and \(\hat\tau\), using the Cramér–Rao inequality [4]. The elements \(I_{ij}\) of the information matrix are defined as follows:

\[ I_{ij}=\int_{0}^{T}\frac{\partial s(t,\theta)}{\partial\theta_i} \frac{\partial g(t,\theta)}{\partial\theta_j}\,dt,\qquad \theta=\{A,\tau\}. \tag{19} \]

We obtain approximate expressions for computing \(\overline{(\tau-\hat\tau)^2}\) and \(\overline{(A-\hat A)^2}\) for the following cases:

a) for known signal width \((k=1)\)

\[ \overline{(\tau-\hat\tau)^2}\geq \frac{\gamma a_0^7}{5}\frac{\sigma^2}{A^2}, \tag{20} \]

\[ \overline{(A-\hat A)^2}\geq \frac{\gamma\sigma^2a_0^5}{44}; \tag{21} \]

b) for the greatest uncertainty in the signal duration \((k\gg1)\)

\[ \overline{(\tau-\hat\tau)^2}\geq 3k\gamma a_0^7\frac{\sigma^2}{A^2}, \tag{22} \]

\[ \overline{(A-\hat A)^2}\geq \frac{\gamma\sigma^2a_0^5}{44}. \tag{23} \]

It follows from expressions (20)–(23) that, when the signal width is known, the only way to increase the accuracy of determining \(\tau\) is to observe signals with large amplitudes. At the same time, the error in determining the amplitude itself is significant. In the case of considerable uncertainty in the signal width, the error in determining the parameter \(\tau\) will be greater than in the first case, while the error in estimating the amplitude will decrease as this uncertainty increases. It should be noted that our conclusions are valid for weak signals; therefore, when attempting to reduce the estimation errors to a lower level, we encounter a threshold phenomenon. The physical interpretation of this phenomenon is as follows: an erroneous interval is selected on which local estimates of the given parameters are made.

Fig. 1 shows the structural block diagram of the system obtained.

Fig. 1. Structural diagram of the algorithm for processing radiometric information.

Next we shall determine the quality of operation of the developed system. From expressions (13)–(17), the test statistic is reduced to the form

\[ \left(\sum_{n=0}^{m} \nu_n L_n\right)^2 \mathop{\gtrless}_{H_0}^{H_1} \frac{2\nu_1}{3}\,\gamma\sigma^2(\Phi-1), \tag{24} \]

or

\[ \left|\sum_{n=0}^{m} l_n\right|^2 \mathop{\gtrless}_{H_0}^{H_1} C,\qquad l_n=\nu_n L_n . \tag{25} \]

Under hypothesis \(H_0\),

\[ \bar l_n=0;\qquad \overline{l_0^2}=\sigma^2\nu_0^2=\mu_0^2; \]

\[ \overline{l_n^2}= \left(\frac{n}{n+1}\right)^2 \frac{T^{2n}}{2n+1}\,\nu_n^2\sigma^2 =\mu_n^2 . \]

Hence one immediately obtains the expression for calculating the conditional probability of a false indication of the presence of a defect,

\[ \alpha= \frac{1}{(2\pi)^{\frac{m+1}{2}}\mu_0\mu_1\ldots\mu_m} \int\limits_{\substack{(\ldots)\\(D)}}\cdots\int\limits^{(m+1)} \exp\left\{-\sum_{n=0}^{m}\frac{l_n^2}{2\mu_n}\right\} \,dl_0\,dl_1\ldots dl_m , \tag{26} \]

where \(D\) is a Euclidean space of dimension \(m+1\), excluding the region \(\Gamma\) lying between the two planes

\[ l_0 + l_1 + \cdots + l_m - C = 0, \]

\[ l_0 + l_1 + \cdots + l_m + C = 0. \]

Under hypothesis \(H_1\),

\[ \overline{I_n} = \frac{AT^n}{n+1}; \]

\[ \sigma^2[I_0] = \sigma^2 = \mu_0^2;\quad \sigma^2[I_n] = \sigma^2 \nu_n^2 \left(\frac{n}{n+1}\right)^2 \frac{T^{2n}}{2n+1} = \mu_n^2. \]

The conditional probability of missing a defect is

\[ \beta = \frac{1}{(2\pi)^{\frac{m+1}{2}}\mu_0\mu_1\cdots\mu_m} \int\cdots\int_{\substack{(m+1)\\(D)}} \exp\left\{ -\sum_{n=0}^{m} \frac{1}{2\mu_n^2} \left(l_n-\frac{AT^n}{n+1}\right)^2 \right\} \times \]

\[ \times\, dl_0\,dl_1\cdots dl_m . \tag{27} \]

Figure 2 shows the decision space for \(m=1\). Obtaining more closed-form expressions for calculating \(\alpha\) and \(\beta\) presents considerable difficulties; therefore, in specific situations it is advisable to use numerical calculation methods on an electronic digital computer, using expressions (26) and (27).

Figure 3 shows the operating characteristics of the given model, calculated for \(m=4\) and \(k=1\). The quantity \(M\) denotes the signal-to-noise ratio averaged over the signal width. In this figure, the dashed line shows the operating characteristics of a relay flaw detector, calculated for the same inspection conditions. Comparison of these operating characteristics makes it possible to quantitatively estimate the relative advantage of the investigated algorithm for processing radiometric information.

Fig. 2. Decision space for \(m=1\).

Fig. 3. Operating characteristics of the investigated flaw-detector model (solid lines) and of a relay flaw detector: \(1\) — \(M=4\), \(2\) — \(3\), \(3\) — \(2\).

The threshold \(C\) is specified by the system optimality criterion (usually the Neyman–Pearson criterion). The minimum activity of the radiation source is determined by the minimum signal amplitude that must be estimated and detected with the specified reliability.

Conclusion

Obtaining this optimal structure of the system is an essential point in constructing real flaw-detection systems, since a real system approximating the optimal one must reflect the principal features of the optimal structure. The model considered, in the indicated range of variation of the signal and noise parameters, made it possible to obtain qualitatively correct results for defect detection and for measuring its parameters. We note that this model is most valid for weak signals. When strong signals appear, the reliability of our judgments will decrease; moreover, the statistical properties of the observed process \(r(t)\) will change substantially. Analysis of the obtained algorithm for processing the process \(r(t)\) shows that small changes in the parameters of the signal and interference do not lead to significant changes in the reliability of the decisions made.

REFERENCES

1. Van Trees, H. Detection, Estimation, and Modulation Theory. Vol. 1. Moscow, “Sovetskoe radio,” 1972.
2. I. I. Privalov. Integral Equations. Moscow, ONTI, 1935.
3. V. I. Smirnov. A Course of Higher Mathematics. Moscow, Fizmatgiz, 1958.
4. H. Cramér. Mathematical Methods of Statistics. Moscow, IL, 1948.