WO2017205914A1 - Material characterisation system and method - Google Patents

Material characterisation system and method Download PDF

Info

Publication number
WO2017205914A1
WO2017205914A1 PCT/AU2017/050514 AU2017050514W WO2017205914A1 WO 2017205914 A1 WO2017205914 A1 WO 2017205914A1 AU 2017050514 W AU2017050514 W AU 2017050514W WO 2017205914 A1 WO2017205914 A1 WO 2017205914A1
Authority
WO
WIPO (PCT)
Prior art keywords
detector
detectors
energy
pulse
spectrum
Prior art date
Application number
PCT/AU2017/050514
Other languages
French (fr)
Inventor
Paul Scoullar
Christopher Mclean
Shane Tonissen
Syed Khusro SALEEM
Original Assignee
Southern Innovation International Pty Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Priority claimed from AU2016902059A external-priority patent/AU2016902059A0/en
Application filed by Southern Innovation International Pty Ltd filed Critical Southern Innovation International Pty Ltd
Priority to AU2017274079A priority Critical patent/AU2017274079B2/en
Priority to BR112018074796-3A priority patent/BR112018074796B1/en
Priority to US16/305,803 priority patent/US11243327B2/en
Priority to CA3025821A priority patent/CA3025821A1/en
Priority to EP17805394.8A priority patent/EP3465180A4/en
Publication of WO2017205914A1 publication Critical patent/WO2017205914A1/en
Priority to ZA2018/08640A priority patent/ZA201808640B/en

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N23/00Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
    • G01N23/02Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material
    • G01N23/06Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and measuring the absorption
    • G01N23/083Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and measuring the absorption the radiation being X-rays
    • G01N23/087Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and measuring the absorption the radiation being X-rays using polyenergetic X-rays
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V5/00Prospecting or detecting by the use of ionising radiation, e.g. of natural or induced radioactivity
    • G01V5/20Detecting prohibited goods, e.g. weapons, explosives, hazardous substances, contraband or smuggled objects
    • G01V5/22Active interrogation, i.e. by irradiating objects or goods using external radiation sources, e.g. using gamma rays or cosmic rays
    • G01V5/224Multiple energy techniques using one type of radiation, e.g. X-rays of different energies
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01JMEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
    • G01J3/00Spectrometry; Spectrophotometry; Monochromators; Measuring colours
    • G01J3/02Details
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01JMEASUREMENT OF INTENSITY, VELOCITY, SPECTRAL CONTENT, POLARISATION, PHASE OR PULSE CHARACTERISTICS OF INFRARED, VISIBLE OR ULTRAVIOLET LIGHT; COLORIMETRY; RADIATION PYROMETRY
    • G01J3/00Spectrometry; Spectrophotometry; Monochromators; Measuring colours
    • G01J3/28Investigating the spectrum
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/25Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
    • G01N21/255Details, e.g. use of specially adapted sources, lighting or optical systems
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/84Systems specially adapted for particular applications
    • G01N21/85Investigating moving fluids or granular solids
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/84Systems specially adapted for particular applications
    • G01N21/88Investigating the presence of flaws or contamination
    • G01N21/93Detection standards; Calibrating baseline adjustment, drift correction
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N23/00Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
    • G01N23/02Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material
    • G01N23/04Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and forming images of the material
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N23/00Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
    • G01N23/02Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material
    • G01N23/06Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and measuring the absorption
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N23/00Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
    • G01N23/02Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material
    • G01N23/06Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and measuring the absorption
    • G01N23/083Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and measuring the absorption the radiation being X-rays
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N23/00Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
    • G01N23/02Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material
    • G01N23/06Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and measuring the absorption
    • G01N23/10Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and measuring the absorption the material being confined in a container, e.g. in a luggage X-ray scanners
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N23/00Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00
    • G01N23/02Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material
    • G01N23/06Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and measuring the absorption
    • G01N23/12Investigating or analysing materials by the use of wave or particle radiation, e.g. X-rays or neutrons, not covered by groups G01N3/00 – G01N17/00, G01N21/00 or G01N22/00 by transmitting the radiation through the material and measuring the absorption the material being a flowing fluid or a flowing granular solid
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V5/00Prospecting or detecting by the use of ionising radiation, e.g. of natural or induced radioactivity
    • G01V5/20Detecting prohibited goods, e.g. weapons, explosives, hazardous substances, contraband or smuggled objects
    • G01V5/22Active interrogation, i.e. by irradiating objects or goods using external radiation sources, e.g. using gamma rays or cosmic rays
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/25Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
    • G01N2021/258Surface plasmon spectroscopy, e.g. micro- or nanoparticles in suspension
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/84Systems specially adapted for particular applications
    • G01N21/88Investigating the presence of flaws or contamination
    • G01N21/93Detection standards; Calibrating baseline adjustment, drift correction
    • G01N2021/933Adjusting baseline or gain (also for web inspection)
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2223/00Investigating materials by wave or particle radiation
    • G01N2223/30Accessories, mechanical or electrical features
    • G01N2223/345Accessories, mechanical or electrical features mathematical transformations on beams or signals, e.g. Fourier
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2223/00Investigating materials by wave or particle radiation
    • G01N2223/40Imaging
    • G01N2223/423Imaging multispectral imaging-multiple energy imaging
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2223/00Investigating materials by wave or particle radiation
    • G01N2223/50Detectors
    • G01N2223/501Detectors array
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2223/00Investigating materials by wave or particle radiation
    • G01N2223/50Detectors
    • G01N2223/505Detectors scintillation
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2223/00Investigating materials by wave or particle radiation
    • G01N2223/60Specific applications or type of materials
    • G01N2223/643Specific applications or type of materials object on conveyor

