Motion trajectory measurement of explosively-driven flyer by detonator based on multi-point velocity testing

The motion trajectory of the flyer plays a crucial role in evaluating the shock detonation performance in detonator experiments. We used a multi-point Photon Doppler Velocimetry device to measure the flyer’s motion trajectory. When the flyer’s free surface deforms, the optical axis forms an instantaneous angle with the normal velocity at the reflected spot, causing the spot to undergo a moving displacement. The average velocity of this displacement has components both along and perpendicular to the optical axis. The velocity component perpendicular to the optical axis does not affect the beat frequency of the reflected spot along the beam path. However, the optical phase difference changes with the displacement on the surface, representing the instantaneous average velocity of the normal displacement. The tilt angle can be obtained by combining the velocity vector along the beam path with the normal velocity vector. The normal displacement and tilt angle from the detection point, along with the detection azimuth of each laser probe, are used to calculate the multi-point three-dimensional coordinates of the flyer as it changes in space over time. The coordinates affected by noise are corrected using Kalman filtering, and curved surfaces are drawn to obtain the motion trajectory as it changes with time. Based on the flyer’s deformation characteristics and the limitations of the Doppler shift, we analyze and model the signal data obtained from the detonator explosion experiment. This algorithm provides theoretical support for evaluating the shock detonation performance of different micro-charged columns by studying the explosively driven flyer.

1 Introduction

The detonator, as the most representative example of the new high-security detonation class of pyrotechnics, plays a crucial role in detonation dynamics. With the evolution of detonator charges towards passivation and quantification of micro-charges, the motion trajectory of the flyer has become increasingly important for both research and practical applications [1]. However, due to its small size (typically on the millimeter scale) and high-speed, instantaneous nature, conventional measurement methods struggle to accurately capture the flyer’s dynamic parameters [2]. Current studies on flyer shapes mainly rely on numerical simulations [3], which cannot directly measure dynamic changes. Photon Doppler Velocimetry (PDV), however, offers great potential for these studies, providing accurate, non-contact measurements of velocity, displacement, and acceleration on moving surfaces. PDV is particularly well-suited for high-speed temporal analysis and visualization of velocity fields in transient conditions like blasting [4] and other high-speed measurement scenarios [5].

PDV is sensitive to the surface reflectivity characteristics and motion state of the target, and it is necessary to select an appropriate measurement method that considers the specific measurement requirements and environmental conditions. A measurable Doppler shift is generated only when the velocity of the target changes along the optical axis [6]. The measured velocity of the Doppler shift is the projection of the actual velocity onto the incoming and outgoing direction vectors of the optical path. The integral of the measured velocity will only accurately reflect the change in distance if the normal velocity at the reflected spot is exactly parallel to the optical axis [7]. In explosively-driven flyer experiments, the surface shape of the flyer constantly changes, and this change of direction is not entirely parallel to the optical axis [8], causing the reflected spot to undergo displacement on the flyer’s surface. The average velocity of this displacement includes components along and perpendicular to the optical axis, with the perpendicular component not affecting the measurable Doppler shift. However, the time-dependent phase of the Doppler signal contains the optical phase difference generated by the displacement of the reflected spot on the flyer’s surface [9]. Angular information can be obtained from the time-dependent phase of the Doppler signal according to the deformation characteristics of the flyer and the limitation that the Doppler shift only measures the motion of the portion arising from the velocity component along the beam. This angle information represents the bending degree of the flyer’s surface during flight.

A multi-point PDV device can synchronously detect the surface of the flyer. As shown in Figure 1, the three-dimensional motion trajectory of the detection point can be calculated based on the detection azimuth θ and detection radius r of the probe, the surface tilt angle β of the flyer, and the normal displacement D. The flyer’s surface is fitted based on motion trajectories detected by the multi-point. This method provides a valuable reference for evaluating the impact initiation performance of different micro-charged columns by analyzing the explosively driven flyer. It also contributes to a better understanding of the motion characteristics of the flyer during explosion dynamics.

Figure 1

Figure 1. Diagram of the probe detecting the flyer’s surface.

2 Motion trajectory measurement based on multi-point PDV device

2.1 Based on the beat frequency and phase PDV measurement principle

The measured velocity of the PDV is the projection of the actual velocity onto the incoming $\widehat{i}$ and outgoing $\widehat{o}$ direction vectors of the optical path, as shown in Equation 1.

${\mathit{v}}^{*}=\overrightarrow{\mathit{v}}·\frac{\widehat{\mathit{o}}-\widehat{\mathit{i}}}{\mathbf{2}}={\mathit{v}}_{\mathit{n}}\left(\frac{\mathit{cos}{\mathit{\theta }}_{\mathbf{2}}+\mathit{cos}{\mathit{\theta }}_{\mathbf{1}}}{\mathbf{2}}\right)+{\mathit{v}}_{\mathit{t}}\left(\frac{\mathit{sin}{\mathit{\theta }}_{\mathbf{2}}-\mathit{sin}{\mathit{\theta }}_{\mathbf{1}}}{\mathbf{2}}\right)(1)$

Where ${\theta }_1$ is the illumination angle, ${\theta }_2$ is the collection angle on both sides of the surface normal. In specular reflection (${\theta }_1={\theta }_2$), the velocity of the target along the path of the beam does not carry information about transverse motion $v_t$. The apparent velocity is $v_n\mathit{cos}\theta$.

As shown in Figure 2, when the angle between the normal velocity $v_n\left(t\right)$ of the reflected spot and the optical axis is 0°, the relationship between the normal velocity $v_n\left(t\right)$ and the Doppler shift $f_d\left(t\right)$ is shown in Equation 2.

${\mathit{f}}_{\mathit{d}}\left(\mathit{t}\right)=\mathbf{2}{\mathit{v}}_{\mathit{n}}\left(\mathit{t}\right)/{\mathit{\lambda }}_{\mathbf{0}}(2)$

