全球回荡的冲击波

原文链接: https://www.nv5geospatialsoftware.com/Learn/Blogs/Blog-Details/the-shockwave-heard-around-the-world

17400 给这篇文章评分:

4.1

全球回荡的冲击波

Jim Pendleton 2022年5月10日，星期二

洪阿汤加-洪阿哈阿帕伊火山喷发大气波的 GOES-16、GOES-17 和 Meteosat-8 影像动画

继2022年1月15日南太平洋洪阿汤加-洪阿哈阿帕伊火山喷发之后，最初在线发布的一些数据产品和可视化成果给了我启发。特别是，由马萨诸塞大学洛厄尔分校的 Matthew Barlow 在 EOS 科学新闻网站上展示的辐射率差值动画，这些动画使用了 NOAA 的 GOES-17 高级基线成像仪（ABI，由 L3Harris 开发）数据，揭示了一个非常有趣的现象。在动画中，地球大气层因喷发冲击波而像钟一样鸣响，类似于固体地球在地震后的“鸣响”。

观察地面气压传感器的研究人员报告称，在冲击波环绕地球的数天时间里，可以多次探测到大气扰动。我希望在时间和空间范围上扩展最初的中层对流层动画，看看这种持续效应是否在大气更高层也可见。（提示：是的！）

对于此示例，我下载了“全盘 ABI L2 云和湿度影像”产品，并将影像限制为仅使用第9波段。这个水汽波段位于红外波段的6.9微米中心，被称为“中层对流层水汽波段”。在该波长波段中看不到地表特征，因此在最终动画中显然需要大陆轮廓线来帮助观看者定位。

GOES-16和GOES-17卫星分别从西经75.2度和137.2度的地球静止轨道收集多光谱全盘ABI影像数据。（上图显示的是GOES-17数据的一个示例。）

我在美国的NV5 Geospatial同事David Starbuck最近发表了一篇博客，其中包含了一个通过IDL和AWS下载GOES-16和GOES-17数据的工具，所以我将此练习留给读者。

与此同时，EUMETSAT的Meteosat-8卫星上的自旋增强可见光和红外成像仪（SEVIRI）也收集全盘影像。在此次事件发生时，它位于东经41.5度、印度洋上方的地球静止轨道。它的视角覆盖了GOES-16和GOES-17几乎看不到的地球区域。该探测器的第6波段波长范围中心位于7.35微米，其波段范围与GOES ABI仪器的波段范围重叠。

由于洪阿汤加火山（南纬20.545度，西经175.394度）的对跖点大致位于马里、阿尔及利亚和尼日尔共享边界的上方，Meteosat-8的空间覆盖提供了在压力波首次环绕地球的中途观察其干涉图样的机会。

EUMETSAT网站提供了Meteosat-8数据集的访问权限。下载历史数据的选项需要注册，但您可以通过免费的在线帐户申请访问。登录帐户后，提供了多种下载数据的机制。我选择了API端点和产品“EO:EUM:DAT:MSG:MSG15-RSS - 快速扫描高速率SEVIRI Level 1.5影像数据 – MSG”，然后按日期浏览。典型的图像文件下载网页如左图所示。

图像文件采用EUMETSAT的https://api.eumetsat.int/data/browse/collections “MSG”格式，扩展名为“.nat”，如上图所示。（对Meteosat MSG格式文件的支持最近在ENVI® 5.6.2中引入）。这些文件非常大，因此我编写了一些ENVI代码，仅为后续分析提取感兴趣的单个波段。

pro meteosat_extract_band, band_number = band_number

compile_opt strictarr

on_error, 2

band = band_number ne !null ? band_number : 6

e = envi(/current)

if (e eq !null) then e = envi(/headless)

dir = dialog_pickfile(Title = 'Select directory of Meteosat .nat files', $

    /dir, /must_exist)

if (dir eq '') then return

files = file_search(filepath(root = dir, 'MSG1-*.nat'), count = c)

if (c eq 0) then return

copybandmetadata = [ $

    'DATA GAIN VALUES', $

    'DATA OFFSET VALUES', $

    'BAND NAMES', $

    'WAVELENGTH', $

    'FWHM']

foreach file, files, index do begin

  print, (index + 1).tostring() + '/' + c.tostring(), ' ', file

  bandfile = filepath(root = file_dirname(file), $

      file_basename(file, '.nat') + '-Band-' + band.tostring() + '.dat')

  ; don't reprocess if the band file exists
  if (~file_test(bandfile)) then begin
      r1 = e.openraster(file)
      d1 = r1[0].getdata(band = band - 1)
      newraster = ENVIRaster(d1, URI = bandfile, INHERITS_FROM = r1[0])
      metadata = r1[0].metadata
      newmetadata = newraster.metadata
      foreach m, copybandmetadata do begin
          newmetadata[m] = (metadata[m])[band - 1]
      endforeach
      newraster.save
      newraster.close
  endif

  if (r1 ne !null && obj_valid(r1[0])) then foreach r, r1 do r.close
  r1 = !null

