David P.Williams论文系列合成孔径声纳图像快速目标检测:一种新算法及大规模性能分析

郭铁拳

已于 2022-11-25 15:30:16 修改

阅读量1.4k

点赞数 3

文章标签：目标检测人工智能

于 2022-11-22 20:55:57 首次发布

本文链接：https://blog.csdn.net/weixin_44072988/article/details/127844001

版权

摘要

本文提出了一种新的无监督合成孔径声纳(SAS)水下目标探测算法。该方法利用高分辨率的高质量SAS图像，这些图像上的目标有很多像素。该方法的一个特别新颖的组成部分还可以检测沙波纹并估计它们的方向。通过采用级联结构和利用积分图像表示，使整个算法加速。因此，该方法可以通过自主水下航行器(AUV)上的声纳数据实现近实时探测目标。不需要训练数据，因为提出的方法是自适应定制的环境特征的感知数据收集在原地。为了验证和评估所提出的检测算法的性能，对包含各种类水雷目标的SAS图像进行了大规模研究。数据是由MUSCLE AUV在2008年至2012年期间在不同地理位置和不同环境条件下进行的六次大型海上实验中收集的。分析检查了探测性能作为目标类型、相位、范围、图像质量、海底环境和地理位置的函数。据我们所知，这项研究基于近3万张SAS图像，总共覆盖约160公里的海床，涉及1100多个目标探测机会，代表了迄今为止使用SAS数据对目标探测性能进行的最广泛的系统、定量评估。通过分析揭示了对目标探测性能影响最大的变量，即图像质量和海底环境条件。还提出了利用传感器数据进行自适应AUV测量的方法。

Index Iterms——算法，AUV，水雷对抗，性能分析，合成孔径声呐

介绍

在军事行动中，必须打击敌人部署水雷，以确保部队和物资的安全运输。在和平时期的民用部门，冲突后遗留下来的这些地雷可能引爆，有可能破坏全球供应链和航运路线。在这两种情况下，可靠地探测这些隐蔽的战争武器的能力是至关重要的。声纳提供的水下环境的高分辨率成像证明对探测海底的水雷特别有用。由于海洋机器人技术的进步，用于处理这一任务的声纳数据通常由自主水下航行器(AUV)收集。由于这种地雷对抗(MCM)任务的内在危险性和时间敏感性，下一个紧迫的优先事项是在AUV中嵌入智能，以便它能够立即对其收集的数据作出反应。通过就地调整测量路线和有效分配资源，AUV可以为手头的任务收集最丰富的数据，同时减少任务时间，改善“数据泛滥”和操作人员疲劳的问题。然而，这样做需要一种算法，能够在处理能力有限的自动水下航行器上近乎实时地进行可靠的地雷探测。现有的检测算法在研究社区受到一些限制，排除了这种可能性。

高分辨率合成孔径声纳[2]的出现意味着探雷任务现在大大容易了。SAS系统提供的优秀分辨率直接转化为目标上的更多像素，这是基于图像的检测[3]最关键的参数。因此，利用这种改进的图像分辨率的相对简单的检测算法足以实现良好的检测性能。这项工作提出了一种计算速度快、易于调整的探测器设计，旨在并利用高分辨率的高质量SAS图像，其高分辨率允许在目标上有许多像素。

现有的大多数水雷探测算法，最初设计用于低分辨率侧扫描声纳图像，而不是SAS图像，只是搜索水雷[4]-[9]特征的高光阴影模式。(一个显著的例外是[10]，它是专门为双波段、高频和宽带SAS数据开发的。)高光是来自水雷本身的声波回声的结果，而投射的阴影是由于水雷(特别是它高于海底的高度)和传输信号的掠射角度之间的几何形状。不幸的是，大多数现有的水雷探测算法至少有四个主要限制中的一个。首先，所有这些检测算法都隐含地(错误地)假设整个图像的图像质量一致优秀。其次，许多检测算法没有充分利用目标回波和阴影的距离依赖特性，而这些特性是可以根据几何形状预测的。第三，大多数检测算法并没有根据测量发生的环境条件量身定制，而是依赖于来自不同地点的训练数据；除了可能导致的潜在数据不匹配之外，使用训练数据还会限制预期可以检测到的目标类型。第四，一些检测算法，如[8]，依赖于相对缓慢的匹配滤波型方法，采用离散傅里叶变换，排除了近实时检测的可能性。有一个例外相对于速度而言，是[9]中提出的算法也使用了积分图像，尽管与我们提出的方法相比以一种根本不同的方式。[9]中的算法受到前面提到的其他限制，它通过在一种需要访问一组训练数据进行学习的蛮力策略中考虑数千个基于haar的通用特征来合并检测和分类阶段。