Figure 2

Figure 2. The angle is 0° between the optical axis and the normal velocity of the reflected spot.

As shown in Figure 3, due to the deformation of the surface of the flyer, there is an instantaneous angle $\beta$ between the optical axis and the normal velocity of the detection point. The reflected spot undergoes a moving displacement on the surface of the flyer. The average velocity of this displacement includes components along and perpendicular to the optical axis. The displacement component $∆x$ in the direction perpendicular to the optical axis [11, 12] contributes to speckle variation but not the measurable Doppler shift at non-relativistic speeds [10]. According to Equation 3, the reflection point undergoes a moving displacement $∆s$ on the surface, and the optical phase difference $∆\phi$ of the continuous wave signal changes accordingly. Based on physics principles, the integral velocity vector aligns with the instantaneous direction of the displacement.

$∆\mathit{\phi }=\mathbf{2}\mathit{\pi }\frac{∆\mathit{s}}{{\mathit{\lambda }}_{\mathbf{0}}}=\frac{\mathbf{4}\mathit{\pi }}{{\mathit{\lambda }}_{\mathbf{0}}}{\int }_{\mathbf{0}}^{\mathit{t}}{\mathit{v}}_{\mathit{n}}\left(\mathit{t}\right)\mathit{d}\mathit{t}(3)$

Figure 3

Figure 3. The angle is β between the optical axis and the normal velocity of the reflected spot.

Therefore, in the presence of an instantaneous angle $\beta$ between the optical axis and the normal velocity of the detection point, the transverse displacement perpendicular to the optical axis does not affect the measurable Doppler shift. However, the velocity solved by the time-dependent phase is the vector sum of the velocity along and perpendicular to the optical axis. The relationship between the velocity of the reflected spot along the optical axis $v_d\left(t\right)$ and the normal velocity $v_n\left(t\right)$ is shown in Equation 4.

${\mathit{v}}_{\mathit{d}}\left(\mathit{t}\right)={\mathit{v}}_{\mathit{n}}\left(\mathit{t}\right)·\mathit{cos}\mathit{\beta }(4)$

The measurable Doppler shift of the interference signal is proportional to the velocity of the detection point along the beam path, and the phase is proportional to the displacement. The time-frequency analysis method is used to obtain the beat frequency curve over time and then solve for the velocity, which is the velocity analysis mode. The phase is solved directly by analyzing the interference signal, and the displacement is calculated and then differentiated to obtain the velocity, which is the displacement analysis mode.

In the current PDV literature, the effect of signal phase changes due to the motion displacement of reflected spots on target surfaces is seldom discussed. However, in our study, when the detection target is a dynamic surface, the displacement of the reflected spot alters the optical phase difference. The displacement solved by the optical phase difference in this paper does not refer to the forward movement of the flyer along the optical axis, but rather to the vector sum of the reflected spot’s displacement both along and perpendicular to the optical axis. The direction of the displacement vector sum follows the normal of the reflected spot. Therefore, the normal velocity in this paper explicitly refers to the instantaneous average velocity of this displacement.

Furthermore, due to the displacement component of the reflected spot on the flyer’s surface perpendicular to the optical axis, we address the limitations of measurable Doppler shifts by relating the velocity obtained from beat frequency and phase to angular information. This aspect of our work represents a significant departure from existing studies, which have not fully considered these displacement components in their models.

2.2 Principle of measuring the trajectory of the flyer

As shown in Figure 4, the multi-point 3D coordinates of the flyer in space, which change with time, can be calculated based on parameters such as detection radius, normal displacement, detection azimuth, and tilt angle using the PDV device.

Figure 4

Figure 4. Flow chart for measuring the trajectory of the flyer.

First, the detection azimuths $\theta$ and detection radius r of the four probes are designed and arranged in appropriate positions to cover the target area. Subsequently, the PDV device is used to simultaneously collect the Doppler signals of the flyer’s motion from the four probes. During the measurement process, deformation of the flyer is inevitable, causing instantaneous displacements of the normal velocity at the detection points. By performing velocity analysis and displacement analysis on the signals, the velocity along the optical axis and the normal velocity of the detection points are obtained, and the tilt angle $\beta$ is solved through association. Each signal is processed to calculate the tilt angle $\beta$ and the normal displacement D of the corresponding detection point.

Using the data collected by the four probes, a mathematical model is applied to calculate the time-varying 3D coordinates (X, Y, Z) of multiple points on the flyer surface, thereby obtaining the time-varying multi-point 3D positions of the flyer. Finally, Kalman filtering is employed to correct the coordinates affected by noise, and the corrected multi-point 3D coordinates are plotted into a surface to provide a complete representation of the flyer’s motion trajectory over time.

2.2.1 Velocity along the optical axis demodulation based on STFT

The photon Doppler signal after the short-time Fourier transform (STFT) [13], in the absence of high-frequency noise interference, exhibits a single peak spectrum in each window function spectrum. The beat frequency curve of the Doppler signal is obtained by extracting the maximum amplitude spectrum frequency from each window spectrum. Equation 2 is then used to demodulate the velocity $v_d\left(t\right)$ along the optical axis.

2.2.2 Normal velocity demodulation based on phase

The photon Doppler interference signal [6] can be written as

$\mathit{i}\left(\mathit{t}\right)=\mathit{a}\left(\mathit{t}\right)+\mathit{b}\left(\mathit{t}\right)\mathit{cos}\left[\mathit{\phi }\left(\mathit{t}\right)+{\mathit{\phi }}_{\mathbf{0}}\right](5)$

Among them, $a\left(t\right)$ is the DC voltage amplitude, $b\left(t\right)$ is the AC voltage amplitude, $\phi \left(t\right)$ is the optical phase difference, and ${\phi }_0$ is the initial phase. Equation 5 can also be expressed as Equation 6 [14]

