Electromagnetic Property Sensing and Channel Reconstruction Based on Diffusion Schrödinger Bridge in ISAC

In the channel estimation stage, the BS adopts fully digital precoding structure where the number of radio frequency (RF) chains $N_{RF}$ is equal to the number of transmitting antennas $N_{t}$. The central frequency of the signals is denoted by $f_{c}$ with the corresponding wavelength $\lambda_{c}$. The number of subcarriers is denoted by $K$ and the frequency spacing between adjacent subcarriers is denoted by $\Delta_{f}$. The number of the transmitted symbols in each subcarrier is denoted by $I$. We assume a quasi-static environment where the channels remain unchanged throughout the sensing period.

Since only the signals scattered by the target carry the information of its EM property, we may send the pilot signals towards the target by beamforming properly at the transmitter. Let the subscript $k$ represent that the physical quantity is associated with the $k$-th subcarrier. Denote $\mathbf{H}_{k}\in\mathbb{C}^{N_{r}\times N_{t}}$ as the overall echo channel from the transmitter to the receiver. Thus, the received signals can be formulated as

where $\mathbf{w}_{k}\in\mathbb{C}^{N_{t}\times 1}$ is the pilot symbol on the $k$-th subcarrier; $\mathbf{n}_{k}\sim\mathcal{CN}\left(\mathbf{0},\sigma_{k}^{2}\mathbf{I}_{N_{r}}\right)$ is the CSCG noise at the receiver of the $k$-th subcarrier. Denote $\mathbf{W}_{k}=\left[\mathbf{w}_{k,1},\mathbf{w}_{k,2},\cdots,\mathbf{w}_{k,I}\right]\in\mathbb{C}^{N_{t}\times I}$ as the pilot matrix stacked by time, and denote $\mathbf{N}_{k}=\left[\mathbf{n}_{k,1},\mathbf{n}_{k,2},\cdots,\mathbf{n}_{k,I}\right]\in\mathbb{C}^{N_{r}\times I}$. Then the overall received pilot signals $\mathbf{Y}_{k}\in\mathbb{C}^{N_{r}\times I}$ can be formulated in a compact form

In order to extract the EM property of the target, the BS first needs to estimate $\mathbf{H}_{k}$ by applying the least square (LS) method as

where $\hat{\mathbf{H}}_{k}$ is the minimum variance unbiased (MVU) estimated sensing channel. We then stack $\hat{\mathbf{H}}_{k}$ for $k\in\{1,\cdots,K\}$ into a 3rd order tensor $\hat{\mathcal{H}}\in\mathbb{C}^{K\times N_{r}\times N_{t}}$ as

According to [11, 12, 13], the EM property of the target is implicitly encoded in the received echo signals that are transmitted through the sensing channel. Thus, we can leverage $\hat{\mathcal{H}}$ as the prior information to reconstruct the EM property of the target.

Suppose that the EM property of the target exhibits isotropy [25, 26, 27]. The process of sensing the EM characteristic is essentially about determining the contrast function $\chi_{k}(\mathbf{r})$, which represents the discrepancy in the complex relative permittivity of the target as compared to the surrounding air. Considering the relative permittivity and conductivity of air to be nearly $1$ and $0$ Siemens per meter (S/m) respectively, the contrast function can be defined as [28, 29]

where $\epsilon_{r}(\mathbf{r})$ denotes the real relative permittivity at point $\mathbf{r}$, $\sigma(\mathbf{r})$ denotes the conductivity at point $\mathbf{r}$, $\omega_{k}=2\pi f_{k}$ denotes the angular frequency of the EM waves, and $\epsilon_{0}$ is the vacuum permittivity.

