montajMagmapFiltering.pdf

montaj MAGMAP Filtering 2D Frequency Domain Processing of Potential Field Data Extension for Oasis montaj 6.4

TUTORIAL

www.geosoft.com

The software described in this manual is furnished under license and may only be used or copied in accordance with the terms of the license. Manual release date: 07/09/2007. Written by, Nancy Whitehead and Chris Musselman. Please send comments or questions to [email protected] Copyright © Geosoft Inc. 2007. All rights reserved. No part of this publication may be reproduced, stored in a retrieval system or transmitted, in any form, or by any means, electronic, mechanical, photo-copying, reading, or otherwise, without prior consent from Geosoft Inc. Program Copyright© Geosoft Inc. 2007. All rights reserved. Geosoft and Oasis montaj are registered trademarks of Geosoft Inc. GEOSOFT, Oasis are trademarks of Geosoft Inc. Windows®, and Windows NT are either registered trademarks or trademarks of Microsoft Corporation.

Geosoft Incorporated 8th Floor 85 Richmond St. W. Toronto, Ontario M5H 2C9 Canada Tel: (416) 369-0111 Fax: (416) 369-9599 Web Site: www.geosoft.com E-mail: [email protected]

Contents Geosoft License Agreement

1

Finding Help Information

3

Contacting Technical Support

3

Chapter 1: System Capabilities and Concepts

5

Who Should Use This Document

5

How this Document is Organized

6

MAGMAP Menu and Processing System

6

Navigating the MAGMAP Menu

6

Understanding Processing Sequence

7

Preparing Grids

8

Specifying Filter

8

Filtering and Post-processing

9

Chapter 2: Before You Begin

10

Creating a Project

10

Loading the MAGMAP Menu

11

Chapter 3: montaj MAGMAP Filtering Tutorials

13

Tutorial 1: Preparing Data for Processing

14

Tutorial 2: Displaying Grids Prior to Processing and Analysis

14

Tutorial 3: Preparing Grids for Transformation

15

Tutorial 4: Applying the Forward FFT

17

Tutorial 5: Setting Filters

17

Tutorial 6: Applying Filters and the Inverse Transform

20

Tutorial 7: Displaying Filtered Grids for Analysis

21

Tutorial 8: Calculating Radially Averaged Power Spectra

22

Tutorial 9: Displaying Radially Averaged Spectra

23

Tutorial 10: Calculating and Displaying 2D Power Spectra

24

Tutorial 11: Interactive Spectrum Filtering

26

Tutorial 12: Calculating Analytic Signal

34

Tutorial 13: Calculating Tilt Derivative

37

Tutorial 14: Sharing Results

40

Chapter 4: 2D Fast Fourier Transform and MAGMAP Theory 42 The FFT Algorithm

45

The Energy Spectrum

47

Chapter 5: Preparing Grids for the FFT Processing

49

Trend Removal Algorithm

49

Grid Expansion Algorithm

50

Grid Filling Algorithm

50

Minimizing "Ringing" from Grid Filling

50

Applying Maximum Entropy Prediction

51

Controlling Edge Effects

51

Limiting Strong Anomalies Near Grid Edges

51

Limiting Strong Anomalies by Magnitude

52

Setting Trend Removal, Grid Expansion, and Filling Parameters

Chapter 6: Specifying FFT Filters

52

54

Creating Your Own Filter Control File

54

MAGMAP Filters

55

BPAS – Bandpass Filter

56

BTWR – Butterworth Filter

57

CNDN – Downward Continuation

58

CNUP – Upward Continuation

58

COSN – Cosine Roll-off Filter

59

DCOS – Directional Cosine Filter

60

DENS – Apparent Density Calculation

61

DPAS – Directional Pass/Reject Filter

62

DRVX – Derivative in the X Direction

63

DRVY – Derivative in the Y Direction

63

DRVZ – Derivative in the Z Direction

63

GFILT – Gravity Earth Filter

64

GAUS – Gaussian Regional/Residual Filter

64

GNRL – General Radially Symmetric Filter

65

GPSD – Pseudo-Gravity Filter

66

HPAS – High-pass Filter

67

INTG – Vertical Integration

67

LPAS – Low-pass Filter

68

OPTM – Weiner Optimum Filter

68

REDE – Reduce to the Magnetic Equator

70

REDP – Reduce to the Magnetic Pole

70

SUSC – Apparent Susceptibility Calculation

71

TXYZ – Conversion between Field Components

72

Filter Examples

73

Susceptibility Map

73

Second Vertical Derivative

74

De-corrugation

Chapter 7: Applying the Inverse FFT

74

75

Processing Option

75

Selection

75

Result

75

References

76

1

Geosoft License Agreement GEOSOFT agrees to supply the Licensed Program(s) as specified in my purchase order. Geosoft shall grant me a nontransferable, non-exclusive license to use the Licensed Program(s), subject to the Terms and Conditions herein contained. Should there be a separate signed agreement between you and Geosoft, or between your company and Geosoft, pertaining to the licensed use of this software, that agreement shall take precedence over the terms of this agreement. 1.

DEFINITIONS:

In this Agreement: "Licensed Program(s)" means the actual copy of all or any portion of Geosoft’s proprietary software technology, computer software code, components, dynamic link libraries (DLLs) licensed through the Geosoft license server, including any modifications, improvements or updates provided by GEOSOFT. “Effective Date” is the date the Geosoft license is installed. This date is recorded by the Geosoft License server when the Licensed Program(s) is installed. "Services" means the Services described on Section 4. "Termination" means the occurrences contemplated by Section 6 and 7. 2.

LICENSE:

GEOSOFT grants to me a non-transferable and non-exclusive license to use the Licensed Program(s) for my own purposes whereby the Licensed Program(s) are being used only by myself, on one computer, at any one time. Title and all intellectual property rights in and to the License Program(s), including, without limitation, copyright, trade secrets and trade marks, shall remain with GEOSOFT. I agree to refrain from raising any objection or challenge to such intellectual property rights, or from assisting or causing or permitting other(s) to do so, during the term of the Agreement and thereafter I may not assign this Agreement or any part thereof or sub-license the rights granted herein, or lend, rent, time-share, sell or lease the software without the prior written consent of GEOSOFT. I may not attempt to reverse engineer, de-compile or disassemble the software. I may not make any attempt to circumvent the License Manager that controls the access to the software use. 3.

TERM:

The Term of this Agreement shall commence on the Effective Date and shall continue until termination, as described in Section 6. 4.

SERVICES:

(i) According to the terms of my initial purchase, GEOSOFT shall make available to me, without additional fees such corrections and improvements to the Licensed Program(s) as may be generally incorporated into the Licensed Program(s) by GEOSOFT. (Normally this will be for a period of twelve (12) months). (ii) GEOSOFT has a strong commitment to customer service and product support. GEOSOFT offers me, subject to applicable Service Charge(s), continuing support in the form of email or telephone advice and other assistance in problem diagnosis and the correction of errors or faults in the Licensed Program(s) during the life of this License. When a problem occurs which appears to be related to errors or faults in the Licensed Program(s), I may contact GEOSOFT and GEOSOFT will make an honest effort to solve the problem. However, GEOSOFT cannot guarantee service results or represent or warrant that all errors or program defects will be corrected. Also it is to be noted that each Licensed Program is designed to operate on a Windows NT (sp 6 or later), Windows 2000 or Windows XP platform. (iii) Further, if I request service relating to the modification of the Licensed Program(s) to meet a particular need or to conform with a particular operating environment, GEOSOFT may, at its discretion, modify the Licensed Program(s) to meet these particular needs, subject to applicable Services Charge(s). However, all intellectual property or other rights which may arise from such modifications shall reside with GEOSOFT. 5.

PROTECTION AND SECURITY OF LICENSED PROGRAM

I agree that all additions, modifications, revisions, updates and extensions to the Licensed Program(s) shall be subject to all of the terms and conditions in this agreement. I acknowledge that all copies of the Licensed Program(s), provided by GEOSOFT or made by me pursuant to this Agreement, including, without limitation, translations, compilations, partial copies, modifications, derivative materials and/or updated materials, are proprietary, and the property of GEOSOFT, and may not be distributed to any other persons, without GEOSOFT’s prior written consent. I will not provide or otherwise make the Licensed Program(s) available to anyone in any form without GEOSOFT's prior written consent.

2 6.

TERMINATION:

This agreement shall terminate upon the termination date, if any, specified in your purchase agreement with Geosoft. This agreement may be terminated only upon thirty-days prior written notice to GEOSOFT. GEOSOFT may terminate this Agreement upon prior written notice effective immediately if I fail to comply with any of the terms and conditions of this Agreement. This Agreement shall terminate automatically upon the institution, or consenting to the institution of proceedings in insolvency or bankruptcy, or upon a trustee in bankruptcy or receiver being appointed for me/us for all or a substantial portion of my/our assets. 7.

EVENTS UPON TERMINATION:

I shall forthwith discontinue use of the Licensed Program(s), on the day Termination shall occur and agree not to resume such use in the future without written authorization from GEOSOFT. I shall uninstall and remove all software from my computer. Within thirty days after Termination, I shall destroy all physical and digital copies of the Licensed Program(s). This obligation relates, without limitation, to all copies in any form, including translations, compilations, derivatives and updated materials, whether partial or complete, and whether or not modified or merged into other materials as authorized herein. 8.

WARRANTY:

GEOSOFT does not warrant that the functions contained in the Licensed Program will meet my requirements or will operate in the combinations which may be selected for use by me, or that the operation of the Licensed Program will be uninterrupted or error free or that all program defects will be corrected. Each Licensed Program shall be furnished to me in accordance with the terms of this Agreement. No warranties, either express or implied, are made to me regarding the Licensed Program. THE FOREGOING WARRANTIES ARE IN LIEU OF ALL OTHER WARRANTIES, EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OR MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. 9.

LIMITATION OF REMEDIES

I agree to accept responsibility for the use of the programs to achieve my intended results, and for the results obtained from use of said Program(s). I therefore accept complete responsibility for any decision made based on my use of the aforementioned Licensed Program(s). In no event shall GEOSOFT be liable for any damages arising from performance or non-performance of the Licensed Program(s), or for any lost profits, lost savings or other consequential damages, even if GEOSOFT has been advised of the possibility of such damages, or for any claim against me by any other party. 10. GENERAL: I agree that this Agreement is a complete and exclusive statement of the agreement with GEOSOFT. This Agreement supersedes all previous Agreements with respect to the Licensed Programs, with the exception of a current signed Technical Service Agreements. GEOSOFT is not responsible for failure to fulfill its obligations under the Agreement due to causes beyond its control. Should any part of This Agreement for any reason be declared invalid, such declaration shall not affect the remaining portion which shall remain in full force and effect as if this Agreement had been executed without the invalid portion thereof. The relationship between the parties is that of independent contractors. Nothing contained in this Agreement shall be deemed to constitute or create a partnership, association, joint venture or agency. The provision of this Agreement shall be binding upon me and GEOSOFT and my respective successors and permitted assigns. This Agreement will be governed by the laws of the Province of Ontario and applicable laws of Canada. 11. YEAR 2000: The Licensed Programs have been tested to conform to DISC PD2000 1:1998 Year 2000 Conformity Requirements (www.bsi.org.uk/disc/year2000/2000.html), with the exception of clause 3.3.2, paragraph b. Section 3.3.2 paragraph b) requires that inferences for two-digit year dates greater than or equal to 50 imply 19xx, and those with a value equal to or less than 50 imply 20xx. The Licensed Programs will recognize all two digit years as 19xx. This is to prevent errors importing historical data that pre-dates 1950. All dates that follow 1999 must use four digit dates in the Licensed Programs.

3