$\mathit{i}\left(\mathit{t}\right)=\mathit{a}\left(\mathit{t}\right)+\mathit{c}\left(\mathit{t}\right)\mathit{exp}\left(\mathit{i}\mathbf{2}\mathit{\pi }{\mathit{f}}_{\mathbf{0}}\right)+{\mathit{c}}^{*}\mathit{exp}\left(-\mathit{i}\mathbf{2}\mathit{\pi }{\mathit{f}}_{\mathbf{0}}\right)(6)$

in which

$\mathit{c}\left(\mathit{t}\right)=\left(\mathbf{1}/\mathbf{2}\right)\mathit{b}\left(\mathit{t}\right)\mathit{exp}\left[\mathit{i}\mathit{\phi }\left(\mathit{t}\right)\right](7)$

The Fourier transform of $i\left(t\right)$ can be expressed as Equation 8.

$\mathit{I}\left(\mathit{f}\right)=\mathit{A}\left(\mathit{f}\right)+\mathit{C}\left({\mathit{f}}_{\mathbf{0}}-{\mathit{f}}_{\mathit{d}}\right)+{\mathit{C}}^{*}\left({\mathit{f}}_{\mathbf{0}}+{\mathit{f}}_{\mathit{d}}\right)(8)$

The spectrum is divided into three parts, where $A\left(f\right)$ is the zero-frequency component of the DC signal, $C\left(f_0-f_d\right)$ and $C^{*}\left(f_0+f_d\right)$ represent the spectrum including the signal amplitude and phase. The spectrum information of $C\left(f_0-f_d\right)$ and $C^{*}\left(f_0+f_d\right)$ is redundant. Setting $A\left(f\right)$ and $C^{*}\left(f_0+f_d\right)$ to zero, only the spectrum component of $C\left(f_0-f_d\right)$ is retained [15, 16]. Then perform the inverse Fourier transform (IDFT) on $C\left(f_0-f_d\right)$ to obtain $c\left(t\right)\mathit{exp}\left(i2\pi f_0\right)$, and subsequently perform logarithmic transformation, as shown in Equation 9. It can be observed that the phase difference $\phi \left(t\right)$ is completely separated from the representative amplitude change $b\left(t\right)$.

$\mathit{ln}\left[\mathit{c}\left(\mathit{t}\right)\mathit{exp}\left(\mathit{i}\mathbf{2}\mathit{\pi }{\mathit{f}}_{\mathbf{0}}\right)\right]=\mathit{ln}\left[\left(\mathbf{1}/\mathbf{2}\right)\mathit{b}\left(\mathit{t}\right)\right]+\mathit{i}\left[\mathit{\phi }\left(\mathit{t}\right)+{\mathit{\phi }}_{\mathbf{0}}\right](9)$

The phase distribution between -π∼π can be obtained by taking its imaginary part, as shown in Equation 10.

$\mathit{\phi }\left(\mathit{t}\right)=\mathit{I}\mathit{m}\left\{\mathit{ln}\left[\mathit{c}\left(\mathit{t}\right)\mathit{exp}\left(\mathit{i}\mathbf{2}\mathit{\pi }{\mathit{f}}_{\mathbf{0}}\right)\right]\right\}={\mathit{\phi }}_{\mathbf{0}}+\mathit{\phi }\left(\mathit{t}\right)(10)$

By unwrapping $\phi \left(t\right)$ and removing the initial phase ${\phi }_0$, the optical phase difference information $\phi \left(t\right)$ can be obtained. According to the relationship between phase and displacement, the displacement can be expressed as.

$\mathit{D}\left(\mathit{t}\right)=\frac{\mathit{\phi }\left(\mathit{t}\right)\mathit{\lambda }}{\mathbf{4}\mathit{\pi }}(11)$

Among them, $\lambda$ is the wavelength of the light source. By differentiating the displacement, the instantaneous average speed $v\left(t\right)$ of the displacement can be obtained, as shown in Equation 12.

$\mathit{v}\left(\frac{{\mathit{t}}_{\mathbf{1}}+{\mathit{t}}_{\mathbf{2}}}{\mathbf{2}}\right)=\frac{\mathit{D}\left({\mathit{t}}_{\mathbf{2}}\right)-\mathit{D}\left({\mathit{t}}_{\mathbf{2}}\right)}{{\mathit{t}}_{\mathbf{2}}-{\mathit{t}}_{\mathbf{1}}}(12)$

Since the digital differentiator is very sensitive to noise, to reduce the effect of noise in practical applications, the displacement can be processed by piecewise averaging first, followed by the differential operation. Digital low-pass differentiators are typically required to have a simple, low-order structure and enable fast differentiation [17]. Based on the analysis of minimum squared error, optimum cut-off frequency, noise amplification factors, and signal-to-noise ratio for digital differentiators in Ref. [18], we selected the Smooth-I DD digital differentiator. Smooth-I DD consists of a moving average filter of length 2M₁+1 and a simple digital differentiator of length 2M₂+ 1, where M₂ = [M/2], M₁ = M - M₂, and [·] denotes the operator that rounds the nearest integer. When the length of Smooth-I DD is 2M + 1, it is represented as shown in Equation 13.

$\mathit{y}\left(\mathit{n}\right)={\displaystyle \sum _{\mathit{k}=\mathbf{1}}^{\mathit{M}}}\frac{\left[\mathit{x}\left(\mathit{n}+\mathit{k}\right)-\mathit{x}\left(\mathit{n}-\mathit{k}\right)\right]}{\mathit{M}\left(\mathit{M}+\mathbf{1}\right)},\mathit{M}\ge \mathbf{2}(13)$

2.2.3 Motion trajectory tracking based on the Kalman filter