haar特征常用于人脸检测，它是利用积分图像计算的，描述灰度图像的变化。

本文有两个主要贡献。第一个贡献是引入了一种新的算法，用于在声纳图像中探测水下目标，克服了一些限制，使AUV上的近实时探测成为可能。这项工作的第二个主要贡献是进行了大规模的性能评估，以验证所提出的算法，并研究在SAS图像中检测水雷目标的能力。该分析通过使用2008年至2012年在不同地理位置和不同环境条件下进行的6次大型海上实验中收集的真实SAS数据，检验了探测性能作为目标类型、相位、距离、图像质量和海底环境的函数。据我们所知，这项研究基于近3万张SAS图像，总共覆盖约160公里的海床，涉及1100多个目标探测机会，代表了迄今为止使用SAS数据对目标探测性能进行的最广泛的系统、定量评估。分析揭示了对水雷探测性能影响最大的变量，这可以进一步为智能、自适应、完全自主水下探测的发展提供信息，不需要任何“人工参与”。

虽然这里的工作主要集中在水雷探测任务上，但算法也可以用于广泛的其他与目标探测相关的水下应用，包括栖息地测绘[11]、海底分类[12]、考古[13]和管道监测[14]。重要的是，通过简单的修改，检测算法可以用于侧扫描声纳图像，这仍然比SAS图像更常用，假设物体的尺寸相对于声纳分辨率仍然很大。但所提出的方法也适用于为基于图像处理的任务使用其他传感器模式的领域。例如，快速检测框架可以用于地面目标检测[15]，而波纹检测算法将用于需要检测沙丘[16]的石油勘探应用。

本文的其余部分组织如下。第二节详细介绍了所提出的检测算法。第三节解释了真实的、测量的声纳数据，第四节给出了探测结果。第五节给出了算法和结果的简要讨论，第六节给出了结束语。

提出的检测算法

综述

这里详细描述的目标检测算法是一种基于整体图像的方法，能够在配备声纳的AUV上执行近实时检测。该算法的早期初步版本首先在[17]中发表(在会议论文中)，然后，在纹波检测阶段发展之后，在[18]中发表，但随后还进行了进一步的重大修改。在这些变化中，包括加入了有助于降低误报率的侧回声计算，在过滤步骤中使用上下文信息来拒绝共同构成不现实的高密度接触的报警，以及一个新的最终检测置信度评分(基于两个新特性)，更准确地对报警进行排序。新的最终检测分数的实质性价值将在第四节中详细说明。此外，这项工作很重要，而且足够新颖，因为此前从未进行过对检测算法的大规模性能评估。

首先，我们简要概述了用于创建所提议的水雷探测算法的动机和一般原理。总体目标是创建一个健壮的算法，可以可靠地检测构成威胁的水下物体。然而，该算法也必须在处理能力有限的AUV上快速执行。由于不允许人工干预，算法必须足够灵活，能够从现场收集的数据感知和适应不断变化的环境条件。同时，提出的方法应该直接解决和克服困扰现有检测算法的局限性。为了满足这些要求，开发了广泛的特定于领域的问题研究方法。在可能的情况下，重点还放在根据应用程序的基本物理和几何结构调整算法上。结果是严格和灵活性的混合。此外，级联算法体系结构通过在三个主要阶段(阴影检测、波纹检测和回声检测)的每个阶段对图像的越来越小的部分进行操作，从而最小化计算成本。

第II-B-II-E节将描述检测方法的每一步，通过让其他研究人员复制该方法。由于该方法非常适合可视化，因此算法描述将由图来补充说明每个步骤对一个特别具有挑战性的声纳图像的影响，如图1所示。

该SAS图像将被用来作为一个运行的例子来解释检测算法的每一步。图像的左侧以沙波为特征，右侧为平坦的沙质海床。图片中还有7个类似地雷的目标，它们的位置将在后面公布。图像包含2001×7333像素。

阴影检测

