【论文精读系列】之《Turbulence-Induced 2D Correlated Image Distortion》其二

最新推荐文章于 2024-09-11 15:45:55 发布

西园啜茗

最新推荐文章于 2024-09-11 15:45:55 发布

阅读量447

点赞数 3

分类专栏： CV 文章标签：计算机视觉图形学

本文链接：https://blog.csdn.net/qq_39800591/article/details/108151231

版权

CV 专栏收录该内容

7 篇文章 2 订阅

订阅专栏

【论文精读系列】之《Turbulence-Induced 2D Correlated Image Distortion》其二

3 Background（背景）
- 3.1 Turbulence and Variance（湍流和方差）
- 3.2 Pair-wise Distortion Correlation（成对失真的相关性）
4 Designing the Field Autocorrelation（场自相关设计）

论文地址：《Turbulence-Induced 2D Correlated Image Distortion》

3 Background（背景）

3.1 Turbulence and Variance（湍流和方差）

3D湍流的特点是在一定范围内具有独立的流体涡流。一个涡流在最大（外部）空间尺度 $L_0$ 处分解为较小的涡流，而较小的涡流类似上述过程进一步减小。这个过程一直持续到最小的涡流达到湍流的内尺度 $l_0$ 。惯性范围由 $L_0,l_0]$ 定义。在那里，动能被转移到尺度越来越小的漩涡中。在 $l_0$ 以下，流体运动中的能量完全消散为热量。湍流的外部尺度通常被建模为 $L_0=0.4h$ ，其中 $h$ 是地面以上的高度。
湍流伴随着压力和温度的波动，这会影响介质的折射率。在统计上均匀且局部各向同性的湍流中，折射波动的特征用折射率结构常数 $C_n^2$ 刻画，它表示湍流强度，且以 $m^{-2/3}$ 为单位进行测量。 $C_n^2$ 的值越大表示湍流越强烈，而 $C_n^2=0$ 表示没有湍流的介质。在空气中， $C_n^2$ 的范围一般在 $10^{-17}-10^{-13}m^{-2/3}$ ，夏季的白天这个值一般在 $8\times10^{-13}m^{-2/3}$ 的高值。由于湍流强度取决于环境，因此通常 $C_n^2$ 会在成像路径上变化，特别易受低海拔地区的高压和高海拔地区的强风的影响，这两者都会影响天文望远镜的放置和性能。在本文中，我们将重点放在水平视图上，这是远距离观测等应用中的典型视图。
焦距为 $f$ 的相机在距离 $L$ 处观察到一个物点 $o$ 。在没有湍流的情况下，视线由单位矢量 ${\boldsymbol o}$ 表示。如图2所示，由于折射率波动，湍流导致从物点发出的光线发生随机折射。因此，此光线以随机到达角（ $A A$ ）到达相机，该角度是围绕 ${\boldsymbol o}$ 方向的扰动。 $A A$ 扰动会扭曲 $o$ 在像平面的投影。

（图2：在没有湍流的情况下，到达目标点 $o$ 的视线是 ${\boldsymbol o}$ ，然后将该点投影到 ${\boldsymbol p}$ 上。这里 $L$ 是路径长度， $f$ 是焦距， $D$ 是透镜孔径直径。湍流介质中的随机折射产生一个到达角（ $A A$ ）扰动，导致 ${\boldsymbol p}+{\boldsymbol e}({\boldsymbol p})$ 的畸变投影）

球形光的波前的 $A A$ 方差为： $\sigma_{AA}^2=\gamma(q)C_n^2LD^{-1/3} \tag{1}$ 其中 $D$ 是相机镜头的直径， $q$ 是菲涅尔数： $q=\frac{D}{\sqrt{\lambda L}} \tag{2}$ 其中 $\lambda$ 是光线的波长。球面波的系数 $\gamma(q)$ 为： $\gamma(q)=\frac{\sqrt{3}}{16}\Gamma(\frac{1}{6})\Gamma(\frac{8}{3})(\frac{2}{\beta})^{\frac{1}{3}}[1+\frac{6}{17}\frac{8}{3}\cdot Re\{(-\frac{\pi(\beta q)^2}{2j})^\frac{1}{6}[_2F_1(\frac{1}{6},\frac{17}{6};\frac{23}{6};1+\frac{\pi(\beta q)^2}{2j})]\}] \tag{3}$ 其中 $\beta=0.5216$ 且： $\Gamma(\mu)=\int_0^{\infty}\tau^{\mu-1}e^{-\tau}d\tau \tag{4}$ 是Gamma函数。在公式（3）中， $_2F_1()$ 是超几何函数： $_2F_1(a,d;c;\zeta)=\sum_{m=0}^{\infty}\frac{(a)_m(d)_m}{(c)_m}\frac{\zeta^m}{m!} \tag{5}$ 对于 $q\gg 1$ ，公式（3）可以简化为： $\gamma(q)\approx1.09 \tag{6}$