The detection distance and angle between the probes are shown in Figure 5. Probe one is located at the coordinate origin, probes 2, 3, and four are each at a distance r from the origin, and the angle between them is 120°.

Figure 5

Figure 5. Design of detection range and orientation.

Assume that the trajectory coordinates measured by probe 1, probe 2, probe 3, and probe four are (X₁, Y₁, Z₁), (X₂, Y₂, Z₂), (X₃, Y₃, Z₃), and (X₄, Y₄, Z₄), respectively. Then the multi-point 3D coordinate (X, Y, Z) matrix of the flyer motion trajectory is shown in Equation 14.

$\begin{array}{ll}\hfill \left[\begin{array}{ccc}{\mathit{X}}_{\mathbf{1}}& {\mathit{Y}}_{\mathbf{1}}& {\mathit{Z}}_{\mathbf{1}}\\ {\mathit{X}}_{\mathbf{2}}& {\mathit{Y}}_{\mathbf{2}}& {\mathit{Z}}_{\mathbf{2}}\\ {\mathit{X}}_{\mathbf{3}}& {\mathit{Y}}_{\mathbf{3}}& {\mathit{Z}}_{\mathbf{3}}\\ {\mathit{X}}_{\mathbf{4}}& {\mathit{Y}}_{\mathbf{4}}& {\mathit{Z}}_{\mathbf{4}}\end{array}\right]& =\left[\begin{array}{ccc}\mathbf{0}& {\mathit{D}}_{\mathbf{1}}\mathit{sin}\left({\mathit{\beta }}_{\mathbf{1}}\right)& {\mathit{D}}_{\mathbf{1}}\mathit{cos}\left({\mathit{\beta }}_{\mathbf{1}}\right)\\ \mathbf{0}& {\mathit{D}}_{\mathbf{2}}\mathit{sin}\left({\mathit{\beta }}_{\mathbf{2}}\right)& {\mathit{D}}_{\mathbf{2}}\mathit{cos}\left({\mathit{\beta }}_{\mathbf{2}}\right)\\ \mathbf{0}& {\mathit{D}}_{\mathbf{3}}\mathit{sin}\left({\mathit{\beta }}_{\mathbf{3}}\right)& {\mathit{D}}_{\mathbf{3}}\mathit{cos}\left({\mathit{\beta }}_{\mathbf{3}}\right)\\ \mathbf{0}& {\mathit{D}}_{\mathbf{4}}\mathit{sin}\left({\mathit{\beta }}_{\mathbf{4}}\right)& {\mathit{D}}_{\mathbf{4}}\mathit{cos}\left({\mathit{\beta }}_{\mathbf{4}}\right)\end{array}\right]\hfill \\ \hfill & +\left[\begin{array}{ccc}\mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathit{r}& \mathbf{0}\\ -\mathit{r}·\mathit{cos}\mathbf{30}°& -\mathit{r}·\mathit{sin}\mathbf{30}°& \mathbf{0}\\ \mathit{r}·\mathit{cos}\mathbf{30}°& -\mathit{r}·\mathit{sin}\mathbf{30}°& \mathbf{0}\end{array}\right]\hfill \end{array}(14)$

Among them, D is the normal displacement, and β is the angle between the optical axis and the normal direction of the detection point, also regarded as the tilt angle. The displacement and tilt angle determine the changes in the (Y, Z) coordinates of all detection points, and the X coordinate is completely determined by the orientation matrix P, as shown in Equation 15.

$\mathit{P}=\left[\begin{array}{ccc}\mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathit{r}& \mathbf{0}\\ -\mathit{r}·\mathit{cos}\mathbf{30}°& -\mathit{r}·\mathit{sin}\mathbf{30}°& \mathbf{0}\\ \mathit{r}·\mathit{cos}\mathbf{30}°& -\mathit{r}·\mathit{sin}\mathbf{30}°& \mathbf{0}\end{array}\right](15)$

Affected by the vibration of the detection environment, noise and signals are collected together, and the coordinates of each detection point (Y, Z) inevitably produce errors. Therefore, the motion trajectory of the flyer is corrected based on Kalman filtering. The flyer is measured in polar coordinates, with the measurement parameters being the displacement D and the tilt angle β. When the tracking filter is performed directly in the polar coordinate system, even if the target undergoes the simplest uniform linear motion, its form becomes a complex nonlinear equation. Therefore, the Cartesian coordinate system is chosen as the tracking coordinate system to simplify the state equation. After appropriately linearizing the measurement equations, existing tracking filter algorithms can be applied [19]. A commonly used linearization method for measurement equations is the coordinate change method, which converts polar coordinate measurements into equivalent Cartesian coordinate measurements and converts the additive noise in polar coordinate measurements into equivalent additive noise in Cartesian coordinate measurements.

The flyer undergoes variable acceleration motion, with time units on the order of microseconds. The velocity at adjacent moments is close in the case of high-frequency resolution. The acceleration at moment k-1 can be used to estimate the velocity at moment k in the state equation, as shown in Equation 16.

$\mathit{v}\left(\mathit{k}+\mathbf{1}\right)=\mathit{v}\left(\mathit{k}\right)+\mathit{a}\left(\mathit{k}\right){\mathit{T}}_{\mathit{s}}\approx \mathit{v}\left(\mathit{k}\right)\left[\mathbf{1}+\mathit{a}\left(\mathit{k}\right)\mathit{T}/\mathit{v}\left(\mathit{k}-\mathbf{1}\right)\right](16)$

Let the state vector $X=\left[\begin{array}{cccc}y& v_y& z& v_z\end{array}\right]$ of the observation equation, and the state transition at time k+1 is shown in Equation 17.