所提出的检测算法的第一阶段建立图像中阴影的位置。作为图像预处理步骤，对声呐图像进行了简单的中值归一化处理。(具体来说，具体来说，设 $m_{i}^{r}$ 为图像的第 $i$ 行的像素中值，然后将第 $i$ 行中的每个像素值除以 $m_{i}^{r}$ 归一化，并对每一行重复此操作。设 $m_{i}^{c}$ 为图像的 $j$ 列的像素中值，然后将第 $j$ 列中的每个像素值除以 $m_{i}^{c}$ 归一化，并对每一行重复此操作)。然后将中位归一化声纳图像转换为积分图像表示。随后，利用该积分图像快速估计背景图和阴影图，通过对比提供阴影(与地雷一致)的位置，在算法的第二阶段进一步检查。

此处的python代码如下所示：

def img_mean_norm(img):
    h, w = img.shape[:2]#h是行，w是列
    mr = np.mean(img, axis=0)#列均值
    mc = np.mean(img, axis=1)#行均值
    result = np.zeros((h, w), dtype=np.float64)
    for u in range(h):
        result[u,:] = img[u,:]/mc[u]
        for v in range(w):
            result[:,v] = img[:,v]/mr[v]
    return result

效果如下：

1）积分图像：

积分图像[19]是一种图像表示，它允许在恒定时间内在任何尺度或位置上快速计算矩形、类似哈尔的特征(因为计算不依赖于输入的大小)。在算法的后续阶段，这样的矩形特征对于评估物体的某些区分特征(如回声和阴影级别)将是无价的。使用积分图像表示也是使流检测成为可能的关键。(在处理积分图像计算时，应注意确保使用足够的数值精度，以避免意外的截断错误。)

因此，我们不是在声纳图像的基于像素的系统上操作，而是立即转换为一个积分图像系统。从原始声纳图像出发，构建相应的积分图像如下：(需要注意的是，下面的变量是垂直方向运行的，水平方向运行的。)积分图像中某一位置(x,y)的值对应于(x,y)的上方和左侧(包括)的像素的和。这就是

$I(x,y)=\sum _{x^{'}\leq x,y^{'}\leq y}A(x^{'},y^{'})$

利用递归关系可以快速生成积分图像：

$I(x,y)=I(x-1,y)+z(x,y)$

其中 $z(x,y)$ 是原始图像的一行像素的累计和：

$Z(x,y)=Z(x,y-1)+A(x,y)$

图1中声纳图像对应的积分图像表示如图2所示。

这里要算积分图像的话直接用OpenCV的库函数就行，积分图像的处理会改变矩阵的大小，各自从axis=0和axis=1轴上延伸1个单位长度，可以看见这里的第一行和第一列都是0，（看上面的公式用到了递归，出现了x-1，所以要补0）此处的python代码如下所示：
imageIntegral = cv2.integral(np.uint8(img_men), None, cv2.CV_32FC1)
#imageIntegral = cv2.normalize(imageIntegral,None,0,255,cv2.NORM_MINMAX)#这个是归一化处理
plt.imshow(imageIntegral)
plt.colorbar()
plt.show()
#img_men是经过中值化处理的图像矩阵
imageIntegral = imageIntegral[1:,1:]
#删除补0部分
下面是上述图像展示，下面是经过标准化处理的。

2）背景估计：

积分图像的首次应用是在声纳图像背景图 $B$ 的估计中。背景图的目的是建立海床的混响水平，从而确定图像中阴影的位置。

混响水平是海底成分的一个重要函数。例如，软粘土海床的混响水平将低于硬填料沙子海床。事实上，海底的自然变化会导致混响水平在一个特定地点甚至在同一幅图像内发生巨大变化。

因此，在这项工作中，通过使用测量数据中的海底局部特征来获得背景的可靠估计，[20]和[21]中也是这样做的。根据数据本身调整估计所提供的固有灵活性提高了该方法在成功处理不同环境时的鲁棒性。自适应估计还消除了训练数据不匹配[22]，[23]的可能性，因为没有使用训练数据。

背景估计使用分割窗口模板[24]计算，该模板由三个大小相等(堆叠)的水平带(即矩形)组成，其中顶部和底部带贡献背景值，而中间带不贡献背景值。每个带的大小都一样，沿轨道方向 $b_x = 2.5m$ ，距离方向上 $b_y=5m$ ，在选择这些尺寸时，所关心的对象(即水雷)将完全包含在中间波段内，因此不影响背景估计。(对于涉及其他类型对象的不同应用程序，这些参数将相应地进行调整。)

