网络建站一般多少钱,最近10条重大新闻,开源官网,上街区做网站#x1f3af;要点
#x1f3af;时域、时频域以及时间和频率相关联偏振特性分析三种算法 | #x1f3af;时域波参数估计算法 | #x1f3af;机器学习模型波形指纹分析算法 | #x1f3af;色散曲线和频率相关波分析算法 | #x1f3af;动态倾斜校正算法 | #x1f3af;声…要点
时域、时频域以及时间和频率相关联偏振特性分析三种算法 | 时域波参数估计算法 | 机器学习模型波形指纹分析算法 | 色散曲线和频率相关波分析算法 | 动态倾斜校正算法 | 声学地震波像分析
Python,MATLAB和C/C物理信号滤波合成检测模拟用例
Python和MATLAB数字信号波形和模型模拟
Python量化噪声卷积信号和傅里叶时频分析
C数字化声音信号处理和数控振荡合成
Python嵌入式动态用户调制解调响应式射频信号
Python嵌入式片上系统逻辑电路物理信号处理
Arduino生物波反馈和环境检测外套
Python线性代数傅里叶分析和动态系统模拟分析之一
Python线性代数数字图像和小波分析之二
Python数值和符号算法计算及3D视图物理数学波形方程
MATLAB_ESP32有限脉冲响应FIR无限脉冲响应IIR滤波器
语言内容分比 Python波场逆问题
在许多应用中声波或电磁波被用来观察事物水下声学、雷达、声纳、医学超声波、探地雷达、地震成像、全球地震学。在其中一些应用中测量是被动的我们记录调查对象发出的波并试图推断其位置、大小等。在主动测量中会产生波并记录对象的响应。
为了处理逆问题我们必须首先了解正问题。在最基本的形式中波动方程由下式给出 L [ c ] v ( t , x ) q ( t , x ) ( 1 ) L[c] v(t, x)q(t, x)\qquad(1) L[c]v(t,x)q(t,x)(1) 且 L [ c ] [ 1 c ( x ) 2 ∂ 2 ∂ t 2 − ∇ 2 ] L[c]\left[\frac{1}{c(x)^2} \frac{\partial^2}{\partial t^2}-\nabla^2\right] L[c][c(x)21∂t2∂2−∇2] 其中 v : R × R n → R v: R \times R ^n \rightarrow R v:R×Rn→R 表示波场 c : R n → R c: R ^n \rightarrow R c:Rn→R 是传播速度 q : R × R n → R q: R \times R ^n \rightarrow R q:R×Rn→R 是源术语。我们假设源在空间和时间上都有紧支持并且是平方可积的。
为了定义(1)的唯一解我们需要提供边界和初始条件。在上面讨论的应用中通常考虑无界域 ( x ∈ R n ) \left(x \in R ^n\right) (x∈Rn) 和柯西初始条件 v ( 0 , x ) 0 , ∂ v ∂ t ( 0 , x ) 0 v(0, x)0, \frac{\partial v}{\partial t }(0, x)0 v(0,x)0,∂t∂v(0,x)0。
在散射实验中通常根据入射波场和散射波场重写波动方程 v v i v s vv_iv_s vvivs 和散射势 u ( x ) c ( x ) − 2 − c 0 − 2 u(x)c(x)^{-2}- c_0^{-2} u(x)c(x)−2−c0−2 : L [ c 0 ] v i ( t , x ) q ( t , x ) , L [ c 0 ] v s ( t , x ) − u ( x ) ∂ 2 v ∂ t 2 ( t , x ) ( 2 ) L\left[c_0\right] v_i(t, x)q(t, x), \quad L\left[c_0\right] v_s(t, x)-u(x) \frac{\partial^2 v}{\partial t^2}(t, x)\qquad(2) L[c0]vi(t,x)q(t,x),L[c0]vs(t,x)−u(x)∂t2∂2v(t,x)(2) 在弱散射假设下我们可以忽略 u u u和 v s v_s vs的相互作用并将(2)中的 v v v替换为 v i v_i vi。
测量结果通常被视为限制为 [ 0 , T ] × Δ [0, T] \times \Delta [0,T]×Δ 的分散波场其中 Δ ⊂ R n \Delta \subset R ^n Δ⊂Rn 可以是流形或一组点。然后数据可以用 f ( t , r ) f(t,r) f(t,r)表示。在主动实验中通常的做法是考虑对源项集合的测量。源可以是点源在这种情况下 q ( t , x ) w ( t ) δ ( x − s ) q(t, x)w(t) \delta(x-s) q(t,x)w(t)δ(x−s) 且 s ∈ Σ s \in \Sigma s∈Σ。或者事件场 v i v_i vi 可以由 s ∈ Σ s \in \Sigma s∈Σ 给出和参数化。然后我们可以用 f ( t , s , r ) f(t, s, r) f(t,s,r) 表示数据其中 s ∈ Σ 、 r ∈ Δ s \in \Sigma、r \in \Delta s∈Σ、r∈Δ 和 t ∈ [ 0 , T ] t \in[0, T] t∈[0,T]。
基于这个基本设置我们将讨论三个可能的逆问题
逆源问题从已知且恒定 c ( x ) ≡ c 0 c(x) \equiv c_0 c(x)≡c0 的测量中恢复源 q q q。逆散射假设 c 0 c_0 c0 已知从多个已知源的散射场测量中恢复散射势 u ( x ) u(x) u(x)。波形断层扫描从多个已知源的总波场测量中恢复 c ( x ) c(x) c(x)。
下面您将找到一些实践中出现的逆问题的典型示例。
地震定位由源项 w ( t ) q ( x ) w(t) q(x) w(t)q(x) 描述的地震由多个地震仪在位置 Δ { r k } k 1 n 处记录 \Delta\left\{r_k\right\}_{k1}^n 处记录 Δ{rk}k1n处记录。目标是恢复 q q q 以确定地震位置。被动声纳使用数组 Δ { x 0 r p ∣ r ∈ [ − h , h ] } \Delta\left\{x_0r p \mid r \in[-h, h]\right\} Δ{x0rp∣r∈[−h,h]} 记录从不明目标发出的声波其中 $p \in S ^2 $ 表示数组的方向 h h h 表示其宽度。目标是恢复源项 w ( t ) q ( x ) w(t) q(x) w(t)q(x) 以确定源的起源和性质。雷达成像入射平面波由方向 s ∈ Σ ⊂ S 2 s \in \Sigma \subset S ^2 s∈Σ⊂S2 参数化被发送到介质中并由阵列记录其反射响应。目标是恢复散射势。全波形反演在勘探地震学中目标是从表面总波场的测量中恢复 c c c Σ Δ { x ∣ n ⋅ x 0 } \Sigma\Delta\{x \mid n \cdot x0\} ΣΔ{x∣n⋅x0}。超声断层扫描目标是康复 来自物体周围的源和接收器的总波场。
我们研究逆源问题的一个变体其中 q ( t , x ) δ ( t ) u ( x ) q(t, x)\delta(t) u(x) q(t,x)δ(t)u(x) 和 u u u 在 Ω ⊂ R n \Omega \subset R ^n Ω⊂Rn 上得到紧凑支持。常量 c c c 的前向运算符由下式给出 K u ( t , x ) ∫ Ω u ( x ′ ) g ( t , x − x ′ ) d x ′ K u(t, x)\int_{\Omega} u\left(x^{\prime}\right) g\left(t, x-x^{\prime}\right) d x^{\prime} Ku(t,x)∫Ωu(x′)g(t,x−x′)dx′ 在半径为 ρ \rho ρ 的球体上进行测量。解决逆问题的一种流行技术是反向传播它基于将前向算子的伴随应用于数据。在这种情况下伴随运算符由下式给出 K ∗ f ( x ) ∫ Δ ∫ 0 T g ( t ′ , x ′ − x ) f ( t ′ , x ′ ) d x ′ d t ′ K^* f(x)\int_{\Delta} \int_0^T g\left(t^{\prime}, x^{\prime}-x\right) f\left(t^{\prime}, x^{\prime}\right) d x^{\prime} d t^{\prime} K∗f(x)∫Δ∫0Tg(t′,x′−x)f(t′,x′)dx′dt′ 我们看到 p K ∗ f pK^* f pK∗f 可以通过求解下式得到 L [ c ] w ( t , x ) ∫ Δ f ( t , x ′ ) δ ( x − x ′ ) d x ′ L[c] w(t, x)\int_{\Delta} f\left(t, x^{\prime}\right) \delta\left(x-x^{\prime}\right) d x^{\prime} L[c]w(t,x)∫Δf(t,x′)δ(x−x′)dx′ 使用时间反转格林函数并在 t 0 t0 t0 处求值即 p ( x ) w ( 0 , x ) p(x)w(0, x) p(x)w(0,x)。
为了了解其原理我们研究普通运算符 K ∗ K K^* K K∗K。在时域傅立叶域中对于 c 1 c1 c1运算符变为 f ^ ( ω , x ) ∫ Ω u ( x ′ ) exp ( ı ω ∣ x − x ′ ∣ ) ∣ x − x ′ ∣ d x ′ \widehat{f}(\omega, x)\int_{\Omega} u\left(x^{\prime}\right) \frac{\exp \left(\imath \omega\left|x-x^{\prime}\right|\right)}{\left|x-x^{\prime}\right|} d x^{\prime} f (ω,x)∫Ωu(x′)∣x−x′∣exp(ω∣x−x′∣)dx′ 和 u ( x ) ∬ Δ f ^ ( ω , x ′ ) exp ( − ı ω ∣ x ′ − x ∣ ) ∣ x ′ − x ∣ d x ′ d ω u(x)\iint_{\Delta} \widehat{f}\left(\omega, x^{\prime}\right) \frac{\exp \left(-\imath \omega\left|x^{\prime}-x\right|\right)}{\left|x^{\prime}-x\right|} d x^{\prime} d \omega u(x)∬Δf (ω,x′)∣x′−x∣exp(−ω∣x′−x∣)dx′dω 于是 K ∗ K u ( x ) ∬ Δ ∫ Ω f ( x ′ ) exp ( ı ω ∣ x ′ ′ − x ′ ∣ ) ∣ x ′ ′ − x ′ ∣ exp ( − ı ω ∣ x ′ ′ − x ∣ ) ∣ x ′ ′ − x ∣ d x ′ d x ′ ′ d ω K^* K u(x)\iint_{\Delta} \int_{\Omega} f\left(x^{\prime}\right) \frac{\exp \left(\imath \omega\left|x^{\prime \prime}-x^{\prime}\right|\right)}{\left|x^{\prime \prime}-x^{\prime}\right|} \frac{\exp \left(-\imath \omega\left|x^{\prime \prime}-x\right|\right)}{\left|x^{\prime \prime}-x\right|} d x^{\prime} d x^{\prime \prime} d \omega K∗Ku(x)∬Δ∫Ωf(x′)∣x′′−x′∣exp(ω∣x′′−x′∣)∣x′′−x∣exp(−ω∣x′′−x∣)dx′dx′′dω 对于 ∣ x ′ ′ ∣ ≫ ∣ x ∣ \left|x^{\prime \prime}\right| \gg|x| ∣x′′∣≫∣x∣ 我们可以将其近似为 ∣ x ′ ′ − x ∣ ≈ ∣ x ′ ′ ∣ x ⋅ x ′ ′ / ∣ x ′ ′ ∣ \left|x^{\prime \prime}-x\right| \approx\left|x^{\prime \prime}\right|x \cdot x^{\prime \prime} /\left|x^{\prime \prime}\right| ∣x′′−x∣≈∣x′′∣x⋅x′′/∣x′′∣ 同样对于 ∣ x ′ ′ − x ′ ∣ \left |x^{\prime \prime}-x^{\prime}\right| ∣x′′−x′∣。这称为远场近似。引入 ξ ′ ′ x ′ ′ / ∣ x ′ ′ ∣ x ′ ′ / ρ \xi^{\prime \prime}x^{\prime \prime} /\left|x^{\prime \prime}\right|x^{\prime \prime} / \rho ξ′′x′′/∣x′′∣x′′/ρ单位球面我们发现 K ∗ K u ( x ) ρ ∫ Ω u ( x ′ ) ∬ exp ( w ξ ′ ′ ⋅ ( x ′ − x ) ) d ξ ′ ′ d ω d x ′ K^* K u(x)\rho \int_{\Omega} u\left(x^{\prime}\right) \iint \exp \left(w \xi^{\prime \prime} \cdot\left(x^{\prime}-x\right)\right) d \xi^{\prime \prime} d \omega d x^{\prime} K∗Ku(x)ρ∫Ωu(x′)∬exp(wξ′′⋅(x′−x))dξ′′dωdx′ 此式 k ( x − x ′ ) ∬ exp ( u ω ξ ′ ′ ⋅ ( x ′ − x ) ) d ξ ′ ′ d ω k\left(x-x^{\prime}\right)\iint \exp \left(u \omega \xi^{\prime \prime} \cdot\left(x^{\prime}-x\right)\right) d \xi^{\prime \prime} d \omega k(x−x′)∬exp(uωξ′′⋅(x′−x))dξ′′dω
有时称为点扩散函数。
这种过滤也可以通过乘以 ∣ ξ ∣ 2 |\xi|^2 ∣ξ∣2 在傅立叶域中实现。一个例子如下所示。
import numpy as np
import matplotlib.pyplot as plt
import scipy.sparse as sp
from scipy.ndimage import laplacedef solve(q,c,dt,dx,T1.0,L1.0,n1):gamma dt*c/dxnt int(T/dt 1)nx int(2*L/dx 2)u np.zeros((nx**n,nt))for k in range(1,nt-1):u[:,k1] Au[:,k] Lu[:,k] Bu[:,k-1] (dx**2)*Cq[:,k]return udef getMatrices(gamma,nx,n):l (gamma**2)*np.ones((3,nx))l[1,:] -2*(gamma**2)l[1,0] -gammal[2,0] gammal[0,nx-2] gammal[1,nx-1] -gammaif n 1:a 2*np.ones(nx)a[0] 1a[nx-1] 1b -np.ones(nx)b[0] 0b[nx-1] 0c (gamma)**2*np.ones(nx)c[0] 0c[nx-1] 0L sp.diags(l,[-1, 0, 1],shape(nx,nx))else:a 2*np.ones((nx,nx))a[0,:] 1a[nx-1,:] 1a[:,0] 1a[:,nx-1] 1a.resize(nx**2)b -np.ones((nx,nx))b[0,:] 0b[nx-1,:] 0b[:,0] 0b[:,nx-1] 0b.resize(nx**2)c (gamma)**2*np.ones((nx,nx))c[0,:] 0c[nx-1,:] 0c[:,0] 0c[:,nx-1] 0c.resize(nx**2)L sp.kron(sp.diags(l,[-1, 0, 1],shape(nx,nx)),sp.eye(nx)) sp.kron(sp.eye(nx),sp.diags(l,[-1, 0, 1],shape(nx,nx)))return A,B,C,LL 1.0
T 1.0
dx 1e-2
dt .5e-2
nt int(T/dt 1)
nx int(2*L/dx 2)
c 1u np.zeros((nx,nx))
u[nx//2 - 10:nx//210,nx//2 - 10:nx//210] 1
q np.zeros((nx*nx,nt))
q[:,1] u.flatten()w_forward solve(q,c,dt,dx,TT,LL,n2)theta np.linspace(0,2*np.pi,20,endpointFalse)
xs 0.8*np.cos(theta)
ys 0.8*np.sin(theta)
I np.ravel_multi_index(np.array([[xs/dx nx//2],[ys//dx nx//2]],dtypenp.int), (nx,nx))
d w_forward[I,:]r np.zeros((nx*nx,nt))
r[I,:] d
r np.flip(r,axis1)w_adjoint solve(r,c,dt,dx,TT,LL,n2)fig, ax plt.subplots(nrows2, ncols3, figsize(8, 5), shareyTrue)
plt.gray()ax[0,0].plot(xs,ys,r*)
ax[0,0].imshow(w_forward[:,2].reshape((nx,nx)), extent(-L,L,-L,L))
ax[0,0].set_title(t 0)
ax[0,1].plot(xs,ys,r*)
ax[0,1].imshow(w_forward[:,101].reshape((nx,nx)), extent(-L,L,-L,L))
ax[0,1].set_title(t 0.5)
ax[0,2].plot(xs,ys,r*)
ax[0,2].imshow(w_forward[:,200].reshape((nx,nx)), extent(-L,L,-L,L))
ax[0,2].set_title(t 1)ax[1,0].plot(xs,ys,r*)
ax[1,0].imshow(w_adjoint[:,200].reshape((nx,nx)), extent(-L,L,-L,L))
ax[1,1].plot(xs,ys,r*)
ax[1,1].imshow(w_adjoint[:,101].reshape((nx,nx)), extent(-L,L,-L,L))
ax[1,2].plot(xs,ys,r*)
ax[1,2].imshow(w_adjoint[:,2].reshape((nx,nx)), extent(-L,L,-L,L))plt.show()参阅更新计算思维 | 亚图跨际