Definitions

  • This invention relates to a system and method for material characterisation.
  • the invention may be applied to a wide range of materials, including but not limited to discrete items of material, a particulate substance, a non-particulate substance, mineral ore, coal, cement, mineral sand, fluid, hydrocarbon fluid, soil, agricultural material, food or beverage, metals, plastics, glass, timber, pharmaceuticals, contaminated material, hazardous material, waste, refuse, sewage, effluent, recyclable material, construction debris, natural materials, manufactured products, clothing, electronics, medical devices, automotive components, aerospace components, dentistry devices, and dentistry components.
  • the invention has particular application to characterising coal or mineral ore travelling on a conveyor.
  • the invention also has particular application to determining elemental, mineral ogical or density characteristics of a material, for example the percentage of iron in an iron ore, but in some embodiments may also be applied to determination of other physical characteristics, such as particle or grain size, and surface properties such as texture.
  • the present invention relates to a system and method which has advantages in a combination of high throughput, accuracy, resolution, elemental sensitivity, safety and operating cost.
  • a source of incident X-rays configured to irradiate at least part of the material
  • one or more detectors adapted to detect radiation emanating from within or passing through the material as a result of the irradiation by the incident radiation and thereby produce a detection signal
  • one or more digital processors configured to process the detection signal to characterise at least part of the material
  • the one or more detectors and one or more digital processors are configured to characterise at least part of the material by performing energy resolved photon counting X-ray transmission spectroscopy analysis.
  • a method of characterising at least part of a material comprising the steps of: irradiating at least part of the material with incident X-rays from an X-ray source; detecting radiation emanating from within or passing through the material as a result of the irradiation by the incident X-rays;
  • the words “energy resolved photon counting X-ray transmission spectroscopy” are synonymous with “multi energy X-ray transmission spectroscopy” and “multiple X-ray absorptiometry”.
  • photon counting mean processing and analysing individual packets of radiation, in particular, X-ray photons, in the course of processing the detection signal, and does not necessarily imply obtaining an actual enumeration.
  • packets in relation to incident radiation includes individual massless quantum particles of X-rays (X-ray photons).
  • energy spectrum in relation to a particular detector refers to a generation of energy values of the individual packets of radiation emanating from or passing through the part of the items under investigation as detected over a time interval from the particular detector, which energy values can comprise values over a range, typically continuous, and may be represented as a histogram of detection counts versus a plurality of defined energy bins, the number of bins representing the desired or achievable energy resolution and constituting at least 10 bins but preferably more than 50, 100 or 1000 energy bins.
  • the word "material” refers to any material of solid, liquid or gaseous composition, including but not limited to discrete items of material, a particulate substance, a non-particulate substance, mineral ore, coal, cement, mineral sand, fluid, hydrocarbon fluid, soil, agricultural material, food or beverage, metals, plastics, glass, timber, pharmaceuticals, contaminated material, hazardous material, waste, refuse, sewage, effluent, recyclable material, construction debris, natural materials, manufactured products, clothing, electronics, medical devices, automotive components, aerospace components, dentistry devices, and dentistry components.
  • characterisation means measurement, classification, assessment, or analysis of one or more features or properties of said material, including but not limited to physical, chemical, optical, magnetic or other material properties.
  • Figure 1 shows a material characterisation system according to one embodiment of the invention installed on a conveyor adapted to ore characterisation, adapted to high- capacity conveyors and having an X-ray source above the conveyor with detectors below .
  • Figure 2 shows a material characterisation system according to another embodiment of the invention installed on a conveyor adapted to ore characterisation, adapted to high- capacity conveyors and having an X-ray source below the conveyor with detectors above. S.
  • Figure 3 shows the system of Figure 1 or Figure 2 installed on a sample conveyor.
  • Figure 4 shows a material characterisation system similar to Figure 1 or 2 installed to characterise material falling off the end of a conveyor.
  • Figure 5 shows a detector module comprising a plurality of detector arrays according to an embodiment adapted for ease of maintenance and replacement.
  • Figure 6 and 7 show a detector module variant in which the detector arrays are angled for optimum detection.
  • Figure 8 shows a source module with shielding and an IP rated enclosure according to an embodiment.
  • Figure 9 shows a more detailed view of the detection system and processing electronics according to an embodiment.
  • Figure 10 illustrates a flowchart for a method of effective atomic number (effective Z) processing of the full energy spectra computed by the pulse processing electronics according to an embodiment.
  • Figure 11 is a graph illustrating the removal of pileup of two pulses from a spectrum according to an embodiment.
  • Figure 12 is a graph illustrating the removal of pileup of two and three pulses from a spectrum according to another embodiment.
  • Figure 13 is a graph illustrating the partial removal of pileup of two and three pulses from a spectrum when assuming only 2 pulse pileup according to an embodiment.
  • Figure 14 is a graph illustrating the shape of the spectral smoothing filter when using a rectangular window or a raised cosine pulse window according to an embodiment.
  • Figure 15 illustrates how data is arranged and built up into an image of a scanned sample, prior to further post processing and image display according to an embodiment.
  • Figure 16 is a graph illustrating the un-calibrated received spectra from a plurality of detectors according to an embodiment.
  • Figure 17 is a graph illustrating a set of detector gains calculated based on a calibration procedure according to an embodiment.
  • Figure 18 is a graph illustrating the received spectra from a plurality of detectors after setting the digital gain of the detectors based on the detector gains illustrated in Figure 17.
  • Figure 19 illustrates results of the effective Z interpolation process for a 10% transmission case according to an embodiment.
  • Figure 20 illustrates effective Z plotted against intensity (percent transmission) for a range of material samples tested according to an embodiment.
  • Figure 21, 22, 23 illustrate a detector module comprising an array of detectors according to embodiments in which a linear (or curved) array of scintillation crystals is coupled to an array of pulse producing elements in the form of silicon photomultipliers.
  • Figure 24 and 25 illustrate a detector module comprising single scintillation crystals individually coupled to an array of pulse producing elements in the form of silicon photomultipliers by means of an optical coupling layer interposed between the scintillation crystals and silicon photomultipliers according to an embodiment.
  • Figure 26 illustrates a detector subsystem converting photons into voltage pulses for pulse processing according to an embodiment.
  • Figure 27 illustrates an example of the formation of clusters, where single tiles are ignored according to an embodiment.
  • Figure 28 illustrates an example of effective Z processing steps according to an embodiment.
  • Figure 29 illustrates a table relating to an edge mask L(c) indexed on columns according to an embodiment.
  • Figure 30 illustrates the behaviour of the moving average as it transitions over an edge according to an embodiment.
  • Figures 31 to 37 show aspects of a number of embodiments of the invention adapted to lower capacity conveyors.
  • X-rays passing through matter interact with that matter via a number of modalities including: scattering off crystal planes, causing fluorescence X-ray emission from within the electron structure of the elements; and, scattering off nano-scale structures within the material being scanned.
  • modalities including: scattering off crystal planes, causing fluorescence X-ray emission from within the electron structure of the elements; and, scattering off nano-scale structures within the material being scanned.
  • the system of one of the embodiments described below provides for a detection system capable of estimating the energy of the individual X-ray photons received at the detector. This is achieved using a single detector array per X-ray source, with each of the detectors in the array constructed from an appropriate detector material coupled to a photomultiplier, producing an analog signal comprising a series of pulses - one pulse for each detected X-ray, which may or may not be overlapping when received at the detector.
  • a pulse processing system is then used to generate a histogram for each single detector.
  • This histogram comprises of a count of the number of X-rays falling into each histogram bin in a given time interval.
  • the histogram bins represent a range of energies of the received X-rays, and the histogram is therefore the energy spectrum of the received X- ray beam.
  • Effective Z may be used, for example, to distinguish different contents in a material, enabling assessment of mineral ore or coal quality, or waste material characteristics.
  • the system may also extract other information from the energy spectrum, for example material density, in order to aid in characterisation of the material, which may be implemented in software in the digital processors.
  • Figure 1 shows a material characterisation system according to an embodiment of the invention, installed on a conveyor 1 adapted to convey material in the direction of the arrow A.
  • the system comprises a housing 2 encompassing conveyor 1. Coal is seen progressing along conveyor 1 through aperture 3 in housing 2.
  • Housing 2 houses a source of incident X-rays in the form of an X-ray generator 4 positioned above conveyor 1 which emits a fan-shaped beam of X-rays 5 perpendicular to the direction of travel providing an irradiation region.
  • a plurality of modules 6 is positioned within housing 2 beneath conveyor 1 so as to intersect and detect radiation emanating from or transmitting through the coal radiated by beam 5.
  • Each module 6 comprises a plurality of detectors in an array, wherein each detector comprises a detector subsystem, wherein digital processor electronics may be provided within the modules 6 or separately and in communication with modules 6 as described below.
  • a suitable X-ray generator 4 in this embodiment for effective Z determination is obtainable from Industrial Control Machines SA (www.icmxray.com), for example the SitexTM range. Suitable X-ray generators may be selected for their power and physical design according to the application to specific material and integration requirements.
  • Modules 6 are provided in a plug-in modular constructions, allowing easy removal and replacement on-site, as is X-ray generator 4.
  • a control panel and user interface 7 resides on a side of housing 2, and the system is equipped with wireless communication antenna 8 for remote communication control.
  • the user interface 7 may comprise a display screen or indicator lights to provide visual outputs, system status of emitted results.
  • Wireless communication antenna 8 may connect with a tablet or other portable computer an operator standing nearby.
  • wireless communication antenna 8 may be adapted to communicate with an on-site control centre or over a remote communication system such as the Internet to an off-site remote control centre. Alternatively or in addition, a wired network communication system may be employed.
  • the system may be installed on a main conveyor 1 with the aim of characterising substantially all of the coal produced at the mine, handled by a processing facility, or shipped to or from a port or customer.
  • FIG. 2 an alternative configuration is shown with X-ray generator 4 disposed beneath the conveyor and detectors in modules 6 disposed above the conveyor.
  • This configuration may be preferable in some circumstances to avoid interference caused by X-rays intercepting support structures of conveyor 1, by providing a well-positioned source 4 between and avoiding the under conveyor support structures.
  • a sample analysis facility may be provided off-line comprising housing 2 installed on a sample conveyor 9 which feeds coal samples into the system for analysis placed manually on conveyor 9 through hopper 10.
  • Figure 4 shows another embodiment with the system installed to characterise coal falling off the end of a conveyor 11 in the direction of arrow B.
  • Housing 2 is similar or identical to the embodiment of Figure 1 but is installed with the aperture disposed horizontally so as to catch the coal falling of the end of conveyor 11.
  • a second conveyor 12 is disposed beneath aperture 3 to catch and continue transport of the coal.
  • This embodiment may be designed to provide analysis in applications where the material of the conveyor interferes somewhat with analysis of the detected radiation and where it is therefore advantageous that only the coal is interposed between the X-ray generator and the detectors in the irradiation region.
  • the material may be falling off the end of a conveyor onto a storage pile or into a cargo hold of a ship or other long-range means of transport.
  • the system is modularised to facilitate manufacture, maintenance and replacement. Modularity reduces service component weight, minimises on-site work, allows maintenance tasks to be undertaken by unskilled workers, provides logical cost breakdowns for replacement parts and improves reliability by sealing sensitive components together in each simple module.
  • a detector module 20 is shown in a cutaway illustration as comprising several detector boards 21, each detector board 21 comprising a linear array of detector elements as described below disposed behind collimators 26 are constructed from lead which allow only a narrow slit of x-radiation to impinge on the detector elements.
  • the detector boards in 21 are mounted on an aluminium base of the precision manufactured enclosure, allowing reliable alignment of the detector boards and heat transference.
  • the module also comprises aluminium side support rails 24 which may be provided as an extension of the aluminium base or otherwise rigidly and accurately fixed to the aluminium base.
  • a top cover 23 of the enclosure is composed of 5 mm thick aluminium, which is effectively opaque to X-rays, except for a slit composed of thinner aluminium ( ⁇ lmm) or HDPE which is X-ray translucent disposed above collimators 26, to allow the emanating radiation to enter the detector module and pass through the slits of collimators 26.
  • a cooling circuit passes through custom heatsinks 25 on each detector board 21, and liquid coolant may be circulated through coolant connectors 22 .
  • an alternative embodiment of detector module comprises a taller module which allows each detector board 21 to be angled so as to face the fan -shaped X-ray beam at right angles.
  • Detector arrays 27 can also be seen in the crosssectional representation of Figure 7.
  • the aluminium base 24 and side rail are formed from the same sheet. This embodiment is particularly suitable for the variants described in figures 31-38 where the X-ray source is close to the detector.
  • a source module is shown with similar module enclosure 23 and aluminium base 24.
  • the X-ray source is contained within a source enclosure 28 composed of 5 mm thick lead walls to provide X-ray shielding, cast over an internal aluminium pipe (not shown) to provide rigidity.
  • a thicker 10 mm lead region 36 is provided surrounding the beam exit point of the contained X-ray source and incorporates a collimating slit 35 (which may be an HDPE or thin aluminium window) to provide a substantially planar beam of source X-rays exiting the enclosure for irradiation of the material in the irradiation region.
  • Cooling holes 33 provide internal air circulation and cooling.
  • Source enclosure 28 comprises precision machined end caps 32 and guide tabs 34 for accurate positioning within module enclosure 23. Once source enclosure 28 is slid into position, door 30 is closed.
  • a partitioned electronics and power section 31 is provided also within module enclosure 23, with internal connections 29 to the X-ray source. Cooling may also be provided through coolant passing through multipurpose connectors such as 22 above for the detector module (not shown).
  • the source and detector are positioned around conveyor belt 1 with source module 38 above conveyor belt 1 and detector module 39 below conveyor belt 1.
  • Source module 38 and detector module 39 are accurately positioned on a precision engineered first chassis being "tech" chassis 4.
  • Tech chassis 41 is seated on a second chassis being an external or support chassis 42 through vibration damped mounts such as rubber mounts or wire rope isolators (not shown).
  • Support chassis 42 in turn is permanently bolted to conveyor base 40.
  • Text chassis 41 and support chassis 42 are distinguished in the diagram by forward and backward diagonal lines. The use of a separate "floating" tech chassis 41 enables tight dimension tolerances which are required between source module 38 and detector module 39.
  • Source module 38 and detector module 39 may be removed from tech chassis 41 when required.
  • Figure 32 shows a variant of the embodiment of Figure 31, equipped with a C- shaped support chassis 42 and a similar floating tech chassis 41.
  • the entire tech chassis 41 may be slid out in the direction shown.
  • Figure 33 shows another variant, with a similar C- shaped support chassis and a rotatably attached tech chassis 41.
  • source module 38 is shown underneath a conveyor 1 and detector module 39 is shown above.
  • the tech chassis may be rotated in the direction of the arrow shown to provide good access for removal of entire tech chassis or individual modules by crane as required.
  • Figure 34 shows a cross-section of the embodiment of Figure 31, to illustrate passage of the X-ray beam 44 from source module 38 into detector module 39, and to illustrate the optional inclusion of beam shielding elements in the form of lead plates 43 disposed on a shielding chassis structure 45 which also sits on support chassis 42.
  • Figure 35 shows in exploded view the components of shielding chassis structure 45, lead baffle plates 46, lead top plate 48 (passage slit not shown) and side plates 47 for disposing beneath the conveyor 1.
  • the beam shielding may be necessary because of potential exposure of personnel, depending on the geometry of the system and any containing structures.
  • the source and detector modules may be disposed within a weatherproof enclosure 56 with sliding door 50, shown here on a conveyor 1 equipped with a hood 52.
  • a hood 52 For elevator conveyors, and access platform 57 may be provided with user interface 58 located on platform 57.
  • the electronic and X-ray detection components will now be described.
  • One or more digital processors process the received X-ray pulses from the detector arrays within the detector modules.
  • the digital processors may: i. Reside on the same boards as a detector array within the modules. ii. Reside on separate hardware, housed within the module chassis or outside the chassis in auxiliary equipment (not shown)
  • iii Form part of a Host controller responsible for overall operation of the system housed inside or outside the chassis.
  • Figure 9 shows a more detailed view of the detection system and processing. This figure shows the steps for a single detector. The effective Z may utilise, and image post processing will require, access to the spectra from all detectors.
  • an image of the coal is formed from rows of pixels as the coal passes through the detector region.
  • the invention does not require the formation of an image and a system comprising a single detector is within the scope of the broadest aspect.
  • the detection system and processing electronics comprises:
  • a detector subsystem (301) for each individual detector (with N such subsystems for a lxN detector array), the detector subsystem comprising:
  • Detector material for detecting the incident X-rays (300) and converting each detected X-ray to a light pulse
  • a photomultiplier for receiving and amplifying the incident light pulses into an analog signal comprising pulses (312) that may or may not overlap c.
  • Appropriate analog electronics which may include filtering d.
  • High rate pulse processing (305) for each detector subsystem (301), for example the pulse processing systems disclosed in US Patent No. 7383142, US Patent No. 8812268 and WO/2015/085372, wherein the pulse processing comprises:
  • a gate signal source (306) which outputs a gate signal (314) at a regular pre- configured interval.
  • the gate interval is a constant short interval that determines the histogram accumulation period. This gate interval also determines the pixel pitch in the resulting X-ray image.
  • another appropriate means may be used to control and synchronise the timing of energy histogram collection across all detectors. For example, a suitably precise network timing signal may be used instead of the gate signal.
  • Calibration System (307) which receives input from appropriate analog and digital signals and then communicates the desirable calibration parameters back to the various processing blocks. The calibration system performs:
  • Effective Z computation (308) which takes the computed energy spectra in each detector during each gate interval and determines the effective Z of the sample. This in turn leads to the production of an effective Z image.
  • Intensity image generation including:
  • High penetration or high contrast image determined by integration of selected energy bands from the full energy spectrum.
  • Image post processing and display (311) with features that may include one or more of the following:
  • the images produced for display comprise a number of data elements recorded for each of N detector elements (501) and for each gate interval (500).
  • the data obtained for detector i during gate interval j is used in the production of effective Z, intensity and high penetration/high contrast images as shown in Figure 9.
  • a number of elements are recorded in each pixel (502), including one or more of: 1.
  • Figure 15 illustrates how this data is arranged and built up into an image of the scanned sample, prior to further post processing and image display.
  • the detector subsystem utilises a scintillator (such as phosphor) coupled to an array of PIN diodes to convert the transmitted X-ray into light, and subsequently into an electrical signal.
  • a scintillator such as phosphor
  • the detector subsystem of this embodiment comprises: a) A detector material
  • the detector material may be of dimensions X by Y by Z, or some other shape.
  • the photomultiplier may be a silicon photomultiplier (SiPM) and the coupling means may be a form of optical grease or optical coupling material. It may be desirable to use a form of bracket or shroud to hold the detector in position relative to the photomultiplier.
  • the photomultiplier requires appropriate power supply and bias voltage to generate the required level of amplification of the detected signal. In this embodiment, typically each detector dimensions approximately l-3mm in width and depth, with a height of approximately 2- 6mm, and inter-detector pitch of approximately 2-200mm.
  • a large number of single element detector subsystems may be required to produce each detector array. It may be desirable to group these in an appropriate way, depending on the specific X-ray scanner requirements. Individual elements of detector material may be grouped into a short array of M detectors. Small groups of M detector elements may be mounted onto a single detector board, for example 2, 4 or more groups of M onto one board. The full detector array is then made up of the number of detector boards required to achieve the total number N of detector elements per array.
  • detector subsystems can be arranged in a number of different configurations including: curved arrays of 1 x N devices; linear arrays of 1 x N devices; square or rectangular arrays of N x M devices; or L-shaped, staggered, herringbone or interleaved arrays.
  • three curved arrays of 1 x 24 devices are utilised, disposed within the detector ring housing 31, however, any chosen distribution of detectors and detector boards may be used according to the desired geometry of the system and population of detectors around the detector ring housing 31
  • One example of a detection device used to convert incoming radiation photons into and electrical signal, is the combination of a scintillation crystal, coupled to a silicon photomultiplier (SiPM) or multi-pixel photon counter (MPPC).
  • SiPM silicon photomultiplier
  • MPPC multi-pixel photon counter
  • a scintillation crystal such as LSYO (1701) is used to convert the incoming radiation photon (1700) into UV photons (1703).
  • LYSO scintillation material the peak emission of UV photons occurs at 420 nm, other scintillation material such as those listed in Table 1 may have different emission peaks.
  • a multi -pixel photon counter, or silicon photomultiplier (1704) with sensitivity in the UV region may be used to detect these photons and produce an electrical signal.
  • Figure 21 depicts a linear array of LYSO scintillation crystals (1600), indicative of how single detection devices can be joined together to form a linear array.
  • the individual LYSO crystals (1600) have a cross section of 1.8mm and a height of 5mm, the individual LYSO crystals (1600) are wrapped around the sides in a reflective material to assist in collecting all the UV photons.
  • the pitch of this exemplary array is 2.95mm, the length is 79.2mm and the width of the array is 2.5mm.
  • Figures 22 and 23 depict a detector array from a top view and side view respectively, comprising the linear array of LYSO crystals depicted in Figure 21 coupled to an electrical pulse producing element (1604) on substrate (1605).
  • the electrical pulse producing element may comprise a silicon photomultiplier (SiPM).
  • Enhanced specular reflector (ESR) or aluminium or other reflective foil (1601) is disposed around side surfaces of the scintillation crystals to direct the scintillation photons onto the silicon photomultiplier material (1604) and prevent light leakage (cross-talk) between adjacent detection devices.
  • optical coupling (1606) may be interposed between the LYSO crystals and SiPM, and may comprise any number of known suitable materials, for example, a thin layer of optically transparent adhesive.
  • the detector arrays may be curved rather than linear as shown in Figures 21 and 23 to appropriately fit inside different housing shapes.
  • scintillation crystals (1607) may be individually coupled to electrical pulse producing elements (1604), as depicted in figures 23 and 24. Coupling may be achieved by a number of methods, for example interposing an optically transparent adhesive film (1609) or optical coupling material between the scintillation crystals (1607) and electrical pulse producing elements (1604), where the electrical pulse producing elements (1604) may comprise SiPMs or an MPCC. Coupling may be performed by a 'pick and place' assembly machine to individually align and couple components and coupling material. Scintillation crystals may be wrapped in a reflective material such as a foil or ESR material (1608) to aid in the capture of photons.
  • a reflective material such as a foil or ESR material (1608
  • the LSYO crystals (1600, 1607) may typically have a cross-section (width) approximately l -2mm, a depth of approximately l-2mm, and height of approximately 3 -5mm, where the reflective or ESR film (1601, 1608) is approximately 0.05mm - 0.1mm thick.
  • the cross-section is 1.62mm
  • the depth is 1.24mm
  • the height is approximately 4.07mm
  • the ESR film is 0.07mm thick.
  • the cross sectional area of the scintillator material is preferably greater than 1mm square, and may be greater than 2mm square and less than 5mm square.
  • detector subsystem design uses a scintillator which is compact, robust, cost effective and non-hygroscopic
  • other detector subsystems can be considered. These include detector subsystems which use alternate inorganic or inorganic scintillator materials, the characteristics of some such material are provided in Table 1. Other mechanisms for converting radiation photons into electrical signals could also be considered for the detector subsystem.
  • High Purity Germanium HPGe: Achieves 'gold standard' resolution of 120 eV for the Fe55 X-ray line at 5.9 keV, detectors can be made > 10 mm thick thus detect high energy X-rays up to many 100s of keV.
  • SDD Silicon Drift Diode
  • PIN Diodes The detection efficiency for X-rays up to 60 keV is substantially higher than SDD detectors and falls off to approximately 1% for X-ray energies above 150 keV. These detectors can be operated at room temperature, however, resolution improves with cooling, resolution of the 5.9 keV line is ⁇ 180 eV,
  • Cadmium Zinc Telluride Is a room temperature solid state radiation detector used for the direct detection of mid-energy X-ray and Gamma-ray radiation. It has a detection efficiency for 60 keV X-ray very close to 100% and even for X-rays photons with energies of 150eV the detection efficiency remains greater than 50%.
  • Cesium Iodine (CsI(Tl)) This is a scintillation material used for detection of X-rays in medical imaging and diagnostic applications. The scintillation material is used to convert the X-ray into photons of light which are generally then converted into an electrical signal either by a photomultiplier tube. Csl is a cheap and dense material and has good detection efficiency of X-rays and Gamma-ray to many 100s of keV.
  • the material characterisation system comprises a large number of individual detectors. While each detector and associated electronics is ideally designed to have identical response to incident radiation, in practice this will not be possible. These variations between detectors result in detector to detector variation in energy spectrum output. By properly and fully calibrating the detection system, the energy spectra output from the pulse processing digital processors can be appropriately calibrated so they represent received X-ray intensity in known narrow energy bins.
  • Detector pulse calibration is used to identify the pulse characteristics for each detector required by the pulse processing system. The exact parameters required may vary, depending on the detection system. For typical applications using the pulse processing method disclosed in US Patent No. 7383142 and US Patent No. 8812268, the pulse is modelled as an averaged dual exponential of the form:
  • t 0 is the pulse time of arrival
  • T a is the pulse averaging window
  • A is a pulse scaling factor related to the pulse energy.
  • the processing requires the two parameters a and ⁇ , and the pulse form p(t) which can be obtained via an appropriate calibration method, or from knowledge of the design of the detection subsystem.
  • a suitable method for estimating ⁇ , ⁇ and p(t) from received pulses is described below.
  • Each detector subsystem combined with an analog to digital converter, will have slightly different characteristics due to manufacturing variations. As a result of such component variations, the energy spectra will be scaled differently. Variations other than gain scaling are handled within the Baseline Offset Calibration or the Energy Calibration.
  • the objective of the gain calibration is to achieve alignment of the energy spectra output by the pulse processing electronics across all detectors.
  • the need for absolute accuracy may be reduced or eliminated if per detector energy calibration is applied.
  • Gain calibration may be achieved in a number of ways. The following approach may be applied: Setup a known X-ray source.
  • a material with particular characteristics can be inserted into the beam.
  • lead (Pb) has a known absorption edge at 88 keV.
  • the required gain for each detector is then computed as the ratio of median location to the specific detector feature location.
  • the median or other suitable reference e.g. maximum or minimum
  • Median is chosen so some channels are amplified, and some are attenuated, rather than attenuating all channels to the minimum amplitude.
  • the gains are then applied to each detector channel.
  • the gain may be applied as an analog gain, digital gain, or combination of the two, depending on particular system functionality. For best results, at least part of the gain is digital gain, where arbitrarily fine gain variation can be achieved.
  • the detection subsystem there may be a number of detector cards, each with a number of detector subsystems.
  • the total number of detector subsystems may be several thousand or more. Results from one example of such a detector board are presented here.
  • the example board comprises 108 detectors, with LYSO used as the scintillator material in this case. These detectors are packed into linear arrays of 27 detectors. Each detector board then uses 4 x 27 detector arrays to achieve a total of 108 detectors.
  • each detector is placed above a SiPM, and it is the SiPM that detects and amplifies the emitted photons.
  • the detectors are coupled to the SiPM via an optical grease.
  • the gain of each SiPM is determined by the bias voltage applied, and the SiPM breakdown voltage.
  • Each detector subsystem may have a slightly different baseline level, as measured at the output of the analog to digital converter.
  • the baseline is estimated and removed. Any suitable method can be used including, for example:
  • the pulse processing electronics will produce an energy spectrum that is uncalibrated. That is, the output will comprises a number of counts in a set of histogram bins, but the exact energy of those histogram bins is unknown. In order to achieve accurate effective Z results, knowledge of the energy of each bin is required. [0100] This is achieved as follows:
  • a source with known spectrum peaks is a Bal33 source, with spectral peaks at 31, 80, 160, 302 and 360 keV
  • A, B and C are constants determined from the measured Bal33 spectrum. This formula is inverted to define Energy as a function of Histogram Bin expressed in terms of the same A, B and C.
  • baseline offset is expressed as a function of count rate.
  • any suitable method can be used for this calibration, including injecting a known source spectrum of variable count rate, and recording the spectrum shift as count rate increases.
  • the source has a narrow energy band so the shift can be clearly measured, and also variable energy so the offset can be calibrated as a function of energy if desirable.
  • the residual spectrum is measured with a large mass of material in the beam, sufficient to completely block the X-ray beam, such as a large thickness of steel. In practice, a small level of energy still reaches the detector array, whether from scatter or other mechanisms, and this residual spectrum must be measured so it can be removed from the received spectra during normal operation.
  • the residual spectrum is then measured by averaging the received spectra for a number of gate intervals with the blocking mass in the beam.
  • the pileup parameters can be calibrated in several ways, for example:
  • i Use a narrow energy source, where the energy and count rate can be varied. ii. Measure the received spectrum as the source energy and count rate are varied. iii. Directly measure the ratio of received 2-pulse and 3-pulse pileup to the main signal peak.
  • a high rate pulse processing system (305), such as those disclosed in US Patent No. 7383142, US Patent No. 8812268 or WO/2015/085372, is allocated to each detector subsystem, to perform the following operations on the digitised pulse signal output from the analog to digital converter: a) Baseline tracking and removal, or fixed baseline removal.
  • the intensity value is computed from the energy spectrum generated for each detector i at each gate interval j according to:
  • the thresholds may be pre-set or user configurable.
  • intensity images with varying contrast are generated based on integrating the received spectrum across different energy bands.
  • the system can only utilise the broad energy range inherent in the detector material.
  • arbitrary energy ranges can be used to generate associated intensity images in that energy range.
  • Specific energy ranges can then be defined in order to best isolate and display particular material types, with energy ranges tuned, for example, for organic material, inorganic material, or light, medium or heavy metals.
  • the high contrast / high penetration images are generated for each detector i at each gate interval j according to:
  • the energy band may be user defined or pre-configured. One, two or more different energy bands may be configured to enable the user to select between images of interest.
  • the effective Z processing involves the use of full energy spectra computed by the pulse processing electronics, combined with the energy calibration, to compute an estimate of the effective Z of the sample material.
  • the effective Z processing is performed for every detector, and for each detector proceeds as follows (so for a 1 x N detector array, this process is repeated N times). To reduce computational requirement, the effective Z processing is only performed for received detectors i and gate intervals j that are not declared either impenetrable or empty.
  • the FFT is the first stage of the pileup reduction, which is not required if data compression has already been achieved using an FFT.
  • the pileup reduction can be achieved with a suitable algorithm as outlined below.
  • FFT domain phase shift (404) in order to achieve a desired lateral shift of the energy spectrum. This step has been found to be desirable where a count rate specific baseline shift exists. Note: multiplication by a linearly increasing (with FFT bin) phase term in the FFT domain results in a lateral shift after iFFT. The extent of the lateral shift is determined by the slope of the linear increase.
  • a frequency domain window (405). This window can be used to design a desired smoothing of the energy spectrum.
  • the window design process is outlined below. A good window has been designed to achieve a smooth filtering of the energy spectrum. Filtering of the noise in the energy spectrum allows the possibility of using a reduced number of energy bins in the effective Z computation for overall improvement in computational efficiency.
  • the zero padding inserts data that was truncated after the FFT. It is not essential to insert zeros for all truncated bins. For example, padding less zeros can produce a smaller FFT buffer which it is more computationally efficient to compute the IFFT.
  • the elements N+2 to 2N of the FFT output are the complex conjugate of the elements 2 to N.
  • N+l will be one of the elements set to zero by the zero padding.
  • N is the number of gate
  • E is the energy level of the X-rays.
  • the mass attenuation constants for a given effective Z and given energy define the extent to which the given material Z will attenuate X-rays of energy E.
  • the intensity of received energies at a particular energy will be given by:
  • Mass attenuation data is available at a finite (small) number of energies, perhaps every 10, 20 or 50 keV, whereas the energy spectra created by the method disclosed in this disclosure may be generated at energy spacing as little as 1 keV or even less. In practice a finite number of these energy values will be selected for use in the effective Z computation.
  • mass attenuation table In order to achieve a smooth mass attenuation table at all energies in the energy spectrum, data for intermediate energies for each Z are obtained using cubic spline interpolation or other suitable interpolation method. The mass attenuation values as a function of energy are considered sufficiently smooth that a cubic spline is a good interpolation method to apply.
  • C(Z) is the cost function
  • w(E) represents weights for each sum of the squared errors between the received spectrum and the predicted spectrum.
  • Z there is no particular requirement for effective Z to be integer, and in fact the mass attenuation table may contain values for non-integer values of Z representing composite materials. However, it is clearly not possible to represent a continuum of possible Z values in a finite table. In order to compute Z to arbitrary precision, it is possible to interpolate the cost function to the required resolution using an appropriate interpolation algorithm. The value of Z chosen is then the value which minimises the interpolated cost function.
  • the cost function C(Z) is a smooth function, and therefore an actual floating point or continuous value of Z which minimises this smooth function can be reliably predicted from the curve via some form of interpolation.
  • step 3 above indicates the cost function is computed for all available Z values in the mass attenuation table.
  • efficient search methods can be applied to reduce the computational requirements. Such methods include one or more of the following:
  • the cost function form has been chosen so as to be relatively insensitive to noise on the spectrum.
  • the first step in the calibration process is to obtain a reference spectrum I 0 (B) at each histogram bin B, with no material in the X-ray beam for the detector(s) to be calibrated. Histogram bins will now be denoted by B rather than E to denote that there is no requirement to calibrate the bins in terms of their exact energy.
  • step wedge of the material. That is, a sample of the material that comprises a series of steps of known thickness x. The largest step is ideally sufficient to reduce the X-ray beam to a level where it can be considered impenetrable. Note: other material samples can be used, but such a step wedge is a convenient form to calibrate against.
  • the tables of Tx(Z, B, x) and R(Z, x) together form the calibration tables that are used to estimate effective Z at each pixel (detector/gate interval). As previously stated, they may or may not be a function of detector also, depending on the equivalence of data from all detectors.
  • the frequency domain window is still required prior to iFFT.
  • the received spectrum will be denoted 1(B), the intensity in a series of histogram bins B.
  • B differentiates from the use of E for the previous section where the histogram bins are calibrated in terms of their actual energy.
  • This expected received spectrum is an interpolated received spectrum between the nearest two calibration spectra, but based on the different attenuation observed at each bin.
  • Other forms of interpolation between spectra could be used, but the material specific interpolation used here provides a superior interpolation result.
  • Z there is no particular requirement for effective Z to be integer, and in fact the self calibration table may contain values for non-integer values of Z representing composite materials. However, it is clearly not possible to represent a continuum of possible Z values in a finite table. In order to compute Z to arbitrary precision, it is possible to interpolate the cost function to the required resolution using an appropriate interpolation algorithm. The value of Z chosen is then the value which minimises the interpolated cost function.
  • the cost function C(Z) is a smooth function, and therefore an actual floating point or continuous value of Z which minimises this smooth function can be reliably predicted from the curve via some form of interpolation.
  • Periodically a new pulse calibration may be performed, although pulse parameters have been found to remain sufficiently constant over time that at most daily or perhaps weekly or monthly recalibration of pulse parameters may be required.
  • step dimensions are chosen with the following in mind:
  • a width of 30cm is used to ensure a large number of detectors could be calibrated with 1 calibration scan.
  • i. Projection mode is used to effectively increase the number of pixels that could be calibrated in
  • Step heights were determined to achieve a reasonably uniform transmission spacing from approximately ⁇ 0.5% up to 95%.
  • Step length is 50mm.
  • the scan speed is 8mm per second, so approximately 6 seconds of data could be collected from each step. This is necessary to ensure very accurate calibration spectra, given the accuracy required from the eventual effective Z processing.
  • step thickness i. A mapping of step thickness to step intensity. This table is used to interpolate any measured spectrum into an equivalent material thickness.
  • a series of calibration spectra including a reference spectrum. Each spectrum represents a material thickness, from which spectra for intermediate material thicknesses can be interpolated.
  • the spectrum for material Z, at histogram B, at thickness x (which here corresponds to one of the defined transmission levels) is denoted 7(Z, B, x).
  • the new interpolated material spectrum is obtained as follows:
  • continuous/floating point effective Z is able to correctly determine Z values around 6 - essentially avoiding anomalies at the edges, or at least pushing them to Z values outside the range of interest.
  • intersection - this is acceptable for lower Z but ceases to hold true when moving into the metals. It may be necessary to interpolate several spectra to obtain each new interpolated material.
  • the cost function C(Z) is a smooth function, and therefore an actual floating point or continuous value of Z which minimises this smooth function can be predicted from the curve via some form of interpolation.
  • the quadratic model is acceptable where the cost function is well behaved, smooth and relatively noise free. Where very short integration times are used to collect the energy spectrum histograms, the cost function can become noisy and may converge on a local minimum due to noise. In this case, and in general, a more sophisticated interpolation model may be required to smooth the cost function and avoid noise effects. This may involve more than 3 points in the interpolation process. [0150]
  • the quadratic model is just a model to ensure a consistent effective Z is obtained for a particular material. It is not intended to be an accurate functional model of the cost function behaviour, and it is not considered necessary. The principle objective is to obtain an estimate of effective Z that is consistent for a particular material, and enables reliable separation of closely spaced materials. The quadratic model achieves this objective.
  • High transmission is more likely to be an organic or lower Z material b.
  • Low transmission is more likely to be a higher Z material, since very large thicknesses of low Z material would be required.
  • Figure 28 indicates an overview of the various optional processing stages that may be implemented in the present method.
  • the tiling algorithm is effectively a block averaging function.
  • the purpose of the tiling algorithm is to average the floating effective Z image over an area (mm 2 ) that represents the smallest region required to be detected of a constant intensity and material composition.
  • the tiling algorithm generates tiles with 50% overlap to ensure that we always capture the object of interest.
  • the tiling algorithm estimates the mean and standard deviation over rectangular tiles in the floating effective Z image.
  • the tile width and height are defined by the user. Tiles are overlapped by 50% in both vertical and horizontal dimensions. Given an image size N r by N c pixels, and a tile dimension T r by T c pixels, the number of tiles in the vertical dimension is (N r /T r )*2. The tile dimensions must be even valued to ensure 50% overlap.
  • the tiling algorithm executes a loop that indexes into each tile and calculates the mean and standard deviation of all pixels in the tile.
  • tile dimensions essentially comes down to a compromise between: 1. The dimensions of the smallest object that must be detected; and
  • the effective Z resolution required The effective Z variance has been observed to reduce almost linearly with number of effective Z pixels averaged, so larger areas yield better effective Z resolution.
  • the clustering algorithm groups tiles that have a common effective Z and are spatially connected.
  • the purpose of the clustering algorithm is to detect objects that span areas larger than the minimum object size as defined by the tile dimensions. Connectedness is defined along edges. Connected tiles are assigned a common cluster ID.
  • the output of the clustering algorithm is a cluster map and a cluster table.
  • the cluster map is a matrix of connected tiles with associated cluster IDs.
  • the cluster table holds information on each cluster ID including the number of tiles in the cluster, and the vertical and horizontal extent of each cluster.
  • the target material detection algorithm is a nearest neighbour classifier.
  • the algorithm classifies individual tiles. There are two steps in the algorithm, training and classification.
  • the training stage establishes a lookup table mapping normalised intensity to floating effective Z for a range of 'target' materials of interest.
  • the lookup table is approximated as a quadratic fit, for which only the quadratic coefficients are stored.
  • the input is the normalised measured tile intensity (Imeas), the measured tile floating effective Z (Zmeas), and a maximum effective Z classification error (deltaZ).
  • the classifier declares positive classification if abs(Ci(Imeas)-Zmeas) ⁇ deltaZ, where Ci is the quadratic function associated with the 1 th target material.
  • the use of both intensity and effective Z in the target profile is an important aspect of this approach.
  • the effective Z is typically not constant with material thickness, and so including the intensity (related to thickness) provides a two dimensional test with far superior discrimination than effective Z alone.
  • Figure 20 shows the effective Z vs intensity for a range of material samples tested, along with the quadratic interpolation. Here the variation of effective Z with intensity is clear. While the example materials are more relevant for a baggage scanning application, the principle is the same in relation to discrimination between mineral types.
  • the purpose of the edge detection algorithm is to ensure that the moving average window does not straddle material boundaries.
  • the edge detection uses amplitude transitions in the intensity image to declare material edges.
  • the input to the edge detection algorithm is the intensity image. Edges are only detected in the horizontal dimension. The reason for not detecting edges in the vertical dimension is that the moving average window only operates in the horizontal dimension. Edges in the intensity image are computed for each detector.
  • a first order gradient operator is used to detect edges.
  • An edge is declared when abs(G) > g where g is a user defined threshold.
  • the purpose of the moving average algorithm is to filter the intensity histograms for each detector so as to increase the effective signal-to-noise ratio.
  • the algorithm generates a filtered intensity histogram a slice k, for each detector, by averaging the measured intensity histograms over a symmetric window centred on slice k.
  • the edge detector plays an important role in ensuring the moving average window does not straddle different materials. If a window overlaps an edge the average is only calculated up to the edge boundaries. The width of the window can be set by the user. On edges, no averaging is performed.
  • Figure 30 illustrates the behaviour of the moving average as it transitions over an edge.
  • One embodiment that may be more computationally efficient is an adaptive moving average approach:
  • any suitable method of high rate pulse processing can be used within the embodiments described herein.
  • the high X-ray flux present in typical X-ray screening systems results in a high pulse count rate, and a high likelihood of receiving overlapping X-ray pulses.
  • Pulse pile-up has long been a problem to contend within applications of high rate radiation spectroscopy.
  • Traditional approaches to pulse shaping use linear filters to shorten pulse duration which can significantly reduce S R and are therefore limited to output rates of a few hundred kc/s.
  • An alternate approach to processing the data from radiation detectors is based on the idea of mathematically modeling data corrupted by pulse pile-up and solving for the required model parameters. By recovering rather than discarding data corrupted by pulse pile-up this technique enables high throughput, low dead-time pulse processing without the traditional loss in energy resolution.
  • the digitised radiation detector time series (y) is modeled as the sum of an unknown number of radiation events (N), with random times of arrival ( ⁇ ), and amplitudes (a), interacting with a radiation detector, that have an expected pulse shape (h) and with a noise process (co).
  • the digitised detector data can be accurately decomposed into the individual component events and the energy of each event determined.
  • Calibration of the detector is the first stage of the algorithm; it takes as input the detector time series data and determines the unit impulse response of the detector (the expected pulse shape from the detector). Refer to Pulse Parameter Calibration for a more detailed summary of the pulse calibration process.
  • the Pulse Localisation stage determines the number of events in the digitised detector data stream and their TOA relative to each other.
  • the detection of events in the digitised detector waveform is accomplished by fitting an exponential model to a fixed number of data points. After the System Characterisation stage the exponential decay of the pulse tail is well characterised.
  • the detection metric (the signal ultimately used to make a decision as to whether a pulse has arrived or not) is formed by fitting an exponential curve to a specified number of data points. This fixed length 'detection window' is run continuously over the digitised detector data and the sum of the squares of the error is computed (this can also be thought of as the sum of the square of the fit residual). This operation results in three distinct modes of operation:
  • Baseline Operation processing data samples when no signal is present. As the data can be quite accurately modeled by an exponent the sum square of the error is at a minimum and remains quite constant.
  • Tail Operation when processing data on the tail of a radiation event the data points are quite accurately modeled as an exponent. Consequently the sum square of the error returns to the same level as the Baseline Operation mode.
  • the final stage of Pulse Localisation is making a decision on the exact number and time of arrival of each of the radiation events in the detector data stream.
  • One approach would be to apply a simple threshold to the detection metric and declare a pulse arrival at the nearest sample to the threshold crossing.
  • a simple threshold crossing is susceptible to noise and only provides ⁇ 0.5 sample accuracy in determining the pulse arrival time.
  • a quadratic peak detection algorithm can be used. Such an approach fits a quadratic to a sliding window of N samples of the detection metric (N maybe equal to 5).
  • the Pulse Energy Estimation stage determines the energy of all the radiation events in the detector data stream. As its input it uses: the a priori knowledge of the detector unit impulse response; the number of events; and their individual time of arrival data.
  • the digitised detector data of equation (18) may also be written in matrix form as:
  • A is an M x N matrix, the entries of which are given by
  • the columns of matrix A contain multiple versions of the unit impulse response of the detector.
  • the starting point of the signal shape is defined by the signal temporal position. For example, if the signals in the data arrive at positions 2, 40, 78 and 125, column 1 of matrix A will have '0' in the first row, the 1 st data point of the unit impulse response in the second row, the 2 nd data point of the unit impulse response in the 3 rd row, etc.
  • the second column will have '0' up to row 39 followed by the signal form.
  • the third column will have '0' up to row 77; the fourth column will have '0' up to row 124 and then the signal form.
  • the size of matrix A is determined by the number of identified signals (which becomes the number of columns), while the number of rows depends on the number of samples in the 'time series'.
  • the final functional stage of the real-time, signal-processing algorithm is the Validation stage. At this stage all the parameters that have been estimated by previous algorithmic stages (pulse shape, number of events, time of arrival and event energy) are combined to reconstruct a 'noise-free' model of the detector data.
  • the method of WO2010068996 for resolving individual signals in detector output data comprises:
  • the detector output data as a digital series (such as a digital time series or a digitised spectrum);
  • modelling the function output according to a model (such as by modelling the function output as a plurality of sinusoids);
  • the signal form may generally be regarded as characterising the interaction between the detector and the radiation (or other detected input) that was or is being used to collect the data. It may be determined or, if known from earlier measurements, calibrations or the like, obtained from (for example) a database.
  • transforming the digital series according to the mathematical transform comprises forming a model of the digital series and transforming the model of the digital series according to the mathematical transform.
  • the method includes determining a plurality of parameters of the transformed signals, such as frequency and amplitude.
  • the transform is a Fourier transform, such as a fast fourier transform or a discrete fourier transform, or a wavelet transform.
  • the transform may be applied somewhat differently to the signal form and digital series respectively.
  • the mathematical transform is the Fourier transform, but the signal form is transformed with a discrete fourier transform and the digital series is transformed with a fast fourier transform.
  • the transform is a Fourier transform and the function is representable as
  • 'corrupt data' some data (which hence is termed 'corrupt data'), as is described below.
  • ' signal' is interchangeable in this context with 'pulse', as it refers to the output corresponding to individual detection events rather than the overall output signal comprising the sum of individual signals.
  • the temporal position (or timing) of a signal can be measured or expressed in various ways, such as according to the time (or position in the time axis) of the maximum of the signal or the leading edge of the signal. Typically this is described as the arrival time ('time of arrival') or detection time.
  • 'detector data' refers to data that has originated from a detector, whether processed subsequently by associated or other electronics within or outside the detector.
  • the signal form may be determined by a calibration process that involves measuring the detector' s impulse response (such as time domain response or frequency domain response) to one or more single event detections to derive from that data the signal form or impulse response.
  • a functional form of this signal form may then be obtained by interpolating the data with (or fitting to the data) a suitable function such as a polynomial, exponential or spline.
  • a filter (such as an inverse filter) may then be constructed from this detector signal form.
  • An initial estimate of signal parameters may be made by convolution of the output data from the detector with the filter.
  • Signal parameters of particular interest include the number of signals and the temporal position (or time of arrival) of each of the signals.
  • the accuracy of the parameter estimation can be determined or 'validated' by comparing a model of the detector data stream (constructed from the signal parameters and knowledge of the detector impulse response) and the actual detector output. Should this validation process determine that some parameters are insufficiently accurate, these parameters are discarded. In spectroscopic analysis using this method, the energy parameters deemed sufficiently accurate may be represented as a histogram.
  • the data may include signals of different forms. In this case, the method may include determining where possible the signal form of each of the signals.
  • the method includes progressively subtracting from the data those signals that acceptably conform to successive signal forms of a plurality of signal forms, and rejecting those signals that do not acceptably conform to any of the plurality of signal forms.
  • the method may comprise detecting a pulse or pulses in said detector output data by:
  • the one or more functions are functions of time.
  • the one or more functions may not be functions exclusively of time.
  • the method may comprise providing the detector output data in, or converting the detector output data into, digital form before fitting the one or more functions to the detector output data.
  • the one or more functions are of the form:
  • v ⁇ t may be calculated numerically, such as by the formula
  • the one or more functions are of the form:
  • the method includes determining a location and amplitude of the pulse with a method comprising:
  • equation (26) becomes This form is therefore suitable for use in a numerically stable method for a calculating ⁇ .
  • Solving equation (26) can be done numerically, such as with a bisection method, especially since the left hand side is monotonic in ⁇ . Determining the left hand side for different values of ⁇ may be done by any suitable technique, such as with a Taylor series expansion for small ⁇ . (In practice, the value of ⁇ will generally be small because noise will generally preclude accurate characterisation of a pulse that started in the distant past.)
  • the method may comprise constraining ⁇ by requiring that Thus, because the left-hand side of the equation is monotonic in ⁇ , the constraint that is
  • This constraint can be implemented in with the constraints that a and ⁇ are not negative and ⁇ > ⁇ .
  • the method may also comprise constraining the amplitude of the pul se. This can be used, for example, to prevent a fitted pulse from being too small or too large. Indeed, referring to equation (30) above, if ⁇ is constrained to lie between -1 and 0 then A lies between Constraining a therefore constrains the amplitude A.
  • the function f is in the form of a function with three exponentials.
  • the time constants are known and dissimilar (so fewer problems of numerical imprecision
  • time constants are known and
  • the function f is of the form:
  • a and ⁇ are scalar coefficients, and the method comprises determining a and b.
  • determining the location comprises determining a location f * (a, b) where:
  • the function f may be a superposition of a plurality of functions.
  • the present invention may comprise a method and apparatus for estimating the location and amplitude of a sum of pulses from noisy observations of detector output data. It presented the maximum-likelihood estimate as the benchmark (which is equivalent to the minimum mean-squared error estimate since the noise is additive white Gaussian noise).
  • the method may comprise low-pass filtering the data before fitting the one or more functions.
  • the method comprises adapting the one or more functions to allow for a low frequency artefact in the detector output data. This may be done, in one example, by expressing the one or more functions as a linear combination of three exponential functions (such as
  • the method comprises forcing any estimates having the pulse starting within the window to start at a boundary of the window. [0231] In a particular embodiment, the method comprises maximizing window size or varying window size.
  • the method comprises transforming the detector output data with a transform before fitting the one or more functions to the detector output data as transformed.
  • the method may also comprise subsequently applying an inverse transform to the one or more functions, though in some cases it may be possible to obtain the desired information in the transform space.
  • the transform may be a Laplace transform, a Fourier transform or other transform.
  • estimating the location of the peak comprises minimizing an offset between the start of a window and a start of the pulse.
  • the method further comprises detecting a pulse or pulses in the data by:
  • a method for locating a pulse in detector output data comprising:
  • each of the one or more functions is a superposition of a plurality of functions.
  • the method comprises transforming the detector output data to produce stepped data or integral data, detecting at least one event, and estimating a pulse energy associated with the event.
  • detecting the at least one event occurs by fitting an expected pulse shape with a sliding window segment of the transformed pulse shape data.
  • the method further comprises the step of detecting peaks in the signal, wherein a detection metric is applied to the transformed data.
  • the detection metric is compared to a simple threshold - if the metric is less than the threshold, then no pulses are deemed present - if it exceeds the threshold, then one or more pulses may be present.
  • Declaration of significant peaks in the detection metric is conducted, when the slope of the peak changing from positive to negative indicates an event.
  • 'corrupt data' uncharacterised data
  • 'signal' is interchangeable in this context with 'pulse', as it refers to the output corresponding to individual detection events rather than the overall output signal comprising the sum of individual signals.
  • the temporal position (or timing) of a signal can be measured or expressed in various ways, such as according to the time (or position in the time axis) of the maximum of the signal or the leading edge of the signal. Typically this is described as the arrival time ('time of arrival') or detection time.
  • the method optionally comprises deleting samples within a set window around the rising edge to ensure the edge region of each pulse, where the real transformed pulse data differs from an ideally transformed pulse, is excluded from the calculations.
  • the method optionally comprises an assessment of variance of the energy estimations in the data, and validation of the modeled data.
  • the method may include building a model of the data from the processed data output and determining the accuracy of the modeling based on a comparison between the detector output data and the model.
  • the method includes creating a model of the detector output using the signal parameters in combination with the detector impulse response.
  • the method may include performing error detection by comparing the actual detector output with the model of the detector output, such as by using least squares.
  • the method may include discarding energy estimates deemed not sufficiently accurate.
  • the method includes presenting all sufficiently accurate energy estimates in a histogram.
  • the correct solution to take is the "positive” solution. It also relies on taking the correct solution to the complex square root.
  • Figure 14 illustrates the result achieved.
  • the rectangular window if used on its own, results in the measured spectrum being convolved with a sine function, with width at the mid amplitude of approx 10 samples, but significant oscillations - around 22% at the first negative going peak.
  • user window w it is possible to achieve a "time" domain response that is approximately raised cosine in nature, with very little oscillatory nature - around 0.2%.
  • the width at mid amplitude increases to around 20 samples.
  • the pulse parameters may be estimated from a time series capture of the digitised detector signal as follows:
  • Histogram the data block - the histogram bins are integers in the range of the sampled data for 14 bit signed data.
  • nthr noiseMean + 2 x noiseSigma.
  • Factors other than 2 may also be used, depending on the application.
  • Optional step If the count rate estimate is less than half or more than double the nominal count rate, redo the noise and signal threshold calculation with a data buffer length computed from the count rate estimate to get closer to numpO detected pulses.
  • the pulse detection state machined is as follows:
  • nump2 600 has been used, but this is somewhat arbitrary and based on how many pulses were actually in the test data Only half of these pulses will eventually be used in the calibration, so nump2 needs to be double the number of pulses that are required (desired) in the calibration process.
  • optimised results set the final values of alpha, beta, Ta. This can be either median or mean values of the optimised parameters.
  • the pulse form p(t) can be computed directly from the formula and the estimated Ta, a and ⁇ .
  • DC offset or signal baseline, used interchangeably
  • This DC offset can arise from various sources including the bias levels of analog electronics, the analog to digital conversion, and the detector itself.
  • Control theory suggests the DC offset error may be tracked and reduced to zero by generating a feedback signal that is proportional to the integral of the signal - however there is a significant problem in the case of pulse processing. Pulses introduce additional features to the signal that have non-zero mean. This introduces a bias dependent on pulse energy, count rate and pulse shape, which corrupts the feedback signal and prevents standard control loop tracking from successfully removing the DC offset.
  • the detector signal output is digitally processed to remove the pulse shaping effects introduced by analog electronics.
  • this processed signal results in a signal shape that has constant value in the regions between pulse arrivals, and a rapid change in value where pulses arrive. If a residual DC offset is present in the detector signal the processed signal changes linearly with time in the regions between pulse arrivals.
  • An error feedback signal that is proportional to the slope of this signal may be formed by taking the difference between two samples. These samples need not be consecutive, but may be separated by ' ⁇ ' samples in time. By choosing an appropriate value for TNT, a signal with a suitable signal to noise ratio may be found for driving a feedback loop.
  • the baseline tracking loop is not updated when a pulse has arrived between the two samples used to generate the feedback error signal.
  • the impact of bias may be further reduced by preventing the baseline tracking loop from updating when a pulse has arrived within a guard region on either side of the samples used to generate the feedback error signal.
  • Transient from a previous detected pulse has not decayed. Transient can be considered to last M samples after a pulse is detected.
  • the processed signal x(n) is about to reach positive overflow and wrap around to a large negative value. Do not process if x(n) is within a threshold ⁇ of positive or negative overflow.
  • k « 1 is the update gain, and is chosen to achieve the desired balance between fast response and noise on the DC estimate.
  • the same hardware can be used for tracking multiple baseline offsets in multiple channels in a time division multiplexed scheme.
  • the values for the tracking loop variables for each channel are stored/loaded when switching between channels.
  • the baseline tracking loop is prevented from updating until transient effects of the detector channel change has passed.
  • the purpose of the reference calculation is to establish the mean intensity for each detector. This value is used to scale all intensity histograms to unit energy. This is commonly referred to as normalisation.
  • a reference intensity is calculated for each detector. The reference intensity is calculated as the mean intensity over the first N slices in a scan. The intensity is the 1st bin in the FFT or the sum of all complex -valued elements in the FFT vector.
  • the reference is measured during an interval where:
  • X-rays are stabilised, so X-ray flux is not varying and will not vary for the
  • the current implementation uses a duration measured in slices. This can cause problems when the slice rate is slowed below 5 ms for example - the reference collection can run into the sample under test. This needs to be corrected in 2 ways to be fully robust:
  • systems according to the invention can comprise a means for sorting, diverting, removing or selecting at least a part of the material according to the characterisation of at least part of the material, such as may be applicable in a recycling plant.

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Biochemistry (AREA)
  • Analytical Chemistry (AREA)
  • Immunology (AREA)
  • Pathology (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Toxicology (AREA)
  • High Energy & Nuclear Physics (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geophysics (AREA)
  • Analysing Materials By The Use Of Radiation (AREA)
  • Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
  • Radiology & Medical Imaging (AREA)

Abstract

The invention provides a system and method for characterising at least part of a material comprising: a source of incident X-rays (4, 28) configured to irradiate at least part of the material; one or more detectors (300,302,312,1701,1704,1600,1607,1608,1604) adapted to detect radiation emanating from within or passing through the material as a result of the irradiation by the incident radiation (1700) and thereby produce a detection signal (313); and one or more digital processors (304-311,2000-2009) configured to process the detection signal (313) to characterise at least part of the material; wherein the one or more detectors (300,302,312,1701, 1704,1600,1607,1608,1604) and one or more digital processors (304-311,2000-2009) are configured to characterise at least part of the material by performing energy resolved photon counting X-ray transmission spectroscopy analysis.

Description

MATERIAL CHARACTERISATION SYSTEM AND METHOD
[0001] This invention relates to a system and method for material characterisation. In its broadest aspect, the invention may be applied to a wide range of materials, including but not limited to discrete items of material, a particulate substance, a non-particulate substance, mineral ore, coal, cement, mineral sand, fluid, hydrocarbon fluid, soil, agricultural material, food or beverage, metals, plastics, glass, timber, pharmaceuticals, contaminated material, hazardous material, waste, refuse, sewage, effluent, recyclable material, construction debris, natural materials, manufactured products, clothing, electronics, medical devices, automotive components, aerospace components, dentistry devices, and dentistry components. The invention has particular application to characterising coal or mineral ore travelling on a conveyor. The invention also has particular application to determining elemental, mineral ogical or density characteristics of a material, for example the percentage of iron in an iron ore, but in some embodiments may also be applied to determination of other physical characteristics, such as particle or grain size, and surface properties such as texture.
BACKGROUND
[0002] Various types of radiation-based techniques have been applied to characterisation of materials, such as X-ray, neutron activated gamma-ray, laser induced fluorescence, and other electromagnetic techniques. However, existing systems can be slow, intermittent and inaccurate; have limited capability to detect light elements or low concentrations; may produce poor representivity of results; and, in some cases require onerous radiation safety measures. Each existing technique has strengths and weaknesses over a wide range of criteria, including detection speed, maximum tonnage throughput, characterisation accuracy, material penetration depth, spatial resolution, elemental sensitivity, safety and operating cost.
[0003] There is therefore a need for an improved material characterisation system which can overcome at least some of the limitations of conventional material characterisation systems. [0004] The present invention relates to a system and method which has advantages in a combination of high throughput, accuracy, resolution, elemental sensitivity, safety and operating cost.
SUMMARY OF THE INVENTION
[0005] In accordance with a first broad aspect of the invention there is provided a system for characterising at least part of a material comprising:
a source of incident X-rays configured to irradiate at least part of the material;
one or more detectors adapted to detect radiation emanating from within or passing through the material as a result of the irradiation by the incident radiation and thereby produce a detection signal; and
one or more digital processors configured to process the detection signal to characterise at least part of the material;
wherein the one or more detectors and one or more digital processors are configured to characterise at least part of the material by performing energy resolved photon counting X-ray transmission spectroscopy analysis.
[0006] In accordance with a second broad aspect of the invention there is provided a method of characterising at least part of a material, the method comprising the steps of: irradiating at least part of the material with incident X-rays from an X-ray source; detecting radiation emanating from within or passing through the material as a result of the irradiation by the incident X-rays;
producing a detection signal;
processing the detection signal to characterise at least part of the material; and characterising at least part of the material by performing energy resolved photon counting X-ray transmission spectroscopy analysis.
[0007] Features of embodiments of the broad aspects are set out in the accompanying claims and the following description.
[0008] Throughout this specification including the claims, unless the context requires otherwise, the words "energy resolved photon counting X-ray transmission spectroscopy" are synonymous with "multi energy X-ray transmission spectroscopy" and "multiple X-ray absorptiometry". Further, the words "photon counting" mean processing and analysing individual packets of radiation, in particular, X-ray photons, in the course of processing the detection signal, and does not necessarily imply obtaining an actual enumeration.
[0009] Throughout this specification including the claims, unless the context requires otherwise, the term "packets" in relation to incident radiation includes individual massless quantum particles of X-rays (X-ray photons).
[0010] Throughout this specification including the claims, unless the context requires otherwise, the words "energy spectrum" in relation to a particular detector refers to a generation of energy values of the individual packets of radiation emanating from or passing through the part of the items under investigation as detected over a time interval from the particular detector, which energy values can comprise values over a range, typically continuous, and may be represented as a histogram of detection counts versus a plurality of defined energy bins, the number of bins representing the desired or achievable energy resolution and constituting at least 10 bins but preferably more than 50, 100 or 1000 energy bins.
[0011] Throughout this specification including the claims, unless the context requires otherwise, the word "material" refers to any material of solid, liquid or gaseous composition, including but not limited to discrete items of material, a particulate substance, a non-particulate substance, mineral ore, coal, cement, mineral sand, fluid, hydrocarbon fluid, soil, agricultural material, food or beverage, metals, plastics, glass, timber, pharmaceuticals, contaminated material, hazardous material, waste, refuse, sewage, effluent, recyclable material, construction debris, natural materials, manufactured products, clothing, electronics, medical devices, automotive components, aerospace components, dentistry devices, and dentistry components.
[0012] Throughout this specification including the claims, unless the context requires otherwise, the word "characterisation" means measurement, classification, assessment, or analysis of one or more features or properties of said material, including but not limited to physical, chemical, optical, magnetic or other material properties. BRIEF DESCRIPTION OF THE DRAWINGS
[0013] Figure 1 shows a material characterisation system according to one embodiment of the invention installed on a conveyor adapted to ore characterisation, adapted to high- capacity conveyors and having an X-ray source above the conveyor with detectors below .
[0014] Figure 2 shows a material characterisation system according to another embodiment of the invention installed on a conveyor adapted to ore characterisation, adapted to high- capacity conveyors and having an X-ray source below the conveyor with detectors above. S.
[0015] Figure 3 shows the system of Figure 1 or Figure 2 installed on a sample conveyor.
[0016] Figure 4 shows a material characterisation system similar to Figure 1 or 2 installed to characterise material falling off the end of a conveyor.
[0017] Figure 5 shows a detector module comprising a plurality of detector arrays according to an embodiment adapted for ease of maintenance and replacement.
[0018] Figure 6 and 7 show a detector module variant in which the detector arrays are angled for optimum detection.
[0019] Figure 8 shows a source module with shielding and an IP rated enclosure according to an embodiment.
[0020] Figure 9 shows a more detailed view of the detection system and processing electronics according to an embodiment.
[0021] Figure 10 illustrates a flowchart for a method of effective atomic number (effective Z) processing of the full energy spectra computed by the pulse processing electronics according to an embodiment.
[0022] Figure 11 is a graph illustrating the removal of pileup of two pulses from a spectrum according to an embodiment.
[0023] Figure 12 is a graph illustrating the removal of pileup of two and three pulses from a spectrum according to another embodiment. [0024] Figure 13 is a graph illustrating the partial removal of pileup of two and three pulses from a spectrum when assuming only 2 pulse pileup according to an embodiment.
[0025] Figure 14 is a graph illustrating the shape of the spectral smoothing filter when using a rectangular window or a raised cosine pulse window according to an embodiment.
[0026] Figure 15 illustrates how data is arranged and built up into an image of a scanned sample, prior to further post processing and image display according to an embodiment.
[0027] Figure 16 is a graph illustrating the un-calibrated received spectra from a plurality of detectors according to an embodiment.
[0028] Figure 17 is a graph illustrating a set of detector gains calculated based on a calibration procedure according to an embodiment.
[0029] Figure 18 is a graph illustrating the received spectra from a plurality of detectors after setting the digital gain of the detectors based on the detector gains illustrated in Figure 17.
[0030] Figure 19 illustrates results of the effective Z interpolation process for a 10% transmission case according to an embodiment.
[0031] Figure 20 illustrates effective Z plotted against intensity (percent transmission) for a range of material samples tested according to an embodiment.
[0032] Figure 21, 22, 23 illustrate a detector module comprising an array of detectors according to embodiments in which a linear (or curved) array of scintillation crystals is coupled to an array of pulse producing elements in the form of silicon photomultipliers.
[0033] Figure 24 and 25 illustrate a detector module comprising single scintillation crystals individually coupled to an array of pulse producing elements in the form of silicon photomultipliers by means of an optical coupling layer interposed between the scintillation crystals and silicon photomultipliers according to an embodiment.
[0034] Figure 26 illustrates a detector subsystem converting photons into voltage pulses for pulse processing according to an embodiment. [0035] Figure 27 illustrates an example of the formation of clusters, where single tiles are ignored according to an embodiment.
[0036] Figure 28 illustrates an example of effective Z processing steps according to an embodiment.
[0037] Figure 29 illustrates a table relating to an edge mask L(c) indexed on columns according to an embodiment.
[0038] Figure 30 illustrates the behaviour of the moving average as it transitions over an edge according to an embodiment.
[0039] Figures 31 to 37 show aspects of a number of embodiments of the invention adapted to lower capacity conveyors.
DETAILED DESCRIPTION OF EMBODIMENTS
[0040] It is convenient to describe the invention herein in relation to particularly preferred embodiments. However, the invention is applicable to a wide range of methods and systems and it is to be appreciated that other constructions and arrangements are also considered as falling within the scope of the invention. Various modifications, alterations, variations and or additions to the construction and arrangements described herein are also considered as falling within the ambit and scope of the present invention.
[0041] In addition to X-rays being attenuated when transmitted through matter, X-rays passing through matter interact with that matter via a number of modalities including: scattering off crystal planes, causing fluorescence X-ray emission from within the electron structure of the elements; and, scattering off nano-scale structures within the material being scanned. These forms of interaction slightly modify the energy spectrum of the transmitted X-ray beam and by detecting and analyzing this change in energy spectrum it is possible to deduce elemental specific information about the item through which the X-ray beam passed.
[0042] The system of one of the embodiments described below provides for a detection system capable of estimating the energy of the individual X-ray photons received at the detector. This is achieved using a single detector array per X-ray source, with each of the detectors in the array constructed from an appropriate detector material coupled to a photomultiplier, producing an analog signal comprising a series of pulses - one pulse for each detected X-ray, which may or may not be overlapping when received at the detector.
[0043] A pulse processing system is then used to generate a histogram for each single detector. This histogram comprises of a count of the number of X-rays falling into each histogram bin in a given time interval. The histogram bins represent a range of energies of the received X-rays, and the histogram is therefore the energy spectrum of the received X- ray beam. There may be a large number of histogram bins - for example 1024 or more separate energy bands.
[0044] The system of the described embodiment uses this full high resolution energy spectrum to obtain a highly accurate estimate of the material's effective atomic number (effective Z or Zeff), resulting in a vastly superior classification of the material. Effective Z may be used, for example, to distinguish different contents in a material, enabling assessment of mineral ore or coal quality, or waste material characteristics.
[0045] The system may also extract other information from the energy spectrum, for example material density, in order to aid in characterisation of the material, which may be implemented in software in the digital processors.
High Level Overview
[0046] Figure 1 shows a material characterisation system according to an embodiment of the invention, installed on a conveyor 1 adapted to convey material in the direction of the arrow A. The system comprises a housing 2 encompassing conveyor 1. Coal is seen progressing along conveyor 1 through aperture 3 in housing 2. Housing 2 houses a source of incident X-rays in the form of an X-ray generator 4 positioned above conveyor 1 which emits a fan-shaped beam of X-rays 5 perpendicular to the direction of travel providing an irradiation region. A plurality of modules 6 is positioned within housing 2 beneath conveyor 1 so as to intersect and detect radiation emanating from or transmitting through the coal radiated by beam 5. Each module 6 comprises a plurality of detectors in an array, wherein each detector comprises a detector subsystem, wherein digital processor electronics may be provided within the modules 6 or separately and in communication with modules 6 as described below. A suitable X-ray generator 4 in this embodiment for effective Z determination is obtainable from Industrial Control Machines SA (www.icmxray.com), for example the Sitex™ range. Suitable X-ray generators may be selected for their power and physical design according to the application to specific material and integration requirements.
[0047] Modules 6 are provided in a plug-in modular constructions, allowing easy removal and replacement on-site, as is X-ray generator 4.
[0048] A control panel and user interface 7 resides on a side of housing 2, and the system is equipped with wireless communication antenna 8 for remote communication control. The user interface 7 may comprise a display screen or indicator lights to provide visual outputs, system status of emitted results. Wireless communication antenna 8 may connect with a tablet or other portable computer an operator standing nearby. In addition, wireless communication antenna 8 may be adapted to communicate with an on-site control centre or over a remote communication system such as the Internet to an off-site remote control centre. Alternatively or in addition, a wired network communication system may be employed.
[0049] The system may be installed on a main conveyor 1 with the aim of characterising substantially all of the coal produced at the mine, handled by a processing facility, or shipped to or from a port or customer.
[0050] Referring now to Figure 2, an alternative configuration is shown with X-ray generator 4 disposed beneath the conveyor and detectors in modules 6 disposed above the conveyor. This configuration may be preferable in some circumstances to avoid interference caused by X-rays intercepting support structures of conveyor 1, by providing a well-positioned source 4 between and avoiding the under conveyor support structures.
[0051] Referring now to Figure 3, in an alternative configuration a sample analysis facility may be provided off-line comprising housing 2 installed on a sample conveyor 9 which feeds coal samples into the system for analysis placed manually on conveyor 9 through hopper 10.
[0052] Figure 4 shows another embodiment with the system installed to characterise coal falling off the end of a conveyor 11 in the direction of arrow B. Housing 2 is similar or identical to the embodiment of Figure 1 but is installed with the aperture disposed horizontally so as to catch the coal falling of the end of conveyor 11. A second conveyor 12 is disposed beneath aperture 3 to catch and continue transport of the coal. This embodiment may be designed to provide analysis in applications where the material of the conveyor interferes somewhat with analysis of the detected radiation and where it is therefore advantageous that only the coal is interposed between the X-ray generator and the detectors in the irradiation region. In other embodiments, the material may be falling off the end of a conveyor onto a storage pile or into a cargo hold of a ship or other long-range means of transport.
[0053] In an embodiment of the invention, the system is modularised to facilitate manufacture, maintenance and replacement. Modularity reduces service component weight, minimises on-site work, allows maintenance tasks to be undertaken by unskilled workers, provides logical cost breakdowns for replacement parts and improves reliability by sealing sensitive components together in each simple module.
[0054] In Figure 5, a detector module 20 is shown in a cutaway illustration as comprising several detector boards 21, each detector board 21 comprising a linear array of detector elements as described below disposed behind collimators 26 are constructed from lead which allow only a narrow slit of x-radiation to impinge on the detector elements. The detector boards in 21 are mounted on an aluminium base of the precision manufactured enclosure, allowing reliable alignment of the detector boards and heat transference. The module also comprises aluminium side support rails 24 which may be provided as an extension of the aluminium base or otherwise rigidly and accurately fixed to the aluminium base. A top cover 23 of the enclosure is composed of 5 mm thick aluminium, which is effectively opaque to X-rays, except for a slit composed of thinner aluminium (<lmm) or HDPE which is X-ray translucent disposed above collimators 26, to allow the emanating radiation to enter the detector module and pass through the slits of collimators 26. A cooling circuit passes through custom heatsinks 25 on each detector board 21, and liquid coolant may be circulated through coolant connectors 22 . [0055] Referring now to Figure 6 also illustrated in cutaway and Figure 7 in cross-section, an alternative embodiment of detector module comprises a taller module which allows each detector board 21 to be angled so as to face the fan -shaped X-ray beam at right angles. Detector arrays 27 can also be seen in the crosssectional representation of Figure 7. In this embodiment, the aluminium base 24 and side rail are formed from the same sheet. This embodiment is particularly suitable for the variants described in figures 31-38 where the X-ray source is close to the detector.
[0056] Referring now to Figure 8, a source module is shown with similar module enclosure 23 and aluminium base 24. The X-ray source is contained within a source enclosure 28 composed of 5 mm thick lead walls to provide X-ray shielding, cast over an internal aluminium pipe (not shown) to provide rigidity. A thicker 10 mm lead region 36 is provided surrounding the beam exit point of the contained X-ray source and incorporates a collimating slit 35 (which may be an HDPE or thin aluminium window) to provide a substantially planar beam of source X-rays exiting the enclosure for irradiation of the material in the irradiation region. Cooling holes 33 provide internal air circulation and cooling. A corresponding slit 37 on a top face of module enclosure 23 allows exit of the X- ray beam. Source enclosure 28 comprises precision machined end caps 32 and guide tabs 34 for accurate positioning within module enclosure 23. Once source enclosure 28 is slid into position, door 30 is closed. A partitioned electronics and power section 31 is provided also within module enclosure 23, with internal connections 29 to the X-ray source. Cooling may also be provided through coolant passing through multipurpose connectors such as 22 above for the detector module (not shown).
[0057] Before proceeding to discuss in detail the detection system and processing, reference is now made to Figures 31-37 which describe an embodiment of the system suitably adapted to smaller capacity conveyors.
[0058] In the embodiment of Figure 31, the source and detector are positioned around conveyor belt 1 with source module 38 above conveyor belt 1 and detector module 39 below conveyor belt 1. Source module 38 and detector module 39 are accurately positioned on a precision engineered first chassis being "tech" chassis 4. Tech chassis 41 is seated on a second chassis being an external or support chassis 42 through vibration damped mounts such as rubber mounts or wire rope isolators (not shown). Support chassis 42 in turn is permanently bolted to conveyor base 40. Text chassis 41 and support chassis 42 are distinguished in the diagram by forward and backward diagonal lines. The use of a separate "floating" tech chassis 41 enables tight dimension tolerances which are required between source module 38 and detector module 39. Source module 38 and detector module 39 may be removed from tech chassis 41 when required.
[0059] Figure 32 shows a variant of the embodiment of Figure 31, equipped with a C- shaped support chassis 42 and a similar floating tech chassis 41. The entire tech chassis 41 may be slid out in the direction shown.
[0060] Figure 33 shows another variant, with a similar C- shaped support chassis and a rotatably attached tech chassis 41. In this variant, source module 38 is shown underneath a conveyor 1 and detector module 39 is shown above. The tech chassis may be rotated in the direction of the arrow shown to provide good access for removal of entire tech chassis or individual modules by crane as required.
[0061] Figure 34 shows a cross-section of the embodiment of Figure 31, to illustrate passage of the X-ray beam 44 from source module 38 into detector module 39, and to illustrate the optional inclusion of beam shielding elements in the form of lead plates 43 disposed on a shielding chassis structure 45 which also sits on support chassis 42. Figure 35 shows in exploded view the components of shielding chassis structure 45, lead baffle plates 46, lead top plate 48 (passage slit not shown) and side plates 47 for disposing beneath the conveyor 1. The beam shielding may be necessary because of potential exposure of personnel, depending on the geometry of the system and any containing structures.
[0062] Referring now to Figure 36 and 37, the source and detector modules may be disposed within a weatherproof enclosure 56 with sliding door 50, shown here on a conveyor 1 equipped with a hood 52. For elevator conveyors, and access platform 57 may be provided with user interface 58 located on platform 57. [0063] The electronic and X-ray detection components will now be described. One or more digital processors process the received X-ray pulses from the detector arrays within the detector modules. Depending on the implementation architecture, the digital processors may: i. Reside on the same boards as a detector array within the modules. ii. Reside on separate hardware, housed within the module chassis or outside the chassis in auxiliary equipment (not shown)
iii. Form part of a Host controller responsible for overall operation of the system housed inside or outside the chassis.
iv. A combination of the above.
[0064] While the above-described embodiments utilise a variety of convenient housing geometry, in the broadest aspect of the invention many other geometries are possible.
[0065] Further, while the embodiments described below contemplate a large number of detectors each producing signals for analysis, the system in principle is operable with a single detector and in its broadest aspect there are any number of detectors.
[0066] Figure 9 shows a more detailed view of the detection system and processing. This figure shows the steps for a single detector. The effective Z may utilise, and image post processing will require, access to the spectra from all detectors.
[0067] In this embodiment, an image of the coal is formed from rows of pixels as the coal passes through the detector region. In its broadest aspect, the invention does not require the formation of an image and a system comprising a single detector is within the scope of the broadest aspect. The detection system and processing electronics comprises:
1. A detector subsystem (301) for each individual detector (with N such subsystems for a lxN detector array), the detector subsystem comprising:
a. Detector material for detecting the incident X-rays (300) and converting each detected X-ray to a light pulse
b. A photomultiplier for receiving and amplifying the incident light pulses into an analog signal comprising pulses (312) that may or may not overlap c. Appropriate analog electronics, which may include filtering d. An optional variable gain amplifier (302). Fixed analog gain may also be used, or it may not be desirable to use additional gain to the photomultiplier An analog to digital converter (303) for each individual detector, to convert the analog signals into digital values (313).
A variable digital gain (304) for each individual detector to appropriately adjust the digital signal levels prior to processing.
High rate pulse processing (305) for each detector subsystem (301), for example the pulse processing systems disclosed in US Patent No. 7383142, US Patent No. 8812268 and WO/2015/085372, wherein the pulse processing comprises:
a. Baseline tracking and removal, or fixed baseline removal.
b. Detection of incoming pulses.
c. Computation of the energy of each detected pulse.
d. Accumulation of the computed energy values into an energy histogram (energy histogram) (315).
e. Output of accumulated histogram values each time a gate signal is received. f. Reset of the histogram values for the next collection interval.
A gate signal source (306) which outputs a gate signal (314) at a regular pre- configured interval. The gate interval is a constant short interval that determines the histogram accumulation period. This gate interval also determines the pixel pitch in the resulting X-ray image. In the absence of a gate signal source, and gate signal, another appropriate means may be used to control and synchronise the timing of energy histogram collection across all detectors. For example, a suitably precise network timing signal may be used instead of the gate signal. Calibration System (307), which receives input from appropriate analog and digital signals and then communicates the desirable calibration parameters back to the various processing blocks. The calibration system performs:
a. Pulse parameter identification
b. Gain calibration
c. Energy calibration
d. Baseline offset calibration (where fixed baseline is used) e. Count rate dependent baseline shift
8. Effective Z computation (308), which takes the computed energy spectra in each detector during each gate interval and determines the effective Z of the sample. This in turn leads to the production of an effective Z image.
9. Intensity image generation including:
a. Intensity image (309), based on total received energy across the energy spectrum.
b. High penetration or high contrast image (310) determined by integration of selected energy bands from the full energy spectrum.
10. Image post processing and display (311), with features that may include one or more of the following:
a. Image sharpening
b. Edge detection and/or sharpening
c. Image filtering
d. Application of effective Z colour map to colour the image pixels based on identified material.
e. Selection, display and overlay of 2D images for each detector array
i. Effective Z
ii. Intensity
iii. High Penetration / High Contrast images
f. Display of images on an appropriate monitor or other display device.
[0068] As described above, and illustrated in Figure 15, the images produced for display comprise a number of data elements recorded for each of N detector elements (501) and for each gate interval (500).
[0069] The data obtained for detector i during gate interval j is used in the production of effective Z, intensity and high penetration/high contrast images as shown in Figure 9. During the processing, a number of elements are recorded in each pixel (502), including one or more of: 1. The X-ray energy spectra.
2. The computed effective Z value
3. The intensity value (full spectrum summation)
4. High Penetration / High Contrast intensity values computed from integration of one or more energy bands.
[0070] Figure 15 illustrates how this data is arranged and built up into an image of the scanned sample, prior to further post processing and image display.
Detector Subsystem
[0071] The detector subsystem utilises a scintillator (such as phosphor) coupled to an array of PIN diodes to convert the transmitted X-ray into light, and subsequently into an electrical signal.
[0072] When an X-ray impacts the detector it produces an electron charge in the detector proportional to energy of the X-ray, wherein the higher the energy is the more charge is induced in the detector. However, more detailed examination of prior art detector arrays has illustrated that detector systems do not have the resolution to detect individual X-ray photons, and instead they integrate all the charge produced by the detector pixel over a given time period and convert this into a digital value. Where the instantaneous flux of X- rays on the detector pixel is large, a large digital value is produced (a bright pixel in the image) and where few X-rays impact the detector a small digital value is produced (a dark pixel in the image).
[0073] The detector subsystem of this embodiment comprises: a) A detector material
b) A photomultiplier material coupled to the detector material using an appropriate means
c) Analog electronics
[0074] The detector material may be of dimensions X by Y by Z, or some other shape. The photomultiplier may be a silicon photomultiplier (SiPM) and the coupling means may be a form of optical grease or optical coupling material. It may be desirable to use a form of bracket or shroud to hold the detector in position relative to the photomultiplier. The photomultiplier requires appropriate power supply and bias voltage to generate the required level of amplification of the detected signal. In this embodiment, typically each detector dimensions approximately l-3mm in width and depth, with a height of approximately 2- 6mm, and inter-detector pitch of approximately 2-200mm.
[0075] In an X-ray scanning application, a large number of single element detector subsystems may be required to produce each detector array. It may be desirable to group these in an appropriate way, depending on the specific X-ray scanner requirements. Individual elements of detector material may be grouped into a short array of M detectors. Small groups of M detector elements may be mounted onto a single detector board, for example 2, 4 or more groups of M onto one board. The full detector array is then made up of the number of detector boards required to achieve the total number N of detector elements per array.
[0076] In various embodiments, detector subsystems can be arranged in a number of different configurations including: curved arrays of 1 x N devices; linear arrays of 1 x N devices; square or rectangular arrays of N x M devices; or L-shaped, staggered, herringbone or interleaved arrays. In this embodiment, three curved arrays of 1 x 24 devices are utilised, disposed within the detector ring housing 31, however, any chosen distribution of detectors and detector boards may be used according to the desired geometry of the system and population of detectors around the detector ring housing 31
[0077] One example of a detection device, used to convert incoming radiation photons into and electrical signal, is the combination of a scintillation crystal, coupled to a silicon photomultiplier (SiPM) or multi-pixel photon counter (MPPC).
[0078] In such a detection device a scintillation crystal such as LSYO (1701) is used to convert the incoming radiation photon (1700) into UV photons (1703). In the case of LYSO scintillation material the peak emission of UV photons occurs at 420 nm, other scintillation material such as those listed in Table 1 may have different emission peaks. Subsequent to the interaction of the radiation photon (1700) with the scintillation crystal (1701) to produce UV photons (1703) a multi -pixel photon counter, or silicon photomultiplier (1704) with sensitivity in the UV region (such as one with the performance metrics in Table 2) may be used to detect these photons and produce an electrical signal.
[0079] Figure 21 depicts a linear array of LYSO scintillation crystals (1600), indicative of how single detection devices can be joined together to form a linear array. In this indicative example the individual LYSO crystals (1600) have a cross section of 1.8mm and a height of 5mm, the individual LYSO crystals (1600) are wrapped around the sides in a reflective material to assist in collecting all the UV photons. The pitch of this exemplary array is 2.95mm, the length is 79.2mm and the width of the array is 2.5mm.
[0080] Figures 22 and 23 depict a detector array from a top view and side view respectively, comprising the linear array of LYSO crystals depicted in Figure 21 coupled to an electrical pulse producing element (1604) on substrate (1605). The electrical pulse producing element may comprise a silicon photomultiplier (SiPM). Enhanced specular reflector (ESR) or aluminium or other reflective foil (1601) is disposed around side surfaces of the scintillation crystals to direct the scintillation photons onto the silicon photomultiplier material (1604) and prevent light leakage (cross-talk) between adjacent detection devices. Optionally, optical coupling (1606) may be interposed between the LYSO crystals and SiPM, and may comprise any number of known suitable materials, for example, a thin layer of optically transparent adhesive.
[0081] In the current invention, the detector arrays may be curved rather than linear as shown in Figures 21 and 23 to appropriately fit inside different housing shapes.
[0082] In another embodiment, scintillation crystals (1607) may be individually coupled to electrical pulse producing elements (1604), as depicted in figures 23 and 24. Coupling may be achieved by a number of methods, for example interposing an optically transparent adhesive film (1609) or optical coupling material between the scintillation crystals (1607) and electrical pulse producing elements (1604), where the electrical pulse producing elements (1604) may comprise SiPMs or an MPCC. Coupling may be performed by a 'pick and place' assembly machine to individually align and couple components and coupling material. Scintillation crystals may be wrapped in a reflective material such as a foil or ESR material (1608) to aid in the capture of photons. [0083] In any of the embodiments, the LSYO crystals (1600, 1607) may typically have a cross-section (width) approximately l -2mm, a depth of approximately l-2mm, and height of approximately 3 -5mm, where the reflective or ESR film (1601, 1608) is approximately 0.05mm - 0.1mm thick. In a preferred embodiment of the detectors shown in fig. 23 the cross-section is 1.62mm, the depth is 1.24mm, the height is approximately 4.07mm, and the ESR film is 0.07mm thick. The cross sectional area of the scintillator material is preferably greater than 1mm square, and may be greater than 2mm square and less than 5mm square.
[0084] While the exemplar detector subsystem design uses a scintillator which is compact, robust, cost effective and non-hygroscopic, in the broadest aspect of the invention other detector subsystems can be considered. These include detector subsystems which use alternate inorganic or inorganic scintillator materials, the characteristics of some such material are provided in Table 1. Other mechanisms for converting radiation photons into electrical signals could also be considered for the detector subsystem. Some examples of other detector materials options include:
a) High Purity Germanium (HPGe): Achieves 'gold standard' resolution of 120 eV for the Fe55 X-ray line at 5.9 keV, detectors can be made > 10 mm thick thus detect high energy X-rays up to many 100s of keV.
b) Silicon Drift Diode (SDD): SDD detectors measuring relatively low energy radiation. For the same Fe55 line at 5.9 keV SDD detectors have a resolution of approximately 130 eV. Also, these detectors can be operated at higher count rates than HPGe detectors and just below room temperature.
c) PIN Diodes: The detection efficiency for X-rays up to 60 keV is substantially higher than SDD detectors and falls off to approximately 1% for X-ray energies above 150 keV. These detectors can be operated at room temperature, however, resolution improves with cooling, resolution of the 5.9 keV line is ~ 180 eV,
d) Cadmium Zinc Telluride: Is a room temperature solid state radiation detector used for the direct detection of mid-energy X-ray and Gamma-ray radiation. It has a detection efficiency for 60 keV X-ray very close to 100% and even for X-rays photons with energies of 150eV the detection efficiency remains greater than 50%. e) Cesium Iodine (CsI(Tl)): This is a scintillation material used for detection of X-rays in medical imaging and diagnostic applications. The scintillation material is used to convert the X-ray into photons of light which are generally then converted into an electrical signal either by a photomultiplier tube. Csl is a cheap and dense material and has good detection efficiency of X-rays and Gamma-ray to many 100s of keV.
Figure imgf000021_0001
Table 1 - properties of a range of scintillator materials. Geometrical Data
Figure imgf000022_0001
Table 2 - performance data for LYSO scintillators.
Processing Steps
[0085] The following sections outline the steps involved in processing each particular stage of the various algorithms.
1. Calibration
[0086] The material characterisation system comprises a large number of individual detectors. While each detector and associated electronics is ideally designed to have identical response to incident radiation, in practice this will not be possible. These variations between detectors result in detector to detector variation in energy spectrum output. By properly and fully calibrating the detection system, the energy spectra output from the pulse processing digital processors can be appropriately calibrated so they represent received X-ray intensity in known narrow energy bins.
1.1.Detector Pulse Calibration
[0087] Detector pulse calibration is used to identify the pulse characteristics for each detector required by the pulse processing system. The exact parameters required may vary, depending on the detection system. For typical applications using the pulse processing method disclosed in US Patent No. 7383142 and US Patent No. 8812268, the pulse is modelled as an averaged dual exponential of the form:
Figure imgf000023_0001
where a and β are the falling edge and rising edge time constants respectively, t0 is the pulse time of arrival, Ta is the pulse averaging window, and A is a pulse scaling factor related to the pulse energy.
[0088] The processing requires the two parameters a and β, and the pulse form p(t) which can be obtained via an appropriate calibration method, or from knowledge of the design of the detection subsystem. A suitable method for estimating α, β and p(t) from received pulses is described below.
1.2. Detector Gain Calibration
[0089] Each detector subsystem, combined with an analog to digital converter, will have slightly different characteristics due to manufacturing variations. As a result of such component variations, the energy spectra will be scaled differently. Variations other than gain scaling are handled within the Baseline Offset Calibration or the Energy Calibration.
[0090] The objective of the gain calibration is to achieve alignment of the energy spectra output by the pulse processing electronics across all detectors. The need for absolute accuracy may be reduced or eliminated if per detector energy calibration is applied.
[0091] Gain calibration may be achieved in a number of ways. The following approach may be applied: Setup a known X-ray source.
a. A material with particular characteristics can be inserted into the beam. For example, lead (Pb) has a known absorption edge at 88 keV.
b. Make use of the known radiation of the detector material (e.g. LYSO), detected by itself (the self spectrum).
Measure the energy spectrum on each detector, as output by the pulse processing electronics.
Ensure sufficient data is obtained in order to achieve a smooth spectrum with minimal noise.
Select a feature or features on which to perform the alignment. For example,
a. A specific peak in the spectrum
b. An absorption edge (for the case of Pb)
c. The entire spectrum shape (appropriate for LYSO self spectrum)
Compute the histogram bin corresponding to the feature location for each detector. Compute the median of these feature location bins across all detectors.
The required gain for each detector is then computed as the ratio of median location to the specific detector feature location. Note: The median or other suitable reference (e.g. maximum or minimum) is chosen. Median is chosen so some channels are amplified, and some are attenuated, rather than attenuating all channels to the minimum amplitude.
The gains are then applied to each detector channel. The gain may be applied as an analog gain, digital gain, or combination of the two, depending on particular system functionality. For best results, at least part of the gain is digital gain, where arbitrarily fine gain variation can be achieved.
Re-measure the energy spectrum on each detector and confirm that the required alignment has been achieved.
If desirable, compute an updated / refined gain calibration for each detector, and apply the updated calibration to each detector.
Repeat steps 9 and 10 as often as desirable to achieve the required correspondence between spectra from all detectors. [0092] For the methods of effective Z computation outlined in this disclosure, it has been found that spectral alignment to within 1 -2% can be achieved and is desirable for accurate and consistent effective Z results.
[0093] In a practical implementation of the detection subsystem there may be a number of detector cards, each with a number of detector subsystems. The total number of detector subsystems may be several thousand or more. Results from one example of such a detector board are presented here. The example board comprises 108 detectors, with LYSO used as the scintillator material in this case. These detectors are packed into linear arrays of 27 detectors. Each detector board then uses 4 x 27 detector arrays to achieve a total of 108 detectors.
[0094] When X-rays are incident upon a detector, photons are emitted by the LYSO based on the energy of the incident X-ray. Each detector is placed above a SiPM, and it is the SiPM that detects and amplifies the emitted photons. The detectors are coupled to the SiPM via an optical grease. The gain of each SiPM is determined by the bias voltage applied, and the SiPM breakdown voltage. As a result of variations in the LYSO material, quality of coupling between the LYSO and the SiPM, and also variations in the SiPM gain and SiPM material properties, there can be considerable difference in the received pulse energy for a given incident X-ray energy.
[0095] The effect of the variation in detected pulse energy is that the energy spectra from all detectors are not the same. This can be seen in Figure 16, where the uncalibrated received spectra from all 108 detectors are plotted. These energy spectra are measured where a sample of lead (Pb) is in the X-ray beam, and the structure of the Pb spectrum is clearly seen. It can be seen that the tail of the energy spectrum spreads across a range of approximately 150 histogram bins. This means the actual energy per bin is quite different for each detector.
[0096] By following the gain calibration procedure outlined above, a set of detector gains was computed, as shown in Figure 17. From the figure, the calibrated gain value ranges from approximately 0.75 to 1.45. [0097] After setting the digital gain to be equal to the detector gains in Figure 17, the energy spectra from the 108 detectors were re-measured, as shown in Figure 18. It is clear the energy spectra are now well aligned, indicating the success of the gain calibration. The different spectrum amplitude levels reflect the range of factors discussed above that can affect the resulting energy spectrum. In this case, some detectors are capturing an overall greater number of X-rays than others, indicated by the higher spectrum amplitude. Nonetheless, the alignment of the spectral features is very good as required.
1.3. Baseline Offset Calibration
[0098] Each detector subsystem may have a slightly different baseline level, as measured at the output of the analog to digital converter. In order for the pulse processing electronics to accurately estimate the energy of received pulses, the baseline is estimated and removed. Any suitable method can be used including, for example:
1. Offline measurement of baseline offset (with X-rays off):
a. Record and average a series of samples from the detector
b. Use this average as the baseline offset to be subtracted from all data
2. Online baseline offset tracking and adaptation:
a. Use the pulse processing output to estimate and track baseline offset b. Filter the (noisy) tracked baseline values and update the baseline offset register accordingly
c. Use an initial period of convergence with X-rays off, followed by continuous adaptation while X-rays are on
1.4. Energy Calibration
[0099] The pulse processing electronics will produce an energy spectrum that is uncalibrated. That is, the output will comprises a number of counts in a set of histogram bins, but the exact energy of those histogram bins is unknown. In order to achieve accurate effective Z results, knowledge of the energy of each bin is required. [0100] This is achieved as follows:
1. Use a source with known spectrum peaks. One suitable example is a Bal33 source, with spectral peaks at 31, 80, 160, 302 and 360 keV
2. Measure the uncalibrated energy spectrum.
3. Determine the histogram bins corresponding to the known spectrum peaks
[0101] Instead of using a single source with multiple peaks, it is also possible to use a narrow band source with variable (but known) energy, and measure the histogram bin as a function of energy for a range of energies.
[0102] Once a relationship between histogram bins and energy has been measured, it is possible to either:
1. Create a lookup table for the energy of each histogram bin.
2. Estimate parameters of a suitable functional form. For a LYSO/SiPM combination it has been found a quadratic model fits the observed parameters very well. This gives a result of the form:
Figure imgf000027_0001
where A, B and C are constants determined from the measured Bal33 spectrum. This formula is inverted to define Energy as a function of Histogram Bin expressed in terms of the same A, B and C.
[0103] If the variation between detectors is sufficiently small (requiring good component matching and good gain calibration), then a single energy calibration can be applied to all detectors. In this case, averaging the calibration parameters across a number of detectors exposed to the Bal33 source will yield a superior estimate of the Energy Calibration parameters.
[0104] Alternatively, individual calibration table / calibration parameters can be generated for each detector. 1.5. Count rate dependent baseline shift
[0105] Depending on detector / photomultiplier combination, it may be desirable to compensate for a count rate dependent baseline shift. The consequence of this shift is a right shift of the energy spectrum as count rate increases. To properly apply the energy calibration, the spectrum is moved back to the left by a specified number of bins / energy. The calibration required is either: a) A lookup table, defining baseline shift for each count rate, with intermediate results obtained via interpolation.
b) A functional form, where baseline offset is expressed as a function of count rate.
[0106] Any suitable method can be used for this calibration, including injecting a known source spectrum of variable count rate, and recording the spectrum shift as count rate increases. Ideally the source has a narrow energy band so the shift can be clearly measured, and also variable energy so the offset can be calibrated as a function of energy if desirable.
[0107] The need for removal of count rate dependent baseline shift can be diminished or even eliminated if online baseline offset tracking and removal is used.
1.6. Residual Spectrum Calibration
[0108] The residual spectrum is measured with a large mass of material in the beam, sufficient to completely block the X-ray beam, such as a large thickness of steel. In practice, a small level of energy still reaches the detector array, whether from scatter or other mechanisms, and this residual spectrum must be measured so it can be removed from the received spectra during normal operation.
[0109] The residual spectrum is then measured by averaging the received spectra for a number of gate intervals with the blocking mass in the beam.
1.7. Pileup Parameters
[0110] The pileup parameters can be calibrated in several ways, for example:
a) Estimation of pileup parameters from the nature of the received spectra. b) Estimation of pileup parameters from knowledge of the signal, the received pulse count rate, the ADC sampling rate and the pulse detection method.
c) Measurement of the pileup parameters as follows:
i. Use a narrow energy source, where the energy and count rate can be varied. ii. Measure the received spectrum as the source energy and count rate are varied. iii. Directly measure the ratio of received 2-pulse and 3-pulse pileup to the main signal peak.
iv. Form a lookup table of 2-pulse and 3-pulse pileup as a function of count rate and energy.
2. High Rate Pulse Processing
[0111] A high rate pulse processing system (305), such as those disclosed in US Patent No. 7383142, US Patent No. 8812268 or WO/2015/085372, is allocated to each detector subsystem, to perform the following operations on the digitised pulse signal output from the analog to digital converter: a) Baseline tracking and removal, or fixed baseline removal.
b) Detection of incoming pulses.
c) Computation of the energy of each detected pulse.
d) Accumulation of the computed energy values into an energy histogram (energy histogram)
e) Output of the accumulated histogram values each time a gate or other timing signal is received
f) Reset of the histogram values for the next collection interval.
3. Intensity Image
[0112] The intensity value, or more specifically transmission value, is computed from the energy spectrum generated for each detector i at each gate interval j according to:
R i'j = SS (Equation s) where the summations are performed over all histogram bins B (or equivalently, over all energies E), for the received energy spectra (1(B)) and reference energy spectra (I0(B)). [01 13] Elements within the intensity image may be classified as: a) Impenetrable,
Figure imgf000030_0003
b) Empty, or nothing in beam, and set to 1.
Figure imgf000030_0004
The thresholds may be pre-set or user configurable.
Figure imgf000030_0002
4. High Contrast Images
[01 14] Through use of a full energy spectrum, intensity images with varying contrast are generated based on integrating the received spectrum across different energy bands. In existing dual energy X-ray scanners, the system can only utilise the broad energy range inherent in the detector material. When a full energy spectrum is available, arbitrary energy ranges can be used to generate associated intensity images in that energy range. Specific energy ranges can then be defined in order to best isolate and display particular material types, with energy ranges tuned, for example, for organic material, inorganic material, or light, medium or heavy metals.
[01 15] The high contrast / high penetration images are generated for each detector i at each gate interval j according to:
Figure imgf000030_0001
where El and E2 are the lower and upper limits of energy range E12. The energy band may be user defined or pre-configured. One, two or more different energy bands may be configured to enable the user to select between images of interest.
5. Effective Z processing
[01 16] The effective Z processing involves the use of full energy spectra computed by the pulse processing electronics, combined with the energy calibration, to compute an estimate of the effective Z of the sample material. The effective Z processing is performed for every detector, and for each detector proceeds as follows (so for a 1 x N detector array, this process is repeated N times). To reduce computational requirement, the effective Z processing is only performed for received detectors i and gate intervals j that are not declared either impenetrable or empty.
5.1. Preliminary operations.
[0117]
1. With reference to Figure 10, compress the energy spectrum data (400) using an FFT, and discard all but the first N bins (which are selected such that the discarded bins contain little or no signal). Note: this step is optional, but for a system configuration where effective Z is computed on a central processing computer, it enables a significant communication bandwidth reduction. Transfer of 32 complex FFT bins for a 512 point histogram requires only 1/8 of the communication bandwidth.
2. Perform spectrum integration (402), by averaging a number 2S+1 of received FFT'ed energy spectra. This spectrum integration increases the measurement time available for computing the effective Z without reducing the spatial resolution at which the intensity image is computed. The integration is done across gate intervals j-S < j < j + S, so as to perform a moving average centered on gate interval j. If no integration is required, set S = 0.
3. Perform pileup reduction (403). The FFT is the first stage of the pileup reduction, which is not required if data compression has already been achieved using an FFT. The pileup reduction can be achieved with a suitable algorithm as outlined below.
4. If desirable, apply a FFT domain phase shift (404) in order to achieve a desired lateral shift of the energy spectrum. This step has been found to be desirable where a count rate specific baseline shift exists. Note: multiplication by a linearly increasing (with FFT bin) phase term in the FFT domain results in a lateral shift after iFFT. The extent of the lateral shift is determined by the slope of the linear increase.
5. Prior to iFFT, apply a frequency domain window (405). This window can be used to design a desired smoothing of the energy spectrum. The window design process is outlined below. A good window has been designed to achieve a smooth filtering of the energy spectrum. Filtering of the noise in the energy spectrum allows the possibility of using a reduced number of energy bins in the effective Z computation for overall improvement in computational efficiency.
6. Zero pad the FFT data, insert the complex conjugate into the second half of the FFT buffer (406), and apply iFFT (407). At this point a smoothed energy spectrum is obtained in the form of a histogram.
The zero padding inserts data that was truncated after the FFT. It is not essential to insert zeros for all truncated bins. For example, padding less zeros can produce a smaller FFT buffer which it is more computationally efficient to compute the IFFT. For a real vector x, and FFT size 2N, the elements N+2 to 2N of the FFT output are the complex conjugate of the elements 2 to N. Here N+l will be one of the elements set to zero by the zero padding.
7. Subtract the residual spectrum for each detector. As noted previously, this removes any spectrum that would be present even in the presence of a completely blocking material.
8. Apply the energy calibration curve/function (408) to convert the histogram bins into energy values. Note: Alternatively the energy calibration can be applied within the effective Z routine itself. At this stage the output is a smooth calibrated energy spectrum (409).
9. If required, perform spectrum integration across adjacent detectors so integration over 2P+1 energy spectra for detectors i-P < i < i + P. While integration over gate intervals can be performed in the FFT domain, integration over adjacent detectors can only be performed after energy calibration has been applied, since the raw histogram bins of adjacent detectors may not correspond to the same energy. By performing 2D spectrum integration the material identification performance can be improved compared to performing the effective Z processing on a single pixel.
5.2. Reference spectrum measurement.
[0118] In order to compute effective Z (and also the intensity / high contrast images), a reference spectrum is obtained with X-rays on, but before the sample reaches the X-ray beam. Within a given machine design, there will be a delay between the time X-rays are turned on and when the sample reaches the X-ray beam during which the reference spectrum can be collected. The process is as follows: 1. Turn on X-rays.
2. Wait for X-ray beam to stabilise. This can be achieved by a time delay or by filtering X-ray counts until the variation diminishes below a specified threshold.
3. Collect and sum N X-ray energy spectra I0 (E, n) (that is, collect the energy spectrum recorded at the end of N successive gate intervals) at the output of pulse processing electronics.
4. Divide the sum of spectra by N to compute an average reference spectrum so
Figure imgf000033_0001
where is the reference number of counts at energy E, N is the number of gate
Figure imgf000033_0003
intervals, and E is the energy level of the X-rays.
[01 19] If at any time during the reference collection a sample is detected in the X-ray beam, then the accumulation of reference spectra ceases and the average of M collected spectra can be used for the reference, or the measurement terminated if M is insufficient.
5.3. Load or create a table of mass attenuation constants
[0120] The mass attenuation constants for a given effective Z and given energy define the extent to which the given material Z will attenuate X-rays of energy E. In particular, the intensity of received energies at a particular energy will be given by:
Figure imgf000033_0002
where 1(E) is the received number of counts at energy E, I0 (E) is the reference number of counts at energy E, ma(Z,E) is the mass attenuation constant for material with effective atomic number Z at energy E, p is the material density and x is the material thickness relative to the reference thickness used in the creation of the mass attenuation data.
[0121] Mass attenuation data is available at a finite (small) number of energies, perhaps every 10, 20 or 50 keV, whereas the energy spectra created by the method disclosed in this disclosure may be generated at energy spacing as little as 1 keV or even less. In practice a finite number of these energy values will be selected for use in the effective Z computation. [0122] In order to achieve a smooth mass attenuation table at all energies in the energy spectrum, data for intermediate energies for each Z are obtained using cubic spline interpolation or other suitable interpolation method. The mass attenuation values as a function of energy are considered sufficiently smooth that a cubic spline is a good interpolation method to apply.
5.4. Effective Z computation
[0123] The effective Z processing then proceeds as follows:
1. For each detector, and each gate period (a specified detector at a specified gate period defining a pixel in the resultant image), a calibrated energy spectrum will be measured as outlined in the "preliminary operations" section. Effective Z processing is not performed for energy spectra classified as impenetrable or empty.
2. Determine a set of energy bins to be used for effective Z computation.
a. Based on the received spectrum, identify the energy region where sufficient counts are received.
b. These will be the spectrum bins where the counts exceed some predetermined threshold.
c. Alternatively, determine the energies where the transmission (ratio of received to reference spectrum) exceeds a threshold.
3. For each Z value for which mass attenuation data is available at each of the energy bins, perform the following operations:
a. Estimate the material thickness for the assumed Z. One possible method is to estimate the thickness at one energy value E according to
Figure imgf000034_0001
where 1(E) is the received number of counts at energy E, /o (£") is the reference number of counts at energy E, ma(Z,E) is the mass attenuation constant for material with effective atomic number Z at energy E, p is the material density and x is the material thickness relative to the reference thickness used in the creation of the mass attenuation data. An improved thickness estimate can be obtained by averaging the thickness estimate at a number of energies to reduce the impact of noise at the single energy. It is not desirable to estimate x explicitly, the combined parameter px is sufficient.
b. Compute a predicted spectrum for this Z, based on the reference spectrum recorded previously, the thickness parameter and the ma table according to
Figure imgf000035_0001
computed at all selected energies E, where l(Z, E) is the predicted spectrum. c. Compute a cost function for this Z as the sum of the squared errors between the received spectrum and the predicted spectrum under the assumption of material Z
Figure imgf000035_0002
where C(Z) is the cost function, and w(E) represents weights for each sum of the squared errors between the received spectrum and the predicted spectrum.
The weights w(E) can be chosen to be unity, or alternatively w(E) = 1(E) will result in a cost function that gives lower weight to regions of the received spectrum where the number of counts is small, and greater weight to regions where more counts are received.
4. For this pixel (constituting an energy spectrum received from a specific detector during a specific gate period), compute the estimated effective Z as the Z value which minimises the cost function:
Figure imgf000035_0003
[0124] It should be noted that there is no particular requirement for effective Z to be integer, and in fact the mass attenuation table may contain values for non-integer values of Z representing composite materials. However, it is clearly not possible to represent a continuum of possible Z values in a finite table. In order to compute Z to arbitrary precision, it is possible to interpolate the cost function to the required resolution using an appropriate interpolation algorithm. The value of Z chosen is then the value which minimises the interpolated cost function. The cost function C(Z) is a smooth function, and therefore an actual floating point or continuous value of Z which minimises this smooth function can be reliably predicted from the curve via some form of interpolation.
[0125] In addition, it is also noted that step 3 above indicates the cost function is computed for all available Z values in the mass attenuation table. In practice, depending on the behavior of the cost function, efficient search methods can be applied to reduce the computational requirements. Such methods include one or more of the following:
1. Gradient search
2. Best first search
3. Some form of pattern search
[0126] The cost function form has been chosen so as to be relatively insensitive to noise on the spectrum.
6. Effective Z Processing Using Material Calibration
[0127] In practice, due to detector and processing characteristics that can be difficult to characterise, it can be difficult to achieve accurate energy calibration across all detectors, all count rates and all spectrum bins.
[0128] An alternative method has been developed whereby the system is calibrated using varying thickness samples of known materials. The aim is to calibrate the expected received spectra as a function of material, material thickness, and energy histogram bins. This avoids the requirement for absolute energy calibration, and also largely avoids the effect of spectrum shift with count rate (if present). The need for pileup removal may also be eliminated.
6.1. Material (Self) Calibration Process
[0129] Ideally, with good gain calibration, the received spectra from all detectors are consistent with each other, and so it is only desirable to obtain calibration data at one detector for use at all detectors. In practice, it is likely to be desirable to obtain calibration data for groups of adjacent detectors or possibly every detector, depending on the consistency between detectors. [0130] The first step in the calibration process is to obtain a reference spectrum I0 (B) at each histogram bin B, with no material in the X-ray beam for the detector(s) to be calibrated. Histogram bins will now be denoted by B rather than E to denote that there is no requirement to calibrate the bins in terms of their exact energy.
[0131] Then, for each material, to calibrate:
1. Ascertain the effective Z of the material (either by independent measurement, or by specification of material purity).
2. Obtain a "step wedge" of the material. That is, a sample of the material that comprises a series of steps of known thickness x. The largest step is ideally sufficient to reduce the X-ray beam to a level where it can be considered impenetrable. Note: other material samples can be used, but such a step wedge is a convenient form to calibrate against.
3. Scan the step wedge at the required detector location. The result will be a series of uncalibrated energy spectra recorded along each step of the material (the number of spectra will depend on the sample dimensions, the scanning speed and the gate period).
4. Sum the spectra received on each step to minimise the noise in the spectra. These spectra are denoted I(Z, B, x), since they are a function of the material, the histogram bin and the material thickness. Note also that I(Z, B, 0) is just the reference spectrum
Io(B).
5. Compute the transmission characteristic for all materials, histogram bins and thickness as
Figure imgf000037_0002
6. Compute the total transmission as a function of Z and x as
Figure imgf000037_0001
again noting that R(Z, 0)
[0132] The tables of Tx(Z, B, x) and R(Z, x) together form the calibration tables that are used to estimate effective Z at each pixel (detector/gate interval). As previously stated, they may or may not be a function of detector also, depending on the equivalence of data from all detectors.
[0133] Clearly it is desirable to calibrate against samples of all possible materials, however in practice only a subset of the full continuum of materials and mixtures can be sampled. To achieve table entries for intermediate Z values it is desirable to interpolate both Tx and R functions to intermediate values of Z to expand the table coverage.
[0134] Having obtained the calibration tables, it is now possible to estimate effective Z for an unknown material sample as follows.
6.2. Preliminary operations.
[0135] Preliminary operations are substantially the same as described above, with the following comments:
1. It may not be desirable to perform pileup removal.
2. It may not be desirable to perform lateral spectrum shift to compensate for count rate dependent baseline shift.
3. The frequency domain window is still required prior to iFFT.
4. The energy calibration curve is not applied, as there is no requirement for absolute energy calibration with this method, but removal of residual spectrum may still be required.
5. Integration of spectrum can be performed across gate intervals, and across detectors as described below.
6. The received spectrum will be denoted 1(B), the intensity in a series of histogram bins B. The use of B differentiates from the use of E for the previous section where the histogram bins are calibrated in terms of their actual energy.
6.3. Reference spectrum measurement.
[0136] The reference spectrum is obtained in exactly the same manner as described above, but is now denoted I0(B), denoting the use of histogram bins, rather than energy. Effective Z computation. The effective Z processing then proceeds as follows: For each detector, and each gate period (a specified detector at a specified gate period defining a pixel in the resultant image), an uncalibrated energy spectrum 1(B) will be measured as outlined elsewhere. Again, effective Z processing is not performed for energy spectra classified as impenetrable or empty.
Determine a set of histogram bins to be used for effective Z computation:
a. Based on the received spectrum, identify the region where sufficient counts are received.
b. These will be the spectrum bins where the counts exceed some predetermined threshold. Choose B: 1(B) > Imin
c. Alternatively, determine the bins where the transmission (ratio of received to reference spectrum) exceeds a threshold.
d. Alternatively, use all available histogram bins and apply a weighting in the cost function to remove undesired bins from the cost calculation.
e. Note: ultimately reducing the total number of histogram bins processed will achieve improved computational efficiency, so use of every bin is not ideal. Compute the total received X-rays as a ratio to the reference.
Figure imgf000039_0001
For each Z value for which calibration data is available, perform the following operations:
a. Estimate the material thickness from the total received transmission R and the calibration table values of R(Z, x) for this material Z. This achieved via i. Interpolation of the curve of R(Z, x) at the measured value of R to obtain a corresponding x, via for example cubic spline interpolation. ii. From the calibrated R(Z, x) obtain a functional form x = f(R, Z) to compute x as a function of material and transmission. b. From the table of R(Z, x) find xx andx2 such that R(Z, xx) < R < R(Z, x2).
Note that x1 = 0 corresponds to the reference spectrum, and if the received transmission R is smaller than a table entry then use the final 2 entries for xx and x2, and the result will be an extrapolation to a thicker material. c. Now use the calibrated transmission tables Tx(Z, B, x) to determine local mass attenuation coefficients for each histogram bin according to:
Figure imgf000040_0001
d. Then compute an expected received spectrum according to
Figure imgf000040_0002
This expected received spectrum is an interpolated received spectrum between the nearest two calibration spectra, but based on the different attenuation observed at each bin. Other forms of interpolation between spectra could be used, but the material specific interpolation used here provides a superior interpolation result.
e. Compute a cost function C(Z) for this Z as the sum of the squared errors between the received spectrum and the predicted spectrum under the assumption of material Z
Figure imgf000040_0003
The weights w(B) can be chosen to be unity, or alternatively w(B) = 1(B) will result in a cost function that gives lower weight to regions of the received spectrum where the number of counts is small, and greater weight to regions where more counts are received.
5. For this pixel (constituting an energy spectrum received from a specific detector during a specific gate period), compute the estimated effective Z as the Z value which minimises the cost function:
Figure imgf000040_0004
[0138] It should be noted that there is no particular requirement for effective Z to be integer, and in fact the self calibration table may contain values for non-integer values of Z representing composite materials. However, it is clearly not possible to represent a continuum of possible Z values in a finite table. In order to compute Z to arbitrary precision, it is possible to interpolate the cost function to the required resolution using an appropriate interpolation algorithm. The value of Z chosen is then the value which minimises the interpolated cost function. The cost function C(Z) is a smooth function, and therefore an actual floating point or continuous value of Z which minimises this smooth function can be reliably predicted from the curve via some form of interpolation.
[0139] The same form of efficient search methods can be used to reduce computation and avoid an exhaustive search over all materials Z in the calibration table.
6.5. System Adaptation
[0140] Some system parameters will vary over time, so the system adapts in order to maintain calibration over time:
1. Gain calibration update.
a. Calibration spectra are measured during periods when X-rays are off b. Gain is updated according to changes observed in the measured spectra c. Gain is A*old gain + B*new gain, where A+B=l, and B will be small to avoid noise and allow slow adaptation.
2. Pulse calibration update.
a. Periodically a new pulse calibration may be performed, although pulse parameters have been found to remain sufficiently constant over time that at most daily or perhaps weekly or monthly recalibration of pulse parameters may be required.
3. Baseline offset update
a. This may be done during periods when X-rays are off, in the same way as the initial calibration is performed - a short dataset is required for baseline offset.
b. Adapted continuously via a baseline tracking algorithm as described below.
4. Energy calibration, count rate dependent spectrum shift, pileup parameters and residual spectrum may require occasional off-line recalibration. It may also be found that for a given machine these rarely, if ever, require recalibration. 7. Effective Z Processing Example
[0141] The following is an overview of the process used to calibrate detector boards, in particular implementing a self calibrating process and with the option of using 'floating point' effective Z computation:
1. Obtain calibration wedges of known material, ideally pure or close to pure
elements. For the current calibration, 3 materials were used:
a. Carbon (Z=6)
b. Aluminium (Z=13)
c. Stainless steel (Z = approximately 26)
2. The step dimensions are chosen with the following in mind:
a. A width of 30cm is used to ensure a large number of detectors could be calibrated with 1 calibration scan.
i. Projection mode is used to effectively increase the number of pixels that could be calibrated in
ii. With 30cm wedges, 2 calibration heights can cover the 5 detector boards with sufficient overlap to avoid edge effects.
b. Step heights were determined to achieve a reasonably uniform transmission spacing from approximately < 0.5% up to 95%.
i. For carbon, achieving less than 0.5% transmission required
approximately 300mm of material.
ii. For heavier metals, achieving 95% transmission required very thin samples, 0.5mm and less. For metals such as tin (not used here) this was extremely difficult.
c. Step length is 50mm. When scanned at 4% normal speed, the scan speed is 8mm per second, so approximately 6 seconds of data could be collected from each step. This is necessary to ensure very accurate calibration spectra, given the accuracy required from the eventual effective Z processing.
3. The material wedges are scanned, and the resulting data is processed offline in Matlab as follows: a. For each pixel of each scan, determine the start and end locations of the step.
b. Allow some margin, to avoid any effects near the step edges.
c. For each step identified:
i. Extract the binary data corresponding to the measured spectrum at each slice of the step.
ii. Integrate all data, so as to establish a very accurate spectrum (with >5 seconds of data)
iii. Compute the corresponding total intensity (relative to a long term average of the source spectrum measured during the same calibration run with nothing in beam)
d. Create a table for each step (intensity) of each material, containing: i. A mapping of step thickness to step intensity. This table is used to interpolate any measured spectrum into an equivalent material thickness.
ii. A series of calibration spectra, including a reference spectrum. Each spectrum represents a material thickness, from which spectra for intermediate material thicknesses can be interpolated.
[0142] The calibration data for the 3 materials with Z = 6, 13, 26 is adequate for the production of a 3 colour image, classifying the material as either organic (close to Z=6), inorganic / light metal (close to Z=13) or metal (close to Z=26). In order to achieve accurate effective Z estimation to separate materials down to +/- 0.2 Z or better, it is necessary to obtain calibration data from a much larger set of materials, from which a continuous estimate of Z could then be obtained. It is not practical to run calibration scans for all materials from Z=3 to Z=92, so a range of additional calibration data sets were obtained by interpolation. Calibration sets for Z = 3 to Z = 13 were obtained from interpolation/extrapolation of the carbon and aluminium data sets. Calibration sets for effective Z = 13 to Z = 50 may be obtained from interpolation of the aluminium and stainless steel data sets. [0143] For every pixel in the scanner, the procedure for obtaining the additional calibration data sets is as follows:
1. For each of Z = 6, 13, and 26, interpolate the calibration spectra to a new set of intensities. For the current demonstration, the intensities (in percent) used were: 95, 90, 80, 70, 60, 50, 40, 30, 20, 15, 10, 6, 4, 2, 1 ,0.5, 0.2. At this point there is now a calibration table for each material at a common set of intensities. The process now is to create calibration tables for other materials at those same set of common intensities.
2. For the current embodiment the set of required materials is Z = 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 16, 18, 20, 22, 24, 26, 30, and 35. The spectrum for material Z, at histogram B, at thickness x (which here corresponds to one of the defined transmission levels) is denoted 7(Z, B, x). For each required material, for each intensity, the new interpolated material spectrum is obtained as follows:
a. If Z new <= 13
Figure imgf000044_0001
b. If Z new > 13
Figure imgf000044_0002
3. The new tables are then included in the total set of calibration data used for
effective Z processing.
4. All calibration tables are saved in a file format appropriate for input to the
effective Z software.
[0144] There are some important points to note here:
1. For Z < 6, the process is an extrapolation rather than interpolation - one of the coefficients becomes negative, while the other is greater than 1. While it appears to work acceptably well down to Z=3, some care needs to be taken with extrapolation as it may diverge quickly.
2. Similarly for Z > 26, the process is extrapolation. Here it would be better to
include calibration data for tin (Sn) and then also lead (Pb) to fill out the available calibration data. The challenge with these higher Z materials is to get reliable calibration curves at 90% transmission - the material sample must be very thin to achieve this transmission.
3. By including values of Z below 3, the cost function behaves reasonably well in the region of interest around Z=6. This ensures the process for computing
continuous/floating point effective Z is able to correctly determine Z values around 6 - essentially avoiding anomalies at the edges, or at least pushing them to Z values outside the range of interest.
4. More sophisticated interpolation may be required to avoid common point of
intersection - this is acceptable for lower Z but ceases to hold true when moving into the metals. It may be necessary to interpolate several spectra to obtain each new interpolated material.
5. The overall performance is somewhat limited by the use of just 3 actual measured materials. In fact it is quite remarkable that such excellent performance has been achieved given the spacing of the calibration materials.
6. Higher Z materials such as lead have absorption edges, and so some consideration eventually needs to be given to these materials if accurate performance at high effective Z is to be achieved. To this point, absorption edges have not been specifically incorporated into the model.
[0145] The cost function C(Z) is a smooth function, and therefore an actual floating point or continuous value of Z which minimises this smooth function can be predicted from the curve via some form of interpolation.
[0146] The results of the interpolation process are shown in Figure 19 for the 10% transmission case. It can be seen there is a very smooth progression through all of the materials, and this is what results in the discriminating power of the effective Z processing. Any measured material in the effective Z range can be placed somewhere in this set of curves, with the exact displacement from calibration curves used to determine a very accurate estimate of the material effective Z.
[0147] Implementation of floating point effective Z has been based around a quadratic interpolation, using the cost function values at the Z value that minimises the cost function, and the Z value either side of this, with some special consideration at the edges. This approach has yielded effective Z results which (with sufficient spectrum integration) have accurately resolved materials where the known effective Z difference is less than 0.2.
[0148] The process for computation of a continuous / floating point estimate of effective Z is as follows:
1. Compute the cost function C(Z) at every Z value in the calibration table.
2. Find the value of Z for which C(Z) is minimised.
3. Find the values of Z and associated cost function values for the Z values either side of the Z value which minimises C(Z).
4. The coefficients a quadratic model are estimated, where the model in the region of the minimum is:
Figure imgf000046_0005
where n is noise on the cost function. This is in turn modelled in matrix equation form for the 3 Z values Zl, Z2 and Z3 and associated cost function values C I, C2 and C3, with
Figure imgf000046_0003
and where
Figure imgf000046_0001
and the solution is obtained via matrix inversion as:
Figure imgf000046_0004
The general form inv(H'H)H'C is used to accommodate the case where more than 3 Z and C values are used to estimate the quadratic coefficients.
5. Compute the value of Z which minimises the quadratic function. The value of Z at which there is a turning point is simply:
Figure imgf000046_0002
By design, one of the Z values (usually Z2 except at edges) is a minimum, so the resulting optimal Z can be assumed to be a minimum rather than a maximum. At edges, there can be issues, but this needs to be treated separately. In the disclosed embodiment, the following observations are made: By design, one of the Z values (usually Z2 except at edges) is a minimum, so the resulting optimal Z can be assumed to be a minimum rather than a maximum. At edges, there can be issues, but this needs to be treated separately.
If the value of Z which minimises the cost function is at either edge, then it is necessary to use 2 values to one side of the minimum Z value. The floating point calculation may then become an extrapolation to a point outside the range of Z values in the calibration table. When this occurs, the Z estimate can diverge quickly, so care needs to be taken to put a limit on then maximum or minimum Z value estimate (i.e. how much extrapolation is allowed)
The inclusion of Z values below 6 and above 26 is designed to ensure the edge effects described above do not adversely affect the Z estimates, particularly in the region of interest around Z=6.
There are likely more computationally efficient ways of obtaining the quadratic coefficients and associated estimate of the minimum. This has not been explored at this stage.
In fact, for a given set of Z values, all necessary matrices and inverses can be computed offline and stored for more efficient use, since they only depend on the Z values and Z spacing, not the measured spectrum and cost function.
The quadratic model is acceptable where the cost function is well behaved, smooth and relatively noise free. Where very short integration times are used to collect the energy spectrum histograms, the cost function can become noisy and may converge on a local minimum due to noise. In this case, and in general, a more sophisticated interpolation model may be required to smooth the cost function and avoid noise effects. This may involve more than 3 points in the interpolation process. [0150] The quadratic model is just a model to ensure a consistent effective Z is obtained for a particular material. It is not intended to be an accurate functional model of the cost function behaviour, and it is not considered necessary. The principle objective is to obtain an estimate of effective Z that is consistent for a particular material, and enables reliable separation of closely spaced materials. The quadratic model achieves this objective.
[0151] The floating point effective Z algorithm was tested on a range of material samples. There were several observations made about the performance.
1. At high transmission values - corresponding to very thin / low attenuation
samples - the cost function could become noisy, and also the calibration curves for higher Z often poorly interpolated to >90% transmission. As a result, the output tended to over emphasise higher Z values.
2. At very low transmission, and in the vicinity of large changes in transmission
levels such as near the edge of metal blocks, some scatter could be present in the received spectrum resulting in an output biased towards organic, even where the material was known to be metal.
3. A-priori knowledge of likely effective Z gave the following heuristics:
a. High transmission is more likely to be an organic or lower Z material b. Low transmission is more likely to be a higher Z material, since very large thicknesses of low Z material would be required.
[0152] As a result of these observations, a weighting v(Z,I) was introduced in order to tune the cost function as a function of both intensity and Z. These cost function weights are tuned in order to ensure the effective Z output is as required for known test samples. In the disclosed embodiment the implementation has been confined to 3 discrete regions:
1. High transmission, I > high threshold
For high transmission, the output was found to be somewhat biased towards high Z when fairly thin organic materials were scanned. Hence the cost weights were low for organics, and increasing at higher Z.
2. Mid Transmission, low threshold < I < high threshold
In the mid range of transmission, the output was generally consistent with the expected effective Z, hence only a very slight cost function weighting was applied.
3. Low Transmission, I < low threshold
At very low transmission, it was found that higher Z materials were occasionally mis-identified as low Z materials. This was especially true near edges of metal blocks, where scatter could allow an excess of low energy X-rays to reach the detector. As a result, at low transmission, the cost weights were designed to increase the cost of low Z materials to produce the higher Z output more consistently. One side effect of this approach is that very thick organic materials start to be identified as metals at low transmission. This can only really be overcome by removing the underlying source of the mis-identification, being an excess of low energy X-rays in the received spectrum.
8. Effective Z Processing Implementation
[0153] The following sections outline in further detail the individual processing stages and algorithms. Figure 28 indicates an overview of the various optional processing stages that may be implemented in the present method.
8.1. Tiling
[0154] The tiling algorithm is effectively a block averaging function. The purpose of the tiling algorithm is to average the floating effective Z image over an area (mm2) that represents the smallest region required to be detected of a constant intensity and material composition. The tiling algorithm generates tiles with 50% overlap to ensure that we always capture the object of interest. The tiling algorithm estimates the mean and standard deviation over rectangular tiles in the floating effective Z image. The tile width and height are defined by the user. Tiles are overlapped by 50% in both vertical and horizontal dimensions. Given an image size Nr by Nc pixels, and a tile dimension Tr by Tc pixels, the number of tiles in the vertical dimension is (Nr/Tr)*2. The tile dimensions must be even valued to ensure 50% overlap. The tiling algorithm executes a loop that indexes into each tile and calculates the mean and standard deviation of all pixels in the tile.
[0155] The choice of tile dimensions essentially comes down to a compromise between: 1. The dimensions of the smallest object that must be detected; and
2. The effective Z resolution required. The effective Z variance has been observed to reduce almost linearly with number of effective Z pixels averaged, so larger areas yield better effective Z resolution.
[0156] In addition, the idea of tiling and clustering has been used to avoid the need to implement sophisticated image segmentation at this time. It was felt that to get accurate effective Z measurements would in any case require large contiguous blocks of uniform material, so the tiling and clustering approach would be only marginally inferior to full image segmentation. Nonetheless, image segmentation may ultimately prove advantageous for highly irregular shapes, especially where some more sophisticated object recognition approaches may be used in conjunction with effective Z measurement.
8.2. Clustering
[0157] The clustering algorithm groups tiles that have a common effective Z and are spatially connected. The purpose of the clustering algorithm is to detect objects that span areas larger than the minimum object size as defined by the tile dimensions. Connectedness is defined along edges. Connected tiles are assigned a common cluster ID. The output of the clustering algorithm is a cluster map and a cluster table. The cluster map is a matrix of connected tiles with associated cluster IDs. The cluster table holds information on each cluster ID including the number of tiles in the cluster, and the vertical and horizontal extent of each cluster.
[0158] The clustering algorithm performs row-wise scanning of the tiled image. If tile P(r,c) is connected to a tile in the set A = (P(r,c-1), P(r-l,c+l), P(r-l,c), P(r-l,c-l)} then it is assigned the cluster ID. If P(r,c) is not connected to the set A but is connected a tile in the set B = (P(r,c+1), P(r+l,c-l), P(r+l,c), P(r+l,c+l)} then the tile is assigned a new cluster ID. In the case of connectedness with tiles in set A, it is possible for the P(r-l,c+l) to have a different cluster ID to others in the set. In this case a cluster merge is performed. This is achieved by a simply replacing one cluster ID with the other, the specific order is unimportant. The sets A and B are adapted at eight boundary conditions, four along the image edges and four at the image vertices. [0159] Figure 27 depicts the formation of clusters, where single tiles are ignored.
8.3. Target Material Detection
[0160] The target material detection algorithm is a nearest neighbour classifier. The algorithm classifies individual tiles. There are two steps in the algorithm, training and classification. The training stage establishes a lookup table mapping normalised intensity to floating effective Z for a range of 'target' materials of interest. In the current implementation, the lookup table is approximated as a quadratic fit, for which only the quadratic coefficients are stored.
[0161] During the classification stage, the input is the normalised measured tile intensity (Imeas), the measured tile floating effective Z (Zmeas), and a maximum effective Z classification error (deltaZ). For each material in the training set, the classifier declares positive classification if abs(Ci(Imeas)-Zmeas) < deltaZ, where Ci is the quadratic function associated with the 1th target material.
[0162] The use of both intensity and effective Z in the target profile is an important aspect of this approach. The effective Z is typically not constant with material thickness, and so including the intensity (related to thickness) provides a two dimensional test with far superior discrimination than effective Z alone.
[0163] Figure 20 shows the effective Z vs intensity for a range of material samples tested, along with the quadratic interpolation. Here the variation of effective Z with intensity is clear. While the example materials are more relevant for a baggage scanning application, the principle is the same in relation to discrimination between mineral types.
8.4. Edge Detection
[0164] The purpose of the edge detection algorithm is to ensure that the moving average window does not straddle material boundaries. The edge detection uses amplitude transitions in the intensity image to declare material edges. The input to the edge detection algorithm is the intensity image. Edges are only detected in the horizontal dimension. The reason for not detecting edges in the vertical dimension is that the moving average window only operates in the horizontal dimension. Edges in the intensity image are computed for each detector. A first order gradient operator is used to detect edges. The gradient operator mask width, and the gradient threshold, are defined by the user. Given the following edge mask L(c) indexed on columns as depicted in figure 29, the gradient is G = sum(L(c).*Inorm(c)) where Inorm is the normalised intensity. An edge is declared when abs(G) > g where g is a user defined threshold.
8.5. Moving Average
[0165] The purpose of the moving average algorithm is to filter the intensity histograms for each detector so as to increase the effective signal-to-noise ratio. The algorithm generates a filtered intensity histogram a slice k, for each detector, by averaging the measured intensity histograms over a symmetric window centred on slice k. The edge detector plays an important role in ensuring the moving average window does not straddle different materials. If a window overlaps an edge the average is only calculated up to the edge boundaries. The width of the window can be set by the user. On edges, no averaging is performed. Figure 30 illustrates the behaviour of the moving average as it transitions over an edge.
[0166] One embodiment that may be more computationally efficient is an adaptive moving average approach:
1. Compute effective Z on every slice in the presence of edges.
2. Compute effective Z based on 50% overlapping of moving average windows (eg every 5 pixels for 11 pixels MA length).
[0167] This can provide 3-5x improvement in computational speed depending on exact configuration. Additional Details
1. High Rate Pulse Processing
[0168] In principle, any suitable method of high rate pulse processing can be used within the embodiments described herein. However, the high X-ray flux present in typical X-ray screening systems results in a high pulse count rate, and a high likelihood of receiving overlapping X-ray pulses.
[0169] Pulse pile-up has long been a problem to contend within applications of high rate radiation spectroscopy. Traditional approaches to pulse shaping use linear filters to shorten pulse duration which can significantly reduce S R and are therefore limited to output rates of a few hundred kc/s. An alternate approach to processing the data from radiation detectors is based on the idea of mathematically modeling data corrupted by pulse pile-up and solving for the required model parameters. By recovering rather than discarding data corrupted by pulse pile-up this technique enables high throughput, low dead-time pulse processing without the traditional loss in energy resolution.
[0170] The disclosures of international patent publications WO2006029475, WO2009121130, WO2009121131, WO2009121132, WO2010068996, WO2012171059 and WO2015085372 are useful in the current invention achieving high rate pulse processing with reduction in pulse pileup rejection and are all incorporated herein by reference in their entirety as useful in embodiments of the current invention as if repeated here verbatim, and the applicant reserves the right to incorporate any terminology and concepts disclosed in the above international patent publications in future claim language amendments in the current application.
[0171] The following account includes a selection from the techniques disclosed in the above international patent publications adapted to the current invention, but persons skilled in the art will appreciate that all of these techniques are potentially useful and choice among the alternative approaches is guided by satisfaction of various competing performance constraints including processing speed, energy determination accuracy and maximum count rate. 1.1.Model-Based, High-Throughput Pulse Processing - Method 1
[0172] The algorithm briefly described here, and in more detail in WO2006029475 (incorporated by reference), for processing the data from radiation detectors is a model - based, real-time, signal-processing algorithm that characterises the output of the radiation detector y[n] as shown below:
Figure imgf000054_0001
[0173] The digitised radiation detector time series (y) is modeled as the sum of an unknown number of radiation events (N), with random times of arrival (τ), and amplitudes (a), interacting with a radiation detector, that have an expected pulse shape (h) and with a noise process (co).
[0174] Therefore, so as to fully characterise the digitised output of the radiation detector, it is desirable to estimate: the expected impulse response of the detector; the number of events in the digitised detector time series; the time of arrival of each of those radiation events; and the individual energies of each event. Once these parameters have been determined, the digitised detector data can be accurately decomposed into the individual component events and the energy of each event determined.
System Characterisation
[0175] Calibration of the detector is the first stage of the algorithm; it takes as input the detector time series data and determines the unit impulse response of the detector (the expected pulse shape from the detector). Refer to Pulse Parameter Calibration for a more detailed summary of the pulse calibration process.
Pulse Localisation
[0176] After the unit impulse response of the detector has been determined this is used by the Pulse Localisation stage to determine the number of events in the digitised detector data stream and their TOA relative to each other. [0177] The detection of events in the digitised detector waveform is accomplished by fitting an exponential model to a fixed number of data points. After the System Characterisation stage the exponential decay of the pulse tail is well characterised. The detection metric (the signal ultimately used to make a decision as to whether a pulse has arrived or not) is formed by fitting an exponential curve to a specified number of data points. This fixed length 'detection window' is run continuously over the digitised detector data and the sum of the squares of the error is computed (this can also be thought of as the sum of the square of the fit residual). This operation results in three distinct modes of operation:
1. Baseline Operation: processing data samples when no signal is present. As the data can be quite accurately modeled by an exponent the sum square of the error is at a minimum and remains quite constant.
2. Event Detection: when a radiation event enters the detection window the data can no longer be accurately modeled as an exponent (the data could be consider non- differential at T=0 the exact arrival time of the radiation event). Consequently the sum square of the errors will increase. This detection metric will continue to increase until the radiation event is positioned in the middle of the detection window.
3. Tail Operation: when processing data on the tail of a radiation event the data points are quite accurately modeled as an exponent. Consequently the sum square of the error returns to the same level as the Baseline Operation mode.
[0178] Using such an exponent pulse fitting operation on the digitised detector produces an ideal detection metric. It remains low during baseline, rises rapidly in response to an event arrival and decays rapidly once the rising edge of the radiation event has subsided. Furthermore, by increasing the number of ADC samples in the fixed length detection window it is possible to suppress the detector noise and accurately detect very low energy events. However, the width of the detection metric (in samples) varies proportionally with the detection window. Consequently, as the detection window gets wider the ability to distinguish two closely separated pulses is diminished. Quadratic Peak Detection
[0179] The final stage of Pulse Localisation is making a decision on the exact number and time of arrival of each of the radiation events in the detector data stream. One approach would be to apply a simple threshold to the detection metric and declare a pulse arrival at the nearest sample to the threshold crossing. However, a simple threshold crossing is susceptible to noise and only provides ±0.5 sample accuracy in determining the pulse arrival time. To have more accurate pulse arrival time and to be robust against noise (of particular importance when dealing with low energy signals close to the noise floor) a quadratic peak detection algorithm can be used. Such an approach fits a quadratic to a sliding window of N samples of the detection metric (N maybe equal to 5). In order for a peak to be declared we examine the decomposition and declare a peak if the curvature is within a permitted range, the constant is over a threshold, and the linear term has change from positive to negative. The coefficients can also be used to determine sub-sample time of arrival.
Pulse Energy Estimation
[0180] The Pulse Energy Estimation stage determines the energy of all the radiation events in the detector data stream. As its input it uses: the a priori knowledge of the detector unit impulse response; the number of events; and their individual time of arrival data. The digitised detector data of equation (18) may also be written in matrix form as:
Figure imgf000056_0003
Figure imgf000056_0002
where A is an M x N matrix, the entries of which are given by
Figure imgf000056_0001
[0181] Thus, the columns of matrix A contain multiple versions of the unit impulse response of the detector. For each of the individual columns the starting point of the signal shape is defined by the signal temporal position. For example, if the signals in the data arrive at positions 2, 40, 78 and 125, column 1 of matrix A will have '0' in the first row, the 1st data point of the unit impulse response in the second row, the 2nd data point of the unit impulse response in the 3rd row, etc. The second column will have '0' up to row 39 followed by the signal form. The third column will have '0' up to row 77; the fourth column will have '0' up to row 124 and then the signal form. Hence the size of matrix A is determined by the number of identified signals (which becomes the number of columns), while the number of rows depends on the number of samples in the 'time series'.
[0182] Once the system matrix has been created it is possible to solve for the desired energies of each radiation event by calculating the pseudo inverse of matrix A: x = inv(A'A) .A' . y (Equation 21)
Data Validation
[0183] The final functional stage of the real-time, signal-processing algorithm is the Validation stage. At this stage all the parameters that have been estimated by previous algorithmic stages (pulse shape, number of events, time of arrival and event energy) are combined to reconstruct a 'noise-free' model of the detector data.
[0184] By subtracting this model of the detector data from the actual digitised detector time series, the accuracy of the estimated parameters can be determined. Much like examining the residual from a straight line fit of a data set, if the magnitude of the residuals is small, the parameters well describe the data. However, if large residuals are observed, the detector data has been poorly estimated and that portion of the data can be rejected.
1.2. Model-Based, High-Throughput Pulse Processing - Method 2
[0185] The algorithm briefly described here, and in more detail in WO2010068996 (incorporated by reference), for processing the data from radiation detectors is a model - based, real-time, signal-processing algorithm wherein the signal processing is at least in part conducted in a transform space.
[0186] In one embodiment, the method of WO2010068996 for resolving individual signals in detector output data comprises:
obtaining or expressing the detector output data as a digital series (such as a digital time series or a digitised spectrum);
obtaining or determining a signal form (or equivalently the impulse response) of signals present in the data;
forming a transformed signal form by transforming the signal form according to a mathematical transform;
forming a transformed series by transforming the digital series according to the mathematical transform, said transformed series comprising transformed signals;
evaluating a function of at least the transformed series and the transformed signal form (and optionally of at least one parameter of the transformed signals) and thereby providing a function output;
modelling the function output according to a model (such as by modelling the function output as a plurality of sinusoids);
determining at least one parameter of the function output based on the model; and determining a parameter of the signals from the at least one determined parameter of the function output.
[0187] It will be understood by the skilled person that individual signals in detector output data may also be described as individual pulses in a detector output or in a detector output signal (in which case signal form could be referred to as pulse form).
[0188] The signal form may generally be regarded as characterising the interaction between the detector and the radiation (or other detected input) that was or is being used to collect the data. It may be determined or, if known from earlier measurements, calibrations or the like, obtained from (for example) a database.
[0189] In some embodiments, transforming the digital series according to the mathematical transform comprises forming a model of the digital series and transforming the model of the digital series according to the mathematical transform.
[0190] In certain embodiments, the method includes determining a plurality of parameters of the transformed signals, such as frequency and amplitude.
[0191] In certain particular embodiments, the transform is a Fourier transform, such as a fast fourier transform or a discrete fourier transform, or a wavelet transform. Indeed, in certain embodiments the transform may be applied somewhat differently to the signal form and digital series respectively. For example, in one embodiment the mathematical transform is the Fourier transform, but the signal form is transformed with a discrete fourier transform and the digital series is transformed with a fast fourier transform.
[0192] In one embodiment, the transform is a Fourier transform and the function is representable as
Figure imgf000059_0001
where X(k) is the transformed series and H(k) is the transformed signal form.
[0193] Thus, this method endeavours to determine a parameter of the signals and hence of as much of the data as possible, but it will be appreciated that it may not be possible to do so for some data (which hence is termed 'corrupt data'), as is described below. It will be understood that the term ' signal' is interchangeable in this context with 'pulse', as it refers to the output corresponding to individual detection events rather than the overall output signal comprising the sum of individual signals. It will also be appreciated that the temporal position (or timing) of a signal can be measured or expressed in various ways, such as according to the time (or position in the time axis) of the maximum of the signal or the leading edge of the signal. Typically this is described as the arrival time ('time of arrival') or detection time.
[0194] It will also be understood that the term 'detector data' refers to data that has originated from a detector, whether processed subsequently by associated or other electronics within or outside the detector.
[0195] The signal form (or impulse response) may be determined by a calibration process that involves measuring the detector' s impulse response (such as time domain response or frequency domain response) to one or more single event detections to derive from that data the signal form or impulse response. A functional form of this signal form may then be obtained by interpolating the data with (or fitting to the data) a suitable function such as a polynomial, exponential or spline. A filter (such as an inverse filter) may then be constructed from this detector signal form. An initial estimate of signal parameters may be made by convolution of the output data from the detector with the filter. Signal parameters of particular interest include the number of signals and the temporal position (or time of arrival) of each of the signals.
[0196] The particular signal parameters of interest can then be further refined.
[0197] The accuracy of the parameter estimation can be determined or 'validated' by comparing a model of the detector data stream (constructed from the signal parameters and knowledge of the detector impulse response) and the actual detector output. Should this validation process determine that some parameters are insufficiently accurate, these parameters are discarded. In spectroscopic analysis using this method, the energy parameters deemed sufficiently accurate may be represented as a histogram. [0198] The data may include signals of different forms. In this case, the method may include determining where possible the signal form of each of the signals.
[0199] In one embodiment, the method includes progressively subtracting from the data those signals that acceptably conform to successive signal forms of a plurality of signal forms, and rejecting those signals that do not acceptably conform to any of the plurality of signal forms.
1.3. Model-Based, High-Throughput Pulse Processing - Method 3
[0200] The algorithm briefly described here, and in more detail in WO2012171059 (incorporated by reference), for processing the data from radiation detectors is a model - based, real-time, signal-processing algorithm wherein determining a location and amplitude of pulses within the signal is achieved by fitting a function to detector output data.
[0201] The method may comprise detecting a pulse or pulses in said detector output data by:
sliding a window across the data to successive window locations;
identifying possible pulses by performing pulse fitting to the data in the window at each window location;
determining which of the possible pulses have a pulse start falling before and near the start of the respective window location and a peak amplitude exceeding the standard deviation of the noise in the window at the respective window location; and
identifying as pulses, or outputting, those of said possible pulses that have a pulse start falling one, two or three samples before the start of the respective window location and a peak amplitude exceeding the standard deviation of the noise in the window at the respective window location.
[0202] In many embodiments, the one or more functions are functions of time.
[0203] In some of those embodiments, however, the skilled person will appreciate that the one or more functions may not be functions exclusively of time.
[0204] The method may comprise providing the detector output data in, or converting the detector output data into, digital form before fitting the one or more functions to the detector output data. [0205] In one embodiment, the one or more functions are of the form:
Figure imgf000061_0001
[0206] In this embodiment, v{t) may be calculated numerically, such as by the formula
Figure imgf000061_0002
for t = 0, 1, 2, ... (with v(0) = 0).
[0207] Although mathematically the
Figure imgf000061_0003
above formula may be used to evaluate v(t) numerically. Furthermore, the above formula remains correct even when reducing in that instance to v
Figure imgf000061_0014
Figure imgf000061_0004
[0208] In one embodiment, the one or more functions are of the form:
Figure imgf000061_0005
and the method includes determining a location and amplitude of the pulse with a method comprising:
defining a reference pulse pit) as a convolution of
Figure imgf000061_0006
determining the location τ and amplitude with
Figure imgf000061_0007
Figure imgf000061_0008
[0209] The skilled person will appreciate that the present aspect of the invention contemplates different but mathematically equivalent expressions of this approach.
[0210] The sk
Figure imgf000061_0009
illed person will also appreciate that:
Figure imgf000061_0010
[0211] Expanding
Figure imgf000061_0015
gives the two equations:
Figure imgf000061_0011
where In the limit as β becomes equal to a, the constant γ becomes 1, and
Figure imgf000061_0012
equation (26) becomes This form is therefore suitable for use in a numerically
Figure imgf000061_0013
stable method for a calculating τ.
[0212] If \β— α \ is very small, care needs to be taken with the calculation of y. This may be done by summing the first few terms in the Taylor expansion:
Figure imgf000062_0001
[0213] Solving equation (26) can be done numerically, such as with a bisection method, especially since the left hand side is monotonic in τ. Determining the left hand side for different values of τ may be done by any suitable technique, such as with a Taylor series expansion for small τ. (In practice, the value of τ will generally be small because noise will generally preclude accurate characterisation of a pulse that started in the distant past.)
[0214] The linear approximation in τ of equation (26) is and is exact if
Figure imgf000062_0003
Figure imgf000062_0004
The exact, general solution (theoretically) is , the Taylor series
Figure imgf000062_0002
expansion of which is:
Figure imgf000062_0005
which is valid provided |x| < 1.
[0215] The method may comprise constraining τ by requiring that
Figure imgf000062_0008
Thus, because the left-hand side of the equation is monotonic in τ, the constraint that is
Figure imgf000062_0009
equivalent to the constraint on a and b that 0 < b < ca where the scalar c is given by
Figure imgf000062_0006
Indeed, if
Figure imgf000062_0007
Thus, it is possible to provide a constrained optimisation.
This constraint can be implemented in with the constraints that a and β are not negative and α > β.
[0216] The method may also comprise constraining the amplitude of the pul se. This can be used, for example, to prevent a fitted pulse from being too small or too large. Indeed, referring to equation (30) above, if τ is constrained to lie between -1 and 0 then A lies between Constraining a therefore constrains the amplitude A.
Figure imgf000063_0003
[0217] According to another particular embodiment, the function f is in the form of a function with three exponentials. In a certain example of this embodiment, the time constants are known and dissimilar (so fewer problems of numerical imprecision
Figure imgf000063_0004
arise), and the method includes fitting the curve:
Figure imgf000063_0005
[0218] In another example of this embodiment, the time constants are known and
Figure imgf000063_0006
in ascending order such that and fitting the function f includes using basis vectors:
Figure imgf000063_0007
Figure imgf000063_0001
[0219] For reference, if the time-constants differ, then
Figure imgf000063_0002
Figure imgf000063_0008
where
Figure imgf000063_0009
[0220] Note, however, that— unlike the previous 'double-exponential' case, in which there were two unknowns (viz. the location and the amplitude of the pulse) and two equations (coming from the two basis vectors), in this 'three-exponential' case there are two unknowns but three equations. There are therefore many different ways of inverting these equations (thereby recovering the location and the amplitude of the pulse), and generally this will be the strategy that is robust to noise.
[0221] In another particular embodiment, the function f is of the form:
Figure imgf000064_0001
wherein a and β are scalar coefficients, and the method comprises determining a and b.
[0222] This approach may not be suitable in applications in which but in some
Figure imgf000064_0008
applications it may be known that this is unlikely to occur, making this embodiment acceptable.
[0223] In one example of this embodiment, determining the location comprises determining a location f*(a, b) where:
Figure imgf000064_0002
[0224] It will be appreciated that this embodiment, which uses
Figure imgf000064_0005
has the disadvantage that these terms converge as β approaches α (unlike the terms v{t) and
Figure imgf000064_0007
in the above-described embodiment, which remain distinct. Indeed might be said to
Figure imgf000064_0006
correspond to the tail of a pulse that occurred at—∞ (whereas v{t) represents a pulse occurring at time 0).
[0225] The function f may be a superposition of a plurality of functions.
[0226] The method may include determining the pulse amplitude by evaluating f = f(t) at
Figure imgf000064_0003
[0227] Thus, the present invention may comprise a method and apparatus for estimating the location and amplitude of a sum of pulses from noisy observations of detector output data. It presented the maximum-likelihood estimate as the benchmark (which is equivalent to the minimum mean-squared error estimate since the noise is additive white Gaussian noise).
[0228] The method may comprise low-pass filtering the data before fitting the one or more functions.
[0229] In one embodiment, however, the method comprises adapting the one or more functions to allow for a low frequency artefact in the detector output data. This may be done, in one example, by expressing the one or more functions as a linear combination of three exponential functions (such as
Figure imgf000064_0004
[0230] In a certain embodiment, the method comprises forcing any estimates having the pulse starting within the window to start at a boundary of the window. [0231] In a particular embodiment, the method comprises maximizing window size or varying window size.
[0232] In one embodiment, the method comprises transforming the detector output data with a transform before fitting the one or more functions to the detector output data as transformed.
[0233] This approach may be desirable in applications in which the analysis is simplified if conducted in transform space. In such situations, the method may also comprise subsequently applying an inverse transform to the one or more functions, though in some cases it may be possible to obtain the desired information in the transform space.
[0234] The transform may be a Laplace transform, a Fourier transform or other transform.
[0235] In one embodiment, estimating the location of the peak comprises minimizing an offset between the start of a window and a start of the pulse.
[0236] In a particular embodiment, the method further comprises detecting a pulse or pulses in the data by:
sliding a window across the data to successive window locations;
identifying possible pulses by performing pulse fitting to the data in the window at each window location;
determining which of the possible pulses have a pulse start falling before and near the start of the respective window location and a peak amplitude exceeding the standard deviation of the noise in the window at the respective window location; and
identifying as pulses, or outputting, those of the possible pulses that have a pulse start falling one, two or three samples before the start of the respective window location and a peak amplitude exceeding the standard deviation of the noise in the window at the respective window location.
[0237] According to a second class of embodiments a method for locating a pulse in detector output data is implemented, comprising:
fitting a plurality of functions to the data;
determining a function of best fit, being whichever of said functions optimises a chosen metric when modelling said data; and
determining a location and an amplitude of a peak of said pulse from said function of best fit.
[0238] In one embodiment, each of the one or more functions is a superposition of a plurality of functions.
1.4. Model-Based, High-Throughput Pulse Processing - Method 4
[0239] The algorithm briefly described here, and in more detail in WO2015085372 (incorporated by reference), for processing the data from radiation detectors is a model - based, real-time, signal-processing algorithm wherein resolving individual signals in the detector output data comprises transforming detector data to produce stepped data, or using data that is already in a stepped form, and detecting at least one signal and estimating a parameter of the signal based at least partially on the stepped data.
[0240]. The method comprises transforming the detector output data to produce stepped data or integral data, detecting at least one event, and estimating a pulse energy associated with the event.
[0241] In some embodiments detecting the at least one event occurs by fitting an expected pulse shape with a sliding window segment of the transformed pulse shape data.
[0242] In some embodiments the method further comprises the step of detecting peaks in the signal, wherein a detection metric is applied to the transformed data. In some embodiments, the detection metric is compared to a simple threshold - if the metric is less than the threshold, then no pulses are deemed present - if it exceeds the threshold, then one or more pulses may be present. Declaration of significant peaks in the detection metric is conducted, when the slope of the peak changing from positive to negative indicates an event.
[0243] It will be appreciated that it may not be possible to adequately characterise all data (uncharacterised data is termed 'corrupt data'); such corrupt data may optionally be rejected. It will be understood that the term 'signal' is interchangeable in this context with 'pulse', as it refers to the output corresponding to individual detection events rather than the overall output signal comprising the sum of individual signals. It will also be appreciated that the temporal position (or timing) of a signal can be measured or expressed in various ways, such as according to the time (or position in the time axis) of the maximum of the signal or the leading edge of the signal. Typically this is described as the arrival time ('time of arrival') or detection time. [0244] The method optionally comprises deleting samples within a set window around the rising edge to ensure the edge region of each pulse, where the real transformed pulse data differs from an ideally transformed pulse, is excluded from the calculations.
[0245] The method optionally comprises an assessment of variance of the energy estimations in the data, and validation of the modeled data.
[0246] The method may include building a model of the data from the processed data output and determining the accuracy of the modeling based on a comparison between the detector output data and the model.
[0247] In one exemplary embodiment of the method 4, the method includes creating a model of the detector output using the signal parameters in combination with the detector impulse response. In another exemplary embodiment, the method may include performing error detection by comparing the actual detector output with the model of the detector output, such as by using least squares.
[0248] The method may include discarding energy estimates deemed not sufficiently accurate. In one embodiment, the method includes presenting all sufficiently accurate energy estimates in a histogram.
2. Pulse Pileup Reduction
[0249] Even where an appropriate high rate pulse processing method is used, there will still be situations where it is not possible to distinguish between closely spaced pulse arrivals. Such a situation occurs when multiple pulses arrive within the window in which the pulse detection algorithm is able to determine the arrival of distinct pulses. Depending on the ADC sampling rate, pulse arrival statistics, and detector electronics, the total amount of pileup may still be in the order of 5% at 1 Mc/s. Pileup can be the result of detecting 2 pulses as a single pulse, however detecting 3 pulses as 1 pulse is also possible, while detecting 4 or more pulses as 1 pulse is possible but much less likely. 2.1. Problem Solution - two pulse pileup removal
[0250] If the underlying X-ray energy spectrum is denoted x then the spectrum with two pulse pileup is:
Figure imgf000068_0002
where * denotes convolution, and is the pileup coefficient that is estimated from data
Figure imgf000068_0007
observation or computed from theory. In order to estimate the underlying spectrum x, the following process is performed:
1. Take the FFT of each side, where convolution now becomes multiplication, so
Figure imgf000068_0003
2. At each FFT bin n, solve the quadratic equation bearing
Figure imgf000068_0006
in mind both X(n) and Y(n) are complex. The solution is
Figure imgf000068_0001
3. The correct solution to take is the "positive" solution. It also relies on taking the correct solution to the complex square root.
4. Now compute the spectrum without pileup by taking the inverse FFT of X.
Figure imgf000068_0004
Using the correct pileup coefficient, it is shown that the pileup is completely removed.
2.2. Problem Solution - two and three pulse pileup removal
[0251] In practice, for spectra measured on an X-ray scanning system, when only two pulse pileup was removed, it was observed that there was still some residual pileup at higher energies. This indicated there was some unremoved three (or more) pulse pileup at these higher energies. In order to remove some, and hopefully most, of this residual pileup, the model is now extended to include 3 pulse pileup, so the received spectrum is given by:
Figure imgf000068_0005
where * denotes convolution, and are the pileup coefficients for two and three
Figure imgf000069_0008
pulse pileup respectively. In order to estimate the underlying spectrum x, the following process is performed:
1. Take the FFT of each side, where convolution now becomes multiplication, so
Figure imgf000069_0002
2. For each of N bins in the FFT, Solve the cubic equation
Figure imgf000069_0003
bearing in mind both X and Y are complex. Like the quadratic, there is a closed form solution, however the solution to the cubic is considerably more complicated as follows:
a. Divide through by k2 so each equation is now of the form
Figure imgf000069_0004
noting that X(n) and c(n) are complex.
b. Compute
Figure imgf000069_0001
c. Compute
Figure imgf000069_0005
d. Check P, and negate if desirable. If Re(conj(R). P) < 0, then P = -P. This is to ensure the correct cube roots are obtained at the next step.
e. Compute
Figure imgf000069_0006
f. Compute
Figure imgf000069_0007
g. Compute the 3 solutions to the cubic equation:
Figure imgf000070_0001
3. Choose the solution to allocate to X(n) . The correct solution is that with smallest magnitude of r2and r3.
Figure imgf000070_0002
4. Now compute the spectrum without pileup by taking the inverse FFT of X
Figure imgf000070_0003
[0252] Using the correct pileup coefficients, it is shown in Figure 1 1 that the pileup is completely removed. If the same data is processed with the quadratic solver, which assumes only two pulse pileup, it can be seen in Figure 12 there is still residual pileup in the spectrum at higher energy values, and slight distortion of the spectrum at lower energy.
3. Design of optimal spectral smoothing window
[0253] Smoothing of the energy spectrum is particularly useful in X-ray screening systems where the duration for spectrum measurement may be very short in order to achieve a high spatial resolution in the sample image. It has been found that typical energy spectra produced by a broad energy X-ray scanning system tend to have almost exclusively low frequency components. Initially to reduce communication bandwidth, but also to reduce computational requirement and provide the added benefit of smoothing the spectrum, the spectrum data are passed through an FFT.
[0254] After FFT, the majority of the FFT bins are discarded, as it is only necessary to keep approximately 1/8 of the FFT bins in order to accurately reconstruct the energy spectrum. For example if there are 512 histogram bins computed, only 32 complex FFT bins are retained. The last 32 complex FFT bins are just the complex conjugate of these bins, and the remaining 448 bins contain (almost) no information.
[0255] The effect of discarding these FFT bins is to:
1. Provi de noi se rej ecti on .
2. Filter the reconstructed spectrum (after iFFT)
However, if a rectangular FFT window is applied, after iFFT the measured spectrum is substantially convolved with a sine function. This is not desirable due to the long extent of the sine function, and large ringing.
[0256] To improve the FFT window function design, the following approach was adopted:
1. Specify the desired "time domain" window. In this example, a raised cosine pulse is used.
2. Take the FFT of the desired window (made symmetric about 0 to give only real FFT output).
3. Multiply this result by the existing rectangular window only.
4. Further multiply the result by a window that has a slight tapering on the edge to further reduce ringing resulting from multiplying by a rectangular window.
[0257] Figure 14 illustrates the result achieved. The rectangular window if used on its own, results in the measured spectrum being convolved with a sine function, with width at the mid amplitude of approx 10 samples, but significant oscillations - around 22% at the first negative going peak. By careful definition of user window w it is possible to achieve a "time" domain response that is approximately raised cosine in nature, with very little oscillatory nature - around 0.2%. However, the width at mid amplitude increases to around 20 samples.
[0258] While the FFT and data truncation were used to reduce communication and computational burden, the additional benefit of an appropriately designed window function is that the received energy spectra are smoothed before processing, resulting in a significant reduction in noise in the effective Z estimates, and the potential for using less bins in the effective Z estimate while achieving a similar result. 4. Pulse Parameter Calibration
[0259] The following is a suitable method for the calibration of the received pulse parameters a and β for pulses of the form:
Figure imgf000072_0001
where a and β are the falling edge and rising edge time constants respectively, t0 is the pulse time of arrival, Ta is the pulse averaging window, and A is a pulse scaling factor related to the pulse energy.
[0260] The pulse parameters may be estimated from a time series capture of the digitised detector signal as follows:
1. Obtain a number of samples of the digitised detector signal, obtained during a period with X-rays on, and overall pulse rate is low enough that isolated pulses can be extracted. Depending on the pulse parameters, with fast pulses it may be adequate to use approx 500k samples at 100 MHz sample rate, and at a count rate of up to 500k pulses per second.
2. Extract a block of samples of length (numpO x sampleRate / nominalCountRate).
For numpO = 40, sample Rate 100 MS/s, nominal count rate 100 kcs, this is 40,000 samples.
3. Calculate the noise threshold nthr
a. Histogram the data block - the histogram bins are integers in the range of the sampled data for 14 bit signed data.
Figure imgf000072_0002
b. Find the bin with the highest value. This is the estimated noise mean.
c. Find the bin where the level falls to 0.63 of the peak. The difference from the peak is the estimated noise std deviation (sigma)
d. Set the noise threshold at 2 sigma from the mean, nthr = noiseMean + 2 x noiseSigma. Factors other than 2 may also be used, depending on the application.
4. Calculate the signal threshold sthr
a. Filter the data block with the "jump" filter of the form [-1 -1 -1 0 1 1 1] b. Set the detection threshold at nthr, and increment in steps of 4 x noise Sigma. c. Threshold the filtered data, and determine the number of runs where the data exceeds sthr. A "run" is a continuous sequence of samples that all exceed sthr, terminated on each end by a sample below sthr.
d. Continue incrementing the detection threshold until at step k nruns(k) - nruns(k-l) >= -1. That is, until the number of runs stops decreasing. (Note: this stopping criteria might produce a pessimistic threshold at higher count rates).
e. Set sthr to be the current threshold at step k.
Estimate the count rate as nruns(k) / (buffer length / sample Rate).
Optional step: If the count rate estimate is less than half or more than double the nominal count rate, redo the noise and signal threshold calculation with a data buffer length computed from the count rate estimate to get closer to numpO detected pulses. Implement the pulse detection state machine. First, detect numpl = 50 pulses to estimate the pulse length lenp (initially set lenp to 0). Then detect nump2 = 600 pulses for full parameter estimation and optimisation. The pulse detection state machined is as follows:
a. Enter at "seekPulse" state
b. When a value exceeds sthr, enter "detPulse" state
c. In "detPulse" state, look for value below sthr. Enter "seekEndPulse" state d. In "seekEndPulse" state
i. If value > sthr, a new detection has occurred before the end of the pulse. Enter "pulsePileUp" state
ii. If value < nthr and pulseLength > lenp, end of a valid pulse is detected - record pulse start/end/length parameters and re-enter "seekPulse" state
e. In "pulsePileUp" state, look for value below sthr, then enter "seekEndPileup" state
f. In "seekEndPileup", change state
i. If value > sthr, a new detection has occurred before the end of a pileup event, indicating further pileup. Return to "pulsePileUp" state. ii. If value < nthr and pulseLength > lenp, the end of the pileup event has been reached. Record the pulse details and mark as pileup so it is not used in calibration. In practice all details about this pulse event could be discarded as it will not be used in calibration. For the first numpl valid (isolated) pulses, do the following:
a. Compute time of arrival (tO), rising edge exponent (beta), falling edge exponent (alpha), averaging time (Ta), maximum signal (Smax), time of maximum (tmax), pulse energy (E).
b. Some pulses may be rejected at this point if it appears there is actually more than one pulse (undetected pileup) - this is indicated by multiple zero crossings in the derivative of the filtered data (zero derivative = local max/min location).
c. Set the pulse length estimate to 7 / median(alpha). This yields an approximate value for the sample at which the pulse will fall to 0.001 of the peak value. A more accurate computation can be obtained using alpha and beta if required, but the 0.001 threshold is somewhat arbitrary in any case, and in the tail the pulse is slowly converging to zero.
Return to step 8 and obtain nump2 pulses. nump2 = 600 has been used, but this is somewhat arbitrary and based on how many pulses were actually in the test data Only half of these pulses will eventually be used in the calibration, so nump2 needs to be double the number of pulses that are required (desired) in the calibration process.
For each of the nump2 pulses:
a. Compute time of arrival (tO), rising edge exponent (beta), falling edge exponent (alpha), averaging time (Ta), maximum signal (Smax), time of maximum (tmax), pulse energy (E). Again, some pulses may be rejected if they appear to be undetected pileup.
b. Sort the pulses into increasing energy sequence.
c. Compute the upper and lower quartile energy values. Discard pulses in the upper and lower quartiles. This effectively removes outlier energy values from the sample, although in a mixture of pulse energies this may not be the best thing to do. In fact it may be better to sort on alpha, beta, or least squares cost function and discard on this basis. For now the Energy sort seems adequate.
d. For the remaining pulses, now only half of nump2 (so about 300 if nump2 = 600)
i. Compute an estimated pulse shape from the parameters alpha, beta, Ta, tO.
ii. Normalise the actual received pulse by its energy.
iii. Compute the cost function = sum of squared errors between the estimated and actual pulses (both are normalised to be nominally unit energy).
iv. Perform an iterative least squares optimisation to obtain optimal estimates of alpha, beta, Ta, along with final cost function and number of iterations for convergence of the least squares optimiser. Note: An approximate Gauss Newton LS optimiser has been implemented. Instead of full 3x3 Jacobian, a series of ID Jacobians on each dimension are computed. These are numerical derivatives, so could be subject to substantial error. This means the trajectory is not always in the optimal direction, with greater risk of divergence if the function is not well behaved. It is not recommended to use the function in this form, but if an efficient LS optimiser is not available a more robust implementation could be provided.
11. From the least squares optimised results, set the final values of alpha, beta, Ta. This can be either median or mean values of the optimised parameters. The value of tO can be set arbitrarily such that either a) tO = 0 (and the pulse therefore has some signal for sample k < 0) or b) tO = ceil(Ta), in which case the pulse is zero at k=0 and has positive value from k>=l .
12. The pulse form p(t) can be computed directly from the formula and the estimated Ta, a and β.
5. Method for baseline tracking
[0261] In order to correctly determine the energy of the pulses, it is desirable to account for DC offset (or signal baseline, used interchangeably) in the signal from the detector. This DC offset can arise from various sources including the bias levels of analog electronics, the analog to digital conversion, and the detector itself. Control theory suggests the DC offset error may be tracked and reduced to zero by generating a feedback signal that is proportional to the integral of the signal - however there is a significant problem in the case of pulse processing. Pulses introduce additional features to the signal that have non-zero mean. This introduces a bias dependent on pulse energy, count rate and pulse shape, which corrupts the feedback signal and prevents standard control loop tracking from successfully removing the DC offset.
[0262] To overcome this problem, the detector signal output is digitally processed to remove the pulse shaping effects introduced by analog electronics. When no DC offset is present, this processed signal results in a signal shape that has constant value in the regions between pulse arrivals, and a rapid change in value where pulses arrive. If a residual DC offset is present in the detector signal the processed signal changes linearly with time in the regions between pulse arrivals. An error feedback signal that is proportional to the slope of this signal may be formed by taking the difference between two samples. These samples need not be consecutive, but may be separated by 'Ν' samples in time. By choosing an appropriate value for TNT, a signal with a suitable signal to noise ratio may be found for driving a feedback loop.
[0263] In order to reduce the impact of bias introduced by pulse events, the baseline tracking loop is not updated when a pulse has arrived between the two samples used to generate the feedback error signal.
[0264] The impact of bias may be further reduced by preventing the baseline tracking loop from updating when a pulse has arrived within a guard region on either side of the samples used to generate the feedback error signal.
[0265] It should be noted that due to the biasing caused by pulse arrival, the value of the processed detector signal increases whenever a pulse arrives. This eventually causes the internal registers used to store the value of the signal to overflow. The value of the processed signal is monitored, and when overflow is detected, the baseline tracking loop is prevented from updating until the effects of overflow have passed. [0266] Denoting the processed pulse signal at sample n as x(n), the following steps summarise the procedure for computing the update to the DC offset estimate, denoted DC(n) :
1. Compute the difference between signal samples N samples apart
2. Determine if the update is to be applied. Do not apply the DC update if
a. A pulse arrival is detected at a sample between samples n and n + N. b. Transient from a previous detected pulse has not decayed. Transient can be considered to last M samples after a pulse is detected.
c. The processed signal x(n) is about to reach positive overflow and wrap around to a large negative value. Do not process if x(n) is within a threshold Δ of positive or negative overflow.
3. If the DC update is to be applied, then compute the DC offset update as
Figure imgf000077_0001
where k « 1 is the update gain, and is chosen to achieve the desired balance between fast response and noise on the DC estimate.
[0267] Finally, the same hardware can be used for tracking multiple baseline offsets in multiple channels in a time division multiplexed scheme. The values for the tracking loop variables for each channel are stored/loaded when switching between channels. The baseline tracking loop is prevented from updating until transient effects of the detector channel change has passed.
6. Collimation
[0268] Very tight collimation may be used within the scanner in order to minimise the effect of scatter on the measured spectrum. This is particularly important where transitions from high to low or low to high intensity occur. The overall results of the system have shown that scatter has been largely addressed through the inclusion of tight collimation.
7. Reference Calculation
[0269] The purpose of the reference calculation is to establish the mean intensity for each detector. This value is used to scale all intensity histograms to unit energy. This is commonly referred to as normalisation. A reference intensity is calculated for each detector. The reference intensity is calculated as the mean intensity over the first N slices in a scan. The intensity is the 1st bin in the FFT or the sum of all complex -valued elements in the FFT vector.
[0270] There is also a reference histogram computed in the same way - by averaging the measured energy histograms for the first N slices. The reference histogram is used to normalise all measured histograms to ensure any run to run variations in X-ray flux do not impact the effective Z computation.
[0271] The reference is measured during an interval where:
1. X-rays are stabilised, so X-ray flux is not varying and will not vary for the
duration of the scan (in practice the Smiths source does vary - especially when failing - and this can impact results)
2. Before the arrival of the sample under test.
[0272] The current implementation uses a duration measured in slices. This can cause problems when the slice rate is slowed below 5 ms for example - the reference collection can run into the sample under test. This needs to be corrected in 2 ways to be fully robust:
1. Use the configured slice rate and scanning speed to compute a duration for which the reference is collected, not a set number of slices.
2. Incorporate the object detection signal to ensure reference collection is stopped immediately if a sample is detected before the reference duration completes - the user should be warned when this occurs as performance would not be guaranteed.
[0273] More accurate and consistent effective Z may be obtained if a longer reference collection duration was used.
CONCLUSION
[0274] Persons skilled in the art will also appreciate that many variations may be made to the invention without departing from the scope of the invention, which is determined from the broadest scope and claims. [0275] For example, while the invention is described in some embodiments with reference to a material characterisation system which forms an image as the material passes through the detector, the formation of an image is not an essential element of the broadest aspect. The more detectors which are provided, the more accurate will be the material characterisation. The number of detectors needs to be determined on the basis of processing speed, cost and accuracy requirements of the application.
[0276] Further, while the invention is described in the embodiment above in relation to or associated with conveyors, the system and method applies in the broadest aspects to systems without material handling, and to material handling systems which are material translation systems adapted to translate or transport material. Such material translation systems include for example a conveyor, a bench, a moving floor, a bucket elevator, a pipe or a falling stream of the material.
[0277] It is further anticipated that systems according to the invention can comprise a means for sorting, diverting, removing or selecting at least a part of the material according to the characterisation of at least part of the material, such as may be applicable in a recycling plant.
[0278] In the claims which follow and in the preceding description of the invention, except where the context requires otherwise due to express language or necessary implication, the word "comprise" or variations such as "comprises" or "comprising" is used in an inclusive sense, i.e. to specify the presence of the stated features but not to preclude the presence or addition of further features in various embodiments of the invention. Further, any method steps recited in the claims are not necessarily intended to be performed temporally in the sequence written, or to be performed without pause once started, unless the context requires it.

Claims

CLAIMS:
1. A system for characterising at least part of a material comprising:
a source of incident X-rays configured to irradiate at least part of the material;
one or more detectors adapted to detect radiation emanating from within or passing through the material as a result of the irradiation by the incident radiation and thereby produce a detection signal; and
one or more digital processors configured to process the detection signal to characterise at least part of the material;
wherein the one or more detectors and one or more digital processors are configured to characterise at least part of the material by performing energy resolved photon counting X-ray transmission spectroscopy analysis.
2. The system of claim 1, wherein the material, the source of incident X-rays and the one or more detectors are positioned to irradiate at least part of the material at an irradiation location while the material is travelling on, in, through or along a material handling system.
3. The system of claim 2, wherein the material handling system is a material translation system adapted to translate or transport the material.
4. The system of claim 3, wherein the material translation system comprises at the irradiation location a conveyor, a bench, a moving floor, a bucket elevator, or a pipe.
5. The system of claim 3 or claim 4, wherein the source of incident X-rays is positioned on one side of the material translation system and the one or more detectors is positioned on one or more opposite sides of the material translation system.
6. The system of claim 2, wherein the source of incident X-rays is positioned above the material handling system and the one or more detectors are positioned beneath the material translation system.
7. The system of claim 2, wherein the source of incident X-rays is positioned beneath the material handling system and the one or more detectors are positioned above the material handling system.
8. The system of claim 1, wherein the source of incident X-rays and the one or more detectors are positioned to irradiate the material while the material is falling.
9. The system of claim 8, wherein at least the source of incident X-rays and the one or more detectors are disposed within a housing having an aperture through which the material passes and the material is irradiated as the material is falling.
10. The system of claim 2, wherein the source of incident X-rays and the one or more detectors are positioned on a first chassis, the first chassis engaging with a second chassis through vibration damping, the second chassis being attached to the material handling system.
11. The system of claim 10, wherein the source of incident X-rays and the one or more detectors are disposed within one or more modular enclosures adapted for removal and replacement from the first chassis.
12. The system of any one of claims 1 to 11, wherein:
the one or more detectors is or are adapted to detect packets of emitted radiation emanating from within or passing through the material as a result of the irradiation by the incident radiation, each said detector being configured to produce an electrical pulse caused by the detected packets having a characteristic size or shape dependent on an energy of the packets; and
the one or more digital processors is or are configured to process each electrical pulse to determine the characteristic size or shape and to thereby accumulate a detector energy spectrum for each detector of the energies of the packets detected, and to characterise at least part of the material based on the detector energy spectra.
13. The system of claim 12, wherein each packet of emitted radiation is a photon and the one or more detectors comprise one or more scintillation detectors each composed of a scintillation material adapted to produce electromagnetic radiation by scintillation from the photons and a pulse producing element adapted to produce the electrical pulse from the electromagnetic radiation.
14. The system of claim 13, wherein the pulse producing element comprises a photon- sensitive material and the one or more scintillation detectors are a plurality of scintillation detectors arranged side-by-side in one or more detector arrays of individual scintillator elements of the scintillation material each covered with reflective material around sides thereof and disposed adjacent and optically coupled to a photon-sensitive material.
15. The system of claim 14, wherein the scintillation material is lutetium -yttrium oxyorthosilicate (LYSO).
16. The system of any one of claims 13 to 15, wherein the photon-sensitive material is a silicon photomultiplier (SiPM).
17. The system of any one of claims 14 to 16, wherein the individual scintillator elements of one or more of the detector arrays present a cross-sectional area to the emitted radiation of greater than 1 square millimetre.
18. The system of claim 17, wherein the cross-sectional area is greater than 2 square millimetres and less than 5 square millimetres.
19. The system of any one of claims 1 to 18, wherein the one or more digital processors are further configured with a pileup recovery algorithm adapted to determine the energy associated with two or more overlapping pulses.
20. The system of any one of claims 1 to 19, wherein the one or more digital processors is configured to compute an effective atomic number Z for each of at least some of the detectors based at least in part on the corresponding detector energy spectrum.
21. The system of claim 20, wherein the one or more digital processors is configured to compute the effective atomic number Z for each of at least some of the detectors by: determining a predicted energy spectrum for a material with effective atomic number Z having regard to an estimated material thickness deduced from the detector energy spectrum and reference mass attenuation data for effective atomic number Z; and
comparing the predicted energy spectrum with the detector energy spectrum.
22. The system of claim 20, wherein the one or more digital processors is configured to compute the effective atomic number Z for each of at least some of the detectors by: determining a predicted energy spectrum for a material with effective atomic number Z having regard to a calibration table formed by measuring one or more materials of known composition; and
comparing the predicted energy spectrum with the detector energy spectrum.
23. The system of any one of claims 21 to 22, wherein the one or more digital processors is configured to perform the step of comparing by computing a cost function dependent on a difference between the detector energy spectrum and the predicted energy spectrum for a material with effective atomic number Z.
24. The system of any one of claims 1 to 23, wherein a gain calibration is performed on each detector individually to provide consistency of energy determination among the detectors and the one or more digital processors is further configured to calculate the detector energy spectrum for each detector taking into account the gain calibration.
25. The system of any one of claims 1 to 24 wherein a count rate dependent calibration is performed comprising adaptation of the detector energy spectra for a count rate dependent shift.
26. The system of any one of claims 1 to 25 wherein a system parameter dependent calibration is performed on the detector energy spectra comprising adaptation for time, temperature or other system parameters.
27. The system of any one of claims 1 to 26, wherein the one or more digital processors is further configured to reduce a communication bandwidth or memory use associated with processing or storage of the detector energy spectra, by performing a fast Fourier transform of the energy spectra and removing bins of the fast Fourier transform having little or no signal to produce reduced transformed detector energy spectra.
28. The system of claim 27, wherein the one or more digital processors is further configured to apply an inverse fast Fourier transform on the reduced transformed detector energy spectra to provide reconstructed detector energy spectra.
29. The system of claim 27 or 28, wherein the one or more digital processors is further configured with a specific fast Fourier transform window optimised to minimise ringing effects of the fast Fourier transform.
30. The system of any one of claims 1 to 29, wherein the one or more digital processors is further configured with a baseline offset removal algorithm to remove a baseline of a digital signal of electrical pulse prior to further processing.
31. The system of any one of claims 1 to 30, wherein the one or more digital processors is further configured to produce an image of the material composed of pixels representing the characterisation of different parts of the material deduced from the detector energy spectra.
32. The system of any one of claims 1 to 31, wherein the one or more digital processors is further configured to perform one or more of tiling, clustering, edge detection or moving average based on the effective atomic numbers determined for said plurality of detectors.
33. The system of any one of claims 1 to 32, further comprising visual or sound outputs of system status or characterisation results.
34. The system of any one of claims 1 to 33, further comprising communication connections to one or more of a screen, a tablet or other computer operated by an operator nearby, at an on-site control centre, or at a remote off-site control centre.
35. The system of any one of claims 1 to 34, wherein at least one of the source of incident X-rays, the one or more detectors and the one or more digital processors is modular, being adapted to be removed and replaced at a site of installation of the system.
36. The system of any one of claims 1 to 35, wherein the material comprises material selected from the group consisting of: discrete items of material, a particulate substance, a non-particulate substance, ore, coal, cement, mineral sand, fluid, hydrocarbon fluid, soil, agricultural material, food or beverage, metals, plastics, glass, timber pharmaceuticals, contaminated material, hazardous material, waste, refuse, sewage, effluent, recyclable material, construction debris, natural materials, manufactured products, clothing, electronics, medical devices, automotive components, aerospace components, dentistry devices, and dentistry components.
37. The system of any one of claims 1-36, further comprising a means for sorting, diverting, removing or selecting at least a part of the material according to the characterisation of at least part of the material.
38. A method of characterising at least part of a material, the method comprising the steps of:
irradiating at least part of the material with incident X-rays from an X-ray source; detecting radiation emanating from within or passing through the material as a result of the irradiation by the incident X-rays; producing a detection signal;
processing the detection signal to characterise at least part of the material; and characterising at least part of the material by performing energy resolved photon counting X-ray transmission spectroscopy analysis.
39. The method of claim 38, comprising irradiating at least part of the material while the material is travelling on, in, through or along a material handling system.
40. The method of claim 39, wherein the material handling system is a material translation system adapted to translate or transport the material.
41. The method of claim 40, wherein the material translation system comprises at the irradiation location a conveyor, a bench, a moving floor, a bucket elevator, a pipe or a falling stream of the material.
42. The method of claim 40 or claim 41, comprising positioning the source of incident X- rays on one side of the material translation system and positioning the one or more on one or more opposite sides of the material translation system.
43. The method of claim 39, comprising positioning the source of incident X-rays above the material handling system and positioning the one or more detectors beneath the material handling system.
44. The method of claim 39, comprising positioning the source of incident X-rays beneath the material handling system and positioning the one or more detectors above the material handling system.
45. The method of claim 38, comprising irradiating at least part of the material while the material is falling.
46. The method of claim 45, comprising irradiating at least part of the material while the material is falling through an aperture in a housing in which the source of incident X-rays and the one or more detectors are disposed.
47. The system of claim 37, wherein the source of incident X-rays and the one or more detectors are positioned on a first chassis, the first chassis engaging with a second chassis through vibration damping, the second chassis being attached to the material handling system.
48. The system of claim 47, wherein the source of incident X-rays and the one or more detectors are disposed within one or more modular enclosures adapted for removal and replacement from the first chassis.
49. The method of claim 38, comprising:
detecting packets of emitted radiation emanating from within or passing through the material as a result of the irradiation by the incident radiation,
producing an electrical pulse having a characteristic size or shape dependent on an energy of the packets;
processing each electrical pulse to determine the characteristic size or shape and to thereby accumulate a detector energy spectrum for each detector of the energies of the packets detected; and,
characterising at least part of the material based on the one or more detector energy spectra.
50. The method of claim 49, wherein each packet of emitted radiation is a photon and the method comprises detecting said photons using one or more scintillation detectors wherein each scintillation detector comprises:
a scintillation material adapted to produce electromagnetic radiation by scintillation from the photons; and,
a pulse producing element adapted to produce the electrical pulse from the
electromagnetic radiation.
51. The method of claim 50, comprising: producing a pulse using a photon-sensitive material; and,
arranging a plurality of scintillation detectors side-by-side in one or more detector arrays of individual scintillator elements of the scintillation material each covered with reflective material around sides thereof and disposed adjacent and optically coupled to a photon-sensitive material.
52. The method of claim 50, wherein the scintillation material is lutetium-yttrium oxyorthosilicate (LYSO).
53. The method of claim 51, wherein the photon-sensitive material is a silicon
photomultiplier (SiPM).
54. The method of any one of claims 50 to 53, comprising presenting a cross-sectional area of scintillation material to the emitted radiation of greater than 1 square millimetre.
55. The method of claim 54, wherein the cross-sectional area is greater than 2 square millimetres and less than 5 square millimetres.
56. The method of any one of claims 38 to 55, comprising determining the energy associated with two or more overlapping pulses using a pileup recovery algorithm.
57. The method of any one of claims 38 to 56, comprising computing an effective atomic number Z for each of at least some of the detectors based at least in part on the corresponding detector energy spectrum.
58. The method of claim 57, comprising computing the effective atomic number Z for each of at least some of the detectors by:
determining a predicted energy spectrum for a material with effective atomic number Z having regard to an estimated material thickness deduced from the detector energy spectrum and reference mass attenuation data for effective atomic number Z; and
comparing the predicted energy spectrum with the detector energy spectrum.
59. The method of claim 57, comprising computing the effective atomic number Z for each of at least some of the detectors by: determining a predicted energy spectrum for a material with effective atomic number Z having regard to a calibration table formed by measuring one or more materials of known composition; and
comparing the predicted energy spectrum with the detector energy spectrum.
60. The method of any one of claims 58 and 59, comprising performing the step of comparing by computing a cost function dependent on a difference between the detector energy spectrum and the predicted energy spectrum for a material with effective atomic number Z.
61. The method of any one of claims 38 to 60, comprising performing gain calibration on each detector individually to provide consistency of energy determination among the detectors and calculating the detector energy spectrum for each detector taking into account the gain calibration.
62. The method of any one of claims 38 to 61, comprising performing a count rate dependent calibration comprising adaptation of the detector energy spectra for a count rate dependent shift.
63. The method of any one of claims 38 to 61, comprising performing a system parameter dependent calibration on the detector energy spectra comprising adaptation for time, temperature or other system parameters.
64. The method of any one of claims 38 to 61, comprising reducing a communication bandwidth or memory use associated with processing or storage of the detector energy spectra, by performing a fast Fourier transform of the energy spectra and removing bins of the fast Fourier transform having little or no signal to produce reduced transformed detector energy spectra.
65. The method of claim 64, comprising applying an inverse fast Fourier transform on the reduced transformed detector energy spectra to provide reconstructed detector energy spectra.
66. The method of any one of claims 64 to 65, comprising minimising ringing effects of the fast Fourier transform by applying a fast Fourier transform window.
67. The method of any one of claims 38 to 66, comprising applying a baseline offset removal algorithm to remove a baseline of a digital signal of electrical pulse prior to further processing.
68. The method of any one of claims 38 to 67, comprising producing an image of the material composed of pixels representing the characterisation of different parts of the material deduced from the detector energy spectra.
69. The method of any one of claims 38 to 68, comprising performing one or more of tiling, clustering, edge detection or moving average based on the effective atomic numbers determined for said plurality of detectors.
70. The method of any one of claims 38 to 69, comprising generating visual or sound outputs of system status or characterisation results.
71. The method of any one of claims 38 to 70, comprising communicating with one or more of a screen, a tablet or other computer operated by an operator nearby, at an on-site control centre, or at a remote off-site control centre.
72. The method of any one of claims 38 to 71, wherein the material comprises material selected from the group consisting of: discrete items of material, a particulate substance, a non-particulate substance, ore, coal, cement, mineral sand, fluid, hydrocarbon fluid, soil, agricultural material, food or beverage, metals, plastics, glass, timber pharmaceuticals, contaminated material, hazardous material, waste, refuse, sewage, effluent, recyclable material, construction debris, natural materials, manufactured products, clothing, electronics, medical devices, automotive components, aerospace components, dentistry devices, and dentistry components.
73. The method of any one of claims 38 to 72, comprising sorting, diverting, removing or selecting at least a part of the material according to the characterisation of at least part of the material.
PCT/AU2017/050514 2016-05-30 2017-05-30 Material characterisation system and method WO2017205914A1 (en)

Priority Applications (6)

Application Number Priority Date Filing Date Title
AU2017274079A AU2017274079B2 (en) 2016-05-30 2017-05-30 Material characterisation system and method
BR112018074796-3A BR112018074796B1 (en) 2016-05-30 2017-05-30 MATERIAL CHARACTERIZATION SYSTEM AND METHOD
US16/305,803 US11243327B2 (en) 2016-05-30 2017-05-30 System and method for material characterization
CA3025821A CA3025821A1 (en) 2016-05-30 2017-05-30 Material characterisation system and method
EP17805394.8A EP3465180A4 (en) 2016-05-30 2017-05-30 Material characterisation system and method
ZA2018/08640A ZA201808640B (en) 2016-05-30 2018-12-20 Material characterisation system and method

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
AU2016902059A AU2016902059A0 (en) 2016-05-30 Material characterisation system and method
AU2016902059 2016-05-30

Publications (1)

Publication Number Publication Date
WO2017205914A1 true WO2017205914A1 (en) 2017-12-07

Family

ID=60478498

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/AU2017/050514 WO2017205914A1 (en) 2016-05-30 2017-05-30 Material characterisation system and method

Country Status (7)

Country Link
US (1) US11243327B2 (en)
EP (1) EP3465180A4 (en)
AU (1) AU2017274079B2 (en)
BR (1) BR112018074796B1 (en)
CA (1) CA3025821A1 (en)
WO (1) WO2017205914A1 (en)
ZA (1) ZA201808640B (en)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2020070473A1 (en) * 2018-10-01 2020-04-09 Smiths Heimann Sas Correction of images
EP3734260A3 (en) * 2019-05-03 2020-12-02 WIPOTEC GmbH Method and apparatus for x-ray inspection of products, in particular food products
WO2024108257A1 (en) * 2022-11-24 2024-05-30 Southern Innovation International Pty Ltd Density measuring

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP3465177A4 (en) * 2016-05-30 2020-03-04 Southern Innovation International Pty Ltd Material characterisation system and method
US20210354911A1 (en) * 2016-07-15 2021-11-18 Cleanrobotics Technologies, Inc. Waste receptacle with sensing and interactive presentation system
US11268857B2 (en) * 2018-05-25 2022-03-08 Toyo Corporation Spectrum analysis method and spectrum analysis apparatus
FR3087035B1 (en) * 2018-10-09 2020-10-30 Commissariat Energie Atomique PROCESS FOR CORRECTING A SPECTRAL IMAGE
BR112021007205A2 (en) * 2018-10-18 2021-08-10 Security Matters Ltd. system for detecting and identifying at least one predetermined foreign element in a substance, and method for inspecting a substance.
CN111814711B (en) * 2020-07-15 2023-08-08 中国矿业大学 Image feature quick matching method and system applied to mine machine vision
KR20240110614A (en) * 2021-11-16 2024-07-15 루마필드, 인크. X-ray device with shroud

Citations (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2006029475A1 (en) 2004-09-16 2006-03-23 Southern Innovation International Pty Ltd Method and apparatus for resolving individual signals in detector output data.
WO2009121132A1 (en) 2008-03-31 2009-10-08 Southern Innovation International Pty Ltd Radiation imaging method with individual signal resolution
WO2009121130A1 (en) 2008-03-31 2009-10-08 Southern Innovation International Pty Ltd Screening method and apparatus
WO2009121131A1 (en) 2008-03-31 2009-10-08 Southern Innovation International Pty Ltd Method and apparatus for borehole logging
WO2010068996A1 (en) 2008-12-18 2010-06-24 Southern Innovation International Pty Ltd Method and apparatus for resolving piled-up pulses by using a mathematical transform
WO2012171059A1 (en) 2011-06-14 2012-12-20 Southern Innovation International Pty Ltd Method and apparatus for identifying pulses in detector output data
EP2405260B1 (en) * 2010-07-09 2014-05-07 Xnext S.r.l. Method and apparatus for performing non-invasive x-ray inspections of objects
WO2014091444A1 (en) * 2012-12-14 2014-06-19 Koninklijke Philips N.V. Detector unit with pulse shaper
US8855809B2 (en) 2011-09-01 2014-10-07 Spectramet, Llc Material sorting technology
WO2015062618A1 (en) * 2013-10-28 2015-05-07 Orexplore Ab Apparatus and method for x-ray transmission analysis of a mineral or electronic waste sample
WO2015085372A1 (en) 2013-12-11 2015-06-18 Southern Innovation International Pty Ltd Method and apparatus for resolving signals in data
US20150338545A1 (en) * 2014-05-23 2015-11-26 Radiabeam Technologies, Llc System and method for adaptive x-ray cargo inspection
WO2016082006A1 (en) * 2014-11-30 2016-06-02 Southern Innovation International Pty Ltd Device and method for material characterisation

Family Cites Families (68)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5600700A (en) 1995-09-25 1997-02-04 Vivid Technologies, Inc. Detecting explosives or other contraband by employing transmitted and scattered X-rays
US5943388A (en) * 1996-07-30 1999-08-24 Nova R & D, Inc. Radiation detector and non-destructive inspection
US6442233B1 (en) * 1998-06-18 2002-08-27 American Science And Engineering, Inc. Coherent x-ray scatter inspection system with sidescatter and energy-resolved detection
US7248716B2 (en) * 2001-07-06 2007-07-24 Palantyr Research, Llc Imaging system, methodology, and applications employing reciprocal space optical design
US7072440B2 (en) * 2001-10-19 2006-07-04 Control Screening, Llc Tomographic scanning X-ray inspection system using transmitted and Compton scattered radiation
ATE331947T1 (en) * 2001-10-22 2006-07-15 Miconos Technology Ltd X-RAY INSPECTION SYSTEM FOR PRODUCTS, ESPECIALLY FOOD
US7065175B2 (en) * 2003-03-03 2006-06-20 Varian Medical Systems Technologies, Inc. X-ray diffraction-based scanning system
WO2004095060A2 (en) * 2003-04-23 2004-11-04 L-3 Communications Security and Detection Systems Corporation X-ray imaging technique
US7092485B2 (en) * 2003-05-27 2006-08-15 Control Screening, Llc X-ray inspection system for detecting explosives and other contraband
US6987833B2 (en) * 2003-10-16 2006-01-17 General Electric Company Methods and apparatus for identification and imaging of specific materials
US7099433B2 (en) * 2004-03-01 2006-08-29 Spectramet, Llc Method and apparatus for sorting materials according to relative composition
US7564943B2 (en) * 2004-03-01 2009-07-21 Spectramet, Llc Method and apparatus for sorting materials according to relative composition
US7453987B1 (en) * 2004-03-04 2008-11-18 Science Applications International Corporation Method and system for high energy, low radiation power X-ray imaging of the contents of a target
US7634061B1 (en) * 2004-03-26 2009-12-15 Nova R & D, Inc. High resolution imaging system
EA009923B1 (en) * 2004-03-27 2008-04-28 Тексмаг Гмбх Apparatus for detecting joints in rubber sheets
US8233667B2 (en) * 2004-09-07 2012-07-31 Petro-Model Ehf Apparatus and method for analysis of size, form and angularity and for compositional analysis of mineral and rock particles
WO2007109710A2 (en) * 2006-03-21 2007-09-27 Board Of Regents, The University Of Texas System Optical device for detecting live insect infestation
US7769132B1 (en) 2007-03-13 2010-08-03 L-3 Communications Security And Detection Systems, Inc. Material analysis based on imaging effective atomic numbers
US7742568B2 (en) * 2007-06-09 2010-06-22 Spectrum San Diego, Inc. Automobile scanning system
US7869566B2 (en) * 2007-06-29 2011-01-11 Morpho Detection, Inc. Integrated multi-sensor systems for and methods of explosives detection
CN101435783B (en) * 2007-11-15 2011-01-26 同方威视技术股份有限公司 Method and apparatus for recognizing substance
US7613274B2 (en) * 2007-11-16 2009-11-03 General Electric Company Method and system of energy integrating and photon counting using layered photon counting detector
US20090168958A1 (en) * 2008-01-02 2009-07-02 Cristina Francesca Cozzini Apparatus and method for identifying components in a container
EP2296550B1 (en) * 2008-07-07 2018-03-28 Koninklijke Philips N.V. K-edge imaging
US7965816B2 (en) 2008-08-11 2011-06-21 Control Screening, LLC. Scanning X-ray inspection system using scintillation detection with simultaneous counting and integrating modes
FR2953603A1 (en) * 2009-12-09 2011-06-10 Commissariat Energie Atomique METHOD AND DEVICE FOR RECOGNIZING A MATERIAL USING ITS TRANSMISSION FUNCTION
US8750454B2 (en) 2010-02-25 2014-06-10 Rapiscan Systems, Inc. High-energy X-ray-spectroscopy-based inspection system and methods to determine the atomic number of materials
CN102313752B (en) * 2010-06-30 2014-07-23 清华大学 Article detection equipment and method
WO2012009725A1 (en) * 2010-07-16 2012-01-19 Mayo Foundation For Medical Education And Research System and method for improved energy series of images using multi-energy ct
US8422636B2 (en) * 2010-10-12 2013-04-16 Ge Medical Systems Israel, Ltd. Photon counting and energy discriminating detector threshold calibration
WO2012054381A1 (en) * 2010-10-18 2012-04-26 American Science And Engineering, Inc. System and methods for intrapulse multi-energy and adaptive multi-energy x-ray cargo inspection
US9316537B2 (en) * 2011-06-29 2016-04-19 Minesense Technologies Ltd. Sorting materials using a pattern recognition, such as upgrading nickel laterite ores through electromagnetic sensor-based methods
DK2726711T3 (en) * 2011-06-29 2020-07-27 Minesense Tech Ltd Extraction of extracted ore, minerals or other materials using sensor-based sorting
GB201113224D0 (en) * 2011-08-01 2011-09-14 Kromek Ltd Method for the radiological investigation of an object
US9448326B2 (en) * 2011-08-01 2016-09-20 Kromek Limited Detection and/or classification of materials
WO2013017878A1 (en) * 2011-08-01 2013-02-07 Kromek Limited Object monitoring using multi spectral radiation
BR112014002662A8 (en) 2011-08-04 2017-06-20 Tech Resources Pty Ltd mining material processing
JP2014535039A (en) * 2011-09-30 2014-12-25 アナロジック コーポレイション Photon number correction method and apparatus {PHTONCOUNTCOUNTRECTION}
US8929508B1 (en) * 2012-04-17 2015-01-06 Robert E. Alvarez Energy selective X-ray system with low noise
JP5863547B2 (en) * 2012-04-20 2016-02-16 ヤマハ発動機株式会社 Printed circuit board inspection equipment
EP3369488B1 (en) * 2012-05-01 2021-06-23 Minesense Technologies Ltd. High capacity cascade-type mineral sorting method
US8907290B2 (en) 2012-06-08 2014-12-09 General Electric Company Methods and systems for gain calibration of gamma ray detectors
US20140077007A1 (en) * 2012-09-20 2014-03-20 National Recovery Technologies, Llc Methods of Processing Waste Material to Render a Compostable Product
EP2914179A1 (en) * 2012-11-02 2015-09-09 Analogic Corporation Volumetric and projection image generation
JP6073675B2 (en) * 2012-12-28 2017-02-01 東芝メディカルシステムズ株式会社 X-ray CT apparatus and control program
CN103913472B (en) * 2012-12-31 2016-04-20 同方威视技术股份有限公司 Ct imaging system and method
US8958524B2 (en) * 2013-01-31 2015-02-17 Analogic Corporation Correction of projection data in radiation system
WO2015058702A1 (en) * 2013-10-23 2015-04-30 曹红光 Photon count-based radiation imaging system, method, and apparatus
WO2015070008A1 (en) * 2013-11-08 2015-05-14 Schlumberger Canada Limited Spectral analysis with spectrum deconvolution
EP2871496B1 (en) * 2013-11-12 2020-01-01 Samsung Electronics Co., Ltd Radiation detector and computed tomography apparatus using the same
US9488739B2 (en) * 2013-12-18 2016-11-08 The Board Of Trustees Of The Leland Stanford Junior University Spectral imaging system and method
JP6289108B2 (en) * 2014-01-14 2018-03-07 キヤノンメディカルシステムズ株式会社 Photon counting CT apparatus and photon counting CT data processing method
FR3023001A1 (en) * 2014-06-30 2016-01-01 Commissariat Energie Atomique METHOD FOR ANALYZING A TWO-STAGE OBJECT USING TRANSMITTED RADIATION THEN A DIFFUSION SPECTRUM
CN110090812B (en) * 2014-07-21 2021-07-09 感矿科技有限公司 High capacity separation of coarse ore minerals from waste minerals
JP6529858B2 (en) * 2014-08-20 2019-06-12 キヤノンメディカルシステムズ株式会社 X-ray CT apparatus and X-ray detector
US10451568B2 (en) * 2014-08-22 2019-10-22 Canon Medical Systems Corporation Photon counting X-ray CT apparatus
JP6313168B2 (en) * 2014-09-02 2018-04-18 キヤノンメディカルシステムズ株式会社 X-ray CT apparatus, image processing apparatus, and image processing program
US10539687B2 (en) * 2014-10-16 2020-01-21 Analogic Corporation Indirect conversion detector array
WO2016178682A1 (en) * 2015-05-07 2016-11-10 Analogic Corporation Combined image generation of article under examination and image of test item
US10215879B2 (en) * 2015-05-07 2019-02-26 Morpho Detection, Llc System for detecting counterfeit goods and method of operating the same
FR3037401B1 (en) * 2015-06-15 2017-06-23 Commissariat Energie Atomique CHARACTERIZATION OF A SAMPLE BY DECOMPOSITION BASED ON MATERIALS.
FR3038722B1 (en) * 2015-07-08 2019-06-14 Commissariat A L'energie Atomique Et Aux Energies Alternatives METHOD FOR IDENTIFYING A MATERIAL
US9916669B2 (en) * 2015-08-24 2018-03-13 Analogic Corporation Projection data correction and computed tomography value computation
CN107533019B (en) * 2015-10-23 2020-05-05 株式会社蛟簿 X-ray apparatus, data processing apparatus, and data processing method
US10338012B2 (en) * 2016-03-09 2019-07-02 Toshiba Medical Systems Corporation Photon counting detector and X-ray computed tomography (CT) apparatus
EP3223002B1 (en) * 2016-03-23 2021-04-28 Carl Zeiss X-Ray Microscopy, Inc. Method, tool and computer program product for measurement and estimation of an x-ray source spectrum in computed tomography, and x-ray ct system
CN107356615B (en) * 2016-05-10 2020-01-21 清华大学 Method and system for dual-energy X-ray CT
EP3465177A4 (en) * 2016-05-30 2020-03-04 Southern Innovation International Pty Ltd Material characterisation system and method

Patent Citations (16)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2006029475A1 (en) 2004-09-16 2006-03-23 Southern Innovation International Pty Ltd Method and apparatus for resolving individual signals in detector output data.
US7383142B2 (en) 2004-09-16 2008-06-03 Southern Innovation International Pty Ltd. Method and apparatus for resolving individual signals in detector output data
WO2009121132A1 (en) 2008-03-31 2009-10-08 Southern Innovation International Pty Ltd Radiation imaging method with individual signal resolution
WO2009121130A1 (en) 2008-03-31 2009-10-08 Southern Innovation International Pty Ltd Screening method and apparatus
WO2009121131A1 (en) 2008-03-31 2009-10-08 Southern Innovation International Pty Ltd Method and apparatus for borehole logging
WO2010068996A1 (en) 2008-12-18 2010-06-24 Southern Innovation International Pty Ltd Method and apparatus for resolving piled-up pulses by using a mathematical transform
US8812268B2 (en) 2008-12-18 2014-08-19 Southern Innovation International Pty. Ltd. Method and apparatus for resolving piled-up pulses by using a mathematical transform
EP2405260B1 (en) * 2010-07-09 2014-05-07 Xnext S.r.l. Method and apparatus for performing non-invasive x-ray inspections of objects
US20140107975A1 (en) * 2011-06-14 2014-04-17 Southern Innovation International Pty Ltd Method and apparatus for identifying pulses in detector output data
WO2012171059A1 (en) 2011-06-14 2012-12-20 Southern Innovation International Pty Ltd Method and apparatus for identifying pulses in detector output data
US8855809B2 (en) 2011-09-01 2014-10-07 Spectramet, Llc Material sorting technology
WO2014091444A1 (en) * 2012-12-14 2014-06-19 Koninklijke Philips N.V. Detector unit with pulse shaper
WO2015062618A1 (en) * 2013-10-28 2015-05-07 Orexplore Ab Apparatus and method for x-ray transmission analysis of a mineral or electronic waste sample
WO2015085372A1 (en) 2013-12-11 2015-06-18 Southern Innovation International Pty Ltd Method and apparatus for resolving signals in data
US20150338545A1 (en) * 2014-05-23 2015-11-26 Radiabeam Technologies, Llc System and method for adaptive x-ray cargo inspection
WO2016082006A1 (en) * 2014-11-30 2016-06-02 Southern Innovation International Pty Ltd Device and method for material characterisation

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
See also references of EP3465180A4

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2020070473A1 (en) * 2018-10-01 2020-04-09 Smiths Heimann Sas Correction of images
US11983859B2 (en) 2018-10-01 2024-05-14 Smiths Detection France S.A.S. Correction of images
EP3734260A3 (en) * 2019-05-03 2020-12-02 WIPOTEC GmbH Method and apparatus for x-ray inspection of products, in particular food products
US11493457B2 (en) 2019-05-03 2022-11-08 Wipotec Gmbh Method and device for the X-ray inspection of products, in particular foodstuffs
WO2024108257A1 (en) * 2022-11-24 2024-05-30 Southern Innovation International Pty Ltd Density measuring

Also Published As

Publication number Publication date
CA3025821A1 (en) 2017-12-07
AU2017274079A1 (en) 2019-01-17
AU2017274079B2 (en) 2022-06-30
US11243327B2 (en) 2022-02-08
US20190212273A1 (en) 2019-07-11
ZA201808640B (en) 2022-10-26
EP3465180A4 (en) 2020-03-04
EP3465180A1 (en) 2019-04-10
BR112018074796B1 (en) 2023-03-28
BR112018074796A2 (en) 2019-03-12

Similar Documents

Publication Publication Date Title
US10545257B2 (en) Device and method for material characterization
AU2017274079B2 (en) Material characterisation system and method
AU2017274077B2 (en) Material characterisation system and method
US20110135060A1 (en) High Energy X-Ray Inspection System Using a Fan-Shaped Beam and Collimated Backscatter Detectors
JP2011519415A (en) Radiation imaging using individual signal resolution
WO2004095060A2 (en) X-ray imaging technique
US5751000A (en) Prefilter collimator for PET gamma camera
CN104166153A (en) Method and device for measuring two-dimensional angle distribution of radiation dose rate of radioactive substance
Kwong et al. A noise spectroscopy detector array for non-intrusive cargo inspection
CN106124539B (en) Detector and detection system and method for intelligently dividing energy regions
Langeveld et al. Implementation of Noise Spectroscopy using biased large-area photodiodes
Redus et al. Dead time correction in the DP5 digital pulse processor
Smith et al. Refinement of Gamma Spectroscopy Methods for Unattended UF6 Cylinder Verification
Kappen et al. Nonlinearities in solid-state fluorescence detectors
Kim et al. Impact of electric field non-uniformity on large CdZnTe crystals

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 17805394

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 3025821

Country of ref document: CA

NENP Non-entry into the national phase

Ref country code: DE

REG Reference to national code

Ref country code: BR

Ref legal event code: B01A

Ref document number: 112018074796

Country of ref document: BR

ENP Entry into the national phase

Ref document number: 2017805394

Country of ref document: EP

Effective date: 20190102

ENP Entry into the national phase

Ref document number: 2017274079

Country of ref document: AU

Date of ref document: 20170530

Kind code of ref document: A

ENP Entry into the national phase

Ref document number: 112018074796

Country of ref document: BR

Kind code of ref document: A2

Effective date: 20181130