在 $(x,y)$ 处的背景得分，即 $B(x,y)$ ， $B(x,y)$ 是模板以 $(x,y)$ 为中心的两个矩形中的平均像素值。由于积分图像表示，在给定位置上构成背景值的两个矩形的计算可以快速计算，从积分图像中总共只有8个数组引用。具体来说，在 $(x,y)$ 处的背景得分为：

$B(x,y)=a_{b}^{-1}[I(x-3\frac{\delta _x}{2},y-\frac{\delta _y}{2})-I(x-3\frac{\delta _x}{2},y+\frac{\delta _y}{2})\\-I(x-\frac{\delta _x}{2},y-\frac{\delta _y}{2})+I(x-\frac{\delta _x}{2},y+\frac{\delta _y}{2})\\+I(x+\frac{\delta _x}{2},y-\frac{\delta _y}{2})-I(x+\frac{\delta _x}{2},y+\frac{\delta _y}{2})\\ -I(x+3\frac{\delta _x}{2},y-\frac{\delta _y}{2})+I(x+3\frac{\delta _x}{2},y+\frac{\delta _y}{2})]$

其中 $\delta _x$ 和 $\delta _y$ 是对应于 $b_x$ 和 $b_y$ 的像素， $a_{b}$ 是总和中涉及的像素总数，因此结果是两个矩形区域的平均像素值。

两个矩形的面积是一样的，均为 $\delta _x\times \delta _y$ ，所以有 $a_{b}=2(\delta _x\times \delta _y)$ 。

紧围绕一个位置的区域(即“中间带”)不用于背景计算，因为这些像素，如果对象存在，将与对象相关，并对背景估计产生偏差。背景归一化方案[25]同样试图在计算中排除信号能量。使用图2中的积分图像得到的背景图如图3所示，从图3中可以看出，背景水平确实在图像中不同。

上述部分的python代码如下所示，还是要进行补零的操作：

def backgroundscore(img,deta_x,deta_y):
    a_b = 2*deta_x*deta_y
    h, w = img.shape[:2]#h是行w是列
    temp = np.pad(img, ((int(1.5*deta_x), int(1.5*deta_x)),(int(0.5*deta_y), int(0.5*deta_y))), 'constant')
    x, y = temp.shape[:2]
    back_result = np.zeros((x,y), dtype=np.float64)
    for u in range(int(h)):
        for v in range(int(w)):
            back_result[u,v] = (temp[u - int(1.5*deta_x), v - int(0.5*deta_y)] -\
                                temp[u - int(1.5*deta_x), v + int(0.5*deta_y)] -\
                                temp[u - int(0.5*deta_x), v - int(0.5*deta_y)] +\
                                temp[u - int(0.5*deta_x), v + int(0.5*deta_y)] +\
                                temp[u + int(0.5*deta_x), v - int(0.5*deta_y)] -\
                                temp[u + int(0.5*deta_x), v + int(0.5*deta_y)] -\
                                temp[u + int(1.5*deta_x), v - int(0.5*deta_y)] +\
                                temp[u + int(1.5*deta_x), v + int(0.5*deta_y)])/a_b
    return back_result

利用上述结果，可以得到如下图所示：

3）阴影估计：

然后，再次使用积分图像构造阴影图。阴影图的目的是帮助确定图像中可能由感兴趣的物体产生的阴影的位置。

回想一下，探测算法的目标是探测海底的水雷。我们所拥有的广泛的特定领域的知识——即，可能的水雷的大小和形状——使我们能够估计出这种物体可能延伸到海底以上的最小高度。但这一问题的独特几何结构，包括AUV、声纳和海底物体，为探测算法提供了大量可利用的信息。任何物体在(假定为平坦的)海底上方都会投下阴影。此外，投射阴影的长度将是物体和声纳的相对位置的函数。对于在海底以上 $a$ 高度作业的AUV，以及在距离AUV $r$ 范围内的海底高度水雷，简单的几何图形决定了物体投射的影子长度：

$s_{y} = \frac{hr}{a-h}$

$h$ 的值设置为我们感兴趣的探测对象的最小高度，从相关水雷库中提取（这里使用的为 $h=0.25m$ ） $a$ 是直接在AUV上测量的。