Finding Help Information There are several functions included in the basic Oasis montaj help system that may be useful to your work. The entire documentation for the system is available through the online help system. This electronic library of information enables us to constantly update the information and provide you with the most up-to-date information available. The best way to find information in this system is to use the Search tab to perform a full-text search across all help topics. If you still cannot find the information you are looking for, the Online Books help system contains complete Geosoft manuals and tutorials in Adobe PDF format.

Contacting Technical Support The following list provides contact information for Geosoft Technical Support around the world. North America

Europe and North Africa

Geosoft Inc., 85 Richmond St. W., 8th Floor Toronto, Ont., Canada M5H 2C9

Geosoft Europe Ltd. 20/21 Market Place, First Floor Wallingford, Oxfordshire United Kingdom OX10 OAD

Tel +1 (416) 369-0111 Fax +1 (416) 369-9599

Tel: +44 1491 835 231 Fax: +44 1491 835 281

Email: [email protected]


South America

Australia and Southeast Asia

Geosoft Latinoamerica Ltda. Praça Floriano 51 / 19º Andar CEP: 20031-050, Centro Rio de Janeiro, RJ, Brasil

Geosoft Australia Pty. Ltd 14 Prowse Street, West Perth WA 6005 Australia

Tel: 55 (21) 2111-8150 Fax: 55 (21) 2111-8181

Tel +61 (8) 9382 1900 Fax +61 (8) 9382 1911



South and Central Africa Geosoft Africa Ltd. Buren Building, Second Floor Kasteelpark Office Park c/o Nossob & Jochemus Streets Erasmuskloof X3, Pretoria Tel: +27 12 347 4519 Fax: +27 12 347 6936 Email: [email protected]

Chapter 1: System Capabilities and Concepts 5

Chapter 1: System Capabilities and Concepts This document describes the montaj MAGMAP Filtering system, an integrated software package that supports application of common Fourier domain filters to gridded data – with particular emphasis on potential field data. For example, if you want to produce an apparent magnetic susceptibility grid from a total magnetic field grid, you can choose and apply the appropriate filter in MAGMAP. The montaj MAGMAP Filtering system provides three filtering methods designed to help you meet your data filtering requirements: • • •

MAGMAP 1-step filtering Step-By-Step filtering Interactive filtering

The Interactive filtering method enables you to select one of the spectral filters (Bandpass, Butterworth, Cosine Roll-off, and Gaussian Regional/Residual) and interactively select the parameters (using sliders), while seeing the effect on the current power spectrum. The system’s 2D Fast Fourier Transform capabilities enable you to: • •

• • •

Rapidly process gridded datasets by applying a wide range of robust geophysical and mathematical filters. Apply multiple filters in one filtering run (in any order). For example, you can run any combination of geophysical filters (such as the Upward Continuation filter) and/or mathematical filters (such as the Butterworth filter) using a single dialog in the system GUI. Control the filtering process by applying either fast (one-step), expanded (stepby-step), or visual (intractive) filtering. Define and apply your own filters. Interpret grids using spectral analysis products (2D and Radially Averaged Power Spectra).

Who Should Use This Document This document is intended for Earth Scientists who are familiar with methodologies used for acquiring, processing, and presenting Earth Science data, and who want to use the montaj MAGMAP Filtering extension of Oasis montaj to process, analyze, visualize, and interpret potential field data. To use this document effectively, you should: • •

Be familiar with Earth Science data Understand basic methodologies used for processing data, including preparing data for processing, evaluating quality before and after processing, displaying

6 Chapter 1: System Capabilities and Concepts

•

data in profile formats, running computer-based algorithms, and preparing data for final presentation. Understand basic database concepts, such as importing and exporting data, storing data, and applying processes to data.

How this Document is Organized This document includes the following chapters: • • •

•

•

•

•

Chapter 1: System Capabilities and Concepts – Explains how to use this document, and the MAGMAP menu and processing system. Chapter 2: Before You Begin – Describes the procedures used for creating a project and loading the MAGMAP menu in Oasis montaj. Chapter 3: montaj MAGMAP Filtering Tutorials – Provides a set of tutorials that guide you through the system procedures used for filtering, analysing, and visualizing potential field data. Chapter 4: Fast Fourier Transform and MAGMAP Theory – Introduces the theory on which MAGMAP is based and the methodology employed in the system. Chapter 5: Preparing Grids for the FFT Processing – Provides a detailed description of the preparation algorithms (trend removal, grid expansion, and filling) that are implemented in the system. Explains how to set the parameters that control these processes. Tells how to control side effects, such as “ringing”, that can affect the quality of your results. Chapter 6: Specifying FFT Filters – Outlines the three methods of creating a filter control file, and defines the role of the filter control file in the system. Provides a comprehensive description of each filter provided in the system, along with a discussion of the filter’s uses and required parameters. Chapter 7: Applying the Inverse FFT – Explains the options available while applying the Inverse FFT.

The References section directs you to further sources of information related to the FFT theory and interpretation.

MAGMAP Menu and Processing System The montaj MAGMAP Filtering system is designed to provide you with an intuitive interface that guides you through the preparation, processing, analysis, and visualization of potential field data. This section describes the application menu and outlines the major steps in the processing sequence. Navigating the MAGMAP Menu The montaj MAGMAP Filtering menu provides access to the MAGMAP functions. When you select a menu option, the system runs a corresponding Geosoft eXecutable

Chapter 1: System Capabilities and Concepts 7

(GX) – a programmed process that records your input values and implements a specific processing, analysis, or visualization task. The following shows the MAGMAP menu options and the corresponding GXs:

MAGMAP1 GX

FFT2PREP GX FFT2IN GX FFT2CON GX FFT2FLT GX

FFT2PREP GX FFT2IN GX FFT2RSPC GX FFT2SPCFLT GX FFT2FLT GX

FFT2RSPC GX FFT2SMAP GX FFT2PSPC GX

GRIDASIG GX

TILTDRV GX

Understanding Processing Sequence For the sake of mathematical convenience and speed, montaj MAGMAP Filtering applies filters in the wavenumber or Fourier domain. This requires a number of steps, each of which is the responsibility of a separate program in the montaj MAGMAP system. Pre-processing steps involve preparation of the original space domain grid for filtering, after which filters are applied. Post-processing involves returning the filtered data to the same size and shape as the original grid, and replacement of a regional trend. Following are descriptions of each of the processing steps. Note:

You can carry out each of these steps individually, using the MAGMAP menus. The MAGMAP One-step filtering menu option can take you through the entire sequence in one step. This approach will produce adequate results in most situations. However, if you experience ringing problems, edge effects, or any other undesirable side effects that appear to be caused by the pre-

8 Chapter 1: System Capabilities and Concepts

processing steps, you must run the pre-processing steps separately and set parameters to address this issue. Preparing Grids 1. A first-order trend is removed from the data. This is not always necessary, but is

recommended. 2. The grid is expanded to be square, with dimensions that are acceptable to the Fast

Fourier Transform (FFT) used in the MAGMAP system. The system pads the edges of the grid with dummy values. 3. The dummy areas are replaced by reasonably interpolated values so that the grid

becomes smoothly periodic. If you think of the grid as a single square tile where copies of the tiles are laid out edge to edge, the grid pattern should smoothly flow from tile to tile (see sample below).

If the matches are not smooth, an effective 'step' function is introduced at the edge of each grid. This can cause serious side effects in the data when filters are applied in the wavenumber domain. 4. The square and periodic space domain grid is transformed to the wavenumber

domain by the application of FFT. A radially averaged spectrum of the data is also produced for reference and analysis. Specifying Filter 1. The required filters are specified for the wavenumber grid.

Chapter 1: System Capabilities and Concepts 9 2. Using the Interactive filtering menu option, filter parameters can be visualized

and selected interactively. Filtering and Post-processing 1. The selected filters are applied to the wavenumber grid. 2. The filtered wavenumber grid is transformed from the wavenumber domain back

to the space domain. 3. The dummy areas of the original grid are restored to the final filtered grid, and the

grid size is reduced to its original size. 4. If a trend surface has been removed, and if no high-pass filters have been applied

to the data, the filtered trend can be replaced in the grid.

10 Chapter 2: Before You Begin

Chapter 2: Before You Begin This chapter describes how to begin working with the montaj MAGMAP Filtering system in Oasis montaj. The topics discussed in this chapter include: • •

Creating a project Loading the MAGMAP menu

This tutorial uses sample data provided on the Oasis montaj CD and installed in your C:\Program Files\Geosoft\Oasis montaj\data\magmap directory. Before you begin the tutorial, you need to create a working directory to store all your data. The system enables you to access files anywhere. However, it is a good strategy to carefully organize your data (project information and files) before carrying out any processing. To start this tutorial, create a working directory called C:\Tutorial. A general rule to follow when working with Geosoft applications is to avoid working in the Geosoft directory. In these tutorials, you follow this rule by keeping all the working data, found in C:\Program Files\Geosoft\Oasis montaj\data\magmap, in your working directory C:\Tutorial.

Creating a Project Work in Oasis montaj requires an open project. An Oasis montaj "project" encompasses every item in your working directory: the data files in your project (databases, maps, and grids), tools used (including auxiliary tools such as histograms, scatter plots, etc.), and the project setup including the menus you have loaded, map or profile as a processed entity, and the state in which you left this entity the last time you used it. The project also controls your working directory. Projects are saved as (*.gpf) files. If you open an existing project from a directory, the system assumes that all your project files are located in the same directory. To streamline your work, as well as to keep it organized, make sure that your project file is in the same directory as the other files you want to use. We recommend that each project you work on have its own project (*.gpf) file. If you use a number of applications or add-on tools in Oasis montaj that have different menus, you can use the project to display only the menus you require. The Project Explorer tool enables you to browse and open project items. The Project Explorer pane has two tabbed sections. The Data section displays all data files included in the project, and the Tools section organizes and maintains the project tools. To access the Tools section, you click the Tools bar at the bottom of the Project Explorer pane. To return to the Data pane, you click the Data bar at the top of the Project Explorer pane.

Chapter 2: Before You Begin 11 TO

CREATE A PROJECT :

1. Start Oasis montaj. 2. From the File menu, select Project > New.

The New Project dialog is displayed.

Oasis montaj assumes that your data are in the directory containing this project (C:\Tutorial). 3. Specify a name and directory for the project. For example, name the project Magmap and place it in the working directory C:\Tutorial. Note:

4. Click [Save].

The system saves the project and indicates that it is open by adding menus to the menu bar, adding buttons to the toolbar, and by displaying the Project Explorer pane. These are visual clues indicating that you are ready to start working with the system.

Loading the MAGMAP Menu Before you can start working with the montaj MAGMAP Filtering system, you have to load the MAGMAP menu in your project. For more information on setting menus, refer to the Oasis montaj Online Help system (Help > Help Topics).

12 Chapter 2: Before You Begin TO

LOAD THE

MAGMAP

MENU :

1. From the GX menu, select Load Menu or click the Load Menu icon (

) on the

toolbar. The Load Menu dialog is displayed. 2. From the list of files, select “magmap.omn” and click the [Open] button.

The system displays the MAGMAP menu on the main Oasis montaj menu bar.

Chapter 3: montaj MAGMAP Filtering Tutorials 13

