美国投掷 6 枚钻地炸弹袭击伊朗福尔多核设施，钻地弹有多大杀伤力？能摧毁伊朗地下核设施吗？

import numpy as np
import rasterio
import matplotlib.pyplot as plt
from skimage.feature import register_translation
from scipy.ndimage import shift

# 加载卫星图像（前后图）
pre = rasterio.open('fordow_pre.tif').read(1)
post = rasterio.open('fordow_post.tif').read(1)

# 图像对齐：估计位移
shift_yx, error, diffphase = register_translation(pre, post, 100)
post_aligned = shift(post, shift=(shift_yx[0], shift_yx[1]))

# 差分图（变化检测）
diff = post_aligned.astype(float) - pre.astype(float)

# 整合模拟爆炸坑（二值边缘）
y, x = np.indices(diff.shape)
center = np.array(diff.shape) / 2
R = np.sqrt((x - center[1])**2 + (y - center[0])**2)
radius = 5  # km->像素，根据分辨率调整
circle = (np.abs(R - radius) < 1).astype(float)

# 可视化对比叠加
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.title('打击前后变化差值（热度图）')
plt.imshow(diff, cmap='RdBu', vmin=-np.percentile(diff,5), vmax=np.percentile(diff,95))
plt.colorbar(label='强度增减')

plt.subplot(1, 2, 2)
plt.title('模拟坑边界叠加真图')
plt.imshow(post_aligned, cmap='gray')
plt.contour(circle, levels=[0.5], colors='red', linewidths=2)
plt.scatter([center[1]], [center[0]], c='yellow', s=50, label='预计爆心')
plt.legend()
plt.show()

import numpy as np
import matplotlib.pyplot as plt
import rasterio
from rasterio.plot import show
from skimage.feature import register_translation
from scipy.ndimage import shift
import geopandas as gpd

# ──────────────── 1. 图像加载 ────────────────
with rasterio.open('fordow_pre.tif') as src_pre, \
     rasterio.open('fordow_post.tif') as src_post:
    pre = src_pre.read(1)
    post = src_post.read(1)
    meta = src_pre.meta

# ──────────────── 2. 图像对齐 ────────────────
shift_yx, error, diffphase = register_translation(pre, post, 100)
post_aligned = shift(post, shift=(shift_yx[0], shift_yx[1]))

# ──────────────── 3. 差分分析 ────────────────
diff = post_aligned.astype(float) - pre.astype(float)

# ──────────────── 4. InSAR 相干性（示例图） ────────────────
phase = rasterio.open('interferogram_phase.tif').read(1)

# ──────────────── 5. 模拟破坏区域 GIS 輪廓 ────────────────
circle = gpd.GeoDataFrame(
    {'geometry': [gpd.points_from_xy([meta['width']/2], [meta['height']/2])[0].buffer(500)],
     'damage': ['simulated']},
    crs=src_pre.crs
)

# ──────────────── 6. 可视化输出 ────────────────
fig, axs = plt.subplots(2, 2, figsize=(14,12))

# 差分热图
im0 = axs[0,0].imshow(diff, cmap='RdBu',
                      vmin=-np.percentile(diff,5),
                      vmax=np.percentile(diff,95))
axs[0,0].set_title('EO/SAR 差分热图')
plt.colorbar(im0, ax=axs[0,0])

# 投后图 + 破坏轮廓
axs[0,1].imshow(post_aligned, cmap='gray')
circle.boundary.plot(ax=axs[0,1], color='red', linewidth=2, label='模拟破坏区域')
axs[0,1].legend()
axs[0,1].set_title('后期影像 + 模拟破坏叠加')

# InSAR 相干性图
im1 = axs[1,0].imshow(phase, cmap='jet')
axs[1,0].set_title('InSAR 干涉相位图')
plt.colorbar(im1, ax=axs[1,0])

# DEM + 模型叠加（假设有dem.tif）
with rasterio.open('dem.tif') as dem_src:
    dem = dem_src.read(1)
axs[1,1].imshow(dem, cmap='terrain', alpha=0.6)
axs[1,1].imshow(post_aligned, cmap='gray', alpha=0.4)
circle.boundary.plot(ax=axs[1,1], color='red', linewidth=2)
axs[1,1].set_title('DEM + EO + 模拟破坏')

plt.tight_layout()
plt.show()