为了尊重问题的几何不可侵犯性，用于阴影计算的矩形模板的大小是范围相关的，有效地将对简单阴影的搜索精炼为对适当长度的阴影的搜索。这种选择为该方法增加了一层额外的严格性。计算中使用沿着轨道的方向长为 $s_x=0.5m$ ， $s_{y}$ 就是上面的这个式子。

在 $(x,y)$ 上的阴影得分是矩形上的平均像素值[围绕着 $(x,y)$ ]，由于积分图像表示，再次计算迅速为 $S(x,y)=a_s^{-1}[I(x-\frac{\delta _x}{2},y)-I(x-\frac{\delta _x}{2},y+\delta _y)\\-I(x+\frac{\delta _x}{2},y)+I(x+\frac{\delta _x}{2},y+\delta _y)]$

其中 $\delta _x$ 和 $\delta _y$ 是 $s_x$ 和 $s_y$ 对应的像素值， $a_s$ 是涉及的像素总数，结果就是矩形区域的平均像素值。使用图2中的积分图像得到的阴影映射如图4所示。

需要注意的是，由于采用了积分图像表示，模板的范围相关特性不会导致额外的计算复杂性。相比之下，现有的检测算法采用匹配滤波类型的方法，依赖于离散傅里叶变换，不能享受这种范围相关的灵活性。相反，这些方法在本质上是不匹配的(因为问题的几何结构不受尊重)和低劣的；也就是说，要么使用单一的平均范围(这不利于在特定范围内检测目标)，要么针对不同的预期阴影大小重复整个过程多次，这增加了计算时间。

上述部分的代码如下所示，依然需要补零的操作：

def shadowscore(img,sx,h,r,a):
    sx = int( sx )
    sy = int(( h * r ) / ( a - h ))
    a_s = sx * sy
    h, w = img.shape[:2]  # h是行w是列
    temp = np.pad(img, ((int(0.5 * sx), int(0.5 * sx)), (sy, sy)), 'constant')
    x, y = temp.shape[:2]
    shadow_result = np.zeros((x,y), dtype=np.float64)
    for u in range(int(h)):
        for v in range(int(w)):
            shadow_result[u, v] = (temp[u - int(0.5*sx), v] - temp[u - int(0.5*sx), v + sy] 
                             -temp[u + int(0.5*sx), v] + temp[u + int(0.5*sx), v + sy])/a_s
    return shadow_result

可以见得此处的阴影部分比上面那一块跟加明显。

4）解决图像质量差的问题：

具有良好图像质量的声纳图像是指物体和海底特征与声纳系统理想工作条件下观察到的情况大致相同的图像。通常认为所有声纳图像的完整范围都是高质量的，但在海上[26]收集到的真实声纳数据往往不是这样。

特别是在远距离，由于传播和吸收(即衰减)，图像质量往往会显著下降，从而导致信号强度减弱，或多径效应[27]。此外，使用SAS，成功成像的运动校正要求在远距离[28]更加严格。因此，声纳系统最远距离的数据图像质量往往是不够的。

SAS的工作原理是将阵列中重叠元素的接收到的声纳信号(或ping返回)连贯相加。两个接收的离散时间信号 $s_i$ 和 $s_j$ 之间的相互关系为

$k_{ij}[n]=\sum _{m=-\infty}^{\infty}s_i[m]s_j^{*}[n+m]$

其中 $n$ 是时移(在实践中，总和是由信号长度决定的有限界计算的)。连续ping信号之间的归一化峰值相关性，称为相干度[29]，是：

$\rho = \underset{n}{max}\left | \frac{k_{ij[n]}}{\sqrt{k_{ii}[0]k_{jj}[0]}} \right |$

在这项工作中，连续ping返回的平均相干值作为距离的函数，被用作图像质量的替代度量。该量在SAS处理的DCPA(位移相位中心天线)算法[28]中计算，与信噪比[30]直接相关。