$\begin{array}{c}{\mathit{X}}_{\left(\mathit{k}+\mathbf{1}\right)}=\left[\begin{array}{c}{\mathit{y}}_{\left(\mathit{k}+\mathbf{1}\right)}\\ {\mathit{v}}_{\mathit{y}\left(\mathit{k}+\mathbf{1}\right)}\\ {\mathit{z}}_{\left(\mathit{k}+\mathbf{1}\right)}\\ {\mathit{v}}_{\mathit{z}\left(\mathit{k}+\mathbf{1}\right)}\end{array}\right]=\left[\begin{array}{cccc}\mathbf{1}& {\mathit{T}}_{\mathit{s}}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{1}+{\mathit{a}}_{\left(\mathit{k}\right)·}{\mathit{T}}_{\mathit{s}}/{\mathit{v}}_{\mathit{y}\left(\mathit{k}-\mathbf{1}\right)}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{1}& {\mathit{T}}_{\mathit{s}}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{1}+{\mathit{a}}_{\left(\mathit{k}\right)·}{\mathit{T}}_{\mathit{s}}/{\mathit{v}}_{\mathit{z}\left(\mathit{k}-\mathbf{1}\right)}\end{array}\right]\\ \left[\begin{array}{c}{\mathit{y}}_{\left(\mathit{k}\right)}\\ {\mathit{v}}_{\mathit{y}\left(\mathit{k}\right)}\\ {\mathit{z}}_{\left(\mathit{k}\right)}\\ {\mathit{v}}_{\mathit{z}\left(\mathit{k}\right)}\end{array}\right]+{\mathit{W}}_{\mathit{k}-\mathbf{1}}=\mathit{F}{\mathit{X}}_{\mathit{k}}+{\mathit{W}}_{\mathit{k}-\mathbf{1}}\end{array}(17)$

Among them, T_s is the interval between two adjacent moments, X represents the state, F is the state transition matrix, and $W_{k-1}$ is the process noise obeying Gaussian distribution. Assuming that the polar coordinates of the target obtained from the measurement are $Z={\left[\begin{array}{cc}D& \beta \end{array}\right]}^T$, D and β are converted to y and z-axis Cartesian coordinate values by coordinate transformation as shown in Equation 18. Then its measurement equation Z is shown in Equation 19.

$\left[\begin{array}{c}\mathit{y}\\ \mathit{z}\end{array}\right]=\left[\begin{array}{c}\mathit{D}\mathit{s}\mathit{i}\mathit{n}\mathit{\beta }\\ \mathit{D}\mathit{c}\mathit{o}\mathit{s}\mathit{\beta }\end{array}\right](18)$

$\mathit{Z}=\left[\begin{array}{c}\mathit{D}\\ \mathit{\beta }\end{array}\right]+\left[\begin{array}{c}{\mathit{w}}_{\mathit{D}}\\ {\mathit{w}}_{\mathit{\beta }}\end{array}\right]=\left[\begin{array}{c}\sqrt{{\mathit{y}}^{\mathbf{2}}+{\mathit{z}}^{\mathbf{2}}}\\ \mathit{arctan}\left(\mathit{y}/\mathit{z}\right)\end{array}\right]+\left[\begin{array}{c}{\mathit{v}}_{\mathit{D}}\\ {\mathit{v}}_{\mathit{\beta }}\end{array}\right](19)$

The measurement noises $v_D$ and $v_{\beta }$ are mutually independent smooth Gaussian noises, and the variances are ${\sigma }_D^2$ and ${\sigma }_{\beta }^2$ respectively.

The measurement matrix of the Jacobian H is shown in Equation 20.

$\mathit{H}=\frac{\mathit{\partial }\mathit{h}}{\mathit{\partial }\mathit{x}}=\left[\begin{array}{cccc}\frac{\mathit{y}}{\sqrt{{\mathit{y}}^{\mathbf{2}}+{\mathit{z}}^{\mathbf{2}}}}& \mathbf{0}& \frac{\mathit{z}}{\sqrt{{\mathit{y}}^{\mathbf{2}}+{\mathit{z}}^{\mathbf{2}}}}& \mathbf{0}\\ \frac{\mathbf{1}/\mathit{z}}{\mathbf{1}+{\left(\mathit{y}/\mathit{z}\right)}^{\mathbf{2}}}& \mathbf{0}& \frac{-\mathit{y}/{\mathit{z}}^{\mathbf{2}}}{\mathbf{1}+{\left(\mathit{y}/\mathit{z}\right)}^{\mathbf{2}}}& \mathbf{0}\end{array}\right](20)$

When the measurement error in the polar coordinate system is relatively small, the approximate representation of equivalent measurement cross-covariance matrix in the Cartesian coordinate system is obtained through first-order differentiation, as shown in Equations 21, 22.

$\left[\begin{array}{c}\mathit{d}\mathit{y}\\ \mathit{d}\mathit{z}\end{array}\right]=\left[\begin{array}{c}\mathit{d}\mathit{D}·\mathit{sin}\mathit{\beta }+\mathit{d}\mathit{\beta }·\mathit{D}\mathit{c}\mathit{o}\mathit{s}\mathit{\beta }\\ \mathit{d}\mathit{D}·\mathit{cos}\mathit{\beta }-\mathit{d}\mathit{\beta }·\mathit{D}\mathit{s}\mathit{i}\mathit{n}\mathit{\beta }\end{array}\right](21)$