Throughout this document, it is presumed that the electric fields exhibit a harmonic time variation characterized by $e^{-j\omega_{k}t}$ [30]. Let $\lambda_{k}=c/f_{k}$ represent the wavelength and $k_{k}=2\pi/\lambda_{k}$ represent the wave number within the ambient medium. The total electric field and the incident electric field, which propagate through the medium in the $x$, $y$, and $z$ dimensions, are represented by the complex vectors $\mathbf{E}_{k}^{t}({\mathbf{r}})\in\mathbb{C}^{3\times 1}$ and $\mathbf{E}_{k}^{i}({\mathbf{r}})\in\mathbb{C}^{3\times 1}$, respectively. Since the incident electric field is linearly induced by the currents on the transmitting antennas, there is a matrix $\tilde{\mathbf{H}}_{1,k}({\mathbf{r}})\in\mathbb{C}^{3\times N_{t}}$ that linearly maps $\mathbf{w}_{k}$ to $\mathbf{E}_{k}^{i}({\mathbf{r}})$, i.e.,

Upon exposure to the incident field $\mathbf{E}_{k}^{i}({\mathbf{r}})$, the fields $\mathbf{E}_{k}^{i}({\mathbf{r}})$ and $\mathbf{E}_{k}^{t}({\mathbf{r}})$ are governed by the homogeneous wave equation for the incident field and the inhomogeneous wave equation for the total field, respectively [31], i.e.,

Suppose that the BS knows the target is positioned in the region $D$ through prior localization. To address the solutions for (7) and (8), the total electric field within $D$ can be formulated by the 3D Lippmann-Schwinger equation [31, 32, 33]

where $\overline{\overline{\mathbf{G}}}_{k}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\in\mathbb{C}^{3\times 3}$ is the dyadic electric field Green’s function that satisfies

Meanwhile, $\overline{\overline{\mathbf{G}}}_{k}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$ can be formulated as [34]

where $R^{\prime}$ is the distance defined as $R^{\prime}\triangleq\|\mathbf{r}-\mathbf{r}^{\prime}\|_{2}$, $\hat{\mathbf{r}}\in\mathbb{R}^{3\times 1}$ is the unit vector from $\mathbf{r}^{\prime}$ to $\mathbf{r}$, and $g_{k}\left(\mathbf{r},\mathbf{r}^{\prime}\right)$ is the scalar Green’s function defined as $g_{k}\left(\mathbf{r},\mathbf{r}^{\prime}\right)\triangleq\frac{\exp(jk_{k}R^{\prime})}{4\pi R^{\prime}}$ [34].

The echo electric field at the BS’s receiver scattered back from the target can then be formulated as [32, 33]

where $\mathbf{r}_{n}$ denotes the position of the $n$-th receiving antenna. Suppose the receiver can only measure the scalar electric field component in the direction represented by the unit vector $\mathbf{q}\in\mathbb{R}^{3\times 1}$. The received echo signals can also be formulated as

where $\tilde{G}_{r}$ denotes the receiving antenna gain.

We utilize the point cloud representation to concisely and vividly represent the distribution of the target’s EM property. Define the $m$-th normalized 5D point $\mathbf{x}_{m}\in\mathbb{R}^{5\times 1}$ that comprises both the 3D location information and the 2D EM property as

where $x_{m}$, $y_{m}$, and $z_{m}$ denote the coordinates of the $m$-th point in each dimension; $x_{c}$, $y_{c}$, and $z_{c}$ denote the coordinates of the center of the target; $x_{d}$, $y_{d}$, and $z_{d}$ denote the corresponding standard deviations. Here, $\epsilon_{m}$ and $\sigma_{m}$ represent the dielectric constant and conductivity at the $m$-th point, respectively; $\epsilon_{c}$ and $\sigma_{c}$ represent their respective central values; $\epsilon_{d}$ and $\sigma_{d}$ denote their respective standard deviations. Suppose a total of $M$ points $\mathbf{x}_{m}$ in (14) constitute the point cloud that represents the target and is defined as $\mathbf{X}_{\text{data}}\in\mathbb{R}^{M\times 5}$.

The implementation of the point cloud approach presents a more efficient and uncomplicated alternative for discerning the EM property. Point clouds intrinsically facilitate the distinction between the background medium and the target, thereby eliminating the necessity to analyze the known background medium and significantly reducing the computational burden. Additionally, the representation of data through point clouds permits a clear and prompt visualization of the 3D target, which thereby enhances the intuitive interpretation of the inversion outcomes.

