OASIS MONTAJ 7.1 MAGMAP FILTERING 2D Frequency Domain 7.1, 2010, 75 Pag-1 PDF [PDF]

Citation preview

montaj MAGMAP Filtering 2D Frequency Domain Processing of Potential Field Data Extension for Oasis montaj v7.1



TUTORIAL and USER GUIDE



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: 08/02/2010. Written by, Nancy Whitehead and Chris Musselman. Please send comments or questions to [email protected] © 2010 Geosoft Inc. All rights reserved. Geosoft is a registered trademark and Oasis montaj is a registered trademark of Geosoft Inc. Other brand and product names mentioned herein are properties of their respective trademark owners. No part of this publication may be reproduced, stored in a retrieval system or transmitted, in any form, or by any means, electronic, mechanical, photocopying, reading, or otherwise, without prior consent from Geosoft Inc. 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. OMx.t.2010.02 Windows®, and Windows NT are either registered trademarks or trademarks of Microsoft Corporation. Geosoft Incorporated Queen’s Quay Terminal 207 Queen’s Quay West Suite 810, PO Box 131 Toronto, Ontario M5J 1A7 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



4 



Contacting Technical Support



Chapter 1: System Capabilities and Concepts



4 



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



18 



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



25 



Tutorial 12: Calculating Analytic Signal



33 



Tutorial 13: Calculating Tilt Derivative



35 



Tutorial 14: Sharing Results



37 



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



42 



The Energy Spectrum



44 



Chapter 5: Preparing Grids for the FFT Processing



46 



Trend Removal Algorithm



46 



Grid Expansion Algorithm



47 



Grid Filling Algorithm



47 



Minimizing "Ringing" from Grid Filling



47 



Applying Maximum Entropy Prediction



48 



Controlling Edge Effects



48 



Limiting Strong Anomalies Near Grid Edges



48 



Limiting Strong Anomalies by Magnitude



49 



Setting Trend Removal, Grid Expansion, and Filling Parameters



Chapter 6: Specifying FFT Filters



49 



51 



Creating Your Own Filter Control File



51 



MAGMAP Filters



52 



BPAS – Bandpass Filter



53 



BTWR – Butterworth Filter



53 



CNDN – Downward Continuation



54 



CNUP – Upward Continuation



55 



COSN – Cosine Roll-off Filter



56 



DCOS – Directional Cosine Filter



56 



DENS – Apparent Density Calculation



58 



DPAS – Directional Pass/Reject Filter



58 



DRVX – Derivative in the X Direction



59 



DRVY – Derivative in the Y Direction



59 



DRVZ – Derivative in the Z Direction



59 



GFILT – Gravity Earth Filter



60 



GAUS – Gaussian Regional/Residual Filter



60 



GNRL – General Radially Symmetric Filter



61 



GPSD – Pseudo-Gravity Filter



62 



HPAS – High-pass Filter



63 



INTG – Vertical Integration



63 



LPAS – Low-pass Filter



64 



OPTM – Weiner Optimum Filter



64 



REDE – Reduce to the Magnetic Equator



66 



REDP – Reduce to the Magnetic Pole



66 



SUSC – Apparent Susceptibility Calculation



67 



TXYZ – Conversion between Field Components



68 



Filter Examples



69 



Susceptibility Map



69 



Second Vertical Derivative



70 



De-corrugation



Chapter 7: Applying the Inverse FFT



70 



71 



Processing Option



71 



Selection



71 



Result



71 



References



73 



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). TO



LOAD THE



MAGMAP



MENU :



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



toolbar. The Load Menu dialog is displayed.



) on the



12 Chapter 2: Before You Begin 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 18 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 25 Tutorial 12: Calculating Analytic Signal….page 33 Tutorial 13: Calculating Tilt Derivative….page 35 Tutorial 14: Sharing Results ….page 37



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 and Image menu, select Display | 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.



Chapter 3: montaj MAGMAP Filtering Tutorials 15



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), prepare 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 the grid smoothly periodic), and grid filling (to replace all the dummy values with interpolated values).



16 Chapter 3: montaj MAGMAP Filtering Tutorials



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”.



Chapter 3: montaj MAGMAP Filtering Tutorials 17 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.



18 Chapter 3: montaj MAGMAP Filtering Tutorials



The system computes the transformation. It also creates a new Input Transform File, “prep_in_trn.grd” to which you will later apply the Filters and Inverse Transform.



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: 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.