${\mathit{R}}^{\prime }=\left[\begin{array}{cc}{\mathit{\sigma }}_{\mathit{y}}^{\mathbf{2}}& {\mathit{\sigma }}_{\mathit{y}\mathit{z}}\\ {\mathit{\sigma }}_{\mathit{y}\mathit{z}}& {\mathit{\sigma }}_{\mathit{z}}^{\mathbf{2}}\end{array}\right]\approx \left[\begin{array}{cc}{\mathit{\sigma }}_{\mathit{D}}^{\mathbf{2}}{\mathit{sin}}^{\mathbf{2}}\mathit{\beta }+{\mathit{\sigma }}_{\mathit{\beta }}^{\mathbf{2}}{\mathit{D}}^{\mathbf{2}}{\mathit{cos}}^{\mathbf{2}}\mathit{\beta }& \left({\mathit{\sigma }}_{\mathit{D}}^{\mathbf{2}}-{\mathit{\sigma }}_{\mathit{\beta }}^{\mathbf{2}}{\mathit{D}}^{\mathbf{2}}\right)\mathit{sin}\mathit{\beta }\mathit{cos}\mathit{\beta }\\ \left({\mathit{\sigma }}_{\mathit{D}}^{\mathbf{2}}-{\mathit{\sigma }}_{\mathit{\beta }}^{\mathbf{2}}{\mathit{D}}^{\mathbf{2}}\right)\mathit{sin}\mathit{\beta }\mathit{cos}\mathit{\beta }& {\mathit{\sigma }}_{\mathit{D}}^{\mathbf{2}}{\mathit{cos}}^{\mathbf{2}}\mathit{\beta }+{\mathit{\sigma }}_{\mathit{\beta }}^{\mathbf{2}}{\mathit{D}}^{\mathbf{2}}{\mathit{sin}}^{\mathbf{2}}\mathit{\beta }\end{array}\right](22)$

To summarize, the (Y, Z) coordinates of the four detecting points are corrected by Kalman filtering, and the orientation matrix P is added to obtain the multi-point 3D coordinates (X, Y, Z). The motion trajectory of the flyer over time is obtained by plotting a surface using multi-point 3D coordinates.

3 Experiment and data analysis

The multi-point PDV device was used to measure the motion trajectory of the flyer explosively driven by the detonator. The laser wavelength is 1550 nm, and there are four test channels with 25 mW per channel. The detector bandwidth is 10 GHz, the oscilloscope bandwidth is 12.5 GHz, and the acquisition frequency is 10¹¹Hz. During the measurement process, it is necessary to monitor the reflected power collected by the probe. Higher reflected power helps improve the signal-to-noise ratio.

Figure 6 is a schematic sectional view of the detonator detonation experiment using an explosive foil to ignite a micro-charged detonator to generate a shock wave to drive the flyer. Under the action of detonation waves and detonation products, the flyer is sheared and moves at high speed through the acceleration chamber toward the fiber optic probe. It is important to design an appropriate probe fixture to ensure that the probe remains perfectly vertical during the measurement process. Using an optical fiber collimator can help reduce interference from scattered light. In the experiment, the diameter of the optical fiber probe is 1.8 mm, the working distance is 10 mm, and the spot size is 20 μm.

Figure 6

Figure 6. Schematic of the profile of detonator detonation experiment.

Figure 7 shows the physical configuration of the probe clamp, where 1, 2, 3, and four correspond to the oscilloscope channels connected to the laser probes for later signal processing. Probe one is positioned at the coordinate origin, while probes 2, 3, and four are located 3 mm away from the origin, spaced evenly at an angle of 120° from each other.

Figure 7

Figure 7. The physical object of the probe clamp.

As shown in Figure 8, the photon Doppler signals of the explosively-driven flyer are collected by four oscilloscope channels, with the channel numbers corresponding to the probe positions indicated in Figure 7. The PDV signals can have their DC offset removed through preprocessing methods. First, the normal velocity of the photon Doppler signal is demodulated based on the phase. The signal processing steps are then illustrated in detail using the photon Doppler signal collected by oscilloscope channel 3 as an example. The spectrum of the collected Doppler signal often contains a large amount of high-frequency noise due to the vibration of the detection environment. After Fourier transforms the Doppler signal, the inverse Fourier transform is not applied to the entire spectrum in order to reduce the influence of noise. Instead, empirical mode decomposition (EMD) of the photon Doppler signal is performed first. The left column of Figure 9 shows the first five intrinsic mode functions (IMF) after the EMD, and the right column shows the corresponding spectrum of the IMF after the discrete Fourier transform (DFT). Since the beat frequencies are fully contained in the spectrum of the IMF2, only the spectrum of the IMF2 is further processed.

Figure 8

Figure 8. Multi-point PDV device acquires four sets of photon Doppler signals.

Figure 9

Figure 9. Empirical mode decomposition of photon Doppler signal.

The amplitude of the Doppler signal changes periodically with time, with the maximum amplitude spectral frequency being the beat frequency. Since displacement differentiation in phase demodulation is highly sensitive to noise, the phase information of the non-beat frequency must be removed. Any amplitude less than ${10}^{-4}$ is considered noise. The inverse discrete Fourier transform (IDFT) is applied to the signal spectrum after setting the amplitude of the noise spectrum to zero. The wrapped phase in the signal is extracted and then unwrapped. After removing the initial phase, the frequency-shifted phase is set to start from zero. The displacement is calculated based on the relationship between phase and displacement in Equation 11, and the time-displacement curve is obtained by segmented averaging.

The instantaneous average velocity, representing the normal displacement of the reflected spot, is obtained by differentiating the displacement. The velocity component along the optical axis is then derived through velocity demodulation of the beat frequency. Since the time resolutions of the two velocity demodulation algorithms differ, it is necessary to unify the time axis before calculating the tilt angle by combining the two velocity vectors. The initial velocity growth trend is influenced by the algorithmic difference. Therefore, the velocity value from the phase demodulation, which is slightly larger than that from the beat frequency demodulation, is selected as the starting point. The velocity curves from both algorithms are fitted using polynomial regression, and the time axis is unified. The residual sum of squares (RSS) for the fits is one and 0.9996, respectively, confirming the reliability of the data. Figure 10 shows the velocity comparison between the two algorithms after unifying the time axis, and the angle is then calculated by combining the velocity components. The angle increases rapidly in the initial phase and then gradually levels off.

Figure 10

Figure 10. Velocity comparison of the two algorithms.