Let $p_{\text{data}}$ denote the distribution of the EM property and $p_{\text{prior}}$ denote the distribution of the latent extracted from the sensing channel. The latent shares the same dimensions as $\mathbf{X}_{\text{data}}$ and is denoted as $\mathbf{X}_{\text{prior}}\in\mathbb{R}^{M\times 5}$. To estimate the EM property from the sensing channel or to reconstruct the sensing channel from the EM property, we need to link $p_{\text{data}}$ and $p_{\text{prior}}$ using DSB. Specifically, DSB establishes the link through a bidirectional process: the forward process gradually transforms $p_{\text{data}}$ into $p_{\text{prior}}$, while the reverse process maps $p_{\text{prior}}$ back to $p_{\text{data}}$. Both processes can be represented via Markov chains.

Denote $\mathbf{X}_{0:N}$ as the set of $\{\mathbf{X}_{0},\cdots,\mathbf{X}_{N}\}$ that are sequentially generated from $\mathbf{X}_{0}=\mathbf{X}_{\text{data}}$ to $\mathbf{X}_{N}=\mathbf{X}_{\text{prior}}$ over a sequence of $N-1$ intermediate states. The forward transition $p_{i+1\mid i}\left(\mathbf{X}_{i+1}\mid\mathbf{X}_{i}\right)$ is constructed to progressively guide the distribution from $p_{0}=p_{\text{data}}$ to approximate $p_{N}=p_{\text{prior}}$. The joint probability density of $\mathbf{X}_{0:N}$ is

Similarly, the reverse process can be formulated as a Markovian sequence, with the reverse joint density being

where the conditional probability $p_{i\mid i+1}\left(\mathbf{X}_{i}\mid\mathbf{X}_{i+1}\right)$ represents the probability of transitioning from state $\mathbf{X}_{i+1}$ at time $i+1$ to state $\mathbf{X}_{i}$ at time $i$. We can decompose the conditional probability $p_{i\mid i+1}\left(\mathbf{X}_{i}\mid\mathbf{X}_{i+1}\right)$ using Bayesian theorem as

where $p_{i+1\mid i}\left(\mathbf{X}_{i+1}\mid\mathbf{X}_{i}\right)$ is the conditional probability of the forward process, $p_{i}\left(\mathbf{X}_{i}\right)$ is the marginal distribution of state $\mathbf{X}_{i}$, and $p_{i+1}\left(\mathbf{X}_{i+1}\right)$ is the marginal distribution of state $\mathbf{X}_{i+1}$.

However, directly computing the conditional probability $p_{i\mid i+1}\left(\mathbf{X}_{i}\mid\mathbf{X}_{i+1}\right)$ using (17) is generally quite challenging due to the complexity of the involved distributions and the recursive nature of the computation. In practice, score-based generative models (SGMs) adopt a more tractable approach to handle the forward process, and thus could simplify the forward process by modeling it as the gradual addition of Gaussian noise to the states over time. The forward process can be represented as

where $\gamma_{i+1}$ is the noise level parameter at time step $i+1$, and $f_{i}\left(\mathbf{X}_{i}\right)$ is the forward drift term that governs the deterministic part of the state evolution.

For a sufficiently large number of time steps $N+1$, the distribution of the state at the final time step will converge to $p_{\text{prior}}$, which serves as the starting point for the reverse-time generative process.

The forward process (15) and the reverse process (16) can also be described in a continuous-time framework using stochastic differential equations (SDEs). Specifically, the forward process can be modeled by an SDE as

where $t\in[0,T]$, $f_{t}(\mathbf{X}_{t}):\mathbb{R}^{M\times 5}\rightarrow\mathbb{R}^{M\times 5}$ is the drift term function, $g_{t}$ represents the diffusion coefficient, and $\mathbf{B}_{t}$ denotes the standard Brownian motion. The reverse process, on the other hand, involves solving the time-reversed version of the SDE in (19), and is given by

The DSB is an extension of the Schrödinger bridge (SB) problem, which incorporates diffusion model to model uncertainty and variability in dynamic systems. In the context of SB problem, we aim to find an optimal distribution $p^{*}\in\mathscr{P}_{N+1}$ that minimizes the KL divergence from a reference path measure $p^{\text{ref}}\in\mathscr{P}_{N+1}$. The optimization problem is defined as

