常见的Abel逆变换求解方法及其python实现

其中AJI我是写了一个函数，就按照式（7）来写就可以，这里我就不贴上来了。这个离散法很简单，上面的代码我直接写的，没测试，不保证完全正确哈。上一篇论文中总结了一些常见的Abel逆变换求解方法。，应粉丝要求，这篇文章写一下对应的算法实现。详细推导见上一篇文章，结果是。其余就是很简单的离散求和。

Fourier猴

1207人浏览 · 2023-07-13 20:27:52

Fourier猴 · 2023-07-13 20:27:52 发布

上一篇论文中总结了一些常见的Abel逆变换求解方法Abel逆变换及其求解方法，应粉丝要求，这篇文章写一下对应的算法实现。

先简单回顾一下，

1. Abel变换（Abel transform）是一种积分变换，常用于球对称或轴对称函数的分析，其数学表达式为：

$F(y)=2\int_y ^\infty\frac{f(r)rdr}{\sqrt{r^2-y^2}}$ （1）

2.Abel逆变换（ inverse Abel transform），其数学表达式为：

$f(r)=-\frac{1}{\pi}\int_r ^\infty \frac{dF}{dy} \frac{dy}{\sqrt{y^2-r^2}}$ 。（2）

求解方法

1.直接离散法

$f(r_j)=- \frac{1}{\pi} \sum_{i=j}^{N-1}\frac{F(y_{i+1})-F(y_i)}{\sqrt{(y_i+\frac{\Delta y}{2})^2-r_j^2}}$ 。（3）

代码实现：

def IAbel_DD(Fy, N, DeltaY, DeltaR):
    r"""
    Inverse Abel Transform: direct-discretization，DD
    Args:
        | Fy: The integrated intensities Fy
        | N: N+1 is the number of discrete distribution datapoint of Fy
        | DeltaY: spatial separation between data
        | DeltaR: spatial separation between data, DeltaY = DeltaR
    Returns:
        | frDD: fr at different radius, solved by DD method
    """

    frDD = np.zeros(N + 1)
    r = np.zeros(N + 1)
    for j in range(0, N + 1):
        r[j] = DeltaR * j
        frDD[j] = 0
        y = np.zeros(N)
        for i in range(j, N):
            y[i] = DeltaY * i
            Deno = np.sqrt(y[i] + DeltaY/2)**2 - r[j]**2)
            Nume = Fy[i+1] - Fy[i]
            Quo = -Nume/Deno/math.pi
            frDD[j] = frDD[j] + Quo
    return frDD

这个离散法很简单，上面的代码我直接写的，没测试，不保证完全正确哈。

2.Hankel-Fourier变换法：

详细推导见上一篇文章，结果是

$f(r_j)=\frac{1}{2\pi[(2N+1)\Delta x]^2}\sum_{i=-N}^{N}F(y_i)\sum_{k=0}^{N}cos[(\frac{i}{2N+1})k]kJ_0( \frac{jk}{2N+1})$ （4）

其中，J0是0阶Bessel function，可以直接用：

J0 = jv(0, j * k / (2 * N + 1)

其余就是很简单的离散求和。

3.Nestor-Olsen法：

该算法由Nestor和Olsen提出，f(r)与F(r)之间得关系由下式给出：

$f(r_j)=-\frac{2}{\pi \Delta y}\sum_{i=j}^{N-1}F(y_i)B_{j,i}$ (5)

其中，Bj,i满足：

$\left\{\begin{matrix} B_{j,i}=A_{j,i-1}-A_{j,i},i\geqslant j+1 \\ B_{j,i}=-A_{j,i},i=j \end{matrix}\right.$ （6）

其中， $A_{j,i}=\frac{[i^2-(j-1)^2]^{\frac{1}{2}}-[(i-1)^2-(j-1)^2]^{\frac{1}{2}}}{2i-1}$ 。（7）

代码实现：

def IAbel_NO(Fy, N, DeltaY):
    """
    Inverse Abel Transform: Nestor–Olsen method (NO method)
    :param Iy: The integrated intensities Fy
    :param N: N + 1 is the number of discrete distribution datapoint of Iy
    :param DeltaY: spatial separation between data

    :return: frNO: The fr at different radius, by NO method
    """
    Ce1 = -2 / math.pi / DeltaY
    frNO = np.zeros(N + 1)
    for j in range(0, N+1):
        sum = 0
        for i in range(j, N):
            if i == j:
                Aji = AJI(i, j)
                Bji = -Aji
            else:
                Aji = AJI(i, j)
                Aji1 = AJI(i - 1, j)
                Bji = Aji1 - Aji
            sum = sum + Iy[i] * Bji
        frNO[j] = Ce1 * sum
    return frNO