具体而言， $r$ 距离的相干性估计如下。当拥有声速的估计值时，可以确定与距离中心 $r$ 的范围窗口 $[r-\Delta r,r+\Delta r]$ 相对应的信号返回的部分。这些信号部分用于上面的相干性计算，结果为构成图像的 $P$ 对连续ping的第 $i$ 个的 $\rho _i(r)$ 。然后计算相干性 $\rho (r)=1(1/P)\sum _{i=1}^{P}\rho _{i}(r)$ 的平均值，并在ping对之间执行平均值，以提高鲁棒性，但要牺牲顺轨粒度。如果均值相干性 $\rho (r)$ 小于预设阈值，与范围 $[r-\Delta r,r+\Delta r]$ 相关的图像索引 $(x,y)$ 被认为图像质量较差。在每张图像中，这一过程是对距离中心 $r$ 的离散集合执行的，通常为 $\Delta r = 5m$ 。

根据我们对声纳图像的大量视觉检查经验，我们认为在阈值 $\tau _{\rho}=2/3$ (对应3-dB信噪比)以下的范围内的图像质量较差，包含不可靠的阴影信息。(在[29]中采用一个类似的临时阈值 $\tau _{\rho}=3/4$ 来声明图像质量较差。)因此，只要选择一个阈值，就可以很容易地获得二值图像质量地图。

5）感兴趣区域的确定（阶段1）：

在构建背景图 $B$ 和阴影图 $S$ 后，确定图像中可能包含感兴趣目标的区域。任何阴影地图值足够低于相应背景地图值的像素将被声明为将接受进一步调查的感兴趣区域。具体地说，如果 $S(x,y)<\beta B(x,y)$ ， $(x,y)$ 被认为是影子的一部分，因此，这是一个值得进一步调查的区域。缩放因子使得对阴影的要求更加严格(这里， $\beta =0.5$ )。本质上，这定义了当一个像素的值小于周围背景的平均值一定数量时，它就是阴影区域的一部分。背景和阴影图对比的二值化结果，结合上面4）中的内容的二值化图像质量图，如图5所示，其中白色像素对应潜在报警区域，需要进一步研究。远距离的白色垂直条纹(图像质量较差的地方)是图像质量地图的贡献。通过使用背景地图搜索阴影，利用了当地的上下文信息。算法的鲁棒性也通过关注问题的内在方面而得到增强:感兴趣的物体必须高出海底一定的高度，这意味着会投射出特定长度的影子。(该算法不适用于完全埋在地下的水雷。)

背景和阴影图比较是检测级联的第一个数据缩减阶段。它实现了必须进一步检查的像素的大规模减少(在实践中，超过90%的SAS图像像素通常在这一步中被删除)，从而大大减少了所有后续阶段的计算负载。(回想一下，最初的SAS图像包含超过1400万像素。)由于在声纳图像中自然出现的高光(不是由海床上方的物体引起的)是常见的，相反，首先搜索高信号强度的位置不会使数据减少如此显著。

此处的python代码如下：

def interest_area_stag1(img_back,img_shadow,beta):
    h, w = img.shape[:2]  # h是行w是列
    stag1_result = np.zeros((h, w), dtype=np.float64)
    img_back_beta = beta * img_back
    for u in range(int(h)):
        for v in range(int(w)):
            if img_shadow[u, v] < img_back_beta[u, v]:
                stag1_result[u, v] = 255
            else:
                stag1_result[u, v] = 0
    return stag1_result

第一张是 $\beta =0.5$ ，第二张是 $\beta =0.05$ 的效果图，我们可以发现，当 $\beta$ 减小时，周围的点也是相应减小的。

波纹检测

所提出的检测算法的第二阶段检测声纳图像中的波纹。一个特征提取步骤将每个阴影映射到一个新空间中的一个点，在这个新空间中可以有效地检测到波纹，其特征是定向相似的高密度拉长阴影。在算法的第三阶段进一步检查不属于波纹的位置。

1）动机：

沙纹是水下水雷探测中最容易混淆的因素之一，也是误报的频繁来源，因为沙纹会产生高光阴影图案——其大小取决于波纹的振幅、周期和范围——这是类雷物体的特征。此外，当地雷位于波纹场时，地雷和波纹的回声可以混合在一起，它们各自的阴影也可以混合在一起，这使得探测具有挑战性。

由于这里的总体目标是探测地雷，所以我们提出的纹波检测算法更准确地称为麻烦纹波检测算法。该算法不是试图检测图像中所有的沙波纹，而是只尝试检测那些影响我们探测地雷能力的波纹。

2）概述：