Chapter 3: montaj MAGMAP Filtering Tutorials The Fast Fourier Transform (FFT) filtering methods have a wide variety of applications in the Earth Sciences investigations. Typically applied to potential field data derived from geophysical surveys, the FFT filters can be used to remove geologic and cultural noise (i.e., to "clean" data), perform regional/residual separations for interpretation purposes, and to estimate certain physical parameters, such as magnetic susceptibility. In addition, FFT filtering enables geoscientists to evaluate and interpret frequency-dependent relationships in the transformed data via power spectra and other forms of advanced analysis. In this tutorial, you will use Geosoft 2D-FFT system to apply the FFT algorithm and accomplish a subset of these tasks – first filtering gridded data to remove geologic noise, and then extracting vertical and horizontal derivatives for comparison with the original grid. You will also learn how to create maps, visualize grids on the screen, and print maps for interpretation. montaj MAGMAP Filtering enables you to apply a "long" (extended) FFT process, a one-step FFT process, or an interactive FFT process. The long process enables you to control each part of the sequence (prepare grids, apply forward FFT, set filters, and apply inverse FFT), whereas the one-step process performs the grid preparation, forward and inverse FFT steps for you. The interactive process enables you to visualize the filtering parameters and interactively select those parameters that best apply to your data. For completeness, the tutorial provides an example of the expanded (step-by-step) process, including a description of how to calculate power spectra, and an example of the interactive method. To help you learn how to use the MAGMAP system, we provide the following tutorials: • • • • • • • • • • • • • •

Tutorial 1: Preparing Data for Processing….page 14 Tutorial 2: Displaying Grids Prior to Processing and Analysis….page 14 Tutorial 3: Preparing Grids for Transformation….page 15 Tutorial 4: Applying the Forward Fast Fourier Transform….page 17 Tutorial 5: Setting Filters….page 17 Tutorial 6: Applying Filters and the Inverse Transform….page 20 Tutorial 7: Displaying Filtered Grids….page 21 Tutorial 8: Calculating Radially Averaged Power Spectra….page 22 Tutorial 9: Displaying Radially Averaged Power Spectra ….page 23 Tutorial 10: Calculating and Displaying 2D Power Spectra….page 24 Tutorial 11: Performing Interactive Spectrum Filtering….page 26 Tutorial 12: Calculating Analytic Signal….page 34 Tutorial 13: Calculating Tilt Derivative….page 37 Tutorial 14: Sharing Results ….page 40

14 Chapter 3: montaj MAGMAP Filtering Tutorials

Tutorial 1: Preparing Data for Processing The montaj MAGMAP Filtering system is designed to work on processed, gridded data. This tutorial does not explicitly describe the procedure used for preparing your data. In real life, you may be required to prepare the data yourself prior to using the montaj MAGMAP Filtering system. For more information about additional tools you may require for filtering, gridding, and performing other basic processing tasks, contact your Geosoft representative.

Tutorial 2: Displaying Grids Prior to Processing and Analysis Before starting, you may want to display the sample grid we have provided for processing (mag_in.grd). TO

DISPLAY A GRID :

1. From the Grid menu, select Display Grid > Colour-Shaded Grid.

The Color-shaded grid image dialog is displayed.

2. Using the […] button, select the Grid name as “mag_in.grd”. 3. Accept the default values for the rest of the parameters, and click the [New Map]

button.


The system creates a colour-shaded image from the mag_in.grd file, and places that image in a new map window called mag_in.map.

The above grid, especially in its central and upper parts, displays “noisy” magnetic field, which is typically produced by shallow sources. The goal of the filtering process for this grid is minimizing the effects of these shallow magnetic sources to enhance the signature of deeper objects.

Tutorial 3: Preparing Grids for Transformation Before transforming a grid to the wavenumber domain (applying the forward FFT), the grid file so that the grid: • • •

Has dimensions that are acceptable to the FFT. Includes no dummy values. Is periodic on its edges. In other words, a point on the left edge of the grid must match the corresponding point on the right edge, and a point on the top edge must match the corresponding point on the bottom edge, in both value and slope (derivative)

When you are using the extended FFT process in the montaj MAGMAP Filtering system, the first processing step is to prepare your grid using the Prepare Grids menu option. When you are using the one-step FFT process, grid preparation is handled automatically.


Basic grid preparation steps include trend removal, grid expansion (to make grid smoothly periodic), and grid filling (to replace all the dummy values with interpolated values). In this tutorial, we only provide a functional (how-to) description of the grid preparation process. If you would like to learn more about trend removal, grid expansions, and filling algorithms, refer to Chapter 6: Preparing Grids for the FFT Processing. This chapter provides a complete description of the grid preparation process, and tells you how to set trend removal, grid filling, and expansion parameters. TO

PREPARE A GRID FOR THE

FFT

TRANSFORMATION :

1. From the MAGMAP menu, select Step-By-Step Filtering > Prepare Grid.

The FFT2 grid pre-processing dialog is displayed.

2. Using the […] button, enter the Name of Input (Original) Grid File as

“mag_in.grd”. 3. Using the […] button, enter a name for the pre-processed grid in the Name of

Output (Pre-processed) Grid File field as “prep_in”. 4. Leaving the intelligent default values for the rest of the parameters, as shown

above, click the [Start] button.


The system prepares the grid and displays it in your current project as a temporary map file.

The above grid has new dimensions that are acceptable to FFT, includes NO dummy values, and is periodic on its edges.

Tutorial 4: Applying the Forward FFT After the grid preparation you apply the forward FFT. TO

APPLY THE FORWARD

FFT:

1. From the MAGMAP menu, select Step-By-Step Filtering > Forward FFT.

The FFT2IN dialog is displayed.

2. Using the […] button, enter the Name of Input Pre-processed Grid File as

“prep_in.grd”. Note that this should be the default value. 3. Click the [OK] button.

The system computes the transformation.

Tutorial 5: Setting Filters After computing the transform file, you are ready to apply FFT filters. You can apply up to six filters in any order. Following is the list of available filters:

18 Chapter 3: montaj MAGMAP Filtering Tutorials BPAS

Bandpass Filter

GAUS

Gaussian regional/residual filter

BTWR

Butterworth Filter

GNRL

General radially symmetric filter

CNDN

Downward continuation

GPSD

Pseudo-Gravity filter

CNUP

Upward continuation

HPAS

High-pass filter

COSN

Cosine roll-off filter

INTG

Vertical integration

DCOS

Directional cosine filter

LPAS

Low-pass filter

DENS

Apparent density calculation

OPTM

Weiner optimum filter

DPAS

Directional pass/reject filter

REDE

Reduce to the magnetic equator

DRVX

Derivative in the X direction

REDP

Reduce to the magnetic pole

DRVY

Derivative in the Y direction

SUSC

Apparent susceptibility calculation

DRVZ

Derivative in the Z direction

TXYZ

Conversion between Field Components

GFILT

Gravity Earth Filter

For the purpose of this tutorial, we will apply the Butterworth and Upward Continuation filters. For a complete description of each of these filters, refer to Chapter 6: Specifying FFT Filters. TO

SET

FFT

FILTERS :

1. From the MAGMAP menu, select Step-By-Step Filtering > Define Filters.

The MAGMAP FILTER DESIGN dialog is displayed.

2. You can specify a new control file name. For this tutorial, accept the default

control file name “magmap.con”. Click the [OK] button.


The MAGMAP FILTER DESIGN dialog is displayed.

3. From the First filter to apply drop-down list, select “Butterworth Filter”. 4. From the 2nd filter (optional) drop-down list, “Upward Continuation”. 5. Click the [OK] button.

The Butterworth Filter (Filter1) dialog is displayed.

Tip:

The Butterworth filter is excellent for applying straightforward high-pass or low-pass filters to data because you can easily control the degree of filter roll-off while leaving the central wavenumber fixed. If ringing is observed, the degree can be reduced until the results are acceptable.

6. Enter the Cutoff Wavenumber (in ground unit) as “4”. 7. Enter the Filter Order as “8”. 8. From the High or Low Pass drop-down list, select “Low-Pass”. 9. Click the [OK] button.

The Upward Continuation dialog is displayed. Note:

The dialogs that appear here correspond to the type of filters you selected to apply. If you only selected one filter, the system displays a single dialog for


that filter. If you selected more than one filter, multiple dialogs are displayed. You must set the parameters for all the selected filters.

Upward continuation is considered a “clean” filter because it produces almost no side effects that may require application of other filters or processes to correct. Because of this, this filter is often used to remove or minimize the effects of shallow sources and noise in grids. Also, upward continued data may be interpreted numerically, with modeling programs. This is not the case for many other filter processes. 10. Type “200” in the Distance to upward continue (in ground units) field. Tip:

11. Click the [OK] button.

The system sets the filter parameters. You are now ready to apply the filter and inverse transform to obtain an output grid.

Tutorial 6: Applying Filters and the Inverse Transform At this point, you are ready to apply the filters and the inverse FFT process. There are several processing options that you can choose. For description of these options, see the Applying the Inverse FFT chapter, page 75. In this tutorial, you select a filtered grid with post-processing. This option applies the grid logic to restore the original grid dimensions, and also replaces the trend that was removed in the initial grid preparation stage. T O A PPLY

FILTERS AND INVERSE

FFT:

1. From the MAGMAP menu, select Step-By-Step Filtering > Apply Filter.

The FFT2FLT dialog is displayed.

Chapter 3: montaj MAGMAP Filtering Tutorials 21 2. Using the […] button, enter “prep_in_trn.grd” as a value for the Name of Input

Transform (*_trn.grd) File field. 3. Using the […] button, enter “mag_out” as a value for the Name of Output Grid

File field. 4. Using the […] button, enter “magmap.con” as a value for the Name of Filter

Control File field. 5. Using the […] button, enter “mag_in.grd” as a value for the Name of Reference

(Original) Grid File field. 6. Click the [OK] button to select the post-processing option.

The system computes a post-processed, filtered, space-domain grid, and displays it in a temporary map.

Tutorial 7: Displaying Filtered Grids for Analysis After generating the filtered grid, you may want to display it on your screen for comparison with the original total field grid. You have a choice of displaying grids either as standard or shaded grid images. TO

DISPLAY GRIDS :

1. Display the original mag_in and the newly created mag_out as colour-shaded

grids on temporary maps.

22 Chapter 3: montaj MAGMAP Filtering Tutorials 2. Use the right-click pop-up menu to zoom and pan the maps until you are satisfied

with the size and location of each map. At this point, your maps should look similar to the following:

Tutorial 8: Calculating Radially Averaged Power Spectra If required, you can calculate a radially averaged version of your spectrum. TO

CALCUALATE THE RADIALLY AVERAGED SPECTRUM :

1. From the MAGMAP menu, select Spectrum Calculation and Display >Radial

average spectrum. The FFTRSPC dialog is displayed.

2. Using the [Browse] button, enter the Name of Input Transform (grid) File as

“prep_in_trn.grd”. 3. Using the [Browse] button, enter the Name of Output Spectrum File as “radial”. 4. Click the [OK] button.


The system computes the spectrum.

Tutorial 9: Displaying Radially Averaged Spectra You can view the spectra calculated in the previous step in your project. The spectra are automatically formatted to fit on a page when printed. To examine a sample plot, see energy spectrum section in Chapter 4: 2D Fast Fourier Transform and MAGMAP Theory page 42. TO

DISPLAY RADIALLY AVERAGED SPECTRUM :

1. From the MAGMAP menu, select Spectrum Calculation and Display >Display

Spectrum. The Create a spectrum map dialog is displayed.

2. Using the […] button, enter the Input spectrum file name as “radial.SPC”. 3.

Click the [OK] button.


The radial.map file is displayed.

Tutorial 10: Calculating and Displaying 2D Power Spectra If you are performing a more complex interpretation, you may be interested in computing and displaying the 2D power spectrum for your grid. TO

CALCULATE AND DISPLAY

2D

POWER SPECTRA;

1. From the MAGMAP menu, select Spectrum Calculation and Display > 2-D

Power Spectrum.


The FFTPSPC dialog is displayed.

2. Using the […] button, enter the Name of Input Transform (grid) File as

“prep_in_trn.grd”. 3. Using the […] button, enter the Name of Output 2D Spectrum (grid) File as

“FFT_power”. 4. Click the [OK] button.

The system computes the spectrum grid and displays it in a temporary map file.