Finally, the corresponding instantaneous displacement value is found according to the velocity of phase demodulation. In summary, the normal displacement and tilt angle under the same time coordinate axis can be obtained by processing the photon Doppler signal from one of the detection points. A total of four sets of normal displacement and tilt angle information are obtained after applying the same data processing to the remaining three channels. The motion trajectory of the (Y, Z) coordinates of each channel can be calculated. The orientation matrix is calculated based on the distance and orientation between the four probes to determine the magnitude of the X-axis coordinate. As shown in Figure 11, the motion trajectory of the multi-point 3D coordinates (X, Y, Z) is obtained. The x-axis and y-axis are primarily determined by the detection positions of the probes, as shown in Figure 7, while the z-axis represents the displacement of the flyer as it moves toward and impacts the probe. The color bars represent displacement and share the same unit as the z-axis.

Figure 11

Figure 11. Multi-point 3D coordinates.

The measurement signal is inevitably affected by noise due to the vibration of the detection environment, resulting in errors in the solution of the (Y, Z) coordinates. Therefore, Kalman filtering is used to correct each measurement signal’s (Y, Z) coordinates. Kalman filtering may make measured values smoother and more stable, especially in high-noise environments. As shown in Figure 12, the Kalman filter performs optimal estimation based on the measured values, resulting in stable estimated values that align with the growth trend of the measured values.

Figure 12

Figure 12. (Y, Z) coordinate correction based on Kalman filtering.

Figure 13 shows the motion trajectory of the corrected multi-point 3D coordinates. Compared with Figure 11, the transition of the corrected 3D coordinates on the same plane is clearer. Surface fitting can present better results.

Figure 13

Figure 13. Corrected multi-point 3D coordinates.

As shown in Figure 14, curved surfaces are plotted to determine the flight attitude of the flyer based on the motion trajectory of multi-point 3D coordinates (X, Y, Z). The duration of the flyer’s entire motion trajectory change process is 0.2μs. A higher signal acquisition frequency results in higher time resolution, meaning more signal points are collected per unit of time, thus allowing for a more detailed capture of the flyer’s motion trajectory.

Figure 14

Figure 14. Motion trajectory of flyer in space.

The motion trajectory process of an explosively-driven flyer by a detonator is as follows [3]. In the initial stage, the shock wave reaches the center of the flyer, accelerating it initially. As the flyer shears, some of the shock wave energy is consumed at the edges, causing a decrease in the driving force at the outer parts. This results in the edges lagging behind, exhibiting a hemispherical flight attitude. In the mid-stage, the shock wave impacts the side wall of the acceleration chamber and is reflected, reinforcing the driving force near the side wall. Most regions of the flyer accelerate to approximately the same velocity as the center, with only the edges having slightly lower velocities, resulting in a flattened attitude of the flyer. Finally, after undergoing multiple reflections and stabilizing, the shock wave continues to accelerate the flyer, making its flight attitude increasingly spherical.

Figure 15 illustrates the typical flight attitude of the flyer at different times, providing a clear observation of the trajectory transformation process. During the initial and final stages, the flight attitude exhibits a pronounced spherical shape under the shock wave’s impact. However, in the middle stage, the flight attitude appears relatively flat due to the reflection of the shock wave. The detonator’s impact detonation performance can be assessed by quantifying the degree of surface bending exhibited by the flyer. The flight attitude of the flyer is consistent with observations based on X-ray flash photography [2] and aligns with the shape obtained through numerical simulation of the detonator parameters [3].

Figure 15

Figure 15. Typical flight attitude of the flyer at different times.

The measurement errors of the motion trajectory can be analyzed from the following aspects. Firstly, the wavelength stability of the laser source is considered. The wavelength jitter range of the 1550 nm laser source is $\pm 0.4$ nm, resulting in a velocity error range of $\pm 0.6$ m/s. From the perspective of signal processing algorithm accuracy, the sampling frequency is 10¹¹Hz, with the number of sampling points typically set to 2048. The velocity resolution between spectral lines is 37.84 m/s. The oscilloscope bandwidth is 12.5 GHz, and the photodetector bandwidth is 10 GHz. The measurement system is capable of tracking high-speed moving targets with a maximum velocity of 7000 m/s. These errors are within an acceptable range for engineering applications.

To improve velocity measurement accuracy, lasers with narrow linewidth, stable frequencies, and stable output power should be selected, along with photodetectors that have high responsivity, short response time, and low noise-equivalent power. The oscilloscope’s bandwidth and sampling rate should be appropriately chosen based on the target’s velocity range, and signal processing algorithms should be optimized according to the characteristics of the signal. In experiments, optical fiber collimators with low insertion loss should be used to reduce the interference from scattered light. Additionally, adjusting the probe’s position to be perpendicular to the target surface or enhancing the reflectivity of the target surface can further improve the probe’s ability to collect reflected light power.

However, this motion trajectory measurement method still has certain limitations, primarily depending on the dynamic characteristics of the target surface. During the measurement process, the target surface needs to be in a dynamically changing state, with the change direction typically along the optical axis. Additionally, detection points on the dynamic surface will generate displacement components perpendicular to the optical axis. Targets driven by shock waves typically exhibit these characteristics. In the future, by combining with other advanced measurement technologies, it will be possible to compensate for the limitation of Doppler shift, which can only measure the velocity component along the beam direction, thus achieving higher precision displacement measurements and more comprehensive target trajectory analysis.

4 Conclusion