where the marginal distributions $p_{0}$ and $p_{N}$ correspond to the distributions at the start and end points of the process, respectively. Typically, the reference measure path $p^{\text{ref}}$ is generated using the same form of the forward SDE as in (19).

After determining the optimal solution $p^{*}$, we can sample $\mathbf{X}_{0}\sim p_{\text{data}}$ by first generating $\mathbf{X}_{N}\sim p_{\text{prior}}$ and then iteratively applying the backward transition $p_{i\mid i+1}\left(\mathbf{X}_{i}\mid\mathbf{X}_{i+1}\right)$. Alternatively, we can sample $\mathbf{X}_{N}\sim p_{\text{prior}}$ by initially drawing $\mathbf{X}_{0}\sim p_{\text{data}}$ and subsequently applying the forward transition $p_{i+1\mid i}\left(\mathbf{X}_{i+1}\mid\mathbf{X}_{i}\right)$. The SB formulation enables the bidirectional transitions without relying on closed-form expression of $p_{\text{data}}$ and $p_{\text{prior}}$.

Although the SB problem does not have a closed-form solution, it can be tackled using iterative proportional fitting (IPF), which iteratively solves the following optimization problems:

where the superscript $n$ denotes the number of iterations, and the initialization $p^{0}$ is set as $p^{\text{ref}}$. However, implementing IPF in real-world scenarios can be computationally prohibitive, as it requires the calculation and optimization of joint densities.

DSB can be viewed as an approximate method for IPF, which simplifies the optimization of the joint density by decomposing it into a sequence of conditional density optimization tasks. Specifically, the distribution $p$ is split into the forward and backward conditional distributions $p_{i+1\mid i}$ and $p_{i\mid i+1}$, respectively:

It can be shown that optimizing conditional distributions $p_{i+1\mid i}$ and $p_{i\mid i+1}$ in (24) and (25) leads to the optimization of the joint distribution $p$ in (22) and (23) [21]. We assume that the conditional distributions $p_{i+1\mid i}$ and $p_{i\mid i+1}$ are Gaussian distributions, following the assumption commonly used in SGMs and allowing DSB to analytically handle the time reversal process. Consequently, we employ two separate neural networks to model the forward and backward dynamics.

The forward process, which governs the transition of the state from one time step to the next, is mathematically expressed as

where $p_{i+1\mid i}^{n}(\mathbf{X}_{i+1}\mid\mathbf{X}_{i})$ denotes the conditional probability distribution of the state $\mathbf{X}_{i+1}$ given the state $\mathbf{X}_{i}$ at the previous time step. In the multivariate Gaussian distribution $\mathcal{N}\left(\mathbf{X}_{i}+\gamma_{i+1}f_{i}^{n}(\mathbf{X}_{i}),2\gamma_{i+1}\mathbf{I}\right)$, $\gamma_{i+1}$ represents the diffusion coefficient that controls the spread of the distribution, and $f_{i}^{n}(\mathbf{X}_{i})$ represents the drift term that accounts for the deterministic part of the state transition in the forward process. For brevity, we define the forward estimation $F_{i}^{n}(\mathbf{X}_{i})$ as

Conversely, the backward process, which describes the reverse-time evolution of (26), is given by

where $p_{i\mid i+1}^{n}(\mathbf{X}_{i}\mid\mathbf{X}_{i+1})$ represents the conditional probability distribution of the state $\mathbf{X}_{i}$ given the state $\mathbf{X}_{i+1}$ at the subsequent time step. Similar to the forward process, the backward process is modeled as a multivariate Gaussian distribution $\mathcal{N}\left(\mathbf{X}_{i+1}+\gamma_{i+1}b_{i+1}^{n}(\mathbf{X}_{i+1}),2\gamma_{i+1}\mathbf{I}\right)$. The term $b_{i+1}^{n}(\mathbf{X}_{i+1})$ in (28) is the drift term specific to the backward process, which plays an analogous role to $f_{i}^{n}(\mathbf{X}_{i})$ in the forward process and can be computed according to (26) as