Tutorial 11: Interactive Spectrum Filtering The interactive spectrum filtering process enables you to display a radially averaged power spectrum of an FFT2 Transform image, a selected filter profile, and the resultant (filtered) power spectrum profile, in an interactive window. The filters can be selected from a drop-down list or from a predefined control file. Filter parameters can be interactively modified to obtain the best results. The Interactive Spectral Filter dialog includes a profile window, a filter selection drop-down list, and filter-specific parameter fields. The profile pane displays: • • •

The radially averaged power spectrum profile from the original input spectrum file, in black The filter profile in blue The resultant (filtered) profile in red

The values of the filter parameter fields can be entered interactively, by moving the associated slider, or by typing in the field boxes. Interactive spectrum filtering can be used when working with the following filters: Bandpass, Butterworth, Cosine Roll-off, Gaussian Regional/Residual, Upward Continuation, Downward Continuation, Derivative in the Z Direction, and Vertical Integration. For more information on these filters, see Chapter 6: Specifying FFT Filters page 54. The Interactive Filtering process includes the following procedures: • • • • •

Preparing grids for transformation Applying the forward Fast Fourier Transform Calculating the radial average spectrum Defining interactive spectrum filters Appying filters and the inverse transform TO

PREPARE A GRID FOR

FFT

TRANSFORMATION :

1. From the MAGMAP menu, select Interactive Filtering > Prepare Grid.


The FFT2 grid pre-processing dialog is displayed.

2. Using the […] button, enter the Name of Input (Original) Grid File as

“mag_in.grd”. 3.

Using the […] button, enter the Name of Output (Pre-processed) Grid File as “Interactive_prep_grid”.

4. Leaving the the intelligent default values for the rest of the parameters, as shown

above, click the [Start] button.


The system prepares the grid and displays it in your current project in a temporary map file.

TO

APPLY THE FORWARD

F AST F OURIER T RANSFORMATION :

1. From the MAGMAP menu, select Interactive Filtering > Forward FFT.

The FFT2IN dialog is displayed.

2. Using the […] button, enter the Name of Input Pre-processed Grid File as

“Interactive_prep_grid.grd”. This should be the default value. 3. Click the [OK] button.

The system computes the transform. TO

CALCULATE THE RADIALLY AVERAGED SPECTRUM :

1. From the MAGMAP menu, select Interactive Filtering > Radial Average

Spectrum.


The FFTRSPC dialog is displayed.

2. Using the […] button, enter the Name of Input Transform (grid) File as

“Interactive_prep_grid_trn.grd”. 3.

Using the […] button, enter the Name of Output Spectrum File as “Interactive_Radial”.


The system computes the spectrum. TO

SPECIFY INTERACTIVE SPECTRUM FILTERS :

1. From the MAGMAP menu, select Interactive Filtering > Interactive Spectrum

Filters. The Interactive FFT2 radially averaged power spectrum filter dialog is displayed.

2. Using the […] button, enter the Spectrum File Name as “Interactive_Radial.SPC”. 3.

Using the […] button, enter the Control File Name as “Interactive_magmap.con”.



The Interactive Spectral Filter window is displayed.

5. Select the spectral filter from the Filter Name drop-down list.


The selected filter parameters are displayed in the Filters section. For example, for the Butterworth filter, the Filters section looks as follows:

6. Modify the filter parameters by moving the parameter sliders or by typing values

in fields and then pressing [Enter] on the keyboard. The filter profile and filtered spectrum profile are updated and re-displayed accordingly. 7. Optionally, to preview the impact of filtering on your grid, click the [Preview]

button.


The Preview section is appended at the right-hand side of the window.

The Original Grid image shows the unfiltered grid. The Filtered Grid image shows the grid filtered with the currently selected filter parameters. The Filtered Grid image changes in real time when you change the filter parameters in the Filters section. 8. When satisfied with the filtered spectrum, click the [OK] button to save the

selected filter along with its parameter values to the output control file (Interactive_magmap.con). TO

APPLY FILTERS AND INVERSE

FFT:

1. From the MAGMAP menu, select Interactive Filtering > Apply Filter.


The FFT2FLT dialog is displayed.

2. Using the […] button, enter the Name of Input Transform (*_trn.grd) File as

“Interactive_prep_grid_trn.grd”. 3. Using the […] button, enter the Name of Output Grid File as

“Interactive_mag_out”. 4. Using the […] button, enter the Name of Filter Control file as

“Interactive_magmap.con”. 5. Using the […] button, enter the Name of Reference (Original) Grid File as

“mag_in.grd”. Note: Depending on the post-processing option you would like to apply to your data, you can click the [Flt-Inv Only], [Filter Only], or the [OK] button. For more information on these options, see the Applying the Inverse FFT chapter. 6. Click the [OK] button to choose the post-processing option.


The system computes a post-processed, filtered, space-domain grid, and displays that grid in your project.

Tutorial 12: Calculating Analytic Signal Analytic signal is the square root of the sum of squares of the derivatives in the x, y, and z directions: asig = sqrt(dx*dx + dy*dy + dz*dz)

The analytic signal is useful in locating the edges of magnetic source bodies, particularly where remanence and/or low magnetic latitude complicate interpretation. The default Z-derivative method is FFT. However, for very large grids (over 4000 x 4000 cells), using the Convolution method saves a lot of processing memory and time. TO

CALCULATE THE ANALYTIC SIGNAL FOR A GRID :

1. From the MAGMAP menu, select Analytic Signal.


The Calculate the analytic signal for a grid dialog is displayed.

2. Use the […] button next to the Input grid field to select an input grid. 3. In the Output analytic signal grid field, enter an output grid name. 4. From the Z-derivative method drop-down list, select “FFT” or “Convolution”. 5. From the Retain derivative grids drop-down list, select:

• Yes – To keep the intermediate derivative grids (dx.grd, dy.grd, and dz.grd) • No – To delete the intermediate derivative grids on exit 6. Click [OK].


The analytic signal grid is created at the defined location. The following screenshots show the input grid and the output grid, respectively.


Tutorial 13: Calculating Tilt Derivative The Tilt Derivative option calculates the tilt derivative of a grid and, optionally, the total horizontal derivative of the tilt derivative grid. The tilt derivative and its total horizontal derivative are used for mapping shallow basement structures and mineral exploration targets. The default Z-derivative method is FFT. However, for very large grids (over 4000 x 4000 cells), using the Convolution method saves a lot of processing memory and time. The tilt derivative is defined as: TDR = arctan(VDR/THDR)

where VDR and THDR are the first vertical and total horizontal derivatives, respectively, of the total magnetic intensity T VDR = dT/dz THDR = sqrt((dT/dx)^2 + (dT/dy)^2)


The total horizontal derivative of the tilt derivative is defined as: HD_TDR = sqrt((dTDR/dx)^2 + (dTDR/dy)^2) TO

CALCULATE TILT DERIVATIVE FOR A GRID :

1. From the MAGMAP menu, select Analytic Signal.

The Calculate a tilt derivative grid dialog is displayed.

2. Use the […] button next to the Input grid field to select an input grid. 3. In the Output tilt derivative grid (TDR) field, enter an output grid name. 4. Optionally, to create the horizontal derivative grid, in the Output horiz. deriv. of

TDR grid (HD_TDR) field, enter a grid name. If you leave this field blank, the HD TDR grid will not be created. 5. From the Z-derivative method drop-down list, select “FFT” or “Convolution”. 6. Click [OK].


The tilt derivative grid (and, optionally, horizontal derivative grid) is (are) created at the defined location(s), and is (are) displayed in the workspace.


Tutorial 14: Sharing Results You can share your maps as: • •

Soft copies – By exporting the maps to a variety of formats Hard copies – By plotting the maps

Geosoft supports the following export formats for maps: • • • • • • • • • • • • •

Enhanced Metafile (*.emf) ArcView Shapefile (*.shp) CGM Plot (*.cgm) Bitmap (*.bmp) PCX (*pcx) J2K JPEG 2000 (*.j2k) PNG (*.png) TIFF Compressed (*.tif) Encapsulated Post Script (*.eps) Geosoft COLOR Grid (*.grd) ER Mapper RGB (*.ers) ER Mapper ECW Compressor (*.ecw) SVG Export (*.svg)

• • • • • • • • • • • •

Geosoft Plot (*.plt) DXF AutoCAD (*.dxf) DXF AutoCAD v12 (*.dxf) JPEG (*.jpg) JPEG High Quality (*.jpg) JP2 JPEG 2000 (*.jp2) TIFF (*.tif) TIFF JPEG 2000 Compressed (*.tif) GeoTIFF (*.tif) MapInfo TIFF (*.tif) ArcView TIFF (*.tif) MapInfo TAB (*.tab)


For information about exporting data, refer to the Exporting and Archiving topic in the Oasis montaj Online Help system, or to the relevant chapter of the Oasis montaj Tutorial. Oasis montaj uses your installed Windows system drivers to create printer or plotter output. For information about plotting maps, refer to the Oasis montaj Online Help system, or to the relevant chapter of the Oasis montaj Tutorial.

42 Chapter 4: 2D Fast Fourier Transform and MAGMAP Theory

Chapter 4: 2D Fast Fourier Transform and MAGMAP Theory For mathematical convenience, MAGMAP applies filters in the Fourier, or wavenumber, domain. This document assumes that you are familiar with the application of filters to two dimensional data using the Fourier Domain techniques. This chapter provides a short review of some of the basic concepts of the Fourier domain processing. Mathematically, the Fourier transform of a space domain function f(x,y) is defined as: f ( μ, ν ) =

∞ ∞

∫ ∫ f ( x, y ) ⋅ e

− i ( μx +νy )

dxdy

−∞ −∞

The reciprocal relation is: f ( x, y ) =

∞ ∞

1 4π 2

∫ ∫ f ( μ,ν ) ⋅ e

i ( μx + νy )

dμ dν

−∞−∞