This paper utilizes a multi-point PDV device to measure the motion trajectory of an explosively driven flyer by a detonator. Initially, the detection orientation is determined based on the distance and angle between the probes, and then the photon Doppler signal of each detection point is obtained. The velocity of each detection point is calculated separately using both the classic velocity analysis mode and the displacement analysis mode. As the flyer’s surface deforms, an instantaneous angle is generated between the optical axis and the normal velocity of each detection point. The velocity obtained by the STFT represents the velocity along the optical axis, while the velocity obtained by the IDFT represents the instantaneous average velocity of the normal displacement of the detection point. The tilt angle can be derived by combining the velocity vector along the optical axis with the normal velocity vector. The multi-point 3D coordinates of the flyer, which change over time, are obtained from the normal displacement, azimuth, and tilt angle of each detection point. After correcting the motion trajectory using Kalman filtering to mitigate noise and errors, curved surfaces are plotted to visualize the motion trajectory. The reconstructed motion trajectory of the flyer accurately reflects the changes driven by the detonation, providing fundamental test data and theoretical support for evaluating the performance of micro-charged columns by analyzing the explosively driven flyer.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

YNS: Conceptualization, Methodology, Validation, Writing–original draft. SJW: Project administration, Supervision, Writing–review and editing. YKN: Writing review and editing, Data curation, Investigation. KXW: Writing review and editing, Formal Analysis, Resources.

Funding

The author(s) declare that financial support was received for the research, authorship, and/or publication of this article. This work was supported by State Key Laboratory of Transient Chemical Effects and Control, China (Grant No. WDYX23614260206).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declare that no Generative AI was used in the creation of this manuscript.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Yadav HS, Gupta NK. Flyer plate motion and its deformation during flight. Int J Impact Eng (1988) 7(1):71–83. doi:10.1016/0734-743x(88)90013-9

CrossRef Full Text | Google Scholar

2. Willey TM, Champley K, Hodgin R, Lauderbach L, Bagge-Hansen M, May C, et al. X-ray imaging and 3D reconstruction of in-flight exploding foil initiator flyers. J Appl Phys (2016) 119(23):235901. doi:10.1063/1.4953681

CrossRef Full Text | Google Scholar

3. Chen QC, Ma T, Li Y. Numerical simulation of velocity and shape of the flyer driven by HNS⁃Ⅳ explosive. Chin J Energetic Mater (2018) 26(10):814–9. doi:10.11943/CJEM2018054

CrossRef Full Text | Google Scholar

4. Pachman JKM, Nemec O, Majzlik J. A comparison of methods for detonation pressure measurement. Shock Waves (2018) 28(2):217–25. doi:10.1007/s00193-017-0761-5

CrossRef Full Text | Google Scholar

5. Pei H, Huang W, Zhang X, Zheng X. Measuring detonation wave profiles in plastic-bonded explosives using PDV. AIP Adv (2019) 9(1):015306. doi:10.1063/1.5057879

CrossRef Full Text | Google Scholar

6. Dolan DH. What does “velocity” interferometry really measure? AIP Conf Proc (2009) 1195(1):589–94. doi:10.1063/1.3295207

CrossRef Full Text | Google Scholar

7. Briggs ME, Moro EA, Shinas MA, McGrane S, Knierim D. Optical distance measurements to recover the material approach missed by optical velocimetry. J Phys Conf Ser (2014) 500:142005. doi:10.1088/1742-6596/500/14/142005

CrossRef Full Text | Google Scholar

8. Lim S, Baldovi P. Observation of the velocity variation of an explosively-driven flat flyer depending on the flyer width. Appl Sci (2018) 9(1):97. doi:10.3390/app9010097

CrossRef Full Text | Google Scholar

9. Rainey K, Gilbertson S, Kalb D, Beery T. Modulation based ranging for direct displacement measurements of a dynamic surface. Opt Express (2021) 29(14):21174–89. doi:10.1364/oe.422055

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Dolan DH. Extreme measurements with photonic Doppler velocimetry (PDV). Rev Sci Instrum (2020) 91(5):051501. doi:10.1063/5.0004363

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Chen GH, Wang DT, Liu J, Meng JH, Liu SX, Yang QG. A novel photonic Doppler velocimetry for transverse velocity measurement. Rev Sci Instrum (2013) 84(1):013101. doi:10.1063/1.4776186

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Mallick DD, Zhao M, Bosworth BT, Schuster BE, Foster MA, Ramesh KT. A simple dual-beam time-multiplexed photon Doppler velocimeter for pressure-shear plate impact experiments. Exp Mech (2019) 59:41–9. doi:10.1007/s11340-018-0435-y

CrossRef Full Text | Google Scholar

13. Dolan DH, Lemke RW, McBride RD, Martin MR, Harding E, Dalton DG, et al. Tracking an imploding cylinder with photonic Doppler velocimetry. Rev Sci Instrum (2013) 84(5):055102. doi:10.1063/1.4803074

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Takeda M, Ina H, Kobayashi S. Fourier-transform method of fringe-pattern analysis for computer-based topography and interferometry. J Opt Soc Am (1982) 72(1):156–0. doi:10.1364/josa.72.000156

CrossRef Full Text | Google Scholar

15. Liu JB, Ronney PD. Modified Fourier transform method for interferogram fringe pattern analysis. Appl Opt (1997) 36(25):6231–421. doi:10.1364/ao.36.006231

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Karaaliog ̃lu C, Skarlatos Y. Fourier transform method for measurement of thin film thickness by speckle interferometry. Opt Eng (2003) 42(6):1694–8. doi:10.1117/1.1572498

CrossRef Full Text | Google Scholar

17. Usui S, Amidror I. Digital low-pass differentiation for biological signal processing. IEEE Trans bio-medical Eng (1982) 29(10):686–93. doi:10.1109/tbme.1982.324861

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Luo J, Bai J, He P, Ying K. Axial strain calculation using a low-pass digital differentiator in ultrasound elastography. IEEE Trans Ultrason Ferroelectrics, Frequency Control (2004) 51(9):1119–27. doi:10.1109/tuffc.2004.1334844

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Ding LF, Geng FL. Radar principles. 5th ed. Beijing: Xidain University Press (2014). p. 185.

Google Scholar