这里开发的原则性纹波检测算法既适合定义沙纹的基本特征，也适合创建它们的基础物理。当水流或波浪向特定方向流动时，就会形成沙纹。由于这种创建波纹的基本机制，波纹往往存在于大型字段中，而不是作为单个实例。但是波纹场，根据定义，必须包含一定密度的波纹。

虽然波纹的特征是一个主要的方向(由流动的方向决定)，但也会存在相当大的变化。也就是说，波纹不会在一个固定的方向上显现出来。相反，他们会围绕某个主导的方向跨越一系列取向。此外，波纹的振幅和周期也可能不同。因此，要能够检测到这种可变性，就需要一种具有内在灵活性的方法。此外，声纳和构成波纹的沙堆之间的相对几何关系意味着波纹投下的阴影将沿着波纹波峰(或波谷)[34]的方向拉长。

在图1中可以观察到的这三个波纹的基本物理特性，将在后续的波纹检测算法阶段中加以利用。具体地说，所提出的纹波检测算法搜索方向相似的高密度拉长阴影。为了高效地执行这种搜索，从背景-阴影图比较中获得的每个离散阴影区域(参见图5)被映射到一个新的特征空间中，在这个空间中很容易检测到这些特征。然后采用了一种新颖的、非常规的积分图像公式的密度估计用途。

3）特征提取：