where μ and ν are wavenumbers in the x and y directions, respectively, measured in radians per metre if x and y are in given metres. These are related to spatial "frequencies" fx and fy, in cycles per metre. A grid (in the space domain) is transformed to and from the wavenumber domain using Fast Fourier Transform (FFT). The equivalent data set in the wavenumber domain is commonly called “Transform”. A Transform of a grid is composed of wavenumbers, which have units of cycles/metre, and have a real and imaginary component. Just as a grid samples a space domain function at even distance increments, the Transform samples the Fourier domain function at even increments of 1/(grid size) (cycles/metre) between 0 and the Nyquist wavenumber (1/[2*cell size]). A given potential field function in the space domain has a single and unique wavenumber domain function, and vice versa. The addition of two functions (anomalies) in the space domain is equivalent to the addition of their Transforms. The energy spectrum is a 2D function of the energy relative to wavenumber and direction. The radially averaged energy spectrum is a function of wavenumber alone, and is calculated by averaging the energy for all directions for the same wavenumber. The Fourier transform of the potential field produced by a prismatic body has a broad spectrum whose peak location is a function of the depth to the prism’s top and bottom surfaces, and whose magnitude is determined by the prism’s density or magnetization. The peak wavenumber (ω') can be determined by the following expression:

Chapter 4: 2D Fast Fourier Transform and MAGMAP Theory 43

ln (hb ht ) hb − ht

ω′ =

where:

ω' is the peak wavenumber in (radians / metre) ht is the depth to the top hb is the depth to the bottom The spectrum of a bottomless prism peaks at the zero wavenumber according to the following expression (Bhatacharia, 1966): f ( μ , ν ) = e − hr r=

μ2 + ν 2

where h is the depth to the top of the prism. The spectrum for a prism with top and bottom surfaces is: f ( μ ,ν ) = e − ht r − e − hb r

where ht and hb are the depths to the top and bottom surfaces, respectively. As the prism bottom is brought up, the peak moves to higher wavenumbers as illustrated in the following figure. 1 no bottom

top = 4

bottom depth 36 20 12 8

0 0

wavenumber

1


Considering the spectrum of a fixed-size prism, as the prism depth is increased, the peak of the spectrum is shifted to lower wavenumbers (the anomaly becomes broader), and the magnitude of the spectrum is reduced. 1 thickness = 4

top = 4

8 16

0

0

wavenumber

1

An important fact to note in the above figure is that the spectrum of a deep prism does not exceed the magnitude of the same prism at a lesser depth at any wavenumber – only the peak is shifted to lower wavenumbers. Because of this, there is no way to separate the effect of deep sources from shallow sources of the same type using wavenumber filters. This is only possible if the deep sources are of stronger magnitude, or if the shallow sources have a lesser depth extent. When considering a grid that is large enough to include many sources, the log spectrum of this data can be interpreted to determine the statistical depth to the tops of the sources using the following relationship: log E ( r ) = 4πhr

The depth of an “ensemble” of sources is easily determined by measuring the slope of the energy (power) spectrum and dividing it by 4π. A typical energy spectrum for magnetic data may exhibit three parts – a deep source component, a shallow source component, and a noise component.


The following figure illustrates the division of an energy spectrum into these three components.

E

d

N s

n

w MAGMAP is commonly used to enhance information of interest in a given 2D data set, either by removing features considered as “noise”, or by enhancing the features of interest. For example, if you are interested in shallow features in a magnetic map, you might apply a first or second vertical derivative filter to the data in order to enhance shallow features at the expense of anomalies caused by deeper sources. MAGMAP takes advantage of the fact that potential field data, by its nature, is very “broad-band”, so that a single measurement includes the effects due to all the physical (geological) sources. Resolution of the different sources depends on the noise level of the measuring system, and the ability to resolve overlapping signals.

The FFT Algorithm In MAGMAP, Fast Fourier Transform (FFT) is used to convert the space domain grid data to the Fourier domain. The system applies Fast Fourier Transform (FFT) to a space domain grid to produce a folded 2D transform as output. As part of the process, MAGMAP also calculates and saves a radially averaged energy spectrum of the transform. The system creates a Fourier domain grid, which is called a Transform. It has the same name as the input grid, but with the .TRN extension. The transform grid contains a folded discrete Fourier transform of the input grid.


The size of the transform grid element is 4 bytes; each pair of elements represents the real and imaginary component of a complex number. The transform is stored in the same way as the input grid, so that each transform vector (row) represents a vector in the storage direction of the input grid (X for kx=1, Y for kx=-1). The following table illustrates the logical storage of the transform: V (n) -(1/nv) r

ir

ir

i•

r

i

V (n-1)-(2/nv) r

ir

ir

i•

r

i

•

•

•

•

V (n/2+2)

-(1/2v - 1/nv) r

ir

ir

i•

r

i

V (n/2+1)

1/2v (nyq.)

r

ir

ir

i•

r

i

r

ir

ir

i•

r

•

•

•

V (3) 2/nv

r

ir

ir

i•

r

i

V (2) 1/nv

r

ir

ir

i•

r

i

V (1) 0

r

ir

ir

i•

r

i

V (n/2) 1/2v - 1/nv

⇑ Vectors Elements ⇒

0

1/ne

i

•

2/ne

E (1,2) E (3,4) E (5,6)

1/2e (nyq.) E(n+1,n+2)

where: r, i are real and imaginary components of each transform element e, v are element and vector separations (cell size) n is the original grid dimension in cells The transform element separation (1/ne) and vector separation (1/nv) is 1 / (grid dimension) cycles/metre. Since both the grid and the grid cell are square, 1/ne = 1/nv. The Nyquist wavenumber is the largest wavenumber that has been sampled by the grid, and is defined as one over twice the grid cell size (1/2e and 1/2v, which are also equal). Looking at the above table, you can note that each transform vector (row) represents a discrete Fourier row in the direction of the input grid vectors. The Fourier elements within each row start at 0 cycles/metre and extend to the Nyquist wavenumber in 1/ne increments.


As a result, the transform grid has (n/2 + 1) elements per vector, where n is the number of elements per vector in the original grid. The transform is folded at the Nyquist wavenumber in the direction of the grid vectors, so the transform grid has n vectors.

The Energy Spectrum In addition to producing the Fourier transform, the Forward FFT option (FFT2IN.GX) also produces a file containing the radially averaged energy spectrum (2-D Power Spectrum) in a format similar to the following example: / / / / / / /

2-D RADIALLY AVERAGED POWER SPECTRUM WAVENUMBER INTERVAL AVERAGE SPECTRAL DENSITY

DWE = 1.428571E-01 LOG(ETOT) = 1.940511E+01

CYC/KM #_SAMP LOG_P 3_DEPTH ———————————————————————————————0.000000E+00 1 6.953915E+00 2.017633E-01 1.428571E-01 8 6.591711E+00 3.892244E-01 2.857143E-01 12 5.556448E+00 5.608587E-01 4.285714E-01 16 4.578010E+00 3.591071E-01 5.714285E-01 32 4.267114E+00 1.735434E-01 7.142857E-01 28 3.954922E+00 2.310963E-01 8.571429E-01 40 3.437389E+00 3.316897E-01 1.000000E+00 40 2.764027E+00 3.829211E-01 . . . 9.000000E+00 364 -6.672572E+00 1.778083E-01 9.142857E+00 440 -7.611941E+00 2.138515E-01 9.285714E+00 400 -7.440382E+00 -1.844887E-02 9.428572E+00 424 -7.545702E+00 1.084567E-01 9.571428E+00 420 -7.829783E+00 -3.100928E-02 9.714286E+00 416 -7.434367E+00 -3.970939E-02 9.857142E+00 448 -7.687212E+00 2.171140E-01 1.000000E+01 394 -8.213891E+00 2.933830E-01 1.014286E+01 340 -1.024793E+01 * 1.028571E+01 348 -1.104301E+01 * 1.042857E+01 272 -1.107746E+01 *

5_DEPTH * * 4.363967E-01 3.645031E-01 2.545823E-01 2.454431E-01 3.152357E-01 2.908998E-01

5.273030E-02 1.244036E-01 1.012864E-01 1.966618E-02 1.257935E-02 4.879844E-02 * * * * *

The radially averaged energy listed in the third column represents the spectral density (energy) averaged for all grid elements at the wavenumber in the first column. The second column indicates the number of elements that were used to determine the average. The energy is normalized by subtracting the log of the average spectral density. The 3-DEPTH and 5-DEPTH columns are ensemble magnetic depth estimates based on 3 and 5 point averages of the slope of the energy spectrum (Spector and Grant, 1970). The depth to a statistical ensemble of sources is determined by the following expression:


h = -s/ 4π where: h is depth s is the slope of the log (energy) spectrum The above estimates can be used as a rough guide to the depth of magnetic source populations. The system enables you to create and view a radially averaged spectrum automatically. The plot format is shown below.

The above plot illustrates the typical reduction in energy with increasing wavenumber. The depth estimate is a plot of the 5-point depth data from the spectrum file.

Chapter 5: Preparing Grids for the FFT Processing 49

Chapter 5: Preparing Grids for the FFT Processing Grid preparation includes the following basic processes: 1. Removing a first order trend from the grid. The removed trend is stored in the

user area of the grid header, and is filtered together with the zero wavenumber. 2. Expanding the grid dimensions by adding dummy areas to the grid edges to

produce a square grid. By default, the grid size is increased by a minimum of 10%, and then the next largest acceptable dimension is selected. The system uses the FFT algorithm for dimensions up to 2520 x 2520 cells. Beyond that, the algorithm switches to the power of 2 FFT methods. 3. Replacing all dummy values in the grid with the interpolated values from the

valid parts of the grid. The following diagram illustrates the effect of the above processes in one dimension. Because the pre-processed grid must be periodic, it is important to remove a first order trend before expanding and filling. Original Data Trend Removal (GRIDTRND)

Expand (GRIDXPND) and Fill to Periodic (GRIDFILL)

If a significant trend is left in the data, the expansion and filling processes are forced to introduce a large step function in order to make the data periodic. In the Fourier domain, this step function predominates, and might cause ringing problems.

Trend Removal Algorithm The system first removes a first order trend from the data. Because the data is made periodic when filled, this procedure removes the amount of step function that may be required to connect the grid at opposite edges. The trend surface is (by default) calculated using the edge points of the data, so that no strong anomalies within the grid affect the trend. The trend coefficients (first order only) are stored in the grid file header and, during filtering, anything applied to the zero wavenumber of the data is also applied to the trend coefficients.

50 Chapter 5: Preparing Grids for the FFT Processing

Grid Expansion Algorithm The grid must be expanded in size to have dimensions that are acceptable to the FFT: 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 24, 28, 30, 36, 40, 42, 48, 56, 60, 70, 72, 80, 84, 90, 112, 120, 126, 140, 144, 168, 180, 210, 240, 252, 280, 315, 336, 360, 420, 504, 560, 630, 720, 840, 1008, 1260, 1280, or 2520 cells. Also, the expansion provides an area for extending the grid to be smoothly periodic. Grid expansion simply pads the edges of the grid with dummy values. The system allows you to specify a minimum percentage for expansion. The default is 10%, after which the next largest acceptable dimension is selected. In most cases, the default size is acceptable. However, if the wavelength of the anomalies of interest is more than approximately two times the default padded size, you must increase the expanded size by increasing the expansion percentage. If the expansion is too small, any step in the expanded area can adversely affect the anomalies within the data.

Grid Filling Algorithm The system replaces all the dummy values within a grid with the interpolated values. It starts by interpolating the blank areas of the grid by: 1. Replacing dummies within each grid row so that the grid lines are periodic 2. Replacing dummies within each grid column so that the grid columns are periodic 3. Averaging the results from the row and column filling. The averaging procedure

uses inverse distance weighting to the nearest data within each row and column. This procedure effectively fills holes in the data and accounts for irregular edges of the grid. When filling, you must avoid introducing abrupt step functions in the filled grid. You can suspect a step problem if you observe “ringing” in a filtered output grid. Ringing is a symptom of Gibb's Phenomena, which is observed when you modify the Fourier Spectrum of a step function. Filters, by definition, must modify Fourier Spectrum; therefore, filtered maps are susceptible to Gibb's Phenomena if you are not careful when filling the grid. Ringing can be identified as a wave pattern that extends away from, or around, a strong anomaly. The wavelength of the pattern is normally near the size of the strong feature in the data. Because the system interpolates data beyond the edges of the grid, it can introduce step functions that cause ringing to spread into a filtered grid.

Minimizing "Ringing" from Grid Filling MAGMAP provides the following parameters that can be set to minimize ringing problems when they occur. For a complete summary of the trend removal, grid

Chapter 5: Preparing Grids for the FFT Processing 51

expansion, and filling parameters, see the Setting Trend Removal, Grid Expansion, and Filling Parameters section at the end of this chapter. Applying Maximum Entropy Prediction If required, you can use Maximum Entropy Prediction (MEP) to interpolate the data. MEP samples the original data near the grid edges to determine its spectral content. It then predicts a data function that would have the same spectral signature as the original data. This means that if the original data is smooth, the predicted data is smooth, and if the original data is noisy, the predicted data is noisy. As a result, the predicted data will not significantly alter the energy spectrum that would result from the original data alone. Also, this method allows noisy data on one edge of a grid to be gradually interpolated into smooth data on the opposite edge of a grid. Controlling Edge Effects Sometimes, the prediction function can produce large ridges that extend away from the edges of the data. If you suspect ringing caused by an edge effect, look at the preprocessed grid to see if edge filling has produced strong ridges in the filled areas. This typically happens in the originally rectangular grids for which a large expansion is required in one direction in order to make the grid square. If you suspect such an edge problem, you can set the distance at which to roll the data to zero. Ensure that the roll to zero distance is at least as large as the longest anomalies of interest along the edges of the data, otherwise the edge anomalies may be distorted. Limiting Strong Anomalies Near Grid Edges Very strong anomalies that are truncated at the edge of a grid can cause problems in filling because their magnitude is extended into the filled borders of the grid. Step functions parallel to the grid edges result, and ringing can manifest itself as a streak in a filtered grid that extends away from the strong anomaly, or appears on the opposite side of the grid (remember that the grid is considered periodic, so an anomaly on one edge can effect the opposite edge of the grid). If required, you can smoothly limit the magnitude of strong anomalies within a certain distance of the edge of a grid. You can specify the maximum edge magnitude and distance. All anomalies of greater magnitude than the edge limit will be smoothly attenuated starting at half the magnitude limit. In extreme cases, it may be necessary to totally limit edge anomalies by specifying the maximum magnitude of 0. This is similar to applying the Hamming window to the grid, which is the more conventional approach used to handle the edges of data in

52 Chapter 5: Preparing Grids for the FFT Processing

the Fourier processing. However, this method produces a pronounced edge ring around most data grids, which is usually unacceptable. Limiting Strong Anomalies by Magnitude Very strong anomalies within the original data area of a grid can also cause ringing problems in filtered maps. For example, iron formation anomalies in magnetic data can be many orders of magnitude greater than the surrounding anomalies. The magnitude of these anomalies is so great that they dominate the Fourier spectrum, and even small changes to the spectrum result in ringing. Such ringing appears as waves that surround the large anomaly. If required, you can smoothly limit any anomalies that exceed a specified magnitude. Data less than half the limiting magnitude is not changed. Above half the limiting magnitude, data is smoothly attenuated so that it does not exceed the limit. If limited anomalies are wide, this can cause the attenuated anomalies to have flat tops, a fact that you should be aware of when interpreting the resulting filtered maps. If flat-topped anomalies are not wanted, another option is to clip high magnitude anomalies using Geosoft's Grid Windowing GX (for more information, contact your Geosoft representative). Clipped areas must be set to dummy values, and the resulting processed grids will have holes where the anomalies have been clipped.

Setting Trend Removal, Grid Expansion, and Filling Parameters As discussed previously in this chapter, Geosoft provides a variety of capabilities for preparing grids prior to applying the forward FFT. These capabilities are controlled through the Prepare Grid menu option. When you select this option, the system displays the FFT2 grid pre-processing dialog. The following table summarizes the parameters that can be set in this dialog and provides guidelines for setting them. Trend Removal Parameters Type of trend surface to remove

Select order of trend to remove (the default is the first order). Options are: remove mean, first order, second order, third order.

trend based on

Select or type either 'edge points' or 'all points'. The trend surface to remove can be calculated either by using all the valid points in the grid, or by using only the points along the valid edge of the grid. Using the edge points is often better, especially if there are any large-magnitude anomalies within the grid.

Chapter 5: Preparing Grids for the FFT Processing 53 Grid Expansion Parameters % expansion

Type the expansion size (grid is expanded in size by at least this distance as a percentage of the smallest grid dimension). The expansion must be about half the size of the broadest features of interest in the grid. If you are tapering the data to 0, the expansion need not be larger than the taper distance.

Square or rectangular expansion

Select or type either 'square' or 'rectangular'. If the grid is small, or if the wavenumbers of interest approaches the size of the grid, we recommend square expansion because it minimises side effects that result from having different wavenumber samples in the X and Y directions. Rectangular grids can save significant processing time and disk space when working with large grids.

Grid Filling and Ringing Parameters grid fill method

Select or type a fill method. When filling the dummy areas, the new values are determined by extrapolation from the nearest valid parts of the grid. This extrapolation may be based on inverse distance weighting or maximum entropy prediction. Maximum entropy is slower but it creates a filled area more similar in character to the actual data.

roll-off to 0 at a distance of (cells)

Specify the number of cells beyond the valid area at which to roll off to zero. By default, no roll-off is applied. For some grids, the prediction function can become unreasonable at large distances from the valid parts of the grid. In these cases, the data can be forced to zero at a specified distance. This option should only be used on trend-removed grids.

limit all magnitudes to be less than

High-magnitude anomalies can cause problems in filtering systems such as MAGMAP. With this option, anomalies that exceed half the specified limit are smoothly attenuated. The attenuation is started at half the limit, with no values allowed to exceed the limit. This option must be used only on trend-removed grids.

edge magnitude limit

High-magnitude anomalies on the edges of the valid area can produce oscillations in the extrapolated areas. With this option, a limit may be placed on anomalies along the edges of the grid.

54 Chapter 6: Specifying FFT Filters

Chapter 6: Specifying FFT Filters All MAGMAP processes are carried out by the application of filters in the Fourier (wavenumber) domain. Filters are simply multiplied by the transform of the grid on an element-by-element basis. Once a Fourier transform has been created, the application of filters is quite straightforward. When using the Step-By-Step method, you select the Define Filters menu option to either create a new control file or select an existing control file. With the Interactive filtering method, you select the Interactive Spectrum Filters menu option to create a control file. When you are ready to proceed, you select the Apply Filters menu option to apply filters according to instructions defined in your new or existing MAGMAP control file.

Creating Your Own Filter Control File The MAGMAP control file used by the Apply Filter option is an ASCII text file that may be created: • •

•

Using a text file editor By copying, renaming, and editing the mapplot.con file in the C:\Program Files\Geosoft\Oasis montaj\etc directory. The mapplot.con file is a blank control file that contains a brief description of the parameters and filters available in MAGMAP. Using the MAGMAP|Step-By-Step|Define Filtesr menu option or the MAGMAP|Interactive Filtering|Interactive Spectrum Filters menu option.

Following is an example of a MAGMAP control file: 70 /geomagnetic inclination 0 /geomagnetic declination BPAS 0.0001 0.003 1 / BTWR 0.0002 4 0 / TXYZ 0 3 / CNDN 50 / CNUP 200 / COSN 0.001 0.003 2 1 /

The control file must contain 6 or more lines with the following information: line 1:

A title line for reference only.

line 2:

The nominal height of the magnetic sensor above the ground (normally, the flying height), or above the level of magnetic sources. This information is used as the default height for some of the filters, and is not always necessary or relevant.

line 3:

The magnetic inclination (negative in the Southern hemisphere). The inclination is only used by the reduction to the pole (REDP), reduction to the equator (REDE), susceptibility (SUSC), and Weiner optimum (OPTM) filters.

Chapter 6: Specifying FFT Filters 55 line 4:

The magnetic declination in degrees of azimuth relative to true North. The inclination is only used by the reduction to the pole (REDP), reduction to the equator (REDE), susceptibility (SUSC) and Weiner optimum (OPTM) filters. Note that grid files contain the direction of the grid’s Y axis in the grid header as the rotation parameter. Filter takes this direction into account so that the Line 4 parameter can be the true declination. The grid header must be correct. Very often, the grid rotation is reported as 0 in the grid header, in which case the declination specified here must be the declination of magnetic North relative to the grid’s Y axis.

line 5:

The nominal total magnetic field strength in nano-Tesla (gammas). This value is only used by the apparent susceptibility map filter (SUSC) in order to derive susceptibility from magnetization.

line 6+

One or more lines specifying the filters to be applied and their parameters. The following section documents the available filters.

The forward slash character (/) must terminate a line (with the exception of the title line), and user comments may follow the slash. After the fifth line, all lines must start either with a forward slash or a filter name. Each filter option occupies one or more lines and consists of a four-letter mnemonic followed by the optional parameter settings. Parameters (if provided) must be separated by a space. Any number of filters may be applied in a single filtering run. However, only one output transform is produced. Note that because multiplication is commutative, the order in which filters are applied is not relevant.

MAGMAP Filters This section describes the available filters. In each filter-specific subsection, filter options are listed in alphabetical order. Each description shows the mathematical expression of the filter followed by a figure if appropriate, then the control file parameters, and usage notes. The filter expressions use the following basic expressions: μ

X wavenumber (complex, radians/ground_unit)

ν

Y wavenumber (complex, radians/ground_unit)

r =

μ2 + ν2

Wavenumber (radians/ground_unit)

θ = tan −1 ( μ ν )

Wavenumber direction

N

Nyquist wavenumber [1 / (2 * cell size)]

k

Wavenumber in cycles/ground_unit (r=2πk)

The horizontal axis of the figures represents wavenumbers between 0 and the Nyquist frequency. All distance references are multiples of the grid cell size. For example, referring to the filter response drawing for upward continuation (CNUP) filter, if the


grid cell size is 25 ground_units, the Nyquist wavenumber is 0.02 cycles/ground_unit, the filter response curves would represent upward continuation distances of 50, 100, 200 and 400 ground_units. BPAS – Bandpass Filter

L( k ) = 0, for k < k 0 L( k ) = 1, for k 0 ≤ k ≤ k1 L( k ) = 0, for k > k1 1.0

L(k)

0.5

reject

pass

reject

0.0

k0

k1

Wavenumber (cycles/ ground_unit)

Parameters: k0

The low wavenumber cut-off in cycles/ground_unit.

k1

The high wavenumber cut-off in cycles/ground_unit.

0/1

If 1, pass the defined band; if 0, reject the defined band. The default is to pass the band.

Ground_unit The survey ground units used in your grid (e.g., metre, foot, etc.). If your “ground units” are in metre, the low and high wavelength cutoff is in cycles/metre. BPAS can be used to pass or reject a range of wavenumbers from the data. Applying such a simple cut-off filter to an energy spectrum almost invariably introduces a significant amount of ringing (otherwise known as Gibb's Phenomena). We recommend that you use a smoother filter, such as the Butterworth filter (BTWR).

Chapter 6: Specifying FFT Filters 57

BTWR – Butterworth Filter

L(k ) =

1 ⎡ ⎛ k ⎞n ⎤ ⎢1 + ⎜ k ⎟ ⎥ ⎣ ⎝ 0⎠ ⎦ 1.0 n =2

L(k)

4

8

16

0.5

0.0

0

ko

N


Parameters: k0

The central wavenumber of the filter.

n

The degree of the Butterworth filter function. By default, 8.

0/1

A flag (0 or 1). Specifies if a residual (0) high pass or a regional (1) low pass is required. By default, a regional filter is applied.

Ground_unit The survey ground units used in your grid (e.g., metre, foot, etc.). If your “ground units” are in metre, the low and high wavelength cutoff is in cycles/metre. The Butterworth filter is excellent for applying straight forward high-pass and lowpass filters to data because you can easily control the degree of filter roll-off while leaving the central wavenumber fixed. If ringing is observed, the degree can be reduced until the result acceptable. A common but more complicated alternative is the Cosine filter (COSN).


CNDN – Downward Continuation

L( r ) = e hr 20.0

h = 16

h=8

L(r)

h=4 h=2 1.0 0

N


Parameters: h

The distance, in ground units, to continue down relative to the plane of observation.

r

Wavenumber (radians/ground_unit). Note: r = 2πk, where k is cycles/ground_unit

Ground_unit The survey ground units used in your grid (e.g., metre, foot, etc.). Downward continuation is used to enhance the responses from sources at a depth by effectively bringing the plane of measurement closer to the sources. Note that it is not theoretically possible to continue through a potential field source. Since shortwavelength signal can appear to be from shallow sources, it must be removed to prevent high magnitude and short wavelength noise in the processed data. To do this, you usually apply some type of low-pass filter, such as the Butterworth or Weiner Optimum filter. You should use a low-pass filter to remove the short wavelength noise (as determined by the radially averaged energy spectrum) before applying the downward continuation filter. The energy spectrum is also a good guide for determining the depth to which the data can be continued downward. CNUP – Upward Continuation

L( r ) = e − hr 1.0

h=2

L(r) h=8

h=4

h = 16 0.0

0

N



Parameters: h

The distance, in ground units, to continue up relative to the plane of observation.

r

Wavenumber (radians/ground_unit). Note: r = 2πk, where k is cycles/ground_unit.

Ground_unit The survey ground units used in your grid (e.g., metre, foot, etc.). Upward continuation is considered a “clean” filter because it produces almost no side effects that may require application of other filters or processes to correct. Because of this, this filter is often used to remove or minimize the effects of shallow sources and noise in grids. Also, upward continued data may be interpreted numerically and with modeling programs. This is not the case for many other filter processes. COSN – Cosine Roll-off Filter L( k ) = 1, for k < k 0 ⎡π ⎛ k − k0 ⎞ ⎤ L( k ) = cos n ⎢ ⎜ ⎟ ⎥, for k 0 ≤ k ≤ k1 ⎣ 2 ⎝ k1 − k 0 ⎠ ⎦ L( k ) = 0, for k > k1 1.0

L(k)

0.5

0.0

n=2

0

k0

1

0.5

k1

N

Wavenumber cycles/ ground_unit)

Parameters: k0

Low wavenumber starting point of the filter (cut-off wavenumber for high pass or start of roll off for low pass.

k1

High wavenumber end point of the filter (start of roll off for high pass or cut-off wavenumber for low pass.


n

The degree of the cosine function. The default is a degree of 2 for a cosine squared roll off.

0/1

0 for residual (high-pass) filter; 1 for regional (low-pass) filter. The default is a low-pass filter.

Ground_unit The survey ground units used in your grid (e.g., metre, foot, etc.). Because this filter has a smooth shape, and it does not alter the energy spectrum below the start of roll off (or after the end of roll off in high-pass mode), it is commonly used for simple high-pass or low-pass operations. To reduce ringing, the separation between r1 and r0 can be increased. DCOS – Directional Cosine Filter π n L ( θ)) = cos (α − θ + 2 ) , to reject direction α

L( θ) = 1 − cos n (α − θ +

π 2

),

to pass direction α

1.0

L( θ ) 0.5

0.0 α−π/2

n=2

1

0.5

α

θ (direction)

α+π/2


n=2

n=1

α = 70

N

1

N

1

L (u,ν) ν

0 N

0 0

-N

L (u,ν) ν

0 N

0

0

u -N

-N

u

-N

n = 0.5

N

1

L (u,ν) ν

0 N

0

0

-N

u

-N

Parameters:

α

Direction of the filter in degrees (0-360 relative to North).

n

The degree of the cosine function. By default, a degree of 2 is used to give a cosine squared function.

0/1

If 1, apply the filter to pass the specified direction; if 0, apply the filter to reject the specified direction. By default, the direction is rejected.

The directional cosine filter is very good for removing directional features from a grid. The cosine function makes the filter smooth, so directional ringing effects are usually not a problem. The rejection (or pass) notch can be narrowed or widened by setting the degree of the cosine function, so that highly directional features can be isolated. De-corrugation of poorly levelled magnetic data is a common application for this filter (see examples). DENS – Apparent Density Calculation

L( r ) =

r 2π G(1- e -tr )

where: G

Gravitational constant.


r

Wavenumber (radians/ground_unit). Note: r = 2πk, where k is cycles/unit.

Parameters: t

Thickness, in ground units, of the earth model.

d

Background density in g/cm3, to be added to the density contrast map. The default is 0, so the density map is produced relative to the average density.

Ground_unit The survey ground units used in your grid (e.g., metre, foot, etc.). The data must also be downward continued using the CNDN filter to be close to the top surface of the source model. Apparent density mapping assumes that an observed gravity field can be explained by a simple model layer of fixed thickness and varying density. This is a poor model in most cases. DPAS – Directional Pass/Reject Filter θ

0

= 50 θ

N

1

= 70

1

L (u,ν) ν

0 N

0 0

-N

u

-N

Parameters:

θ0

The low cut-off angle in degrees of azimuth (0-360 from North).

θ1

The high cut-off angle in degrees of azimuth (0-360 from North).

0/1

If 1, pass the defined band; if 0, reject the defined band. The default is to pass the band.

As with the band-pass filter, the directional-pass often suffers from Gibb's Phenomena (ringing) because the spectrum is cut quite abruptly. We recommend using the directional cosine filter (DCOS) instead.


DRVX – Derivative in the X Direction L ( μ ) = ( μi )

n

Parameter: n

Order of differentiation.

u

The X component of the wavenumber.

i

i=

−1

The horizontal derivative can be used for creating shaded images, and is required for some modeling algorithms, such as Euler de-convolution. DRVY – Derivative in the Y Direction L(ν ) = (νi )n

Parameter: n


v

The Y component of the wavenumber.

i

i=

−1

The horizontal derivative can be used for creating shaded images, and is required for some modeling algorithms, such as Euler de-convolution. DRVZ – Derivative in the Z Direction

L( r ) = r n Parameter: n


r

Wavenumber (radians/unit). Note: r = 2πk, where k is cycles/unit.

The vertical derivative is commonly applied to total magnetic field data to enhance the shallowest geologic sources in the data. As with other filters that enhance the high-wavenumber components of the spectrum, you must often also apply low-pass filters to remove high-wavenumber noise.


GFILT – Gravity Earth Filter L(r ) = 2πG

(e − rz 1 − e − rz 2 ) r

When r=0:

L(r ) = 2πG( z 2 − z 1 ) Where: G

Gravitational constant

r

Wavenumber (radians/ground_unit). Note: r = 2πk, where k is cycles/unit.

Parameters: z1

Depth to the top of the density layer, in ground units. Must be positive for layers above calculation level.

z2

Depth to the bottom of the density layer, in ground units. Must be positive for layers above calculation level.

Ground_unit The survey ground units used in your grid (e.g., metre, foot, etc.). GAUS – Gaussian Regional/Residual Filter

L (k ) = 1 − e

−k 2 2 k02

1.0

L(k)

0.5

0.0

0

K0

2K 0

3K 0

N


Parameters: k0

Standard deviation of the Gaussian function in cycles/ground_unit (similar to a cut-off point, except that the function magnitude at this point is only 0.39).


0/1

If 0, the residual component is produced; if 1, the regional component is produced. The default is 0.

Ground_unit The survey ground units used in your grid (e.g., metre, foot, etc.). The Gaussian filter is another smooth filter that is often used for low-pass or highpass applications. GNRL – General Radially Symmetric Filter

L( 0) = l0

L( dk ) = l1

L( 2dk ) = l2 M

L( n ⋅ dk ) = l n L( k ) = l n , for ( k ≥ n ⋅ dk ) Parameters: dk

The wavenumber increment (cycles/ground_unit), starting from zero wavenumber, at which the following filter coefficients are applied.

li

The coefficients of the filter function at each wavenumber increment, starting at zero wavenumber. The last value given is used for all higher wavenumbers. More than one line can be used to give the coefficients. A slash character (/) must follow the last coefficient.

Ground_unit The survey ground units used in your grid (e.g., metre, foot, etc.). The general filter is used when a special-purpose filter must be applied to the data. Normally, filter coefficients are between 0 and 1. For example, the following defines a low-pass filter that starts at a wavenumber of 0.003 and rolls off to remove all wavenumbers above 0.007: L(0) = 1 L(0.001) = 1 L(0.002) = 1 L(0.003) = 1 L(0.004) = 0.8 L(0.005) = 0.5 L(0.006) = 0.2 L(0.007) = 0

To define this filter in a control file, you would type” gnrl 0.001 1. 1. 1. 1. .8 .5 .2 0


GPSD – Pseudo-Gravity Filter

L(θ ) =

G∗d / J I =I < [sin(I a )+ i cos(I )∗ cos (D− θ )] 2∗ r , if ( | Ia | |I |), a

Where: I

Geomagnetic inclination

Ia

Inclination for magnitude correction (never less than I). The default is ± 20 degrees. If |Ia| is specified to be less than |I|, it is set to I.

d

Density contrast in g/cm^3.

G

Gravitational constant, 6.670E-8.

J

Magnetization in Gauss.

D

Geomagnetic declination.

Parameters:

θ

Direction of wavenumber in degrees of azimuth.

r


Ground_unit The survey ground units used in your grid (e.g., metre, foot, etc.). Output of Pseudo-Gravity is in mgal if ground_unit is metre, or in 0.3mgal if groud_unit is foot.


HPAS – High-pass Filter L( k ) = 0 , for k < k0 L( k ) = 1 , for k ≥ k0

1.0 reject

L(k)

pass

0.5

0.0

k0

Wavenumber (cycles/ground_unit) Parameters: k0

The cut-off wavenumber in cycles/ground_unit. All wavenumbers below this value are removed.

Ground_unit The survey ground units used in your grid (e.g., metre, foot, etc.). As with the band-pass filter, the high-pass filter is seldom used because the results usually suffer from Gibb's Phenomena (ringing). INTG – Vertical Integration

L(r) = r -1 Parameters: r


Ground_unit The survey ground units used in your grid (e.g., metre, foot, etc.). This filter calculates the vertical integral of the input transform. This is the inverse of the vertical derivative. The zero wavenumber is set to 0.


LPAS – Low-pass Filter

L(k) = 1, for k ≤ k0 L(k) = 0, for k > k0 1.0 pass

L(k)

reject

0.5

0.0

k0

Wavenumber (cycles/ground_unit) Parameters: k0

The cut-off wavenumber in cycles/ground_unit. All wavenumbers above this value are removed.

Ground_unit The survey ground units used in your grid (e.g., metre, foot, etc.). As with the band-pass filter, the low-pass filter is seldom used because the results usually suffer from Gibb's Phenomena (ringing). OPTM – Weiner Optimum Filter

The Weiner optimum filter is intended to remove the effect of white noise from the magnetic data. White noise is high-wavenumber background noise present in the data. Because magnetic signal is stronger in the direction of the inducing field, the signalto-noise ratio varies as a function of both inclination and declination. The Weiner optimum filter takes the variation of signal-to-noise ratio into account when applying the filter. ⎡ φ s (k ,θ ) ⎤ L(k ,θ ) = ⎢ ⎥ , for k < k 0 ⎣φ s (k ,θ ) + φ 0 ⎦ ⎡ φ s (k ,θ ) ⎤ 2 π ⎛ k − k0 ⎞ ⎜⎜ ⎟⎟ , for k 0 ≤ k ≤ k1 L(k ,θ ) = ⎢ ⎥ ⋅ cos ( ) + − k , 2 k k φ θ φ 0 ⎦ 0 ⎠ ⎣ s ⎝ 1 L(k ,θ ) = 0 , for k > k1

[sin

]

I + cos 2 I ⋅ cos 2 (D + θ ) φ s (k ,θ ) = ⋅ (φ k − φ 0 ) sin 4 I + sin 2 I ⋅ cos 2 I + 0.375 cos 4 I 2

2


Where: I


D

Geomagnetic declination

φr

Average radial spectral density

φ0

Average radial spectral density of noise

k0

Wavenumber of start of the noise

k1

Wavenumber at the end of the applied roll-off filter

Parameters: h

The depth at which to interpret an optimum filter from the observed energy spectrum. By default, the depth is taken as the flying height in the control file, or the continuation depth if specified by the Downward Continuation or Apparent Susceptibility filter options (in ground_units).

k0

The wavenumber (cycles/ground_unit) at which to start the highwavenumber roll off. This parameter must be given together with φ0. By default, this point is the point at which the slope of the observed energy spectrum rises above the slope defined by the depth (1/4πh).

k1

The wavenumber (cycles/ground_unit) at which to end the high wavenumber roll off. By default, this point is set to be two times k0.

φ0

Spectral density estimate of noise to be removed by the Weiner filter. This is in terms of the log of spectral density as reported in the second column of the energy spectrum. By default, this is calculated as the average of the spectral density between k0 and k1.

Ground_unit The survey ground units used in your grid (e.g., metre, foot, etc.). The optimum filter is most often used to remove the theoretical effect of all sources that lie above a specified depth. The filter parameters can be specified or calculated automatically based on analysis of the energy spectrum. Although this filter can calculate the parameters of the filter, we recommend that you confirm that the calculated parameters are reasonable. When the energy spectrum is not smooth, the filter can choose the wrong point at which to start the noise calculation. Most often, this point is chosen to be too low, and the resulting maps appear too smooth.


The optimum filter can be quite complex to use and understand. A good alternative is using the Butterworth filter as a low-pass filter. Determine the wavenumber at which sources appear too shallow by interpreting the depth estimate in the energy-spectrum plot. REDE – Reduce to the Magnetic Equator L (θ ) =

[sin

[sin( I ) − i ⋅ cos( I ) ⋅ cos( D − θ )]2 × (− cos 2 ( D − θ )) 2

] [

( Ia ) + cos 2 ( Ia ) ⋅ cos 2 ( D − θ ) × sin 2 ( I ) + cos 2 ( I ) ⋅ cos 2 ( D − θ )

] , if (| Ia |<| I |), Ia = I

where: I


D


L(θ) Direction of the wavenumber vector in degrees of azimuth No parameters. Reduction to the equator is used in low magnetic latitudes to centre the peaks of magnetic anomalies over their sources. This can make the data easier to interpret while not losing any geophysical meaning. Reducing the data to the pole (REDP) does much the same thing, but at low latitudes. A separate magnitude correction is usually required to prevent North-South signal in the data from dominating the results. As a result, reduced to the pole data may present a less “honest” view of the data REDP – Reduce to the Magnetic Pole L(θ ) =

[sin

[sin( I ) − i ⋅ cos( I ) ⋅ cos( D − θ )]2 2

][

( Ia ) + cos 2 ( Ia ) ⋅ cos 2 ( D − θ ) ⋅ sin 2 ( I ) + cos 2 ( I ) ⋅ cos 2 ( D − θ )

], if (| Ia |<| I |), Ia = I

where: I


Ia

Inclination for magnitude correction (never less than I)

D


Parameters: Ia

Inclination to be used for the magnitude correction. The default is ± 20 degrees. If |Ia| is specified to be less than |I|, it is set to I.


Reduction to the pole has a magnitude component (the sin(I) term) and a phase component (the i cos(I) cos(D-θ) term). When reducing to the pole from equatorial latitudes, North-South features can be exaggerated due to the strong magnitude correction (the sin(I) component) that is applied when (D-θ) is π/2 (i.e., a magnetic East-West wavenumber). By specifying higher latitude for the magnitude correction alone, this problem can be reduced or eliminated at the expense of under-correcting the magnitudes of North-South features. A magnitude inclination of 90 causes only the phase component to be applied to the data (no magnitude correction), and a value of 0 (zero) causes phase and magnitude components to be applied over the entire range. SUSC – Apparent Susceptibility Calculation

The susceptibility filter is, in fact, a compound filter that performs a reduction to the pole, downward continuation to the source depth, correction for the geometric effect of a vertical square-ended prism, and division by the total magnetic field to yield susceptibility. L(r ,θ ) =

1 2πF ⋅ Η (r )⋅ Γ(θ )⋅ Κ (r ,θ )

Η (r ) = e − hr Γ(θ ) = [sin Ia + i cos I ⋅ cos ( D − θ

)] 2 ⎛ sin ( ar cosθ ) ⋅ sin( ar sin θ) ⎞ Κ (r ,θ ) = ⎜ ⎟ ⎝

ar cos θ ⋅ ar sin θ

⎠

where: Η(r)

Downward continuation to h

Γ(θ)

Reduction to the pole

Κ(r,θ) Geometric factor of a vertical prism (a⋅a⋅∞ in dimension) I


Ia

Pole reduction magnitude inclination

D


F

Total geomagnetic field strength

r


72 Chapter 6: Specifying FFT Filters θ

Wavenumber direction (azimuth)

a

Half the cell size

Parameters: h

Depth in ground units, relative to the observation level at which to calculate the susceptibility. By default, the flying altitude reported in line 2 of the MAGMAP Control File is used.

Ia

Inclination to which to use the phase component only in the reduction to the pole. The default is ± 20. If |Ia| is specified to be less then |I|, it is set to I.

Ground_unit The survey ground units used in your grid (e.g., metre, foot, etc.). The susceptibility filter calculates the apparent magnetic susceptibility of the magnetic sources using the following assumptions: • • • •

The IGRF is removed from the magnetic field . There is no remnant magnetization All magnetic responses are caused by a collection of vertical, square-ended prisms of infinite depth extent The result is in the e.m.u. units

The validity of the results depends on how well the actual observed field conforms to these assumptions. TXYZ – Conversion between Field Components

TXYZ converts from the total field (T) or vertical field (Z) component of the magnetic field to any other (X,Y, Z, or T) component. The following table gives the filter expressions for all possible component field conversions, To X

To Y

To Z

To T

From X

1

v/u

r / iu

P / iu

From Y

u/v

1

r / iv

P / iv

From Z

iu / r

iv / r

1

P/r

From T

iu / P

iv / P

r/P

1


Where: I

Geomagnetic inclination in degrees from the horizon.

D

Geomagnetic declination in degrees azimuth.

i

−1

u

The X wavenumber

v

The Y wavenumber

r

The radial wavenumber =

u 2 +v 2

α

Directional cosine of the Total magnetic field of X axis.

β

Directional cosine of the Total magnetic field of Y axis.

γ

Directional cosine of the Total magnetic field of Z axis.

P

i ∗ (α ∗ u + β ∗ v) + γ ∗ r

Theoretically, any one directional component can be used to calculate any other component. In reality, only the T and Z fields work well as inputs for this conversion because the formulas for conversion from X or Y contain singularities, i.e., u can go to zero while v is large, or vice versa, causing the filter coefficient to become infinite.

Filter Examples This section provides examples of filter application. Susceptibility Map

The following control file (susc.con) applies an optimum depth filter and creates a magnetic susceptibility map calculated at a depth of 200 metres below the survey elevation. susceptibility map 120 72 -12 55000 SUSC 200 OPTM

/sensor elevation /magnetic inclination /magnetic declination /total field strength /susceptibility map /optimum depth filter

The susceptibility map filter includes a downward continuation to the source depth. Because a downward continuation filter magnitude increases with wavenumber, it tends to also amplify high-wavenumber noise in the data. To prevent this noise from entering the final map, you often need to remove high wavenumbers that are considered noise in the data. This is the function of the Weiner Optimum filter applied by the OPTM option.


Second Vertical Derivative

The following control file creates a second vertical derivative. An upward continuation of 20 metres (equivalent to the cell size) is also applied to reduce the effect of high wavenumber noise. The Hanning filter applied after processing can also serve to reduce noise. second vertical derivative 120 72 -12 55000 DRV2 CNUP 20

/sensor elevation /magnetic inclination /magnetic declination /total field strength /second vertical derivative /continue up 20 metres.

De-corrugation

Given a grid surveyed with a nominal line spacing of 150 metres, which has a line-toline levelling problem, the following control file produces a grid that contains the levelling error only. Following a rule of thumb, the Butterworth high-pass filter is set to four times the line separation in order to only pass frequencies on the order of the line separation. The directional cosine filter is set to pass wavelengths only in the direction of the lines (note that North-South line levelling error produces wavenumbers in the East-West direction, hence DCOS 90). Because levelling error is very directional, you can tighten the directional cosine function to an energy value of 0.5. de-corr., 150m line separation, N-S lines 120 /sensor elevation 72 /magnetic inclination -12 /magnetic declination 55000 /total field strength BTWR .00167 8 0 /high-pass butterworth DCOS 90 0.5 1 /directional cosine.

The resulting grid can be subtracted from the original grid to remove the levelling error. Some tuning of both the center wavenumber of the Butterworth filter and the energy of the directional cosine may be required. To remove the more directional signal, increase the energy of the cosine filter. To remove wider features, decrease the Butterworth cut-off point. If you see ringing in the data, decrease the order of the Butterworth filter.

Chapter 7: Applying the Inverse FFT 75

Chapter 7: Applying the Inverse FFT As described earlier, in the Step-By-Step method, you use the Define Filters menu option to select the filters and define their parameters and the Apply Filters menu option to apply the filters to the transform file. Depending on the processing option you select ([OK], [Flt-Inv Only], or [Filter Only]) on the FFT2FLT dialog, you can (i) apply the selected filter, (ii) apply the Inverse FFT (converts the transform file from the wavenumber domain back to the original space domain) and (iii) apply postprocessing (returns the grid to the original dimensions with the trend information restored). The three options for controlling the final results you obtain following the Apply Filters option are detailed below: Processing Option

Selection

Result

1

Filtered grid with post-processing

[OK]

Filtered, space-domain GRD file with original grid dimensions and trend information restored (i.e., post-processing is added to output file).

2

Filtered grid with no post-processing

[Flt-Inv Only].

Filtered, space-domain GRD file as output file with grid-filled dimensions and trend information NOT restored (i.e., no post-processing is added to output file).

3

Filter only – no inverse transform

[Filter Only]

Filtered transform (*.TRN) file as output file.

When you select the first processing option ([OK]) the transform grid is filtered, the Inverse FFT is applied (returning the transform file from the wavenumber domain back to the space domain) and then post-processing restores the square and periodic grid to its original dimensions and restores the trend information. When you select the second processing option ([Flt-Inv Only]) the transform grid is filtered and the Inverse FFT is applied (returning the transform file from the wavenumber domain back to the space domain). No post-processing is applied. When you select the third processing option ([Filter Only]) the transform grid is filtered. The Inverse FFT and post-processing is not applied.

76 References

References Bhattacharya, B. K., 1966, Continuous spectrum of the total magnetic field anomaly due to a rectangular prismatic body. Geophysics, Vol. 31, p.97-121. Spector, A. and Grant, F. S., 1970, Statistical models for interpreting aeromagnetic data. Geophysics, Vol. 35, No.2, p.293-302. Wiener, N., 1949, Extrapolation, Interpolation, and Smoothing of Stationary Time Series. Cambridge, M.I.T. Press. Burg, J. P., 1975, Maximum Entropy Special Analysis. Unpublished doctoral dissertation. Standford University. 168p. Gupta, V. K., and Grant, F. S., 1985, Mineral exploration aspects of gravity and aeromagnetic survey in Sudbury-Cobalt area, Ontario. SEG; The Utility of Regional Gravity and Magnetic Anomaly Maps, W. J. Hinze (Editor) p.392-411. Winograd, S. On Computing the Discrete Fourier Transform, Mathematics of Computation, Vol.32, N0.141, pp.175-199, Jan. 1978. McClellan, J. H. and Nawab H., Complex General-N Winograd Fourier Transform Algorithm (WFTA), Programs for Digital Signal Processing, IEEE Press, pp. 1.7-1 - 1.7-10, 1979. Walter R. Roest, Jacob Verhoef, and Mark Pilkington, 1992, Magnetic interpretation using the 3D analytic signal. Geophysics, Vol. 57, No.1, p.116-125. Ian N. MacLeod, Keith Jones, and Ting-Fan Dai, 1993, 3D Analytic Signal in the Interpretation of Total Magnetic Field Data at Low magnetic Latitudes. Exploration Geophysics, Vol.24, p. 679-688. Bruno Verduzco, J. Derek Fairhead, Chris M. Green, and Chris MacKenzie, Feb. 2004, New insights into magnetic derivatives for structural mapping. The Leading Edge, p. 116-119.

montajMagmapFiltering.pdf

Recommend Documents