endforeach

end

在任何单一时间的单个图像中，压力波并不容易显现。然而，通过对两个相邻时间步长的图像制作差值影像，可以增强这一特征。辐射率差值影像实际上就是一个时间点的辐射率与前一个时间步长辐射率的差值。

pro meteosat_delta_images

compile_opt strictarr

on_error, 2

e = envi(/current)

if (e eq !null) then e = envi(/headless)

 dir = dialog_pickfile($

    Title = 'Select directory of Meteosat band files for time differencing', $

    /dir, /must_exist)

if (dir eq '') then return

files = file_search(filepath(root = dir, 'MSG1-*-Band-*.dat'), count = c)

files = files.sort() ; ensure they're in time order

if (c eq 0) then return

r0 = e.openraster(files[0])

d0 = r0[0].getdata()

gain = (r0[0].metadata)['DATA GAIN VALUES']

offset = (r0[0].metadata)['DATA OFFSET VALUES']

mask = where(d0 eq 0)

d0 = temporary(d0)*gain + offset

foreach r, r0 do r.close

for i = 1, c - 1 do begin

    print, (i + 1).tostring() + '/' + c.tostring(), ' ', files[i]

    r1 = e.openraster(files[i])

    d1 = r1[0].getdata()

    mask = where(d1 eq 0)

    d1 = temporary(d1)*gain + offset

    delta = d1 - d0

    delta[mask] = 0

    fn = filepath(root = file_dirname(files[i]), $

        file_basename(files[i], '.dat') + '_diff.dat')

    if ~file_test(fn) then begin

        nr = ENVIRaster(delta, URI = fn, INHERITS_FROM = r1[0])

        nr.save

        nr.close

    endif

    foreach r, r1 do r.close

    d0 = d1

endfor

end

对GOES ABI数据集也重复了此过程。对差值影像进行了视觉检查，以从序列中移除任何“噪声”影像。此外，我下载的图像集在时间上跨越两个或多个传感器重叠，因此另一个步骤涉及寻找“最佳”图像序列，以在所有三个传感器上制作连贯的动画。我为最终视频选择了 GOES-17 – GOES-16 – Meteosat-8 – GOES-17 的序列。

为这些数据集生成图像动画的方法有很多。在此示例中，我使用了经典的 IDL®，结合直接图形和对象图形，生成图像叠加和注释。虽然原始影像每边可能超过5,000像素，但对于GOES ABI数据，所有三个数据集的动画输入将标准化为每边仅1,024像素。

以下代码片段获取了针对Meteosat-8数据、在东经41.5度的地球静止卫星视角下，大陆轮廓线在1024 x 1024缓冲区中的像素位置。

window, /free, /pixmap, xsize = 1024, ysize = 1024

map_set, sat_p=[35786023.0,0,0], /satellite, 0., 41.5, /noborder ; East

map_continents, /hi

mask = tvrd(channel = 0)

wdelete, !d.window

xy = array_indices(mask, where(mask ne 0))

（对于GOES数据，只需根据每颗卫星的位置适当修改经度。）

接下来，对于每个传感器，我创建了一个IDL SAVE格式文件，其中包含所有带有灰度图像数据以及大陆叠加层的帧。它们存储在一个有序哈希中，其键设置为从文件名派生的图像时间戳。（请注意，NOAA和METEOSAT的文件命名约定在时间戳方面使用略有不同的格式。）

我制作了两个数据集，一个包含绝对值钳制在10.0的辐射率差值，另一个包含钳制在1.0的差值。正如您将在动画中看到的，10.0的值比1.0的差值数据产生“更平滑”、噪声更少的图像，但它不能那么显著或长时间地显示压力波。

pro generate_meteosat_frames

compile_opt strictarr

on_error, 2

dir = dialog_pickfile($

    Title = 'Select directory of Meteosat difference images', $

    /dir, /must_exist)

if (dir eq '') then return

files = file_search(filepath(root = dir, 'MSG1-*Band-6_diff.dat'), count = c)

if (c eq 0) then return

window, /free, /pixmap, xsize = 1024, ysize = 1024

map_set, sat_p=[35786023.0,0,0], /satellite, 0., 41.5, /noborder ; East

map_continents, /hi

mask = tvrd(channel = 0)

wdelete, !d.window

xy = array_indices(mask, where(mask ne 0))

deltas = ['10.', '1.']

window, /free, xsize = 1024, ysize = 1024, /pixmap

e = envi(/current)

if (e eq !null) then e = envi(/headless)

