实验二 边缘识别实验报告(正文)
边缘化提取实验报告
一、实验目的1. 理解图像边缘检测的基本原理和过程。
2. 掌握常用的边缘检测算法，如Roberts算子、Sobel算子、Prewitt算子、Laplacian算子和Canny算子。
3. 通过实验验证不同边缘检测算法的效果，并分析其优缺点。
4. 了解特征提取的基本原理和方法，对图像边缘进行特征提取。
二、实验原理图像边缘是图像中灰度值或颜色值发生突变的地方，是图像分割和特征提取的基础。
边缘检测的目的是找到图像中灰度值变化明显的区域，即边缘。
边缘检测算法可以分为两类：基于微分算子的边缘检测算法和基于二值化的边缘检测算法。
1. 基于微分算子的边缘检测算法：- 利用一阶导数或二阶导数检测图像边缘。
- 常见的算子有Roberts算子、Sobel算子、Prewitt算子、Laplacian算子等。
2. 基于二值化的边缘检测算法：- 利用图像的二值化处理，将图像分为前景和背景两部分。
- 常见的算法有Otsu算法、Sauvola算法等。
三、实验内容1. 实验材料：- OpenCV库- Python编程环境2. 实验步骤：（1）读取图像：使用OpenCV库读取待检测的图像。
（2）灰度化：将图像转换为灰度图像，以便进行边缘检测。
（3）边缘检测：- 使用Roberts算子检测边缘。
- 使用Sobel算子检测边缘。
- 使用Prewitt算子检测边缘。
- 使用Laplacian算子检测边缘。
- 使用Canny算子检测边缘。
（4）特征提取：对检测到的边缘进行特征提取，如计算边缘长度、宽度、方向等。
（5）结果展示：将检测到的边缘和提取的特征进行可视化展示。
四、实验结果与分析1. Roberts算子：- 效果：Roberts算子对图像噪声敏感，边缘检测效果较差。
- 分析：Roberts算子对图像局部区域进行检测，容易受到噪声的影响。
2. Sobel算子：- 效果：Sobel算子对图像噪声有一定的抑制能力，边缘检测效果较好。
- 分析：Sobel算子使用高斯滤波器对图像进行平滑处理，然后计算图像的一阶导数。
边缘检测实验报告
边缘检测实验报告边缘检测实验报告引言：边缘检测是图像处理中的一项重要任务，它能够有效地提取图像中物体的边界信息，为后续的图像分割、物体识别等任务提供基础。
本实验旨在探究不同的边缘检测算法在不同场景下的表现，并对其进行评估和比较。
一、实验背景边缘检测是图像处理领域的经典问题，早期的边缘检测算法主要基于梯度的计算，如Sobel、Prewitt等。
随着深度学习的发展，基于卷积神经网络的边缘检测方法也取得了显著的进展。
本实验将选择传统的Sobel算子和基于深度学习的Canny算法进行对比。
二、实验步骤1. 数据准备：选择一组包含不同场景、不同复杂度的图像作为实验数据集，确保数据集的多样性和代表性。
2. 算法实现：使用Python编程语言，利用OpenCV库实现Sobel算子和Canny 算法。
对于Sobel算子，我们将尝试不同的卷积核大小和阈值设置。
对于Canny算法，我们将调整高低阈值的取值范围。
3. 实验评估：使用评估指标来衡量不同算法的性能，如准确率、召回率、F1值等。
同时，我们还可以通过可视化的方式来比较不同算法的边缘检测效果。
三、实验结果在实验中，我们选择了10张不同类型的图像进行边缘检测，并使用Sobel算子和Canny算法进行处理。
通过对比实验结果，我们得出以下结论：1. Sobel算子：- 当卷积核大小较小（如3x3）时，Sobel算子能够较好地检测到图像中的细节边缘，但对于噪声较多的图像效果较差。
- 当卷积核大小较大（如7x7）时，Sobel算子能够更好地抑制噪声，但会导致边缘检测结果的模糊。
- 阈值的设置对Sobel算子的效果也有较大影响，较低的阈值可以提高边缘检测的敏感性，但也容易引入噪声。
2. Canny算法：- Canny算法基于梯度的计算和非极大值抑制，能够有效地检测到图像中的边缘，并且对噪声有较好的鲁棒性。
- 高低阈值的设置对Canny算法的效果影响较大，合适的阈值范围可以提高边缘检测的准确性和稳定性。
实验二 边缘识别实验报告(正文)
一、基本原理 (3)1.1 位场边缘识别方法研究进展概述 (3)1.2 数值计算类边缘识别方法的研究及计算 (3)1.2.1 垂向导数（VDR） (3)1.2.2 总水平导数（THDR） (3)1.2.3 解析信号振幅（ASM） (3)1.2.4 倾斜角（TA） (4)1.2.5 二阶导数类边缘识别 (4)二、输入/输出数据格式设计 (4)2.1 主要变量设计 (4)2.2位场边缘识别程序数据格式设计 (5)三、总体设计 (6)四、测试结果 (6)4.1 测试用例 (6)4.2 测试结果 (7)4.3 结果分析 (9)五、结论与建议 (10)5.1 结论 (10)5.2 建议 (10)附录：边缘识别程序源代码 (10)一、基本原理1.1 位场边缘识别方法研究进展概述地质目标体边缘是指断裂构造线、不同地质体的边界线，实际上是具有一定密度或磁性差异的地质体的边界线．在地质体的边缘附近，重、磁异常变化率较大，故所有的边缘识别方法均利用这一特点进行设计。
现有重、磁位场边缘识别方法分为数理统计、数值计算和其他三大类。
1.2 数值计算类边缘识别方法的研究及计算数值类边缘识别方法均利用极大值位置或零值位置确定地质体的边缘位置，其依据的理论基础是二度体铅垂台阶模型的重力异常特征．在该模型边缘处重力异常总水平导数和解析信号振幅达到极大值、垂向导数达到零值．故可以利用这些特征位置来确定二度体铅垂台阶的边缘位置．确定倾斜二度体、不规则二度体及三度体边缘位置的理论均为二度体铅垂台阶模型理论的推广，但确定的边缘位置与真实位置有一定偏差．该偏差随着地质体边界形状、埋深、水平尺寸及物性差异等的变化而变化．因此，边缘识别结果是一种定性或半定量解释结果，与定量解释结果有一定区别，识别结果可作为边缘位置定量反演的初值。
1.2.1 垂向导数（VDR）垂向导数方法利用零值位置确定地质体的边缘位置，重力异常可直接使用，对磁力异常必须转换成磁源重力异常（假重力异常）或化极磁力异常才可以使用。
实验二图像分割与边缘检测
实验二、图像分割与边缘检测一、实验目的依据边缘检测的理论，实现灰度图像一阶和二阶边缘检测方法，启发学生依据边缘特征进行图像分析与识别，提高学生图像处理与分析能力和实际动手能力。
二、实验内容1 编程实现一阶差分边缘检测算法，包括Roberts梯度算子、Prewitt算子和Sobel算子。
2 编程实现二阶差分Laplace边缘检测算法。
3 分析与比较各种边缘检测算法的性能。
LOG算子（高斯拉普拉斯算子）三、实验仪器1计算机；2 MA TLAB程序；四、实验报告内容1叙述实验过程；2提交实验的原始图像和结果图像。
五、思考题1．利用一阶导数和二阶导数特性检测图像边缘，分析检测结果的主要区别。
2．在分析几种典型的边缘检测方法的基础上，设计一种通用的边缘特征检测方法，使之具有上述算子的功能。
六、代码I = imread('coins.png');BW1 = edge(I,'roberts');BW2 = edge(I,'prewitt');BW3 = edge(I,'sobel');BW4 = edge(I,'log');BW5 = edge(I,'canny');figuresubplot(3,2,1);imshow(I),title('原图像');subplot(3,2,2);imshow(BW1),title('robert图像');subplot(3,2,3);imshow(BW2);title('prewitt图像');subplot(3,2,4);imshow(BW3);title('sobel图像');subplot(3,2,5);imshow(BW4);title('log图像');subplot(3,2,6);imshow(BW5);title('canny图像');差分边缘检测方法是最原始、基本的方法。
边缘检测matlab实验报告
边缘检测matlab实验报告引言边缘检测在图像处理领域中是一项十分重要的任务。
它可以帮助我们从图像中提取出物体的边缘信息，对于图像分割、目标识别等任务都具有重要意义。
本实验旨在通过利用MATLAB中提供的边缘检测函数，实现对图像中边缘的提取，并对实验结果进行分析和探讨。
实验步骤1. 导入图像首先，我们需要从MATLAB工作环境中导入需要进行边缘检测的图像。
我们可以使用`imread`函数将图像读入到MATLAB的内存中。
matlabimage = imread('example.jpg');2. 灰度化灰度化是边缘检测的前提条件，它可以将一幅彩色图像转化为灰度图像，使得后续的操作更加简便。
我们可以使用`rgb2gray`函数将彩色图像转化为灰度图像。
matlabgray_image = rgb2gray(image);3. 边缘检测接下来，我们可以使用MATLAB中提供的边缘检测函数进行实际的边缘检测操作。
MATLAB中有许多边缘检测算法可供选择，例如Sobel算子、Canny算子等。
本实验我们选择使用Canny算子进行边缘检测。
matlabedge_image = edge(gray_image, 'Canny');4. 结果显示最后，我们可以使用`imshow`函数将原始图像和边缘检测结果显示出来，以便于观察和分析。
matlabsubplot(1, 2, 1);imshow(gray_image);title('原始图像');subplot(1, 2, 2);imshow(edge_image);title('边缘检测结果');5. 结果分析通过以上步骤，我们可以得到原始图像和边缘检测结果。
我们可以观察边缘检测结果，进一步分析图像中的边缘信息。
同时，我们还可以对不同的边缘检测算法进行对比实验，以评估它们的性能和适用性。
实验结果下图展示了使用Canny算子进行边缘检测的实验结果。
详细的图像分割之边缘检测实验报告
边缘检测实验报告一、实验目的通过课堂的学习，已经对图像分割的相关理论知识已经有了全面的了解，知道了许多图像分割的算法及算子，了解到不同的算子算法有着不同的优缺点，为了更好更直观地对图像分割进行深入理解，达到理论联系实际的目的，特制定如下的实验。
二、实验原理：图像处理有两大类目的：1．改善像质（增强、恢复）；2．图像分析：对图像内容作出描述；其一般的图像处理过程如下：图像分割的算法有：（1）阈值分割原理：(,)(,)(,)EBLf x y Tg x y L f x y T≥⎧=⎨<⎩（2）边缘检测：梯度对应一阶导数，对于一个连续图像函数f(x,y):梯度矢量定义:梯度的幅度:梯度的方向:a) Roberts 算子b) Sobel 算子Roberts 算子[]TTyxy f x f G G y x f ⎦⎤⎢⎣⎡∂∂∂∂==∇),(122)()),((),(y x G G y x f mag y x f +=∇=∇)arctan(),(x y G y x =φ()()()[]()()[]{}21221,,11,1,,+-++++-=j i f j i f j i f j i f j i gc) Prewitt 算子d) Kirsch 算子由K 0~K 7八个方向模板组成，将K0~K7的模板算法分别与图像中的3×3区域乘，选最大一个值，作为中央像素的边缘强度（3）区域分割1 区域生长法 算法描述先对每个需要分割的区域找一个种子像素作为生长的起点，然后将种子像素周围邻域中与种子像素有相似性质的像素合并到种子像素所在的区域中。
将这些新像素当作新的种子像素继续进行上面的过程，直到再没有满足条件的像素可被包括进来。
2 分裂合并法实际中常先把图像分成任意大小且不重叠的区域，然后再合并或分裂这些区域以满足分割的要求,即分裂合并法.一致性测度可以选择基于灰度统计特征(如同质区域中的方差),假设阈值为T ,则算法步骤为:① 对于任一Ri ，如果 ，则将其分裂成互不重叠的四等分； ② 对相邻区域Ri 和Rj ，如果 ，则将二者合并； ③ 如果进一步的分裂或合并都不可能了，则终止算法。
图像的边缘检测实验报告
图像的边缘检测实验报告
《图像的边缘检测实验报告》
图像的边缘检测是计算机视觉领域中的重要技术之一，它可以帮助我们识别图
像中物体的边缘和轮廓，从而实现图像分割、特征提取和目标识别等应用。
在
本次实验中，我们将对几种常用的边缘检测算法进行比较和分析，以评估它们
在不同场景下的性能和适用性。
首先，我们使用了Sobel算子进行边缘检测。
Sobel算子是一种基于梯度的边缘检测方法，它通过对图像进行卷积操作来寻找像素值变化最大的地方，从而找
到图像中的边缘。
实验结果显示，Sobel算子在一些简单场景下表现良好，但
在复杂背景和噪声干扰较大的情况下效果不佳。
接着，我们尝试了Canny边缘检测算法。
Canny算法是一种多阶段的边缘检测
方法，它通过对图像进行高斯滤波、计算梯度、非极大值抑制和双阈值处理等
步骤来检测图像中的边缘。
实验结果显示，Canny算法在复杂场景下表现出色，能够有效地抑制噪声并找到图像中的真实边缘。
最后，我们还尝试了Laplacian算子和Prewitt算子等其他边缘检测算法，并对
它们的性能进行了比较和分析。
实验结果显示，不同的边缘检测算法在不同场
景下表现出各自的优势和劣势，需要根据具体的应用需求来选择合适的算法。
总的来说，本次实验对图像的边缘检测算法进行了全面的比较和分析，为我们
进一步深入理解和应用这些算法提供了重要的参考和指导。
希望通过这些实验
结果，我们能够更好地利用边缘检测技术来解决实际的图像处理问题，为计算
机视觉领域的发展做出更大的贡献。
边缘检测实验报告
一、实验目的通过课堂的学习，已经对图像分割的相关理论知识已经有了全面的了解，知道了许多图像分割的算法及算子，了解到不同的算子算法有着不同的优缺点，为了更好更直观地对图像分割进行深入理解，达到理论联系实际的目的，特制定如下的实验。
二、实验原理检测图像边缘信息，可以把图像看做曲面，边缘就是图像的变化最剧烈的位置。
这里所讲的边缘信息包含两个方面：一是边缘的具体位置，即像素的坐标；而是边缘的方向。
微分算子有两个重要性质：定域性(或局部性)、敏感性(或无界性)。
敏感性就是说，它对局部的函数值变化很敏感，但是因其对变化过于敏感又有了天然的缺陷——不能抵抗噪声。
局部性意思是指，每一点的导数只与函数在该点邻近的信息有关。
主要有两大类基于微分算子的边缘检测技术：一阶微分算子边缘检测与二阶微分算子边缘检测。
这些检测技术采用以下的基本步骤：(1) 将相应的微分算子简化为离散的差分格式，进而简化为模板(记为T)。
(2)利用模板对图像f(m,n)进行运算，获得模板作用后的结果Tf(m,n)。
(3) 提出阈值h,在采用一阶微分算子情形记录下高于某个阈值h 的位置坐标}),(|),{(h n m Tf n m S h ≥=(而采用二阶微分算子情形，一般是对某个阈值0>ε确立}),(|),{(ε≥=n m Tf n m S h )(4) 对集合h S 进行整理，同时调整阈值h 。
Roberts 算子Roberts 算子是一种利用局部差分算子寻找边缘的算子，两个模板分别为⎥⎦⎤⎢⎣⎡-=1001x R ⎥⎦⎤⎢⎣⎡-=0110y R 则，),(j i f R x =)1,1(),(++-j i f j i f),(j i f R y =)1,(),1(+-+j i f j i f算法的步骤为：(1) 首先用两个模板分别对图像作用得到f R x 和f R y ；(2) 对22),(y x R R j i Tf +=，进行阈值判决，若),(j i Tf 大于阈值则相应的点 位于便于边缘处。
详细的图像分割之边缘检测实验报告
边缘检测实验报告一、实验目的通过课堂的学习，已经对图像分割的相关理论知识已经有了全面的了解，知道了许多图像分割的算法及算子，了解到不同的算子算法有着不同的优缺点，为了更好更直观地对图像分割进行深入理解，达到理论联系实际的目的，特制定如下的实验。
二、实验原理：图像处理有两大类目的：1．改善像质（增强、恢复）；2．图像分析：对图像内容作出描述；其一般的图像处理过程如下：图像分割的算法有：（1）阈值分割原理：(,)(,)(,)EBL f x y T g x y L f x y T≥⎧=⎨<⎩（2）边缘检测：梯度对应一阶导数，对于一个连续图像函数f(x,y):梯度矢量定义:梯度的幅度:梯度的方向:a) Roberts 算子b) Sobel 算子Roberts 算子[]TTyxy f x f G G y x f ⎥⎦⎤⎢⎣⎡∂∂∂∂==∇),(122)()),((),(y x G G y x f mag y x f +=∇=∇)arctan(),(x y G y x =φ()()()[]()()[]{}21221,,11,1,,+-++++-=j i f j i f j i f j i f j i gc) Prewitt 算子d) Kirsch 算子由K 0~K 7八个方向模板组成，将K0~K7的模板算法分别与图像中的3×3区域乘，选最大一个值，作为中央像素的边缘强度（3）区域分割1 区域生长法 算法描述先对每个需要分割的区域找一个种子像素作为生长的起点，然后将种子像素周围邻域中与种子像素有相似性质的像素合并到种子像素所在的区域中。
将这些新像素当作新的种子像素继续进行上面的过程，直到再没有满足条件的像素可被包括进来。
2 分裂合并法实际中常先把图像分成任意大小且不重叠的区域，然后再合并或分裂这些区域以满足分割的要求,即分裂合并法.一致性测度可以选择基于灰度统计特征(如同质区域中的方差),假设阈值为T ,则算法步骤为:① 对于任一Ri ，如果 ，则将其分裂成互不重叠的四等分； ② 对相邻区域Ri 和Rj ，如果 ，则将二者合并； ③ 如果进一步的分裂或合并都不可能了，则终止算法。
图像的边缘检测实验处理报告
数字视频图像处理与通信实验实验项目：图像的边缘检测指导老师：***班级：姓名：学号：图像的边缘检测实验报告一;实验目的：1.掌握图像边缘检测的基本概念以及边缘检测的基本方法；2.通过matlab 实验的具体操作来具体掌握空间图像边缘检测的方法；3.通过matlab 实验来验证所学知识，达到学以致用；4.通过matlab 实验来理解roberts 、sobel 、canny 、log 几种算子的原理以及各个算法的优缺点，并加以比较。
二;实验原理：图像的边缘是图像最基本的特征之一。
所谓边缘（或边沿）是指周围像素灰度有阶跃性变化或“屋顶”变化的那些像素的集合。
边缘广泛存在于物体与背景之间、物体与物体之间、基元与基元之间，因此它是图像分割依赖的重要特征。
图像边缘对图像识别和计算机分析十分有用，边缘能勾划出目标物体，使观察者一目了然；边缘蕴含了丰富的内在信息（如方向、阶跃性质、形状等）。
从本质上说，图像边缘是图像局部特性不连续性（灰度突变、颜色突变、纹理结构突变等）的反应，它标志着一个区域的终结和另一个区域的开始。
边缘检测技术是所有基于边界分割的图像分析方法的第一步，首先检测出图像局部特性的不连续性，再将它们连成边界，这些边界把图像分成不同的区域，检测出边缘的图像就可以进行特征提取和形状分析，但各算子有自己的优缺点和适用领域。
Roberts 算子Roberts 算子是一种利用局部差分算子寻找边缘的算子，由下式给出： g(x,y)={[y x f ,(-)1,1(++y x f ]2+[y x f ,(- )1,1(++y x f ]2}21 ，其中f(x,y)是具有整数像素坐标的输入图像，平方根运算使该处理类似于在人类视觉系统中发生的过程。
Roberts 算子边缘定位准，但是对噪声敏感。
适用于边缘明显而且噪声较少的图像分割，在应用中经常用Roberts 算子来提取道路。
Prewitt 边缘算子Prewitt 边缘算子的卷积和如图所示，图像中的每个像素都用这两个核做卷积，取最大值作为输出，也产生一幅边缘幅度图像。
实验报告二
实验报告实验人：熊涛学号：101250182 成绩：一题目：编制一个通用的边缘提取函数。
通过输入不同的参数，能够实现Sobel 算子、Prewitt算子、Roberts算子、Marr算子和Canny边缘检测。
二实验图象：lena.bmp；三完成时间：2005.4.10四基本数学公式：1 二维离散卷积g(x,y)=T*f(x,y)=ΣΣT(i,j)f(x-i+(m-1)/2,y-j+(m-1)/2)2 一维高斯函数D(x)=exp（-x*x）五程序流程图：（附在最后，可选中后放大观看）六说明：1 程序入口函数：Sobel算子 EdgeCheck('lena.bmp','sobel');Prewitt算子 EdgeCheck ('lena.bmp','prewwit');Roberts算子 EdgeCheck ('lena.bmp','roberts');Marr算子 EdgeCheck ('lena.bmp','marr');Canny算子 EdgeCheck ('lena.bmp','canny');2 边缘检测：原始图像（图像滤波）-〉平滑图像（梯度算子）-〉梯度图像（边缘检测）-〉阈值分割（边缘定位）-〉边缘的二值图像七其他存在问题：对于这5种算子，各有优劣，应该根据不同的要求选择不同的算法主要考虑：1 假边缘概率2 丢失边缘概率3 边缘方向角估计误差4 边缘的估计值到真边缘的距离的平方均值5 畸变边缘误差范围八教师评语：N。
图像的边缘检测(实验报告)
数字信号处理实验图像的边缘检测图像的边缘检测一，原理本实验主要是对图像的边缘进行提取，通过对边缘的分析来分析图像的特征。
首先，了解一些术语的定义：边缘点：图像中具有坐标[i,j]且处在强度显著变化的位置上的点。
边缘段：对应于边缘点坐标[i,j]及其方位 ,边缘的方位可能是梯度角。
边缘检测器：从图像中提取边缘（边缘点和边缘段）集合的算法。
轮廓：边缘列表，或者是一条表示边缘列表的拟合曲线。
边缘连接：从无序边缘表形成有序边缘表的过程，习惯上，边缘表的表示采用顺时针方向来排序。
边缘跟踪：一个用来确定轮廓的图像（指滤波后的图像）搜索过程。
边缘就是图像中包含的对象的边界所对应的位置。
物体的边缘以图像局部特性的不连续性的形式出现的，例如，灰度值的突变，颜色的突变，纹理结构的突变等。
从本质上说，边缘就意味着一个区域的终结和另外一个区域的开始。
图像边缘信息在图像分析和人的视觉中十分重要，是图像识别中提取图像特征的一个重要属性。
边缘检测（edge detection）在图像处理和对象识别领域中都是一个重要的基本问题。
由于边缘的灰度不连续性，可以使用求导数的方法检测到。
最早的边缘检测方法都是基于像素的数值导数的运算。
本实验主要是对图像依次进行Sobel算子，Prewitt算子,Roberts算子，Laplace算子和Canny算子运算，比较处理结果。
边缘检测有三个共性准则，1，好的检测结果，或者说对边缘的误测率尽可能低，就是在图像边缘出现的地方检测结果中不应该没有；另一方面不要出现虚假的边缘。
2，对边缘的定位要准确，也就是我们标记出的边缘位置要和图像上真正边缘的中心位置充分接近。
3，对同一边缘要有尽可能低的响应次数，也就是检测响应最好是单像素的。
二，对图像进行各种算子运算本实验中主要是对图像依次进行Sobel算子，Prewitt算子,Roberts算子，Laplace算子和Canny 算子运算。
由于MATLAB对彩色图像不能进行分析。
图像边缘检测实验报告
图像边缘检测实验报告图像边缘检测实验报告引言：图像边缘检测是计算机视觉领域中一项重要的任务，它在许多应用中都起到关键作用。
边缘是图像中不同区域之间的分界线，它们包含了图像中物体的轮廓和形状信息。
因此，准确地检测和提取图像边缘对于目标识别、图像分割和特征提取等任务至关重要。
实验目的：本实验旨在通过实践探索和理解常用的图像边缘检测算法，并对其性能进行评估。
我们将使用不同的算法对一组测试图像进行边缘检测，并比较它们的结果，以了解它们的优缺点和适用场景。
实验方法：1. 数据准备：我们从公开的图像数据库中选择了一组具有不同特征和复杂度的测试图像。
这些图像包括自然风景、人物肖像和建筑物等多种场景，以覆盖不同的应用场景。
2. 算法选择：我们选择了三种常用的图像边缘检测算法进行实验：Sobel算子、Canny算子和Laplacian算子。
这三种算法在实践中被广泛应用，并且具有不同的特点和适用范围。
3. 实验步骤：a) Sobel算子：我们首先将测试图像转换为灰度图像，然后使用Sobel算子对其进行边缘检测。
Sobel算子是一种基于梯度的算法，它通过计算图像中每个像素点的梯度值来检测边缘。
b) Canny算子：接下来，我们使用Canny算子对同一组测试图像进行边缘检测。
Canny算子是一种基于多阶段处理的算法，它首先使用高斯滤波器对图像进行平滑处理，然后计算梯度和非最大抑制，最后进行边缘连接和阈值处理。
c) Laplacian算子：最后，我们使用Laplacian算子对测试图像进行边缘检测。
Laplacian算子是一种基于二阶导数的算法，它通过计算图像中每个像素点的二阶导数值来检测边缘。
实验结果：通过对实验图像的边缘检测，我们得到了以下结果：1. Sobel算子产生了较为明显的边缘线，但在一些复杂场景下容易产生噪声，并且边缘线有时会断裂。
2. Canny算子在平滑处理后能够准确地检测到图像中的边缘，并且能够消除噪声和断裂的边缘线。
边缘检测实验报告
一、实验目的通过课堂的学习，已经对图像分割的相关理论知识已经有了全面的了解，知道了许多图像分割的算法及算子，了解到不同的算子算法有着不同的优缺点，为了更好更直观地对图像分割进行深入理解，达到理论联系实际的目的，特制定如下的实验。
二、实验原理检测图像边缘信息，可以把图像看做曲面，边缘就是图像的变化最剧烈的位置。
这里所讲的边缘信息包含两个方面：一是边缘的具体位置，即像素的坐标；而是边缘的方向。
微分算子有两个重要性质：定域性(或局部性)、敏感性(或无界性)。
敏感性就是说，它对局部的函数值变化很敏感，但是因其对变化过于敏感又有了天然的缺陷一一不能抵抗噪声。
局部性意思是指，每一点的导数只与函数在该点邻近的信息有关。
主要有两大类基于微分算子的边缘检测技术：一阶微分算子边缘检测与二阶微分算子边缘检测。
这些检测技术采用以下的基本步骤：(1)将相应的微分算子简化为离散的差分格式，进而简化为模板(记为T)。
(2)利用模板对图像f(m,n)进行运算，获得模板作用后的结果Tf(m,n)。
(3)提出阈值h,在采用一阶微分算子情形记录下高于某个阈值h的位置坐标S h {(m, n)|Tf(m,n) h}(而采用二阶微分算子情形，一般是对某个阈值0确立S h {(m, n)|Tf(m,n) })(4) 对集合S h 进行整理，同时调整阈值h 。
Roberts 算子Roberts 算子是一种利用局部差分算子寻找边缘的算子，两个模板分别为1 0R xR y0 1 y0 11 0则，RxW, j) = f(i, j) f (i 1, j1)R y f(i,j)=f(i 1,j) f(i, j1)算法的步骤为：(1)首先用两个模板分别对图像作用得到R x f 和R y f ；(2)对Tf(i,j) J|R x |2 |R y 「，进行阈值判决，若Tf (i, j)大于阈值则相应的点 位于便于边缘处。
对于阈值选取的说明：由于微分算子的检测性能受阈值的影响较大， 为此， 针对具体图像我们采用以下阈值的选取方法，对处理后的图像统计大于某一阈 值的点，对这些数据求平均值，以下每个程序均采用此方法，不再做说明。
边缘检测实验报告
图像边缘提取实验报告一、实验目的通过课堂的学习，已经对图像分割的相关理论知识已经有了全面的了解，知道了许多图像分割的算法及算子，了解到不同的算子算法有着不同的优缺点，为了更好更直观地对图像分割进行深入理解，达到理论联系实际的目的，特制定如下的实验。
二、实验原理检测图像边缘信息，可以把图像看做曲面，边缘就是图像的变化最剧烈的位置。
这里所讲的边缘信息包含两个方面：一是边缘的具体位置，即像素的坐标；而是边缘的方向。
微分算子有两个重要性质：定域性(或局部性)、敏感性(或无界性)。
敏感性就是说，它对局部的函数值变化很敏感，但是因其对变化过于敏感又有了天然的缺陷——不能抵抗噪声。
局部性意思是指，每一点的导数只与函数在该点邻近的信息有关。
主要有两大类基于微分算子的边缘检测技术：一阶微分算子边缘检测与二阶微分算子边缘检测。
这些检测技术采用以下的基本步骤：(1) 将相应的微分算子简化为离散的差分格式，进而简化为模板(记为T)。
(2)利用模板对图像f(m,n)进行运算，获得模板作用后的结果Tf(m,n)。
(3) 提出阈值h,在采用一阶微分算子情形记录下高于某个阈值h 的位置坐标}),(|),{(h n m Tf n m S h ≥=(而采用二阶微分算子情形，一般是对某个阈值0>ε确立}),(|),{(ε≥=n m Tf n m S h )(4) 对集合h S 进行整理，同时调整阈值h 。
Roberts 算子Roberts 算子是一种利用局部差分算子寻找边缘的算子，两个模板分别为⎥⎦⎤⎢⎣⎡-=1001x R ⎥⎦⎤⎢⎣⎡-=0110y R 则，),(j i f R x =)1,1(),(++-j i f j i f),(j i f R y =)1,(),1(+-+j i f j i f算法的步骤为：(1) 首先用两个模板分别对图像作用得到f R x 和f R y ；(2) 对22),(y x R R j i Tf +=，进行阈值判决，若),(j i Tf 大于阈值则相应的点 位于便于边缘处。
边缘检测实验报告
边缘检测实验报告 Document number：NOCG-YUNOO-BUYTT-UU986-1986UT图像边缘提取实验报告一、实验目的通过课堂的学习，已经对图像分割的相关理论知识已经有了全面的了解，知道了许多图像分割的算法及算子，了解到不同的算子算法有着不同的优缺点，为了更好更直观地对图像分割进行深入理解，达到理论联系实际的目的，特制定如下的实验。
二、实验原理检测图像边缘信息，可以把图像看做曲面，边缘就是图像的变化最剧烈的位置。
这里所讲的边缘信息包含两个方面：一是边缘的具体位置，即像素的坐标；而是边缘的方向。
微分算子有两个重要性质：定域性(或局部性)、敏感性(或无界性)。
敏感性就是说，它对局部的函数值变化很敏感，但是因其对变化过于敏感又有了天然的缺陷——不能抵抗噪声。
局部性意思是指，每一点的导数只与函数在该点邻近的信息有关。
主要有两大类基于微分算子的边缘检测技术：一阶微分算子边缘检测与二阶微分算子边缘检测。
这些检测技术采用以下的基本步骤：(1) 将相应的微分算子简化为离散的差分格式，进而简化为模板(记为T)。
(2) 利用模板对图像f(m,n)进行运算，获得模板作用后的结果Tf(m,n)。
(3) 提出阈值h,在采用一阶微分算子情形记录下高于某个阈值h 的位置坐标}),(|),{(h n m Tf n m S h ≥=(而采用二阶微分算子情形，一般是对某个阈值0>ε确立}),(|),{(ε≥=n m Tf n m S h )(4) 对集合h S 进行整理，同时调整阈值h 。
Roberts 算子Roberts 算子是一种利用局部差分算子寻找边缘的算子，两个模板分别为⎥⎦⎤⎢⎣⎡-=1001x R ⎥⎦⎤⎢⎣⎡-=0110y R 则，),(j i f R x =)1,1(),(++-j i f j i f),(j i f R y =)1,(),1(+-+j i f j i f算法的步骤为：(1) 首先用两个模板分别对图像作用得到f R x 和f R y ；(2) 对22),(y x R R j i Tf +=，进行阈值判决，若),(j i Tf 大于阈值则相应的点 位于便于边缘处。
区域边缘测量实验报告
一、实验目的1. 理解区域边缘测量的基本原理和方法。
2. 掌握使用测量仪器进行实地测量的技能。
3. 学会绘制区域地形图，并准确标注边缘特征。
4. 培养团队协作和数据处理能力。
二、实验原理区域边缘测量是地理信息系统（GIS）和遥感技术中的重要环节，主要用于确定地物边界、地形变化等地理要素。
实验中，我们将使用全站仪、GPS接收机等现代测量仪器，通过实地测量和数据处理，获取区域边缘数据，并绘制成地形图。
三、实验仪器与工具1. 全站仪一台2. GPS接收机一台3. 三脚架三根4. 水准尺一把5. 记录板一块6. 计算器一台7. 地形图一张四、实验步骤1. 前期准备（1）选择实验区域，了解区域地理环境，确定测量范围和目标。
（2）查阅相关资料，了解实验区域的地形、地貌、地物分布等。
（3）制定实验方案，明确测量路线、测量方法、数据采集和处理流程。
2. 实地测量（1）使用全站仪进行水平角测量，确定测站位置。
（2）使用GPS接收机进行高程测量，获取测站点高程。
（3）根据测量数据，使用水准尺进行地物高程测量。
（4）记录测量数据，包括测站位置、地物类型、高程等。
3. 数据处理（1）使用全站仪数据处理软件，将测量数据转换为坐标系统。
（2）使用GPS数据处理软件，将GPS数据转换为坐标系统。
（3）将全站仪和GPS数据合并，生成完整的地形数据。
4. 绘制地形图（1）根据处理后的数据，绘制地形图。
（2）标注测站位置、地物类型、高程等要素。
（3）对地形图进行修饰，使其美观、易懂。
五、实验结果与分析1. 实验成功获取了实验区域的边缘数据，绘制了地形图。
2. 通过实验，掌握了区域边缘测量的基本原理和方法。
3. 提高了使用测量仪器的技能，培养了团队协作和数据处理能力。
4. 发现了实验过程中存在的问题，如测量精度、数据处理等，为今后的实验提供了参考。
六、实验总结本次区域边缘测量实验取得了预期效果，达到了实验目的。
通过实验，我们了解了区域边缘测量的基本原理和方法，掌握了使用测量仪器的技能，培养了团队协作和数据处理能力。
查找边缘实验报告
一、实验目的边缘检测是图像处理中的重要技术，它可以帮助我们提取图像中的关键特征，进而进行图像分析、图像分割、目标识别等后续处理。
本实验旨在通过学习和实践，掌握边缘检测的基本原理和常用算法，并通过对特定图像的边缘检测实验，验证算法的有效性。
二、实验原理边缘检测的基本原理是：在图像中，边缘是图像灰度值变化剧烈的地方，因此可以通过计算图像灰度值的变化率，来判断图像中是否存在边缘。
常用的边缘检测算法有Sobel算子、Prewitt算子、Roberts算子、Laplace算子等。
三、实验步骤1. 数据准备本实验选用一张自然图像作为实验数据，图像尺寸为640×480像素。
2. 边缘检测算法实现（1）Sobel算子Sobel算子是一种广泛应用于边缘检测的算法，它通过计算图像的梯度方向和大小来检测边缘。
本实验中，我们使用以下公式计算Sobel算子：Gx = Σ[-1(x[i-1][j-1]) + (-2)(x[i-1][j]) + (-1)(x[i-1][j+1]) +1(x[i][j-1]) + 2(x[i][j]) + 1(x[i][j+1]) + 1(x[i+1][j-1]) + 2(x[i+1][j]) + 1(x[i+1][j+1])]Gy = Σ[-1(x[i-1][j-1]) + 2(x[i-1][j]) + 1(x[i-1][j+1]) - 1(x[i][j-1]) + 2(x[i][j]) + 1(x[i][j+1]) - 1(x[i+1][j-1]) + 2(x[i+1][j]) +1(x[i+1][j+1])]其中，Gx和Gy分别表示图像在x和y方向的梯度。
然后，通过计算Gx和Gy的绝对值，得到图像的梯度大小：G = sqrt(Gx^2 + Gy^2)最后，将梯度大小大于阈值的部分视为边缘。
（2）Prewitt算子Prewitt算子与Sobel算子类似，也是一种计算图像梯度的方法。
边缘检测试验报告
边缘检测试验一、实验目的1．进一步理解边缘检测的基本原理。
2．掌握对图像边缘检测的基本方法。
3．学习利用Matlab图像工具箱对图像进行边缘检测。
二、实验原理图像边缘检测大幅度地减少了数据量，并且剔除了可以认为不相关的信息，保留了图像重要的结构属性。
有许多方法用于边缘检测，它们的绝大部分可以划分为两类：基于查找一类和基于零穿越的一类。
基于查找的方法通过寻找图像一阶导数中的最大和最小值来检测边界，通常是将边界定位在梯度最大的方向。
三、实验预习在实验前先编码，编码如下：close allclear allI=imread('wamp.jpg'); %读取图像I1=im2double(I); %将彩图序列变成双精度I2=rgb2gray(I1); %将彩色图变成灰色图[thr, sorh, keepapp]=ddencmp('den','wv',I2);I3=wdencmp('gbl',I2,'sym4',2,thr,sorh,keepapp); %小波除噪I4=medfilt2(I3,[9 9]); %中值滤波I5=imresize(I4,0.2,'bicubic'); %图像大小BW1=edge(I5,'sobel'); %sobel图像边缘提取BW2=edge(I5,'roberts'); %roberts图像边缘提取BW3=edge(I5,'prewitt'); %prewitt图像边缘提取BW4=edge(I5,'log'); %log图像边缘提取BW5=edge(I5,'canny'); %canny图像边缘提取h=fspecial('gaussian',5); %高斯滤波BW6=edge(I5,'zerocross',[ ],h); %zerocross图像边缘提取figure;subplot(2,3,1); %图划分为一行三幅图，第一幅图imshow(I2); %绘图subplot(2,3,2);imshow(BW1);title('Sobel算子');subplot(2,3,3);imshow(BW2);title('Roberts算子');subplot(2,3,4);imshow(BW3);title('Prewitt算子');subplot(2,3,5);imshow(BW4);title('log算子');subplot(2,3,6);imshow(BW5);title('canny算子');四、实验步骤1．打开MATLAB软件；2．利用MATLAB图像工具箱中已有函数进行图像的边缘检测；3．显示原图和处理过的图像。
图像的边缘检测实验报告
图像的边缘检测实验报告图像的边缘检测实验报告一、引言图像处理是计算机科学领域中的一个重要研究方向，而边缘检测作为图像处理的基础任务之一，具有广泛的应用价值。
边缘是图像中灰度或颜色变化较为剧烈的地方，通过检测图像中的边缘可以提取出物体的轮廓、形状等重要信息，从而为后续的图像分析和识别提供基础。
二、实验目的本次实验旨在探究不同的边缘检测算法在图像处理中的应用效果，并通过实验结果分析和比较各算法的优缺点，从而为图像处理领域的研究和应用提供参考。
三、实验方法1. 实验环境：使用Python编程语言，结合OpenCV图像处理库进行实验。
2. 实验数据：选择了包含多种物体和复杂背景的图像作为实验数据，以保证实验的可靠性和准确性。
3. 实验步骤：(1) 读取图像数据，并将其转化为灰度图像。
(2) 对图像进行预处理，如降噪、平滑等操作，以提高边缘检测的效果。
(3) 使用不同的边缘检测算法对图像进行处理，如Sobel算子、Canny算法等。
(4) 分析和比较不同算法的实验结果，评估其优缺点。
四、实验结果与分析1. Sobel算子：Sobel算子是一种基于梯度的边缘检测算法，通过对图像进行卷积操作，提取出图像中的边缘信息。
实验结果显示，Sobel算子能够较好地检测出图像中的边缘，但对于噪声较多的图像效果较差。
2. Canny算法：Canny算法是一种经典的边缘检测算法，通过多步骤的处理过程，包括高斯滤波、计算梯度、非极大值抑制和双阈值处理等，最终得到清晰准确的边缘信息。
实验结果显示，Canny算法能够有效地检测出图像中的边缘，并具有较好的抗噪性能。
3. 其他算法：除了Sobel算子和Canny算法外，还有许多其他的边缘检测算法，如拉普拉斯算子、Roberts算子等，它们各自具有不同的特点和适用范围。
在实验中，我们也对这些算法进行了尝试和比较，发现它们在不同的图像场景下有着各自的优势和局限性。
五、实验总结与展望通过本次实验，我们对图像的边缘检测算法进行了探究和比较。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
一、基本原理 (3)1.1 位场边缘识别方法研究进展概述 (3)1.2 数值计算类边缘识别方法的研究及计算 (3)1.2.1 垂向导数（VDR） (3)1.2.2 总水平导数（THDR） (3)1.2.3 解析信号振幅（ASM） (3)1.2.4 倾斜角（TA） (4)1.2.5 二阶导数类边缘识别 (4)二、输入/输出数据格式设计 (4)2.1 主要变量设计 (4)2.2位场边缘识别程序数据格式设计 (5)三、总体设计 (6)四、测试结果 (6)4.1 测试用例 (6)4.2 测试结果 (7)4.3 结果分析 (9)五、结论与建议 (10)5.1 结论 (10)5.2 建议 (10)附录：边缘识别程序源代码 (10)一、基本原理1.1 位场边缘识别方法研究进展概述地质目标体边缘是指断裂构造线、不同地质体的边界线，实际上是具有一定密度或磁性差异的地质体的边界线．在地质体的边缘附近，重、磁异常变化率较大，故所有的边缘识别方法均利用这一特点进行设计。
现有重、磁位场边缘识别方法分为数理统计、数值计算和其他三大类。
1.2 数值计算类边缘识别方法的研究及计算数值类边缘识别方法均利用极大值位置或零值位置确定地质体的边缘位置，其依据的理论基础是二度体铅垂台阶模型的重力异常特征．在该模型边缘处重力异常总水平导数和解析信号振幅达到极大值、垂向导数达到零值．故可以利用这些特征位置来确定二度体铅垂台阶的边缘位置．确定倾斜二度体、不规则二度体及三度体边缘位置的理论均为二度体铅垂台阶模型理论的推广，但确定的边缘位置与真实位置有一定偏差．该偏差随着地质体边界形状、埋深、水平尺寸及物性差异等的变化而变化．因此，边缘识别结果是一种定性或半定量解释结果，与定量解释结果有一定区别，识别结果可作为边缘位置定量反演的初值。
1.2.1 垂向导数（VDR）垂向导数方法利用零值位置确定地质体的边缘位置，重力异常可直接使用，对磁力异常必须转换成磁源重力异常（假重力异常）或化极磁力异常才可以使用。
垂向导数方法研究历史较早，方法较成熟，应用频率较高．通过我们的试验和研究认为：随着地质体埋深的增大，垂向导数所确定的边缘位置偏离实际位置的距离越大。
1.2.2 总水平导数（THDR）总水平导数是利用其极大值位置来确定地质体的边缘位置，适用于重力异常，对磁力异常必须换算成磁源重力异常或化极磁力异常才可以使用。
1.2.3 解析信号振幅（ASM）解析信号振幅也是利用极大值位置来确定地质体的边缘位置，适用于重力异常和磁力异常。
1.2.4 倾斜角（TA）倾斜角实质上是垂向导数和总水平导数的比值。
由于倾斜角为一阶导数的比值，所以能很好地平衡高幅值异常和低幅值异常，起到边缘增强的效果。
1.2.5 二阶导数类边缘识别为了增强边缘分辨能力，和克服当总水平导数等于０时倾斜角存在“解析奇点”，会使得计算结果不稳定。
二、输入/输出数据格式设计依据上述原理，现对上述各种边缘识别方法实现进行程序设计。
2.1 主要变量设计●cmd_file:参数文件名●input_field_filename：输入位场文件名●output_field_filename：输出转换后位场文件名●expan_2D_method：2D扩边方法选择●Field_Deriv_method：边缘识别方法选择●num_area：计算区域大小选择●mpoint,nline：点数，线数●m0,m1,m2,m3：扩边后X方向上各端点●n0,n1,n2,n3：扩边后Y方向上各端点●Xmin,Xmax：X坐标最小，最大值●ymin,Ymax：Y坐标最小，最大值●Ori_Field：原始场值，二维数组，m3*n3●Deriv_Field：转换后后场值，二维数组，m3*n3●Deriv_operator_x：x方向转换因子，二维数组，m3*n3●Deriv_operator_y：y方向转换因子，二维数组，m3*n3●Deriv_operator_z：z方向转换因子，二维数组，m3*n32.2位场边缘识别程序数据格式设计● cmd 文件，输入参数文件，格式如下：（cmd.txt ）● 原始位场文件，输入文件，grd 格式(ASCII) ● 延拓后位场文件，输出文件，grd 格式(ASCII) ● Grd 格式如下：input_field_filename anomaly.grdoutput_field_filename VDR_THDR_anomaly.grd num_area 0 expan_2D__method 1 Field_Deriv_method 7说明：num_area:取0：表示一般扩边区域取1：表示计算区域相对于一般扩边再放大一倍 expan_2D__method:取1：表示余弦扩边，端点值取边界平均值 取2：表示余弦扩边，端点值取全场区平均值 取3：表示余弦扩边，端点值取常数0取4：表示最小曲率扩边，空白处初值取余弦扩边值，端点初值取边界平均值取5：表示最小曲率扩边，空白处初值取余弦扩边值，端点初值取全场区平均值取6：表示最小曲率扩边，空白处初值取余弦扩边值，端点初值取常数0Field_Deriv_method:取1：计算垂向导数 VDR 取2: 计算总水平导数 THDR 取3: 计算解析信号振幅 ASM 取4: 计算倾斜角 TA取5:计算垂向二阶导数 VDR_VDR 取6:计算倾斜角总水平导数 TA_THDR 取7:计算垂向导数总水平导数 VDR_THDRDSAAmpoint nline Xmin Xmax Ymin Ymax Zmin Zmax Z 11 Z 12三、总体设计四、测试结果4.1测试用例（1）观测面上的重力异常存放在“anomaly.grd”中，坐标单位为m。
（2）形体边界存放在“rectangle.bln”中，坐标单位图 3.1 重、磁异常边缘识别程序设计流程图4.2 测试结果-400-20020040060080010000.0140.0160.0180.020.0220.0240.0260.0280.030.0320.0340.0360.0380.040.0420.0440.0460.0480.050.052-1000-800-600-400-2002004006008001000-1000-800-600-400-20020040060080010000.0060.010.0140.0180.0220.0260.030.0340.0380.0420.0460.050.0540.0580.0620.066图 4.2.3 ASM(解析信号振幅)结果图图 4.2.5 VDR_VDR(垂向二阶导)结果图-1000-800-600-400-2002004006008001000-1000-800-600-400-2002004006008001000-0.0011-0.001-0.0009-0.0008-0.0007-0.0006-0.0005-0.0004-0.0003-0.0002-0.000100.00010.00020.0003图 4.2.4 AT(倾斜角)边缘识别结果图-1000-800-600-400-2002004006008001000-1000-800-600-400-20002004006008001000-90-80-70-60-50-40-30-20-100102030405060708090图 4.2.6 VDR_THDR 结果图(垂向导数总水平导数，余弦扩边方法)-1000-800-600-400-2002004006008001000-1000-800-600-400-20020040060080010004E-0058E-0050.000120.000160.00020.000240.000280.000320.000360.00040.000444.3 结果分析1) 从图4.2.1-图4.2.7可以看出，VDR （垂向导数）、VDR_VDR(垂向二阶导)、AT(倾斜角)实验结果零值线与实际模型（白线框）基本重合，与实验原理一致，说明这几种方法程序实现是正确的。
2）同时从图4.2.1-图4.2.7可以也可看出，THDR （总水平导数）、VDR_THDR （垂向导数总水平导数）、AT_THDR(倾斜角总水平导数) 、ASM(解析信号振幅)实验结果极大值线与实际模型（白线框）基本重合，与实验原理一致，说明这几种方法程序实现也是正确的。
3）该实验模型，从以上各方法边缘识别效果来看，ASM(解析信号振幅)分辨能力较低；二阶导数类方法识别效果比一阶导数类识别效果好；就一阶导数类来看AT(倾斜角)由于其他三种方法；二阶导数类方法中AT_THDR(倾斜角总水平导数)识别最精确，但是从（图4.2.7）看出，其稳定性较差，故相比较而言VDR_THDR （垂向导数总水平导数）、VDR_VDR(垂向二阶导)效果是比较理想的。
4）对比图4.2.6和图4.2.8，可以看出在频率域处理选用扩边方法不同对计算结果有一定影响，其中最小曲率扩边计算效果为佳。
-1000-800-600-400-2002004006008001000-1000-800-600-400-200200400600800100000.150.30.450.60.750.91.051.21.351.51.651.81.952.12.252.42.552.7图 4.2.7 TA_THDR(倾斜角总水平导数)结果图-1000-800-600-400-2002004006008001000-1000-800-600-400-200020040060080010004E-0058E-0050.000120.000160.00020.000240.000280.000320.00036图 4.2.8 VDR_THDR 结果图(垂向导数总水平导数，最小曲率扩边)五、结论与建议5.1 结论此方法在一定程度上是有效的，且从实验结果可以垂向导数总水平导数和垂向二阶导识别的效果较好。
5.2 建议没有学好位场边缘识别方法没有可以能提出的建议附录：边缘识别程序源代码!*********************************************************!! PROGRAM: Fre_PoteField_Deriv! PURPOSE: Frequency Potential field Derivative! Author : gesang! Data: 2013-11!**********************************************************PROGRAM Fre_PoteField_Derivimplicit noneCHARACTER *(80) input_field_filenameCHARACTER *(80) output_field_filenameINTEGER expan_2D__method,num_area,Field_Deriv_methodINTEGER mpoint,nlineINTEGER m0,m1,m2,m3INTEGER n0,n1,n2,n3REAL Xmin,XmaxREAL ymin,YmaxREAL dx,dyREAL,ALLOCATABLE:: Ori_Field(:,:)REAL,ALLOCATABLE:: Field_Real(:,:)REAL,ALLOCATABLE:: Field_Image(:,:)REAL,ALLOCATABLE:: Deriv_Field(:,:)REAL,ALLOCATABLE:: Deriv_operator_x(:,:)REAL,ALLOCATABLE:: Deriv_operator_y(:,:)REAL,ALLOCATABLE:: Deriv_operator_z(:,:) !z方向导数因子!INPUT:!读取文件参数CALL read_cmdinfo(input_field_filename,output_field_filename,&expan_2D__method,Field_Deriv_method,num_area)!读取点数线数及测区范围CALL get_mpoint_and_nline(input_field_filename,mpoint,nline,Xmin,Xmax,& Ymin,Ymax)!计算扩边端点CALL get_expan_num(mpoint,m0,m1,m2,m3,num_area)CALL get_expan_num(nline, n0,n1,n2,n3,num_area)ALLOCATE(Ori_Field(m0:m3,n0:n3))ALLOCATE(Field_Real(m0:m3,n0:n3))ALLOCATE(Field_Image(m0:m3,n0:n3))ALLOCATE(Deriv_Field(m0:m3,n0:n3))ALLOCATE(Deriv_operator_x(m0:m3,n0:n3))ALLOCATE(Deriv_operator_y(m0:m3,n0:n3))ALLOCATE(Deriv_operator_z(m0:m3,n0:n3))!从网格化文件读入初始场值CALL Read_field(input_field_filename,m0,m1,m2,m3,n0,n1,n2,n3,Ori_Field)!PROCESS:!扩边dx=(Xmax-Xmin)/(mpoint-1)dy=(Ymax-ymin)/(nline -1)CALL expan_2D(m0,m1,m2,m3,n0,n1,n2,n3,Ori_Field,dx,dy,&expan_2D__method) !expan_2D__method:扩边方法选择!计算导数因子CALL Deriv_operator_sub(dx,dy,m3,n3,Deriv_operator_x,&Deriv_operator_y,Deriv_operator_z)!求导运算Field_Real=Ori_FieldField_Image=0.0CALL field_Deriv(Field_Real,Field_Image,Deriv_operator_x,Deriv_operator_y,& Deriv_operator_z,m0,m1,m2,m3,n0,n1,n2,n3,Field_Deriv_method) Deriv_Field=Field_Real! OUTPUT:!输出求导后位场文件CALL output_grd(output_field_filename,m0,m1,m2,m3,n0,n1,n2,n3,&mpoint,nline,Deriv_Field,Xmin,Xmax,Ymin,Ymax)DEALLOCATE (Ori_Field)DEALLOCATE (Field_Real)DEALLOCATE (Field_Image)DEALLOCATE (Deriv_Field)DEALLOCATE (Deriv_operator_x)DEALLOCATE (Deriv_operator_y)DEALLOCATE (Deriv_operator_z)END PROGRAM Fre_PoteField_Deriv!********************************开始****************! *********得到输入参数子程序开始***********SUBROUTINE read_cmdinfo(input_field_filename,output_field_filename,&expan_2D__method,Field_Deriv_method,num_area) CHARACTER *(*) input_field_filenameCHARACTER *(*) output_field_filenameINTEGER expan_2D__method,Field_Deriv_method,num_areaINTEGER unitCHARACTER*80 sCALL get_unit(unit)OPEN(unit,file='cmd.txt',status='old')READ(unit,*) s ,input_field_filenameREAD(unit,*) s ,output_field_filenameREAD(unit,*) s ,num_areaREAD(unit,*) s ,expan_2D__methodREAD(unit,*) s ,Field_Deriv_methodCLOSE(unit)END SUBROUTINE read_cmdinfo! ******得到输入参数子程序结束******!****读入grd格式文件中点数和线数及范围开始***** SUBROUTINE get_mpoint_and_nline(filename,mpoint,nline,&Xmin,Xmax,Ymin,Ymax)CHARACTER*(*) filenameINTEGER mpoint,nlineREAL Xmin,Xmax,Ymin,YmaxINTEGER unitCALL get_unit(unit)open(unit,file=filename,status='old',err=999)READ(unit,*)READ(unit,*) mpoint,nlineREAD(unit,*) Xmin,XmaxREAD(unit,*) Ymin,Ymaxclose(unit)RETURNPRINT*,'PAUSESTOPENDSUBROUTINE get_mpoint_and_nline!****读入grd格式文件中点数和线数及范围结束********!**********得到扩边端点子程序开始************ SUBROUTINE get_expan_num(mpoint,m0,m1,m2,m3,flag) IMPLICIT NONEINTEGER mtemp,factor_m,muINTEGER mpoint,m0,m1,m2,m3INTEGER flagmtemp=mpointfactor_m=1DO WHILE ((mod(mtemp,2).eq.0).and.(mtemp.ne.0)) mtemp=mtemp/2End doIF (mtemp.eq.1) THENm3=mpoint*2ELSEmu=int(log(float(mpoint))/0.693147+factor_m)m3=2**muif(flag==1) thenm3=m3*2 !计算区域为原来倍elseendifEND IFm1=1+(m3-mpoint)/2m2=m1+mpoint-1m0=1END SUBROUTINE get_expan_num!*******得到扩边端点子程序结束*************!*****读入grd格式文件中场值子程序开始********SUBROUTINE Read_field(filename,m0,m1,m2,m3,n0,n1,n2,n3,G) PARAMETER (vial=1.701411e+38)CHARACTER*(*) filenameINTEGER m0,m1,m2,m3,n0,n1,n2,n3REAL G(m0:m3,n0:n3)INTEGER unitCALL get_unit(unit)G=vialopen(unit,file=filename)DO i=1,5READ(unit,*)ENDDOREAD(unit,*) ((G(i,j),i=m1,m2),j=n1,n2)CLOSE(unit)ENDSUBROUTINE Read_field!*****读入grd格式文件中场值子程序结束********!**************扩边子程序集开始******************************** !*****扩边子程序调用主子程序开始*****SUBROUTINE expan_2D(m0,m1,m2,m3,n0,n1,n2,n3,&Ori_Field,dx,dy,expan_2D__method)INTEGER m0,m1,m2,m3INTEGER n0,n1,n2,n3REAL Ori_Field(m0:m3,n0:n3)INTEGER expan_2D__methodREAL dx,dyIF(expan_2D__method<=3) THENCALL expan_cos2(m0,m1,m2,m3,n0,n1,n2,n3,&Ori_Field,expan_2D__method)ELSECALL MinCurv(Ori_Field,m0,m1,m2,m3,n0,n1,n2,n3,&dx,dy,expan_2D__method)END IFEND SUBROUTINE expan_2D!*****扩边子程序调用主子程序结束*****!*********2D最小曲率扩边子程序开始***********SUBROUTINE MinCurv(Ori_Field,m0,m1,m2,m3,n0,n1,n2,n3,dx,dy,flag) INTEGER m0,m1,m2,m3,n0,n1,n2,n3REAL Ori_Field(m0:m3,n0:n3)REAL dx,dyINTEGER num !空白点点数INTEGER ,ALLOCATABLE::P_POINT(:),P_LINE(:)INTEGER flagREAL eigvaleigval=1.701411e+38CALL Blank_Point_number(Ori_Field,m3,n3,num,eigval)ALLOCATE(P_POINT(1:num))ALLOCATE(P_LINE(1:num))CALL Blank_Point_position(Ori_Field,m3,n3,num,eigval,P_POINT,P_LINE)CALL expan_cos2(m0,m1,m2,m3,n0,n1,n2,n3,Ori_Field,flag-3)!余弦扩边作为空白部分初值及端点值!****2D网格数据无约束最小曲率迭代CALL MinCurV_2D_net_sub(1E-4,10000,dx,dy,m3,n3,&Ori_Field,num,P_POINT,P_LINE)END SUBROUTINE MinCurv!*********2D最小曲率扩边子程序结束***********!****获得数据中的空白点数开始*****************SUBROUTINE Blank_Point_number(FIELD,mpoint,line,bpn,eigval)IMPLICIT NONEINTEGER mpoint,line,bpn,i,jREAL FIELD(1:mpoint,1:line),eigvalbpn=0DO i=1,mpoint,1DO j=1,line,1IF(FIELD(i,j)>=0.9*eigval) bpn=bpn+1ENDDOENDDOENDSUBROUTINE!****获得数据中的空白点数结束*****************!****获得数据中的空白点位置开始********************SUBROUTINEBlank_Point_position(FIELD,mpoint,line,bpn,eigval,P_POINT,P_LINE)IMPLICIT NONEINTEGER mpoint,line,bpn,i,j,kINTEGER P_POINT(1:bpn),P_LINE(1:bpn)REAL FIELD(1:mpoint,1:line),eigvalk=1DO i=1,mpoint,1DO j=1,line,1IF(FIELD(i,j)>=0.9*eigval.AND.j<bpn) THENP_POINT(k)=iP_LINE(k)=jk=k+1ENDIFENDDOENDDOENDSUBROUTINE!****获得数据中的空白点位置结束********************!****2D无约束最小曲率迭代边界处理开始******SUBROUTINE MinCurV_2D_boundary(U,mpoint,line,alph)REAL U(1-2:mpoint+2,1-2:line+2)DO j=1,line,1U(1-1,j)=2*U(1,j)-U(1+1,j)U(mpoint+1,j)=2*U(mpoint,j)-U(mpoint-1,j)ENDDODO i=1,mpoint,1U(i,1-1)=2*U(i,1)-U(i,1+1)U(i,line+1)=2*U(i,line)-U(i,mpoint-1)ENDDOU(1-1,1-1)=U(1+1,1-1)+U(1-1,1+1)-U(1+1,1+1)U(mpoint+1,1-1)=U(mpoint+1,1+1)+U(mpoint-1,1-1)-U(mpoint-1,1+1)U(1-1,line+1)=U(1+1,line+1)+U(1-1,line-1)-U(1+1,line-1)U(mpoint+1,line+1)=U(mpoint+1,line-1)+U(mpoint-1,line+1)-U(mpoint-1,line-1) i=1alph1=alph**2alph2=-2*(1+alph1)DO j=1,linepij=alph1*(U(i+1,j+1)-U(i-1,j+1)+U(i+1,j-1)-U(i-1,j-1))+alph2*(U(i+1,j)-U(i-1,j)) U(i-2,j)=U(i+2,j)+pijENDDOi=mpointDO j=1,linepij=alph1*(U(i+1,j+1)-U(i-1,j+1)+U(i+1,j-1)-U(i-1,j-1))+alph2*(U(i+1,j)-U(i-1,j)) U(i+2,j)=U(i-2,j)-pijENDDOj=1beta=1./alphbeta1=beta*betabeta2=-2*(1+beta1)DO i=1,mpoint,1qij=beta1*(U(i+1,j+1)-U(i+1,j-1)+U(i-1,j+1)-U(i-1,j-1))+beta2*(U(i,j+1)-U(i,j-1)) U(i,j-2)=U(i,j+2)+qijENDDOj=lineDO i=1,mpoint,1qij=beta1*(U(i+1,j+1)-U(i+1,j-1)+U(i-1,j+1)-U(i-1,j-1))+beta2*(U(i,j+1)-U(i,j-1)) U(i,j+2)=U(i,j-2)-qijENDDOENDSUBROUTINE MinCurV_2D_boundary!****2D无约束最小曲率迭代边界处理结束******!****2D网格数据无约束最小曲率迭代开始************** SUBROUTINE MinCurV_2D_net_sub(eps_abs,iteration_max,dx,dy,&mpoint,line,FIELD,number_eigval,M_eigval,N_eigval) REAL eps_abs,dx,dyINTEGER iteration_max,mpoint,line,number_eigvalREAL FIELD(1:mpoint,1:line)INTEGER M_eigval(1:number_eigval),N_eigval(1:number_eigval) REAL,ALLOCATABLE::U(:,:)ALLOCATE(U(1-2:mpoint+2,1-2:line+2))DO j=1,line,1DO i=1,mpoint,1U(i,j)=FIELD(i,j)ENDDOENDDOeps_error=2*ABS(eps_abs)iteration=0alph=dx/dyalph0=-1./(2*(3+4*alph**2+3*alph**2))alph1=alph**4alph2=2*alph**2alph3=-4*(1+alph**2)alph4=-4*(1+alph**2)*alph**2DO WHILE((eps_error>eps_abs).and.(iteration<iteration_max))eps_error=0.CALL MinCurV_2D_boundary(U,mpoint,line,alph)DO k=1,number_eigval,1i=M_eigval(k)j=N_eigval(k)temp=(U(i+2,j)+U(i-2,j))+alph1*(U(i,j+2)+U(i,j-2))+ &alph2*(U(i+1,j+1)+U(i-1,j+1)+U(i+1,j-1)+U(i-1,j-1))+ &alph3*(U(i+1,j)+U(i-1,j))+alph4*(U(i,j+1)+U(i,j-1)) temp=temp*alph0eps_error=MAX(ABS(temp-U(i,j)),eps_error)U(i,j)=tempENDDOiteration=iteration+1ENDDOprint*,iteration,eps_errorDO j=1,line,1DO i=1,mpoint,1FIELD(i,j)=U(i,j)ENDDOENDDOdeALLOCATE(U,STA T=ierr)ENDSUBROUTINE!**2D网格数据无约束最小曲率迭代结束********!***************2D余弦扩边开始*******SUBROUTINE expan_cos2(m1,m2,m3,m4,n1,n2,n3,n4,G,flag) PARAMETER (pi=3.141592654)INTEGER flagINTEGER m1,m2,m3,m4INTEGER n1,n2,n3,n4REAL G(m1:m4,n1:n4)REAL G1(m1:m4,n1:n4)REAL G2(m1:m4,n1:n4)REAL sumINTEGER numsum=0.0num=0IF(flag==1) then!扩边端点值取边界平均值DO i=m2,m3sum=sum+G(i,n2)+G(i,n3)num=num+2END DODO i=n2+1,n3-1sum=sum+G(m2,i)+G(m3,i)num=num+2END DODO i=m1,m4G(i,n1)=sum/numG(i,n4)=sum/numEND DODO i=n1+1,n4-1G(m1,i)=sum/numG(m4,i)=sum/numEND DOELSE IF(flag==2) then!扩边端点值取全部数据平均值DO i=m2,m3DO j=n2,n3sum=sum+G(i,j)num=num+1END DOEND DODO i=m1,m4G(i,n1)=sum/numG(i,n4)=sum/numEND DODO i=n1+1,n4-1G(m1,i)=sum/numG(m4,i)=sum/numEND DOELSE IF(flag==3) then!扩边端点值取常数，DO i=m1,m4G(i,n1)=0G(i,n4)=0END DODO i=n1+1,n4-1G(m1,i)=0G(m4,i)=0END DOEND IF!扩边G1=GG2=GDO j=n2,n3DO i=m1+1,m2-1G1(i,j)=(G1(m2,j)-G1(m1,j))*cos((m2-i)*1.0/(m2-m1)*pi/2.0)+G1(m1,j) END DODO i=m3+1,m4-1G1(i,j)=(G1(m4,j)-G1(m3,j))*cos((m4-i)*1.0/(m4-m3)*pi/2.0)+G1(m3,j) END DOEND DODO i=m1,m4DO j=n1+1,n2-1G1(i,j)=(G1(i,n2)-G1(i,n1))*cos((n2-j)*1.0/(n2-n1)*pi/2.0)+G1(i,n1) END DODO j=n3+1,n4-1G1(i,j)=(G1(i,n4)-G1(i,n3))*cos((n4-j)*1.0/(n4-n3)*pi/2.0)+G1(i,n3) END DOEND DODO i=m2,m3DO j=n1+1,n2-1G2(i,j)=(G2(i,n2)-G2(i,n1))*cos((n2-j)*1.0/(n2-n1)*pi/2.0)+G2(i,n1) END DODO j=n3+1,n4-1G2(i,j)=(G2(i,n4)-G2(i,n3))*cos((n4-j)*1.0/(n4-n3)*pi/2.0)+G2(i,n3) END DOEND DODO j=n1,n4DO i=m1+1,m2-1G2(i,j)=(G2(m2,j)-G2(m1,j))*cos((m2-i)*1.0/(m2-m1)*pi/2.0)+G2(m1,j) END DODO i=m3+1,m4-1G2(i,j)=(G2(m4,j)-G2(m3,j))*cos((m4-i)*1.0/(m4-m3)*pi/2.0)+G2(m3,j) END DOEND DOG=(G1+G2)/2.0END SUBROUTINE expan_cos2!***************2D余弦扩边开始*******!**************扩边子程序集结束*******************************!************计算导数因子开始*****************SUBROUTINEDeriv_operator_sub(dx,dy,m,n,Deriv_operator_x,Deriv_operator_y,Deriv_operator_z) IMPLICIT NONEREAL dx,dyREAL dzINTEGER m,n,i,jREAL W(1:m,1:n),U(1:m),V(1:n)REAL Deriv_operator_Z(1:m,1:n)REAL Deriv_operator_x(1:m,1:n)REAL Deriv_operator_y(1:m,1:n)CALL WA VE2D_sub(U,V,W,m,n,dx,dy)!计算垂向导数因子,实数Deriv_operator_z=W!计算X和y导数因子，虚数do j=1,ndo i=1,mDeriv_operator_x(i,j)=U(i)Deriv_operator_y(i,j)=V(j)end doend doEND SUBROUTINE Deriv_operator_subSUBROUTINE WA VE2D_sub(U,V,W,m,n,dx,dy)REAL W(1:m,1:n),U(1:m),V(1:n)INTEGER m,nREAL dx,dyREAL i,jCALL Wave1D_sub(n,dy,V)CALL Wave1D_sub(m,dx,U)DO j=1,n !计算园频率DO i=1,mW(i,j)=SQRT(U(i)*U(i)+V(j)*V(j))END DOEND DOEND SUBROUTINE WA VE2D_subSubroutine Wave1D_sub(N,dy,V)REAL dyINTEGER NREAL V(0:N-1)REAL deltfdeltf=(2.0*3.141592654)/(N*dy)Do j=0,N/2,1V(j)=j*deltfEnd doDo j=N/2+1,N-1,1V(j)=(j-n)*deltfEnd doEND subroutine Wave1D_sub!************计算导数因子子程序结束*****************!******************** 求导子程序集开始**************************! **************求导子程序总模块开始**********SUBROUTINEfield_Deriv(Field_Real,Field_Image,Deriv_operator_x,Deriv_operator_y,Deriv_operator_z,& m0,m1,m2,m3,n0,n1,n2,n3,Field_Deriv_method)IMPLICIT NONEINTEGER m0,m1,m2,m3,n0,n1,n2,n3REAL Field_Real(m0:m3,n0:n3),Field_Image(m0:m3,n0:n3)REALDeriv_operator_x(m0:m3,n0:n3),Deriv_operator_y(m0:m3,n0:n3),Deriv_operator_z(m0:m3,n 0:n3)INTEGER Field_Deriv_methodIF(Field_Deriv_method==1) THEN!计算垂向导数CALL VDR(Field_Real,Field_Image,m3,n3,Deriv_operator_z) ELSE IF(Field_Deriv_method==2) THEN!计算总水平导数callTHDR(Field_Real,Field_Image,m3,n3,Deriv_operator_x,Deriv_operator_y)ELSE IF(Field_Deriv_method==3) THEN!计算解析信号振幅callASM_meth(Field_Real,Field_Image,m3,n3,Deriv_operator_x,Deriv_operator_y,Deriv_operat or_z)ELSE IF(Field_Deriv_method==4) THEN!计算倾斜角CALLTA(Field_Real,Field_Image,m3,n3,Deriv_operator_x,Deriv_operator_y,Deriv_operator_z) ELSE IF(Field_Deriv_method==5) THEN!计算垂向二阶导数call VDR_VDR(Field_Real,Field_Image,m3,n3,Deriv_operator_z) ELSE IF(Field_Deriv_method==6) THEN!计算倾斜角总水平导数CALLTA_THDR(Field_Real,Field_Image,m3,n3,Deriv_operator_x,Deriv_operator_y,Deriv_operat or_z)ELSE!计算垂向导数总水平导数callVDR_THDR(Field_Real,Field_Image,m3,n3,Deriv_operator_x,Deriv_operator_y,Deriv_operator_z)END IFEND SUBROUTINE field_Deriv! ************延拓子程序总模块结束**********!***************VDR*******************SUBROUTINE VDR(Field_Real,Field_Image,m3,n3,Deriv_operator_z) INTEGER m3,n3REAL Field_Real(1:m3,1:n3),Field_Image(1:m3,1:n3)REAL Deriv_operator_z(1:m3,1:n3)CALL FFT2(Field_Real,Field_Image,m3,n3,2)CALL multiply_sub(Field_Real,Field_Image,m3,n3,Deriv_operator_z)CALL FFT2(Field_Real,Field_Image,m3,n3,1)END SUBROUTINE VDR!*************VDR******************!************ 相乘运算开始**********SUBROUTINE multiply_sub(Field_Real,Field_Image,m3,n3,Deriv_operator) IMPLICIT NONEINTEGER m3,n3,i,jREAL Field_Real(1:m3,1:n3),Field_Image(1:m3,1:n3)REAL Deriv_operator(1:m3,1:n3)DO i=1,m3DO j=1,n3Field_Real(i,j) =Field_Real(i,j) *Deriv_operator(i,j)Field_Image(i,j)=Field_Image(i,j)*Deriv_operator(i,j)END DOEND DOEND SUBROUTINE multiply_sub!********* 相乘运算结束*************!***************THDR**********************SUBROUTINE THDR(Field_Real,Field_Image,m3,n3,&Deriv_operator_x,Deriv_operator_y)INTEGER m3,n3,I,JREAL Deriv_operator_X(1:m3,1:n3),Deriv_operator_y(1:m3,1:n3)REAL Field_Real(1:m3,1:n3)REAL Field_Image(1:m3,1:n3)REAL Field_Real_x(1:m3,1:n3),Field_Image_x(1:m3,1:n3)REAL Field_Real_y(1:m3,1:n3),Field_Image_y(1:m3,1:n3)CALL FFT2(Field_Real,Field_Image,m3,n3,2)do i=1,m3do j=1,n3Field_Real_x(i,j)=-Deriv_operator_X(i,j)*Field_Image(i,j)Field_Image_x(i,j)=Deriv_operator_X(i,j)*Field_Real(i,j)Field_Real_y(i,j)=-Deriv_operator_y(i,j)*Field_Image(i,j)Field_Image_y(i,j)=Deriv_operator_y(i,j)*Field_Real(i,j)end doend doCALL FFT2(Field_Real_x,Field_Image_x,m3,n3,1)CALL FFT2(Field_Real_y,Field_Image_y,m3,n3,1)do i=1,m3do j=1,n3Field_Real(i,j)=sqrt(Field_Real_x(i,j)**2+Field_Real_y(i,j)**2) end doend doEND SUBROUTINE THDR!*****************THDR*********************!************************ASM_meth******************************** SUBROUTINE ASM_meth(Field_Real,Field_Image,m3,n3,Deriv_operator_x,&Deriv_operator_y,Deriv_operator_z)INTEGER m3,n3REAL Field_Real(1:m3,1:n3),Field_Image(1:m3,1:n3)REAL Field_Real_THDR(1:m3,1:n3),Field_Image_THDR(1:m3,1:n3)REAL Field_Real_VDR(1:m3,1:n3),Field_Image_VDR(1:m3,1:n3)REAL Deriv_operator_z(1:m3,1:n3)REAL Deriv_operator_X(1:m3,1:n3),Deriv_operator_y(1:m3,1:n3)Field_Real_THDR=Field_RealField_Image_THDR=Field_ImageCALL THDR(Field_Real_THDR,Field_Image_THDR,m3,n3,&Deriv_operator_x,Deriv_operator_y)Field_Real_VDR=Field_RealField_Image_VDR=Field_ImageCALL VDR(Field_Real_VDR,Field_Image_VDR,m3,n3,Deriv_operator_z)DO i=1,m3do j=1,n3Field_Real(i,j)=sqrt(Field_Real_THDR(i,j)**2+Field_Real_VDR(i,j)**2) end doend doEND SUBROUTINE ASM_meth!**************************ASM_meth********************************** !************************TA****************************************** SUBROUTINE TA(Field_Real,Field_Image,m3,n3,Deriv_operator_x,&Deriv_operator_y,Deriv_operator_z)INTEGER m3,n3REAL Field_Real(1:m3,1:n3),Field_Image(1:m3,1:n3)REAL Deriv_operator_z(1:m3,1:n3)REAL Deriv_operator_X(1:m3,1:n3),Deriv_operator_y(1:m3,1:n3)REAL Field_Real_THDR(1:m3,1:n3),Field_Image_THDR(1:m3,1:n3)REAL Field_Real_VDR(1:m3,1:n3),Field_Image_VDR(1:m3,1:n3)Field_Real_THDR=Field_RealField_Image_THDR=Field_ImageCALL THDR(Field_Real_THDR,Field_Image_THDR,m3,n3,&Deriv_operator_x,Deriv_operator_y)Field_Real_VDR=Field_RealField_Image_VDR=Field_ImageCALL VDR(Field_Real_VDR,Field_Image_VDR,m3,n3,Deriv_operator_z)DO i=1,m3do j=1,n3Field_Real(i,j)=atand(Field_Real_VDR(i,j)/Field_Real_THDR(i,j)) end doend doEND SUBROUTINE TA!************************TA****************************************** !**********************VDR_VDR*************************************** SUBROUTINE VDR_VDR(Field_Real,Field_Image,m3,n3,Deriv_operator_z) INTEGER m3,n3REAL Field_Real(1:m3,1:n3),Field_Image(1:m3,1:n3)REAL Deriv_operator_z(1:m3,1:n3)CALL FFT2(Field_Real,Field_Image,m3,n3,2)CALL multiply_sub(Field_Real,Field_Image,m3,n3,Deriv_operator_z)CALL multiply_sub(Field_Real,Field_Image,m3,n3,Deriv_operator_z)CALL FFT2(Field_Real,Field_Image,m3,n3,1)END SUBROUTINE VDR_VDR!*********************VDR_VDR**************************************** !********************TA_THDR*************************************** SUBROUTINE TA_THDR(Field_Real,Field_Image,m3,n3,Deriv_operator_x,&Deriv_operator_y,Deriv_operator_z)INTEGER m3,n3REAL Field_Real(1:m3,1:n3),Field_Image(1:m3,1:n3)REAL Deriv_operator_z(1:m3,1:n3)REAL Deriv_operator_X(1:m3,1:n3),Deriv_operator_y(1:m3,1:n3)REAL Field_Real_THDR(1:m3,1:n3),Field_Image_THDR(1:m3,1:n3)REAL Field_Real_TA(1:m3,1:n3),Field_Image_TA(1:m3,1:n3)Field_Real_TA=Field_RealField_Image_TA=Field_ImageCALLTA(Field_Real_TA,Field_Image_TA,m3,n3,Deriv_operator_x,Deriv_operator_y,Deriv_operat or_z)Field_Real_THDR=Field_Real_TAField_Image_THDR=Field_ImageCALLTHDR(Field_Real_THDR,Field_Image_THDR,m3,n3,Deriv_operator_x,Deriv_operator_y) Field_Real=Field_Real_THDREND SUBROUTINE TA_THDR!*******************TA_THDR***************************************** !********************VDR_THDR*************************************** SUBROUTINEVDR_THDR(Field_Real,Field_Image,m3,n3,Deriv_operator_x,Deriv_operator_y,Deriv_operator_z)INTEGER m3,n3REAL Field_Real(1:m3,1:n3),Field_Image(1:m3,1:n3)REAL Deriv_operator_z(1:m3,1:n3)REAL Deriv_operator_X(1:m3,1:n3),Deriv_operator_y(1:m3,1:n3)REAL Field_Real_THDR(1:m3,1:n3),Field_Image_THDR(1:m3,1:n3)REAL Field_Real_VDR(1:m3,1:n3),Field_Image_VDR(1:m3,1:n3)Field_Real_VDR=Field_RealField_Image_VDR=Field_ImageCALL VDR(Field_Real_VDR,Field_Image_VDR,m3,n3,Deriv_operator_z)Field_Real_THDR=Field_Real_VDRField_Image_THDR=Field_ImageCALLTHDR(Field_Real_THDR,Field_Image_THDR,m3,n3,Deriv_operator_x,Deriv_operator_y) Field_Real=Field_Real_THDREND SUBROUTINE VDR_THDR!*******************VDR_THDR***************************************** !*******************************************************************!**********************离散Fourier变换子程序集开始*********!c******************************************************************c !c c !c 2-D FAST FOURIER TRANSFORM FOR COMPLEX FUNCTIONc!c------------------------------------------------------------------c!c DREAL: The real part of the function to be transformed!c DIMAGE: The imaginary part of the function to be transformed!c M: The number of points. M must be the nu-th power of 2!c N: The number of points. N must be the nu-th power of 2!c NF: (=1,reverse transform!c (=2, normal transform!c!c Author: gesang!c TREAL -- Real part of the working array!c TIMAG -- Imaginary part of the working array!c------------------------------------------------------------------cSUBROUTINE FFT2(DREAL,DIMAG,M,N,NF)real dreal(m,n),dimag(m,n)real,ALLOCATABLE::treal(:),timag(:)nmmax=max(m,n)allocate(treal(nmmax),timag(nmmax),STA T=ierr)DO i=1,mIF (n.ne.1) THENdo j=1,ntreal(j)=dreal(i,j)timag(j)=dimag(i,j)end docall fft(treal,timag,n,nf)do j=1,ndreal(i,j)=treal(j)dimag(i,j)=timag(j)end doEND IFEND DODO j=1,nIF(m.ne.1) THENdo i=1,mtreal(i)=dreal(i,j)timag(i)=dimag(i,j)end docall fft(treal,timag,m,nf)do i=1,mdreal(i,j)=treal(i)dimag(i,j)=timag(i)end doEND IFEND DOdeallocate(treal,timag,STAT=ierr)END SUBROUTINE!c******************************************************************c !c c !c 1-D FAST FOURIER TRANSFORM FOR COMPLEX FUNCTION c !c------------------------------------------------------------------c!c XREAL -- Real part of the function to be transformed c!c XIMAG -- Imaginary part of the function to be transformed c!c N ------ The number of points. N must be the nu-th power of 2 c!c NF ----- (=1,reverse transform c!c (=2, normal transform c!c (0, 1, ......, n/2-1, n/2, -n/2+1, ......,-1) c!c------------------------------------------------------------------cSUBROUTINE FFT(XREAL,XIMAG,N,NF)real xreal(n),ximag(n)nu=int(log(float(n))/0.693147+0.001)n2=n/2nu1=nu-1f=float((-1)**nf)k=0DO l=1,nuDO while (k.lt.n)do i=1,n2p=ibitr(k/2**nu1,nu)arg=6.2831853*p*f/float(n)c=cos(arg)s=sin(arg)k1=k+1k1n2=k1+n2treal=xreal(k1n2)*c+ximag(k1n2)*stimag=ximag(k1n2)*c-xreal(k1n2)*sxreal(k1n2)=xreal(k1)-trealximag(k1n2)=ximag(k1)-timagxreal(k1)=xreal(k1)+trealximag(k1)=ximag(k1)+timagk=k+1end dok=k+n2END DOk=0nu1=nu1-1n2=n2/2END DODO k=1,ni=ibitr(k-1,nu)+1if(i.gt.k) thentreal=xreal(k)timag=ximag(k)xreal(k)=xreal(i)ximag(k)=ximag(i)xreal(i)=trealximag(i)=timagend ifEND DOIF(nf.ne.1) THENdo i=1,nxreal(i)=xreal(i)/float(n)ximag(i)=ximag(i)/float(n)end doEND IFEND SUBROUTINEFUNCTION IBITR(J,NU)j1=jitt=0do i=1,nuj2=j1/2itt=itt*2+(j1-2*j2)ibitr=ittj1=j2end doend function!******* 离散Fourier变换子程序集结束**********!*******************求导子程序集结束**************************。