首先，为二进制阴影图中的每个离散阴影(即每个连接区域，或“blob”(斑点区域)计算质心和方向。(没有使用形态学运算符来更改阴影图。)阴影的方向(以度为单位) $\phi$ 是通过将椭圆拟合到阴影区域，然后校正不相等的像素大小而获得的。质心被转换为距离特征(以米为单位) $\chi$ ，它只是质心到任务中固定参考位置的距离(对于这里的单幅图像，图像的“第一个”像素：在初始沿轨道位置的最近距离的像素)。这种特征提取过程有效地将每个阴影映射到 $\mathbb{R}^{2}$ 中新特征空间中的特定位置 $(\chi,\phi )$ 。(为了稍后将该空间作为图像处理，这些特征在方向-维上离散为1°，在位置-维上离散为1m。)

不过，在继续之前，有些阴影要从进一步的考虑中去除。如果 $\phi <d_\phi/2$ 或 $\phi > 180^{\circ}-d_\phi/2$ ，表示阴影几乎是水平的，因此忽略阴影。(在本作品中， $d_\phi/2 = 30^{\circ}$ ）之所以做出这样的选择，是因为根据声纳和海床的几何形状，在那些被禁止的近水平方向上的大部分阴影很可能是由类似地雷的物体(而不是波纹)造成的。此外，这种选择消除了在即将到来的积分图像密度估计中，由于不尊重定向的“环绕”性质(从180°到0°)而产生的任何不一致。将图5中的声纳图像空间中的阴影转化为位置-方向空间的特征提取过程的结果如图6所示。

此处没有明确的给出代码的思路，我的思路首先找出整个图像的所有连通区域，这里可以用OpenCV的函数如下，
cv2.findContours()
来进行遍历，得到所有的连通区域的元组数据，后面仍然用OpenCV的函数进行椭圆拟合,
cv2.fitEllipse()
利用上述函数进行遍历，得到椭圆拟合后的结果，它的输出是：

ellipse = [ (x, y) , (a, b), angle ]
                    （x, y）代表椭圆中心点的位置
                    （a, b）代表长短轴长度，应注意a、b为长短轴的直径，而非半径
                      angle 代表了中心旋转的角度
这里我们可以得到椭圆的中心点，还有椭圆的角度（阴影的角度），这里是说将提取的位置特征和方向特征一起映射到了位置-方向空间中的一个点上，所以说在这里每个连通域在新空间一张图上就是一个点。

所以说横轴就是我们得到的所有角度的范围，纵轴就是质心与极点的距离。这样就可以绘制出对应的图像，然后进行后续的操作，具体代码如下所示：
cnts, hierarchy = cv2.findContours(img_interest_stag1_0,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)
#imag = cv2.drawContours(img_interest_stag1_1,img_countours,3,(0,255,0),3)
#img_interest_stag1_0_el = cv2.ellipse(img_interest_stag1_0,ellipse,(0,255,0),2)
def get_ellipse(cnts,d_fai):
    all_ellipse = []
    all_location = []
    all_angle = []
    for x in range(len(cnts)):
        if len(cnts[x]) >= 5:
            ellipse = cv2.fitEllipse(cnts[x])
            #ellipse = [ (x, y) , (a, b), angle ]
            # (x, y)代表椭圆中心点的位置
            # (a, b)代表长短轴长度，应注意a、b为长短轴的直径，而非半径
            #angle 代表了中心旋转的角度
            all_ellipse.append(ellipse)
            if  pd.isna(ellipse[0][1]) == False and (ellipse[2] < (d_fai/2) or ellipse[2] > (180 - (d_fai/2))):#在这里加判决条件
                all_location.append(ellipse[0])
                all_angle.append(round(ellipse[2]))
    return all_ellipse,all_location,all_angle
all_ellipse,all_location,all_angle = get_ellipse(cnts,30)
def get_distance(location,x,y):#(x,y)是极点
    all_distance = []
    for u in range(len(location)):
        d = ((location[u][0] - x)**2 + (location[u][1] -y )**2)**0.5
        all_distance.append(round(d))
    return all_distance
all_distance = get_distance(all_location,0,0)
def get_distance_angle_2_img(angle,distance):
    max_angle, max_distance = max(angle),max(distance)
    result_img = np.ones((max_distance + 1,max_angle + 1), dtype=np.float64)*255
    for (d,a) in zip(distance,angle):
            result_img[int(d), int(a)] = 0
    return result_img
img_2_da = get_distance_angle_2_img(all_angle,all_distance)
plt.imshow(img_2_da)
plt.colorbar()
plt.show()
由此得到的图是这样的：

这里还是要调一下参数的，这里用的是原文中的参数 $d_\phi/2 = 30^{\circ}$ 得到的结果，我是椭圆拟合了5个点以及5个点以上的连通域，得到的结果，我觉得这里应该也要调一下。

我们可以看到90度左右的点特别多，观察原图也可以得到，斜向下的连通域也特别多。

4）密度估计：

通过将该新特征空间视为二维“图像” $A'$ ，可以容易地计算相应的积分图像 $I'$ 。然后可以利用积分图像公式在该空间中快速执行密度估计。具体而言，将给定位置附近矩形区域中离散阴影的平均数量用作阴影密度估计。

计算中使用的矩形大小为 $d_{\chi }=25m$ 在距离维度上，在方向维度上是 $d_{\phi }=30^{\circ}$ ，选择该尺寸是因为它考虑了有助于定义波纹场的潜在机制的物理接近性和方向可变性。在 $(\chi ,\phi )$ 上的密度估计 $D(\chi ,\phi )$ 为取为以 $(\chi ,\phi )$ 为中心的矩形上的平均像素值（尽管与方向无关），计算如下所示：

$D(\chi ,\phi )=d_{\chi }^{-1}[I'(\chi-\frac{d_{\chi }}{2},\phi -\frac{d_{\phi }}{2})-I'(\chi-\frac{d_{\chi }}{2},\phi +\frac{d_{\phi }}{2})\\-I'(\chi+\frac{d_{\chi }}{2},\phi -\frac{d_{\phi }}{2})+I'(\chi+\frac{d_{\chi }}{2},\phi +\frac{d_{\phi }}{2})]$

尽管密度估计中包含了一系列方向，但平均值仅在位置维度上计算，因为我们不希望在类似方向上区分波纹（因为一般波纹场中自然存在微小的方向变化）。根据图6计算得出的阴影密度估计值如图7所示。

此处很简单，代码如下所示：

def dense_img(img,dx,df):
    imIntegral = cv2.integral(np.uint8(img), None, cv2.CV_32FC1)
    imIntegral = imIntegral[1:, 1:]
    h, w = imIntegral.shape[:2]#h是行w是列
    temp = np.pad(imIntegral, ((int(0.5*dx), int(0.5*dx)),(int(0.5*df), int(0.5*df))), 'constant')
    x, y = temp.shape[:2]
    dense_result = np.zeros((x,y), dtype=np.float64)
    for u in range(int(h)):
        for v in range(int(w)):
            dense_result[u, v] = (temp[(u - int(0.5*dx), (v - int(0.5*df)))]-
                                  temp[(u - int(0.5*dx), (v + int(0.5*df)))]-
                                  temp[(u + int(0.5*dx), (v - int(0.5*df)))]+
                                  temp[(u + int(0.5*dx), (v + int(0.5*df)))]) / (dx)
    return dense_result

结果如下：