Chapter 3: montaj MAGMAP Filtering Tutorials 19 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 Wavelength (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.



20 Chapter 3: montaj MAGMAP Filtering Tutorials



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 71. 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.



2. Using the […] button, enter prep_in_trn.grd as a value for the Name of Input



Transform (*_trn.grd) File field.



Chapter 3: montaj MAGMAP Filtering Tutorials 21 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. You will now comparatively display the original mag_in.grd and the newly created mag_out.grd as colour-shaded grids on temporary maps. TO



DISPLAY GRIDS :



1. From the Grid and Image menu, select Display | Colour-Shaded Grid.



The Color-shaded grid image dialog is displayed. 2. Using the […] button, select the Grid name as mag_in.grd.



22 Chapter 3: montaj MAGMAP Filtering Tutorials 3. Accept the default values for the rest of the parameters, and click the [Current



Map] button. 4. Repeat Steps 1-3 for mag_out.grd. 5. 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.



Chapter 3: montaj MAGMAP Filtering Tutorials 23 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 39. 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”.



24 Chapter 3: montaj MAGMAP Filtering Tutorials 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.



Chapter 3: montaj MAGMAP Filtering Tutorials 25 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:



26 Chapter 3: montaj MAGMAP Filtering Tutorials



• • •



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 51. 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 Applying 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.



Chapter 3: montaj MAGMAP Filtering Tutorials 27 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”.



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.grd. This should be the default value. 3. Click the [OK] button.



28 Chapter 3: montaj MAGMAP Filtering Tutorials



The system computes the transform. It also creates a new Input Transform File, “Interactive_prep_trn.grd” to which you will calculate the Radially Averaged Spectrum and later apply the Filters and Inverse 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 _trn.grd. 3.



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



4. Click the [OK] button.



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”.



Chapter 3: montaj MAGMAP Filtering Tutorials 29 4. Click the [OK] button. The Interactive Spectral Filter window is displayed.



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



30 Chapter 3: montaj MAGMAP Filtering Tutorials



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.



Chapter 3: montaj MAGMAP Filtering Tutorials 31



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).



32 Chapter 3: montaj MAGMAP Filtering Tutorials 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_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.



Chapter 3: montaj MAGMAP Filtering Tutorials 33 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.



34 Chapter 3: montaj MAGMAP Filtering Tutorials 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.



Chapter 3: montaj MAGMAP Filtering Tutorials 35



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)



36 Chapter 3: montaj MAGMAP Filtering Tutorials



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.



Chapter 3: montaj MAGMAP Filtering Tutorials 37



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.



38 Chapter 3: montaj MAGMAP Filtering Tutorials



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.



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



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:



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



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 (Bhattacharya, 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



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



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.



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



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.



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



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.



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



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:



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



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.



46 Chapter 5: Preparing Grids for the FFT Processing



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.



Chapter 5: Preparing Grids for the FFT Processing 47



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



48 Chapter 5: Preparing Grids for the FFT Processing



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



Chapter 5: Preparing Grids for the FFT Processing 49



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.



50 Chapter 5: Preparing Grids for the FFT Processing



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.



Chapter 6: Specifying FFT Filters 51



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: line 2:



line 3:



A title line for reference only. 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. 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.



52 Chapter 6: Specifying FFT Filters



line 4:



line 5:



line 6+



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. 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. 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



Chapter 6: Specifying FFT Filters 53



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). BTWR – Butterworth Filter



L(k ) =



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



54 Chapter 6: Specifying FFT Filters 1.0 n =2



L(k)



4



16



8



0.5



0.0



0



ko



N



Wavenumber (cycles/ ground_unit)



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



L(r)



h=8 h=4 h=2



1.0 0



N



Wavenumber (radians/ground_unit)



Parameters:



h



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



Chapter 6: Specifying FFT Filters 55



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



Wavenumber (radians/ground_unit)



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.



56 Chapter 6: Specifying FFT Filters



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



0.5



n=2



L(k) 0.0



1



0.5



N



0



k0



k1



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 π L ( θ)) = cos n (α − θ + 2 ) , to reject direction α



π L( θ) = 1 − cos n (α − θ + 2 ) , to pass direction α



Chapter 6: Specifying FFT Filters 57 1.0



n=2



L( θ ) 0.5



0.5



1



0.0 α−π/2



α



α+π/2



θ (direction) 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).



58 Chapter 6: Specifying FFT Filters



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).



Chapter 6: Specifying FFT Filters 59



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



Order of differentiation.



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



Order of differentiation.



r



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



60 Chapter 6: Specifying FFT Filters



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 − rz1 − 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 k0 2



1.0



L(k)



0.5



0.0



0



K0



2K 0



3K 0



N



Wavenumber (cycles/ ground_unit)



Chapter 6: Specifying FFT Filters 61



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



62 Chapter 6: Specifying FFT Filters



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



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.). Output of Pseudo-Gravity is in mgal if ground_unit is metre, or in 0.3mgal if groud_unit is foot.



Chapter 6: Specifying FFT Filters 63



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



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.). 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.



64 Chapter 6: Specifying FFT Filters



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 − + s 0 1 0 ⎠ ⎝ ⎦ ⎣



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



Chapter 6: Specifying FFT Filters 65



Where: I



Geomagnetic inclination



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.



66 Chapter 6: Specifying FFT Filters



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 |