For brevity, we define the backward estimation $B_{i}^{n}(\mathbf{X}_{i})$ as

According to (26) and (27), we can compute the probability density function $p_{i+1}^{n}(\mathbf{X}_{i+1})$ at time step $i+1$ as

where the integral normalization factor $(4\pi\gamma_{i})^{-\frac{5M}{2}}$ comes from $\mathbf{X}_{i},\mathbf{X}_{i+1}\in\mathbb{R}^{M\times 5}$.

Taking the logarithm of $p_{i+1}^{n}(\mathbf{X}_{i+1})$ and applying the gradient with respect to $\mathbf{X}_{i+1}$, we obtain

where we employ the Bayesian rule and then use (26) as

Substituting (32) into (29) and utilizing (27), we obtain

Now, considering the Bayesian update step and substituting (34) into (30), we can derive the backward estimation $B_{i+1}^{n}(\mathbf{X}_{i+1})$ as

where the expectation is taken over the joint distribution $p_{i,i+1}^{n}(\mathbf{X}_{i},\mathbf{X}_{i+1})$.

To minimize the difference between $B_{i+1}^{n}(\mathbf{X}_{i+1})$ and the variable on the right-hand-side of (35), we define the loss function $\mathcal{L}_{B_{i+1}^{n}}$ as

Similarly, the loss function $\mathcal{L}_{F_{i}^{n+1}}$ for the mapping function $F_{i}^{n+1}$ can be given by

The loss functions (36) and (37) are derived under the assumption that the distribution $p_{i,i+1}^{n}$ accurately models the underlying dynamics of the system and that the error introduced by approximating $B_{i+1}^{n}$ and $F_{i}^{n+1}$ is minimized in the Frobenius norm sense. The schematic diagram of DSB pipeline with forward and backward processes is shown in Fig. 2, where the latent refers to the compressed features extracted from the sensing channel.

In practical applications, the DSB methodology employs two neural networks to approximate the forward estimation (27) and the backward estimation (30). Let $\alpha^{n}$ and $\beta^{n}$ represent the trainable parameters of two neural networks $F_{\alpha^{n}}(i,\mathbf{X}_{i})$ and $B_{\beta^{n}}(i,\mathbf{X}_{i})$ in the $n$-th iteration. Specifically, the neural network $B_{\beta^{n}}(i,\mathbf{X}_{i})$ is used to approximate the backward estimation $B_{i}^{n}(\mathbf{X}_{i})$, and the neural network $F_{\alpha^{n}}(i,\mathbf{X}_{i})$ is used to approximate the forward estimation $F_{i}^{n}(\mathbf{X}_{i})$. The iterative optimization of $F_{\alpha^{n}}(i,\mathbf{X}_{i})$ and $B_{\beta^{n}}(i,\mathbf{X}_{i})$ is crucial for the DSB methodology.

The DSB methodology proceeds by alternately training the backward network $B_{\beta^{n}}(i,\mathbf{X}_{i})$ and the forward network $F_{\alpha^{n}}(i,\mathbf{X}_{i})$ across multiple iterations indexed by $n$. Specifically, during the $(2n+1)$-th epoch, we train the backward network $B_{\beta^{n}}(i,\mathbf{X}_{i})$, while during the $(2n+2)$-th epoch, we train the forward network $F_{\alpha^{n}}(i,\mathbf{X}_{i})$. The alternating optimization of $B_{\beta^{n}}(i,\mathbf{X}_{i})$ and $F_{\alpha^{n}}(i,\mathbf{X}_{i})$ is designed to minimize (24) and (25), which ensures that the forward and backward processes are accurately modeled. Moreover, both $B_{\beta^{n}}(i,\mathbf{X}_{i})$ and $F_{\alpha^{n}}(i,\mathbf{X}_{i})$ are composed of the same concatsquash layers as in [13].

The convergence of IPF has been rigorously proven in the literature, as demonstrated in [21]. The proof underpins the theoretical soundness of the DSB approach, which confirms that the method will, under appropriate conditions, converge to a solution that satisfies the desired characteristics of DSB.