3.2 Pair-wise Distortion Correlation（成对失真的相关性）

（图3：在没有湍流的情况下，到达两目标点 $o$ 和 $o^{'}$ 的视线分别是 ${\boldsymbol o}$ 和 ${\boldsymbol o'}$ ，然后分别投影到 ${\boldsymbol p}$ 和 ${\boldsymbol p}+{\boldsymbol v}$ ）

现在，考虑两个物点 $o$ 和 $o^{'}$ 。在没有湍流的情况下，它们对应的视线由来自相机投影中心的单位向量 ${\boldsymbol o}$ 和 ${\boldsymbol o'}$ 定义，如图3所示。这两条线形成一个平面 ${\boldsymbol \Pi}={\boldsymbol o}\times{\boldsymbol o'}$ 。两线之间的夹角为： $\theta=arccos({\boldsymbol o}\cdot{\boldsymbol o'}) \tag{7}$ 湍流引起的随机折射扰动，会导致这些视线中的每条视线相应地偏离到达角 $A A (o)$ 和 $A A (o^{'})$ 。 $A A (o)$ 和 $A A (o^{'})$ 之间存在相关性，且相关函数为： $b(\theta)=\frac{Cov[AA(o),AA(o')]}{\sigma_{AA}^2} \tag{8}$ 其中Cov是协方差运算。该相关性通常随着 $\theta$ 的增大而减小。此外，相关性取决于 $A A$ 相对于 ${\boldsymbol \Pi}$ 的方向。假设一个 $A A$ 扰动在 ${\boldsymbol \Pi}$ 平面内，该平面包括 ${\boldsymbol o}$ 和 ${\boldsymbol o'}$ ，那么这个扰动被称为平行的，用 $\parallel$ 表示。如果 $A A$ 扰动垂直于该平面，则表示为 $\bot$ 。一般来说，在 $o$ 处的 $A A$ 扰动既有平行分量，也有垂直分量。
注意，对于固定的 $o$ ， $\parallel$ 和 $\bot$ 轴和分量是根据另一个点 $o^{'}$ 来定义的。由于物点 $o$ 在不同的相对方向上与大量其它 $W$ 个点共存，因此对于任何特定 $o$ 都有 $\parallel$ 和 $\bot$ 的 $\mathcal{O}(W)$ 定义。因此，我们处理一个非平凡的领域。
现在，让我们返回到分析两个特定的对象点 $o$ 和 $o^{'}$ 。使用以下无量纲变量。从物体的角度来看，光圈的角度大小为： $\theta_D=\frac{D}{L} \tag{9}$ 物体间的视角可以标准化为： $\eta=\frac{\theta}{\theta_D} \tag{10}$ 湍流的外部尺度可以标准化为： $\rho=\frac{L_0}{D} \tag{11}$ 对于无量纲参数 $0\leq\xi\leq 1$ ，定义函数： $Q(\xi)=(1-\xi)^{5/3},M_{\eta}=2\eta\xi(1-\xi) \tag{12}$ 这两个分量的相关函数是： $b_{\parallel,\bot}(\eta)=\frac{\int_0^1\{\int_0^{\infty}H(\kappa)[J_0(M_{\eta}\kappa)\mp J_2(M_{\eta}\kappa)]d\kappa\}C_n^2(L\xi)Q(\xi)d\xi}{\int_0^1\{\int_0^{\infty}H(\kappa)d\kappa\}C_n^2(L\xi)Q(\xi)d\xi} \tag{13}$ 其中 $\kappa$ 是一个包含在非负实数中的无量纲积分变量， $J_0$ 和 $J_2$ 是第一类贝塞尔函数，且： $H(\kappa)=\frac{[(2\rho\kappa)^2+1]^{-11/6}J_2^2(\kappa)}{\kappa} \tag{14}$ 公式（13）中的负号对应 $b_{\parallel}(\eta)$ 而正号对应 $b_{\bot}(\eta)$ 。注意公式（13）的分母只是在 $\eta=0$ 时计算的分子，因此 $b_{\parallel,\bot}(0)=1$ 。
在 $C_n^2$ 为常数的水平路径中，公式（13）简化为： $b_{\parallel,\bot}(\eta)=\frac{\int_0^1\{\int_0^{\infty}H(\kappa)[J_0(M_{\eta}\kappa)\mp J_2(M_{\eta}\kappa)]d\kappa\}Q(\xi)d\xi}{\int_0^1\{\int_0^{\infty}H(\kappa)d\kappa\}Q(\xi)d\xi} \tag{15}$ 相关函数（15）不再依赖于 $C_n^2$ 或范围 $L$ ，而是依赖于 $\rho$ 和 $\eta$ 。如图4所示，可以容易地计算相关函数（15）。

（图4：一对像素之间的、两个失真分量（ $\parallel,\bot$ ）的相关函数（15），绘制 $\rho=4$ 和 $\rho=10$ ）

4 Designing the Field Autocorrelation（场自相关设计）

4.1 Autocorrelation of a Stationary Vector Field（平稳向量场的自相关）

（图5：图像平面。在没有失真的情况下，两个目标点投影到 ${\boldsymbol p}$ 和 ${\boldsymbol p}+{\boldsymbol v}$ ，湍流引起的投影位移分别为 ${\boldsymbol e}({\boldsymbol p})$ 和 ${\boldsymbol e}({\boldsymbol p}+{\boldsymbol v})$ 。任何位移矢量 ${\boldsymbol e}({\boldsymbol p})$ 都可以分成平行于 ${\boldsymbol v}$ 和垂直于 ${\boldsymbol v}$ 的分量，这些分量表示为 ${\boldsymbol e}_{\parallel}({\boldsymbol p})$ 和 ${\boldsymbol e}_{\bot}({\boldsymbol p})$ ）

在本节中，我们直接在二维图像域中工作。我们定义了失真向量场，并在此域中定义了矩阵值协方差函数。考虑图5，图像像素位置 ${\boldsymbol p}$ 是 ${\boldsymbol o}$ 的无失真投影。在没有湍流的情况下，目标 $o^{'}$ 的视线 ${\boldsymbol o'}$ 投影到2D像素位置 ${\boldsymbol p}+{\boldsymbol v}$ 。湍流通过随机向量 ${\boldsymbol e}({\boldsymbol p})$ 和 ${\boldsymbol e}({\boldsymbol p}+{\boldsymbol v})$ 变换各自的投影。因此， ${\boldsymbol e}$ 是图像平面上的一个随机向量场。所有的向量 ${\boldsymbol p}$ 、 ${\boldsymbol v}$ 和 ${\boldsymbol e}$ 都在图像域的 $(x, y)$ 坐标系中给出。具体地说， ${\boldsymbol v}$ 的笛卡尔分量是 $v_x$ 和 $v_y$ 。所有空间分量以及焦距 $f$ 都以像素为单位。对于窄视角： $\theta\approx\frac{|\boldsymbol{v}|}{f}=\frac{r}{f} \tag{16}$ 其中 $r=|\boldsymbol{v}|$ 。因此，参考公式（10,15）： $b_{\parallel,\bot}(\eta)=b_{\parallel,\bot}(\frac{r}{f\theta_D}) \tag{17}$ 为了模拟湍流失真图像，我们寻找一个具有正确空间相关性的随机场 ${\boldsymbol e}$ ，并且它是一个均值为零的空间平稳向量场。由于我们对该场的自相关函数感兴趣，我们可以假定 ${\boldsymbol e}$ 的方差为1，而不失一般性：方差可在之后根据需要使用公式（1）进行调整。自相关函数等于自协方差函数，由2×2矩阵给出： ${\boldsymbol C}({\boldsymbol v})=E[{\boldsymbol e}({\boldsymbol p}){\boldsymbol e}^T({\boldsymbol p}+{\boldsymbol v})] \tag{18}$ 其中 $E$ 表示期望， $T$ 表示转置。注意 ${\boldsymbol C}({\boldsymbol v})$ 是 $v_x$ 和 $v_y$ 的矩阵值函数。在平稳性假设下，自相关函数只依赖于两点间的差 ${\boldsymbol v}$ ，并且在 ${\boldsymbol C}(-{\boldsymbol v})={\boldsymbol C}^T({\boldsymbol v})$ 的意义上是对称的。
我们希望找到一个相关函数（18），使其推导出的纵向和横向自相关函数与公式（15）相匹配。定义单位向量： $\hat{\boldsymbol v}={\boldsymbol v}/r \tag{19}$ 根据上图，位移矢量场 ${\boldsymbol e}({\boldsymbol p})$ 的纵向分量是 ${\boldsymbol e}({\boldsymbol p})$ 在 $\hat{\boldsymbol v}$ 上的投影： $e_{\parallel}({\boldsymbol p})=\left\langle\hat{\boldsymbol v},{\boldsymbol e}({\boldsymbol p})\right\rangle=\hat{\boldsymbol v}^T{\boldsymbol e}({\boldsymbol p}) \tag{20}$ 其中 $\left \langle \cdot,\cdot \right \rangle$ 是内积操作，令： ${\boldsymbol R}=\begin{pmatrix} 0 & -1 \\ 1 & 0 \\ \end{pmatrix} \tag{21}$ 表示一个旋转矩阵，它在图像平面上将 ${\boldsymbol v}$ 旋转90度。因此： $\hat{\boldsymbol v}^T{\boldsymbol R}\hat{\boldsymbol v}=0 \tag{22}$ 垂直位移 $e_{\bot}({\boldsymbol p})$ 是 ${\boldsymbol e}({\boldsymbol p})$ 在垂直向量 ${\boldsymbol R}\hat{\boldsymbol v}$ 上的投影： $e_{\bot}({\boldsymbol p})=\left \langle {\boldsymbol R}\hat{\boldsymbol v},{\boldsymbol e}({\boldsymbol p}) \right \rangle=\hat{\boldsymbol v}^T{\boldsymbol R}^T{\boldsymbol e}({\boldsymbol p}) \tag{23}$ 根据公式（18），平行位移和垂直位移的自相关标量函数分别为： $C_{\parallel}({\boldsymbol v})=E[e_{\parallel}({\boldsymbol p})e_{\parallel}^T({\boldsymbol p}+{\boldsymbol v})]=E[\hat{\boldsymbol v}^T{\boldsymbol e}({\boldsymbol p}){\boldsymbol e}^T({\boldsymbol p}+{\boldsymbol v})\hat{\boldsymbol v}]=\hat{\boldsymbol v}^T{\boldsymbol C}({\boldsymbol v})\hat{\boldsymbol v} \tag{24}$ $C_{\bot}({\boldsymbol v})=E[e_{\bot}({\boldsymbol p})e_{\bot}^T({\boldsymbol p}+{\boldsymbol v})]=\hat{\boldsymbol v}^T{\boldsymbol R}^T{\boldsymbol C}({\boldsymbol v}){\boldsymbol R}\hat{\boldsymbol v} \tag{25}$

4.2 A Solution to the Autocorrelation Function（自相关函数的一种解法）

在给定平行和垂直自相关标量函数 $C_{\parallel}(\boldsymbol v)$ 和 $C_{\bot}(\boldsymbol v)$ 的情况下，我们寻求一个分别与公式（24）和公式（25）一致的矩阵值自相关函数 ${\boldsymbol C}(\boldsymbol v)$ 。解决这样一个问题的办法并不是唯一的。让我们使用解类： ${\boldsymbol C}(\boldsymbol v)=A(\boldsymbol v){\boldsymbol I}-B(\boldsymbol v)\hat{\boldsymbol v}\hat{\boldsymbol v}^T \tag{26}$ 其中 $A(\boldsymbol v)$ 和 $B(\boldsymbol v)$ 是标量函数，且 ${\boldsymbol I}$ 是2×2的单位阵。在公式（26）定义的类中，平行位移和垂直位移不相关。实际上，使用公式（26）： $E[e_{\parallel}({\boldsymbol p})e_{\bot}^T({\boldsymbol p}+{\boldsymbol v})]=E[\hat{\boldsymbol v}^T{\boldsymbol e}({\boldsymbol p}){\boldsymbol e}^T({\boldsymbol p}+{\boldsymbol v}){\boldsymbol R}\hat{\boldsymbol v}]=\hat{\boldsymbol v}^T{\boldsymbol C}({\boldsymbol v}){\boldsymbol R}\hat{\boldsymbol v}=\hat{\boldsymbol v}^TA({\boldsymbol v}){\boldsymbol R}\hat{\boldsymbol v}-\hat{\boldsymbol v}^TB({\boldsymbol v})\hat{\boldsymbol v}\hat{\boldsymbol v}^T{\boldsymbol R}\hat{\boldsymbol v}=0 \tag{27}$ 其中我们使用了公式（22）和函数 $A(\boldsymbol v)$ 和 $B(\boldsymbol v)$ 的标量值性质。
为了确定所需的函数 $A$ 和 $B$ ，我们将公式（26）代入等式（24）和（25），则得到： $C_{\parallel}(\boldsymbol v)=A(\boldsymbol v)\hat{\boldsymbol v}^T\hat{\boldsymbol v}-B(\boldsymbol v)(\hat{\boldsymbol v}^T\hat{\boldsymbol v})^2=A(\boldsymbol v)-B(\boldsymbol v) \tag{28}$ $C_{\bot}(\boldsymbol v)=A(\boldsymbol v)\hat{\boldsymbol v}^T{\boldsymbol R}^T{\boldsymbol R}\hat{\boldsymbol v}-B(\boldsymbol v)\hat{\boldsymbol v}^T{\boldsymbol R}^T\hat{\boldsymbol v}\hat{\boldsymbol v}^T{\boldsymbol R}\hat{\boldsymbol v}=A(\boldsymbol v) \tag{29}$ 其中我们使用了 ${\boldsymbol R}^T{\boldsymbol R}={\boldsymbol I}$ 和公式（22,26,28,29）推导出： ${\boldsymbol C}(\boldsymbol v)=C_{\bot}(\boldsymbol v){\boldsymbol I}-[C_{\bot}(\boldsymbol v)-C_{\parallel}(\boldsymbol v)]\hat{\boldsymbol v}\hat{\boldsymbol v}^T \tag{30}$

4.3 Consistency with Turbulence Correlation（与湍流相关的一致性）

给定公式（10,15,17），平行和垂直标量自相关函数必须设置为满足： $C_{\parallel}(\boldsymbol v)=b_{\parallel}(\frac{r}{f\theta_D}),C_{\bot}(\boldsymbol v)=b_{\bot}(\frac{r}{f\theta_D}) \tag{31}$ 因此，理想的自相关函数可以写成： ${\boldsymbol C}(\boldsymbol v)=A(r){\boldsymbol I}-B(r)\hat{\boldsymbol v}\hat{\boldsymbol v}^T \tag{32}$ 使用正函数： $A(r)=b_{\bot}(\frac{r}{f\theta_D}),B(r)=b_{\bot}(\frac{r}{f\theta_D})-b_{\parallel}(\frac{r}{f\theta_D}) \tag{33}$ 在尺寸为 $U\times U$ 的2D空间图像域中，自相关域 $\boldsymbol v$ 的大小为 $2U\times 2U$ 。

（图6：基于4.3详述的示例参数得到的2×2矩阵值 ${\boldsymbol C}(\boldsymbol v)$ 。在每个子图中， $\boldsymbol v$ 域的空间尺寸为2048×2048像素，对应于1024×1024的图像。右上角的比例尺对应于 $C_{xy}$ 和 $C_{yx}$ ，具有负值和正值。右下角的比例尺对应于 $C_{xx}$ 和 $C_{yy}$ ，具有非负值）

我们给出了一个使用以下参数的示例： $\lambda=550nm$ ，像素大小 $4\mu m$ ， $U = 1024$ 像素，Nikon AF-S Nikkor 镜头， $f = 300 m m$ ， $f\#=5.6$ ，高度 $h = 4 m$ 。目标距离 $L = 2 k m$ 且 $C_n^2=3.6\times10^{-13}m^{-2/3}$ 。然后使用公式（9-15）计算 $A A$ 方差和相关性。如图6和7，绘制2×2矩阵 ${\boldsymbol C}(\boldsymbol v)$ 的每个元素，作为 $\boldsymbol v$ 的函数。使用相同的固定成像参数，在本文中（如下所述）进一步创建图8,9,10,13时，遵循此示例。

（图7：图6的10倍放大图）

# -*- coding: UTF-8 -*-
# Theme: Turbulence-Induced 2D Correlated Image Distortion
# Author: Tianqi Shen
# Date: 2020/08/31

import numpy as np
import scipy.special as sc
from scipy.integrate import quad
import matplotlib.pyplot as plt

class AA_variance:
    def __init__(self):
        self.lam    = 550*1e-9
        self.f      = 300*1e-3
        self.f_     = 5.6      # f#=f_=f/D
        self.D      = self.f/self.f_
        self.L      = 2*1e+3
        self.Cn2    = 3.6*1e-13
        self.beta   = 0.5216
        
    # 公式3
    def spherical_wave_coefficient(self, q):
        temp0 = (np.sqrt(3)/16)*sc.gamma(1/6)*sc.gamma(8/3)*np.power(2/beta,1/3)
        temp1 = (np.pi*(self.beta*q)**2)/(2j)
        temp2 = sc.hyp2f1(1/6,17/6,23/6,(1+temp1))
        temp3 = 1+(6/17)*(8/3)*(np.power(-temp1,1/6)*temp2).real
        return temp0*temp3
    # 公式1
    def variance(self):
        # 公式2
        temp0 = self.D / np.sqrt(self.lam * self.L)
        temp1 = self.spherical_wave_coefficient(temp0)
        temp2 = temp1*self.Cn2*self.L*np.power(self.D,-1/3)
        return temp2

class autocorrelation_function:
    def __init__(self):
        self.lam    = 550*1e-9
        self.pix    = 4*1e-6
        self.f      = 300*1e-3
        self.f_     = 5.6      # f#=f_=f/D
        self.h      = 4
        self.L      = 2*1e+3
        self.D      = self.f/self.f_
        self.thetaD = self.D/self.L
        self.L0     = 0.4*self.h
        self.Rho    = self.L0/self.D

    # 公式12
    def Q(self, xi):
        return np.power((1-xi), 5/3)
    # 公式12
    def M(self, xi, eta):
        return 2 * eta * xi * (1-xi)
    # 公式14
    def H(self, kappa, rho):
        temp0 = np.power((2*rho*kappa)**2+1, -11/6)
        temp1 = np.power(sc.jv(2, kappa), 2)
        return temp0 * temp1 / kappa
    # 公式15分子，里面那个积分的被积函数（平行）
    def integrand_function_para(self, kappa, xi, rho, eta):
        temp0 = sc.jv(0, self.M(xi, eta) * kappa)
        temp1 = sc.jv(2, self.M(xi, eta) * kappa)
        return self.H(kappa, rho) * (temp0 - temp1)
    # 公式15分子，里面那个积分的被积函数（垂直）
    def integrand_function_perp(self, kappa, xi, rho, eta):
        temp0 = sc.jv(0, self.M(xi, eta) * kappa)
        temp1 = sc.jv(2, self.M(xi, eta) * kappa)
        return self.H(kappa, rho) * (temp0 + temp1)
    # 公式15分子，外面那个积分的被积函数（平行）
    def integrate_numerator_para(self, xi, rho, eta):
        temp = quad(self.integrand_function_para, 0, np.inf, args=(xi, rho, eta))
        return temp[0] * self.Q(xi)
    # 公式15分子，外面那个积分的被积函数（垂直）
    def integrate_numerator_perp(self, xi, rho, eta):
        temp = quad(self.integrand_function_perp, 0, np.inf, args=(xi, rho, eta))
        return temp[0] * self.Q(xi)
    # 公式15分子（平行）
    def b_numerator_para(self, rho, eta):
        return quad(self.integrate_numerator_para, 0, 1, args=(rho, eta))[0]
    # 公式15分子（垂直）
    def b_numerator_perp(self, rho, eta):
        return quad(self.integrate_numerator_perp, 0, 1, args=(rho, eta))[0]
    # 公式15分母的被积函数
    def integrate_denominator(self, xi, rho):
        temp = quad(self.H, 0, np.inf, args=rho)
        return temp[0] * self.Q(xi)
    # 公式15分母
    def b_denominator(self, rho):
        return quad(self.integrate_denominator, 0, 1, args=rho)[0]
    # 公式15（kind=0平行，kind=1垂直）
    def correlation_function_b(self, rho, eta, kind):
        if kind == 0:
            return self.b_numerator_para(rho, eta)/self.b_denominator(rho)
        elif kind == 1:
            return self.b_numerator_perp(rho, eta)/self.b_denominator(rho)
        else:
            print("There is no such kind! Please choose 0 or 1!")
    # 公式33
    def A(self, vx, vy):
        temp = np.sqrt(vx**2 + vy**2) / (self.f * self.thetaD)
        return self.correlation_function_b(self.Rho, temp, 1)
    # 公式33
    def B(self, vx, vy):
        temp0 = np.sqrt(vx**2 + vy**2) / (self.f * self.thetaD)
        temp1 = self.correlation_function_b(self.Rho, temp0, 1)
        temp2 = self.correlation_function_b(self.Rho, temp0, 0)
        return temp1 - temp2
    # 公式32
    def C(self, vx, vy):
        temp = vx**2 + vy**2
        Cxx = self.A(vx, vy) - self.B(vx, vy) * vx * vx / temp
        Cxy = -self.B(vx, vy) * vx * vy / temp
        Cyx = -self.B(vx, vy) * vy * vx / temp
        Cyy = self.A(vx, vy) - self.B(vx, vy) * vy * vy / temp
        return np.array([Cxx, Cxy, Cyx, Cyy])
    # 计算自相关域v
    def autocorrelation_domain_v(self, halfU):
        cxx = np.zeros((2*halfU, 2*halfU))
        cxy = np.zeros((2*halfU, 2*halfU))
        cyx = np.zeros((2*halfU, 2*halfU))
        cyy = np.zeros((2*halfU, 2*halfU))
        for i in range(0, 2*halfU):
            for j in range(0, 2*halfU):
                vx = (j+1-0.5-halfU) * self.pix
                vy = (halfU+0.5-i-1) * self.pix
                '''
                print("i:", i)
                print("j:", j)
                print("vx:", vx)
                print("vy:", vy)
                print("temp:", vx**2+vy**2)
                print("B(vx,vy):", B(vx, vy))
                print("###########################")
                '''
                temp = self.C(vx, vy)
                cxx[i][j] = temp[0]
                cxy[i][j] = temp[1]
                cyx[i][j] = temp[2]
                cyy[i][j] = temp[3]
        return cxx, cxy, cyx, cyy

class PlotAll:
    def __init__(self):
        
        acf = autocorrelation_function()
        
        plt.figure('b')
        eta_array = np.linspace(0, 400, 401)
        b = np.vectorize(acf.correlation_function_b)
        b_10_perp = b(10, eta_array, 1)
        b_10_para = b(10, eta_array, 0)
        b_4_perp  = b(4, eta_array, 1)
        b_4_para  = b(4, eta_array, 0)
        plt.plot(eta_array, b_10_perp, color='red', linestyle = '-', zorder=2)
        plt.plot(eta_array, b_10_para, color='red', linestyle = '--', zorder=1)
        plt.plot(eta_array, b_4_perp, color='green', linestyle = '-', zorder=4)
        plt.plot(eta_array, b_4_para, color='green', linestyle = '--', zorder=3)
        plt.legend([r'$b_{\perp}(\rho=10)$', r'$b_{\parallel}(\rho=10)$', r'$b_{\perp}(\rho=4)$', r'$b_{\parallel}(\rho=4)$'], fontsize=17)
        plt.xlabel(r'$\eta=\theta/\theta_D$', fontsize=17)
        plt.ylabel('b', fontsize=17)
    
        cxx, cxy, cyx, cyy = acf.autocorrelation_domain_v(512)
        '''
        print("cxx:", cxx)
        print("cxy:", cxy)
        print("cyx:", cyx)
        print("cyy:", cyy)
        '''
        np.savetxt("cxx.txt", cxx)
        np.savetxt("cxy.txt", cxy)
        np.savetxt("cyx.txt", cyx)
        np.savetxt("cyy.txt", cyy)
        plt.figure("cxx")
        plt.imshow(cxx, cmap='jet')
        plt.figure("cxy")
        plt.imshow(cxy, cmap='jet')
        plt.figure("cyx")
        plt.imshow(cyx, cmap='jet')
        plt.figure("cyy")
        plt.imshow(cyy, cmap='jet')
    
        plt.show()

if __name__ == '__main__':
    plotall = PlotAll()

大家可以试一试上面的代码，可以发现跑的效率非常低！
慢的原因有两个：
① 积分（公式15）较复杂；
② 对每个像素点都进行操作；
建议改小两处的值先试试看部分效果：

eta_array = np.linspace(0, 10, 11)

cxx, cxy, cyx, cyy = acf.autocorrelation_domain_v(6)

代码还有待优化及加速，尽请期待！

西园啜茗

关注

3
点赞
踩
2

收藏

觉得还不错? 一键收藏
0
评论
复制链接

分享到 QQ

分享到新浪微博

扫一扫

专栏目录