foreach delta, deltas do begin

    frames = orderedhash()

    for i = 0, c - 1 do begin

        print, (i + 1).tostring() + '/' + c.tostring(), ' ', files[i]

        r1 = e.openraster(files[i])

        d1 = r1[0].getdata()

        t = ((file_basename(files[i], '.dat')).split('-'))[5]

        tvscl, congrid($

            reverse((d1 < float(delta)) > (-float(delta)), 2), 1024, 1024)

        plots, xy[0, *], xy[1, *], psym = 3, color = 'ff0000'x, /device

        frames['c' + (t.split('\.'))[0]] = tvrd(/true)

        foreach r, r1 do r.close

    endfor

    save, frames, /compress, file = filepath(root=file_dirname(files[0]), $

        'Meteosat-8 SEVIRI Band 6_Delta_' + delta + '_Tonga.sav')

endforeach

wdelete, !d.window

end

对两个GOES数据集重复了这些步骤，此练习留给读者。经过阈值处理和缩放的辐射率差值图像，加上大陆叠加层，然后按时间戳在所有三个数据集上进行排序，并为每个图像添加适当的文本注释。

deltas = ['10.', '1.']

foreach delta, deltas do begin

  restore, 'GOES16_ABI_Band_9_Delta_' + delta + '_Tonga.sav'

  allframes = frames

  allsources = hash()

  foreach k, frames.keys() do begin

    allsources[k] = ‘16’

  endforeach

  restore, 'GOES17 + '_ABI_Band_9_Delta_' + delta + ‘ _Tonga.sav'

  foreach k, frames.keys() do begin

    allsources[k] = ‘17’

    AllFrames[k] = frames[k]

  endforeach

  restore, 'Meteosat-8 SEVIRI Band 6_Delta_' + delta + '_degree_Tonga.sav'

  foreach k, frames.keys() do begin

    allsources[k] = ‘M’

    AllFrames[k] = frames[k]

  endforeach

  outframe = 1

  ks = (AllFrames.keys()).sort()

  time0 = 15*24+4+10/60. ; approximate eruption time

  foreach k, ks, index do begin

    source = allsources[k]

    day = strmid(k, 6, 2)

    hr = strmid(k, 8, 2)

    mi = strmid(k, 10, 2)

    t1 = day*24+hr+mi/60.

    deltat = t1 - time0

    hh = long(deltat)

    mm = abs(round((deltat - hh)*60.))

    time = 'Time from eruption (HH:MM) = ' + $

      (deltat lt 0 ? '-' : '') + hh.tostring() + ':' + string(mm, format = '(I2.2)')

    if (index eq 0) then begin

      i = image(AllFrames[k], dim=[1024, 1024], margin = [0, 0, 0, 0])

      t0  = text(.01, .95, ["Hunga Tonga-Hunga Ha'apai eruption", $

        'Sources: NOAA-NASA GOES-16 & -17 ABI L2 Band 9', $

        'Meteosat-8 SEVIRI Band ' + SEVIBand, $

        'Radiance Difference (+/- ' + delta + ')'], $

        /normal, color = 'Yellow')

      t  = text(.01, .92, time, /normal, color = 'Yellow')

    endif else begin

      i.setdata, AllFrames[k]

      t.string = time

    endelse

    write_png, filepath('frame' + string(outframe++, format = '(I3.3)') + '.png'), $

      i.copywindow()

  endforeach

  i.close

endforeach

在此示例中，我首先将帧单独写入不同目录中的PNG输出文件以便检查。然而，接下来的视频创建步骤可以集成到该逻辑中。

pro tongamovie

compile_opt strictarr

on_error, 2

d1 = dialog_pickfile(/must_exist, /dir, $

    title = 'Select directory of 1-delta PNG images')

if (d1 eq '') then return

d2 = dialog_pickfile(/must_exist, /dir, $

    title = 'Select directory of 10-delta PNG images')

if (d2 eq '') then return

n = file_search(d1 + '*.png', count = c)

if (c eq 0) then return

; write the MP4 video to the same directory as where this routine lives

oVid = IDLffVideoWrite(filepath(root=routine_dir(), 'tonga.mp4'), FORMAT='mp4')

fps = 3

; put the 10 and 1 radiance difference images side-by-side

newx = Round(1920/2*.75)

stream = oVid.AddVideoStream(newx*2, newx, fps)

fin = bytarr(3, 2*newx, newx)

for i = 1, c do begin

    f = 'frame' + string(i, format = '(i3.3)') + '.png'

    i1 = congrid(read_png(filepath(root = d1, f)), 3, newx, newx)

    i2 = congrid(read_png(filepath(root = d2, f)), 3, newx, newx)

    fin[0, newx, 0] = i1

    fin[0, 0, 0] = i2

    r = ovid.put(stream, fin)

endfor

end

请注意，在动画中，压力波在完全传播回起源点后仍然清晰可见，这一特征在任何单个静态图像中可能都很难用肉眼识别出来。我尚未进一步跟进该特征在红外波段可见的时间有多长。我猜想，随着波的能量缓慢消散，需要将辐射率钳制阈值（在此示例中当前设置为1.0）设置得越来越低，信号才会完全消失在背景噪声中。

可以安全地假设，在未来的岁月里，许多高级学位的论文都将基于这一独特事件而诞生。