WO2021038621A1 - 画像処理装置および画像処理方法 - Google Patents

画像処理装置および画像処理方法 Download PDF

Info

Publication number
WO2021038621A1
WO2021038621A1 PCT/JP2019/032984 JP2019032984W WO2021038621A1 WO 2021038621 A1 WO2021038621 A1 WO 2021038621A1 JP 2019032984 W JP2019032984 W JP 2019032984W WO 2021038621 A1 WO2021038621 A1 WO 2021038621A1
Authority
WO
WIPO (PCT)
Prior art keywords
atmospheric
aerosol
radiance
visibility
unit
Prior art date
Application number
PCT/JP2019/032984
Other languages
English (en)
French (fr)
Inventor
真梨子 佐藤
中野　貴敬
将敬 鈴木
Original Assignee
三菱電機株式会社
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
Application filed by 三菱電機株式会社 filed Critical 三菱電機株式会社
Priority to PCT/JP2019/032984 priority Critical patent/WO2021038621A1/ja
Priority to JP2021541757A priority patent/JP6964834B2/ja
Publication of WO2021038621A1 publication Critical patent/WO2021038621A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis

Definitions

  • the present invention relates to an image processing apparatus and an image processing method for processing an observed image having image information of a plurality of observed wavelength bands (hereinafter, referred to as a multispectral image).
  • An optical sensor called a multispectral sensor or a hyperspectral sensor detects electromagnetic waves in a plurality of observation wavelength bands among electromagnetic waves radiated from an observation target to generate an observation image having image information for each observation wavelength band. It is possible to do.
  • this type of optical sensor is mounted on an artificial satellite or an aircraft and used for surface observation.
  • the electromagnetic waves detected by the optical sensor include the components of sunlight scattered in the atmosphere and the components of electromagnetic waves that have passed through the atmosphere and reached the sensor. Is done. Therefore, the spectral characteristics of the observed image generated by the optical sensor are different from the spectral characteristics observed on the ground surface. Therefore, when analyzing the spectral characteristics of the earth's surface using observation images, it is necessary to remove the components affected by the atmosphere from the components of the observation images.
  • Non-Patent Document 1 describes a technique for excluding components affected by the atmosphere from components of an observation image generated by MODES using an optical sensor called Moderate Resolution Imaging Spectroradiometer (MODErate Resolution Imaging Spectroradiometer).
  • MODES is an optical sensor having an observation wavelength band of 36 bands, and can detect electromagnetic waves in a wide wavelength range from the visible region to the thermal infrared region.
  • Non-Patent Document 1 uses a 20-band image from the visible region to the thermal infrared region among the 36-band images observable by MODES, and the observation conditions of the sensor (aerosol, water vapor amount, The amount of ozone, cloud region, etc.) is estimated, and the estimated data is used to calculate the spectral reflectance of a 7-band image from the visible region to the infrared region with the influence of the atmosphere reduced.
  • the spectral reflectance that reduces the influence of the atmosphere it is easy to identify and classify subjects with spectral characteristics observed on the ground surface (vegetation or soil, man-made objects, water areas, etc.), or machine learning etc. By doing so, it is possible to grasp the feature amount of the subject that does not depend on the observation conditions.
  • Non-Patent Document 1 Since the conventional technique described in Non-Patent Document 1 uses a 20-band image for estimating the influence of the atmosphere in addition to the 7-band image used for estimating the spectral reflectance of the ground surface, at least An optical sensor capable of observing 27-band images was needed. An optical sensor that observes a wide wavelength range in many bands, such as MODES, can allocate an appropriate band for each observation condition. However, since the number of observable bands is large, the scale of the device including the optical sensor is increased, and when the number of bands used for image processing is large, the load of image processing is also increased.
  • the present invention solves the above problems, and an object of the present invention is to obtain an image processing apparatus and an image processing method capable of reducing the number of images used for image processing and the processing load.
  • the image processing apparatus inputs a plurality of multispectral images from an optical sensor that observes an imaged region on the ground surface and generates a plurality of multispectral images corresponding to each of the plurality of observation wavelength bands.
  • the radiance component of the atmospheric scattered light is the atmospheric scattered radiance Is detected from a multi-spectral image from a plurality of aerosol models set by a scattering characteristic detector that detects the above from a multi-spectral image and a field corresponding to the aerosol and spectral characteristic data of atmospheric scattered radiance and atmospheric transmittance.
  • An atmospheric characteristic extraction unit that extracts rate data, an incident illuminance calculation unit that calculates the incident illuminance of the area to be photographed using the data extracted by the atmospheric characteristic extraction unit, and the radiance value and atmospheric characteristics of the multispectral image. It includes a ground surface radiance calculation unit that calculates the spectral radiance on the ground surface in the area to be photographed using the data extracted by the extraction unit and the incident radiance calculated by the incident radiance calculation unit.
  • the atmospheric scattering radiance detected from the multispectral image is obtained from a plurality of aerosol models in which the field corresponding to the aerosol and the spectral characteristic data of the atmospheric scattering radiance and the atmospheric transmittance are set in association with each other.
  • FIG. 6A is a graph showing the relationship between the visibility corresponding to the aerosol and the atmospheric scattering radiance
  • FIG. 6B is a graph showing the spectral characteristic data of the atmospheric scattering radiance corresponding to the visibility
  • FIG. 6C is the aerosol.
  • FIG. 8A is a functional block diagram showing a hardware configuration that realizes the function of the image processing device according to the first embodiment
  • FIG. 8B is an execution of software that realizes the function of the image processing device according to the first embodiment.
  • It is a functional block diagram which shows the hardware configuration to perform.
  • It is a functional block diagram which shows the structure of the image processing apparatus which concerns on Embodiment 2.
  • FIG. It is a flowchart which shows the image processing method which concerns on Embodiment 2.
  • FIG. 1 is a functional block diagram showing the configuration of the image processing device 1 according to the first embodiment.
  • the image processing device 1 inputs the observed image from the optical sensor that observes the area to be photographed on the ground surface.
  • Optical sensors are sensors called multispectral sensors or hyperspectral sensors that are mounted on platforms such as artificial satellites, aircraft, or space probes so that they can observe the area to be photographed on the surface of a planet with an atmosphere. There is.
  • the image processing device 1 can receive observation image data from the platform by communicating with the platform on which the optical sensor is mounted.
  • FIG. 2 is a diagram showing the configuration of the observed image MSI.
  • Observed image MSI is, N pieces of band (observation wavelength band) B 1, B 2 having a different central wavelength ⁇ each other, ..., multispectral image IB 1 of N frames respectively corresponding to the B N, IB 2, ⁇ ..., and it is configured from the IB N.
  • FIG. 2 shows a case where the number of bands N is 3 or more, the image processing apparatus 1 is not limited to this, and the number of bands N may be 2.
  • the multispectral images IB 1 to IB N all represent the same imaged area. Further, each pixel of the multispectral images IB 1 to IB N can be specified by a set of two-dimensional coordinate values x and y.
  • FIG. 3 is a diagram schematically showing how the ground surface is observed by an optical sensor mounted on the artificial satellite 10.
  • radiance L sensor band B i, i (unit; W / sr / m 2 / ⁇ m) Is detected.
  • the number i is any one of 1 to N.
  • the radiance L sensor, i is expressed by the following equation (1).
  • L sensor, i ⁇ i ⁇ L i + L scat, i ... (1)
  • ⁇ i ⁇ L i is the observation component of the electromagnetic wave reaching the optical sensor passes through the air after being emitted from the surface.
  • tau i is the atmospheric transmittance which depends on the band B i
  • L i is the surface radiance.
  • L scat, i in the above formula (1) is a component in which light scattered in the atmosphere is incident on the optical sensor.
  • the image processing device 1 includes a data storage unit 17 in which an aerosol model database 17A and a skylight database 17B are stored.
  • the aerosol model database 17A stores a plurality of aerosol models in which the field corresponding to the aerosol in the atmosphere is set in association with the spectral characteristic data of the atmospheric scattering radiance L scatt, i and the atmospheric transmittance ⁇ i. ing.
  • the aerosol is suspended fine particles in the atmosphere. Visibility is the distance until the intensity of light propagating in the atmosphere is reduced to a certain percentage (eg, 2%) by being scattered or absorbed by the aerosol in the atmosphere. Visibility corresponding to aerosols is a physical quantity that represents the characteristics of the atmosphere.
  • the skylight database 17B is a database in which the relationship between the aerosol model and the visibility and the illuminance Esky and i of the skylight component is set.
  • Skylight is sunlight that reaches the surface of the earth after being scattered in the atmosphere, or sunlight that reaches the surface of the earth after being reflected by clouds, in addition to the direct light from the sun.
  • i of the skylight component for example, a value obtained by a simulation calculation using MODTRAN or a value actually measured under various environments can be used.
  • the image processing device 1 includes a radiance conversion unit 11, a scattering characteristic detection unit 12, an aerosol / visibility estimation unit 13, an atmospheric characteristic extraction unit 14, an incident illuminance calculation unit 15, and a ground surface reflectance calculation unit 16.
  • Scattering characteristic detector 12 based on the value of the radiance L: sensor, i multispectral image IB i, air scattered radiation from the multispectral image IB i luminance L 'scatt, detects the i.
  • the aerosol / visibility estimation unit 13 uses the atmospheric scattering radiance L' scat, i detected from the multispectral image IB i by the scattering characteristic detection unit 12 from among a plurality of aerosol models stored in the aerosol model database 17A. Estimate the corresponding aerosol model and visibility. Air characteristic extraction unit 14, the aerosol model estimated by aerosol visibility estimator 13, the estimated air scattered radiance L Scatt corresponding to visibility, to extract the spectral characteristic data of i and atmospheric transmittance tau i.
  • Incident illumination calculation unit 15 calculates the incident illuminance E i of the imaging area.
  • Incident illuminance E i is the illuminance E direct of the direct light component of sunlight, i and Skylight component illuminance E sky, the sum of the i, depends on the band B i.
  • the incident illuminance calculation unit 15 uses the aerosol model used when the atmospheric characteristic extraction unit 14 extracts the atmospheric scattering radiance L scat, i and the atmospheric transmittance ⁇ i from the skylight database 17B, and the skylight component corresponding to the visibility.
  • the illuminance E sky, i of the above is acquired, and the incident illuminance E i is calculated.
  • the ground surface reflectance calculation unit 16 uses the radiance L sensor, i , the atmospheric scattering radiance L scatt, i , the atmospheric transmittance ⁇ i , and the incident illuminance E i , and the spectral reflectance on the ground surface in the area to be photographed. Is calculated.
  • the spectral reflectance on the ground surface will be referred to as the ground surface reflectance ⁇ i .
  • FIG. 4 is a flowchart showing the image processing method according to the first embodiment, and shows the procedure of image processing by the image processing device 1.
  • the radiance conversion unit 11 receives the observed image MSI from the optical sensor mounted on the platform by communicating with the platform on which the optical sensor is mounted. Further, instead of inputting the observed image MSI from the optical sensor mounted on the platform, the radiation brightness conversion unit 11 accesses the server on the ground where the data of the observed image MSI is stored and inputs the observed image MSI from the server. You can also do it.
  • the radiance conversion unit 11 converts the pixel values of each pixel of the multispectral images IB 1 to IB N constituting the observed image MSI into radiance L sensor, i (step ST1).
  • the pixel value of each pixel constituting the multispectral images IB 1 to IB N is a numerical value called DN (digital number).
  • DN digital value
  • the radiance L sensor, i of the multispectral image IB i can be calculated by using the following equation (2).
  • ⁇ (band) is a radiance conversion coefficient according to the sensitivity of the luminance value in the optical sensor
  • ⁇ (band) is a radiance conversion coefficient according to the offset of the luminance value in the optical sensor.
  • the radiance conversion unit 11 acquires ⁇ (band) and ⁇ (band) attached to the observed image MSI, and uses the acquired ⁇ (band) and ⁇ (band) to perform a multi according to the following equation (2).
  • radiance L sensor multispectral image IB i, based on the i, from multispectral image IB i, detects the atmospheric scattering radiance L scatt, i (step ST2).
  • the subject having the lowest reflectance is called a "low reflector”.
  • the low reflector include a water area such as the sea surface or a shaded area of an object.
  • the scattering characteristic detection unit 12 detects the radiance of the region corresponding to the low reflector from the radiance L sensor, i of each pixel of the multispectral image IB i.
  • FIG. 5 is a schematic diagram showing an outline of a process for detecting atmospheric scattered radiance L scat, i from the observed image MSI.
  • the scattering characteristic detection unit 12 generates a graph on the right side in FIG. 5 by performing a histogram analysis of the radiance L sensor, i of the multispectral image IB i. Then, the scattering characteristic detection unit 12 detects the radiance of the region corresponding to the low reflector among the pixels constituting the multispectral image IB i.
  • the aerosol / visibility estimation unit 13 uses the aerosol corresponding to the atmospheric scattering radiance L' scat, i detected from the multispectral image IB i from among the plurality of aerosol models stored in the aerosol model database 17A. Estimate the model and visibility (step ST3).
  • the visibility corresponding to the aerosol under various environments is set in association with the spectral characteristic data of the atmospheric scattering radiance L scatt, i and the atmospheric transmittance ⁇ i.
  • FIG. 6A is a graph showing the relationship between the visibility corresponding to the aerosol and the atmospheric scattering radiance L scat, i.
  • the aerosol model for example, data shown in FIG. 6A is set for each band B i.
  • the aerosol / visibility estimation unit 13 refers to a plurality of aerosol models stored in the aerosol model database 17A by using the atmospheric scattering radiance L'scat, i detected by the scattering characteristic detection unit 12. Then, based on the relationship shown in FIG. 6A, the aerosol / visibility estimation unit 13 selects the aerosol model corresponding to the atmospheric scattering radiance L'scat, i and the atmospheric scattering radiance L' scat, from among the plurality of aerosol models. The visibility v1 corresponding to i is estimated.
  • the atmospheric characteristic extraction unit 14 uses the aerosol model estimated by the aerosol / visibility estimation unit 13 to obtain spectral characteristic data of atmospheric scattered radiation brightness L scatt, i and atmospheric transmittance ⁇ i corresponding to the estimated visibility v1.
  • Is extracted step ST4.
  • 6B is atmospheric scattered radiance corresponding to visibility L Scatt, a graph showing the spectral characteristic data of the i, shows an air scattering radiance L scatt, i in each band B1 ⁇ B 8.
  • FIG. 6C is a graph showing the relationship between the visibility v1 corresponding to the aerosol and the atmospheric transmittance ⁇ i.
  • the atmospheric characteristic extraction unit 14 extracts the spectral characteristic data of the atmospheric scattered radiance L scatt, i corresponding to the visibility v1 by referring to the relationship shown in FIG. 6B in the aerosol model estimated by the aerosol / visibility estimation unit 13. To do.
  • the atmospheric scattering radiance L' scat, i estimated from the image and the atmospheric scattering radiance (model value) L extracted from the aerosol model database 17A The scatt and i are regression-analyzed, and the atmospheric scattering radiance value of the image based on the atmospheric scattering radiance L scatt and i extracted from the aerosol model database 17A is calculated.
  • the atmospheric characteristic extraction unit 14 by referring to the relationship shown in Figure 6C in the estimated aerosol model, extracting the atmospheric transmittance tau i corresponding to visibility v1. Spectral characteristic data of atmospheric transmittance ⁇ i is extracted using the visibility estimated for each band B i.
  • the image processing device 1 estimates the data (atmospheric scattered radiance and atmospheric transmittance) used for estimating the influence of the atmosphere from the aerosol model, so that a multispectral image for estimating the influence of the atmosphere can be obtained. Not needed. Therefore, the number of images used for image processing for calculating the ground surface reflectance ⁇ i from the observed image MSI can be reduced, and the number of image processings is reduced accordingly, so that the processing load can also be reduced.
  • the incident illuminance calculation unit 15 calculates the incident illuminance E i of the imaging area (step ST5).
  • FIG. 7 is a diagram showing components of sunlight incident on the surface of the earth. As shown in FIG. 7, the skylight reaches the surface of the earth in addition to the direct light from the sun.
  • the incident illuminance calculation unit 15 calculates the incident illuminance E i according to the following equation (3).
  • E direct, i is the intensity of the directly light component dependent on band B i
  • E sky i is the intensity of the skylight component dependent on band B i.
  • the illuminance E direct, i can be calculated from the solar illuminance E (TOA) at the upper end of the atmosphere (Top Of Atmosphere) and the solar zenith angle ⁇ by using the following equation (4).
  • E i E direct, i + E sky, i ... (3)
  • E direct, i ⁇ i ⁇ E (TOA) cos ⁇ ⁇ ⁇ ⁇ (4)
  • the incident illuminance calculation unit 15 refers to the skylight database 17B and acquires the aerosol model estimated by the aerosol / visibility estimation unit 13 and the illuminance E sky, i of the skylight component corresponding to the visibility.
  • the incident illuminance calculation unit 15 acquires the solar illuminance E (TOA), which is a default value associated with the observed image MSI, and obtains the solar illuminance E (TOA), the illuminance E sky, i of the skylight component, and the atmospheric characteristic extraction unit 14.
  • TOA solar illuminance E
  • the ground surface reflectance calculation unit 16 uses the radiance L sensor, i , the atmospheric scattering radiance L scatt, i , the atmospheric transmittance ⁇ i, and the incident illuminance E i, and according to the following equation (5), the ground surface reflectance ⁇ i. Is calculated (step ST6).
  • the influence of the atmosphere is reduced so that the surface reflectance ⁇ i is close to the original spectral reflectance observed on the surface of the earth.
  • ⁇ i (L sensor, i ⁇ L scat, i ) / ( ⁇ i ⁇ E i / ⁇ ) ⁇ ⁇ ⁇ (5)
  • the image processing device 1 includes a processing circuit that executes the processes from step ST1 to step ST6 in FIG.
  • the processing circuit may be dedicated hardware or a CPU (Central Processing Unit) that executes a program stored in the memory.
  • FIG. 8A is a functional block diagram showing a hardware configuration that realizes the functions of the image processing device 1
  • FIG. 8B is a functional block diagram showing a hardware configuration that executes software that realizes the functions of the image processing device 1. is there.
  • the input / output interface 100 relays the observed image MSI output from the optical sensor (or platform) to the image processing apparatus 1 and the data associated therewith. Further, the input / output interface 100 relays the ground surface reflectance data output from the image processing device 1 to a subsequent device (not shown).
  • the storage device 101 functions as the data storage unit 17 of FIG.
  • the processing circuit 102 may be, for example, a single circuit, a composite circuit, a programmed processor, a parallel programmed processor, or an ASIC (Application Specific Integrated Circuitd). Circuit), FPGA (Field-Programmable Gate Array), or a combination of these is applicable.
  • the processing circuit is the processor 103 shown in FIG. 8B, the radiance conversion unit 11, the scattering characteristic detection unit 12, the aerosol / visibility estimation unit 13, the atmospheric characteristic extraction unit 14, the incident illuminance calculation unit 15, and the incident illuminance calculation unit 15 in the image processing apparatus 1
  • the function of the surface reflectance calculation unit 16 is realized by software, firmware, or a combination of software and firmware.
  • the software or firmware is written as a program and stored in the memory 104.
  • the processor 103 By reading and executing the program stored in the memory 104, the processor 103 reads and executes the radiance conversion unit 11, the scattering characteristic detection unit 12, the aerosol / visibility estimation unit 13, the atmospheric characteristic extraction unit 14, and the incident in the image processing device 1.
  • the functions of the illuminance calculation unit 15 and the surface reflectance calculation unit 16 are realized.
  • the image processing device 1 includes a memory 104 for storing a program in which the processes of steps ST1 to ST6 in the flowchart shown in FIG. 4 are executed as a result when executed by the processor 103.
  • the memory 104 is a program for causing the computer to function as a radiance conversion unit 11, a scattering characteristic detection unit 12, an aerosol / visibility estimation unit 13, an atmospheric characteristic extraction unit 14, an incident illuminance calculation unit 15, and a ground surface reflectance calculation unit 16. May be a computer-readable storage medium in which is stored.
  • the memory 104 is, for example, a RAM (Random Access Memory), a ROM (Read Only Memory), a flash memory, an EPROM (Erasable Programmable Read Only Memory), an EEPROM (Electrically-volatile) semiconductor, an EPROM (Electrically-volatile), or the like.
  • the radiance conversion unit 11, the scattering characteristic detection unit 12, the aerosol / visibility estimation unit 13, and the atmospheric characteristic extraction unit 14 realize their functions by the processing circuit 102, which is dedicated hardware, and the incident illuminance calculation unit 15 and the ground surface.
  • the reflectance calculation unit 16 realizes the function by the processor 103 reading and executing the program stored in the memory 104.
  • the processing circuit can realize the above-mentioned functions by hardware, software, firmware or a combination thereof.
  • the image processing apparatus 1 includes the radiance conversion unit 11 that converts the pixel values of the multispectral image IB i constituting the observed image MSI into radiance, and the multispectral image IB i .
  • the scattering characteristic detection unit 12 that detects the atmospheric scattered radiance from the multispectral image IB i based on the radiance value, and the atmospheric scattered radiance detected from the multi-spectral image IB i from among a plurality of aerosol models.
  • An atmospheric characteristic extraction unit that extracts spectral characteristic data of atmospheric scattering radiance and atmospheric permeability corresponding to the estimated visual range from the corresponding aerosol model and the aerosol / visibility estimation unit 13 that estimates the visibility, and the estimated aerosol model.
  • Atmospheric scattered radiance and atmospheric transmittance are estimated from the aerosol model, so no multispectral imagery and image processing is required to estimate these data. This makes it possible to reduce the number of images used for image processing and the processing load.
  • FIG. 9 is a functional block diagram showing the configuration of the image processing device 1A according to the second embodiment.
  • the image processing device 1A includes a data storage unit 17 in which a skylight database 17B, a ground surface albedo database 17C, and an aerosol model database 17D are stored.
  • the image processing device 1A includes a radiance conversion unit 11, a scattering characteristic detection unit 12, an aerosol / visibility estimation unit 13A, an atmospheric characteristic extraction unit 14, an incident illuminance calculation unit 15, a ground surface reflectance calculation unit 16, and a ground surface albedo calculation unit 18. I have.
  • the functions and configurations of the skylight database 17B, the radiance conversion unit 11, the scattering characteristic detection unit 12, the atmospheric characteristic extraction unit 14, the incident illuminance calculation unit 15, and the ground surface reflectance calculation unit 16 are described in the first embodiment. As explained in.
  • the ground surface albedo database 17C is a database in which the position information of the ground surface and the albedo of the ground surface at that position are set in association with each other.
  • the albedo on the ground surface is an index value indicating the rate at which sunlight is reflected on the ground surface of the earth, and corresponds to the reflectance.
  • the albedo on the ground surface will be referred to as "surface albedo".
  • surface albedo When the surface albedo is high, the proportion of sunlight reflected on the ground surface and scattered in the atmosphere increases, so that the scattered light in the atmosphere increases. For example, if the surface of the earth is concrete or a desert area, the surface albedo is high, and the scattered light in the atmosphere is also increased accordingly.
  • the surface albedo is low, so that the scattered light in the atmosphere is also small.
  • a value obtained by a simulation calculation using MODTRAN or a value actually measured under various environments can be used.
  • each of the plurality of aerosol models is associated with the surface albedo.
  • the visibility corresponding to the aerosol and the surface albedo is set in association with the spectral characteristic data of the atmospheric scattering radiance and the atmospheric transmittance.
  • the ground surface albedo calculation unit 18 calculates the ground surface albedo in the area to be photographed by using the ground surface albedo stored in the ground surface albedo database 17C.
  • the light reflected on the ground surface captured by the observation image MSI is multiple scattered in the atmosphere. Therefore, when the surface albedo calculation unit 18 acquires all the surface albedos corresponding to the imaged area of the observed image MSI from the surface albedo database 17C, the surface albedo calculation unit 18 calculates the average value of the acquired surface albedo.
  • the ground surface albedo calculation unit 18 outputs the average value of the ground surface albedo to the aerosol / visibility estimation unit 13A as the ground surface albedo in the area to be photographed.
  • the aerosol / visibility estimation unit 13A estimates the aerosol model and visibility corresponding to the surface albedo and the atmospheric scattering radiance among the plurality of aerosol models in the aerosol model database 17D. For example, the aerosol / visibility estimation unit 13A inputs the ground surface albedo calculated by the ground surface albedo calculation unit 18, and the scattering characteristic detection unit 12 inputs the atmospheric scattering radiation brightness detected from the multispectral image IB i and inputs it. The aerosol model and visibility corresponding to the surface albedo and atmospheric scattering radiation brightness are estimated from the aerosol model database 17D.
  • FIG. 10 is a flowchart showing the image processing method according to the second embodiment, and shows the procedure of image processing by the image processing device 1A. Since step ST1a and step ST2a of FIG. 10 are the same processes as step ST1 and step ST2 of FIG. 4, and the processes from step ST5a to step ST7a are the same processes as steps ST4 to ST6 of FIG. , The description of these steps will be omitted.
  • the ground surface albedo calculation unit 18 calculates the ground surface albedo of the area to be photographed (step ST3a).
  • the surface albedo calculated by the surface albedo calculation unit 18 is, for example, an average value of the surface albedo acquired from the surface albedo database 17C based on the position information of the imaged area of the observed image MSI.
  • the aerosol / visibility estimation unit 13A is an aerosol model corresponding to the atmospheric scattering radiation brightness detected by the surface albedo calculation unit 18 and the scattering characteristic detection unit 12 among the plurality of aerosol models in the aerosol model database 17D. And the visibility is estimated (step ST4a). The aerosol model and visibility estimated by the aerosol / visibility estimation unit 13A are output to the atmospheric characteristic extraction unit 14.
  • Each function of is realized by a processing circuit. That is, the image processing device 1A includes a processing circuit for executing the processes from step ST1a to step ST7a shown in FIG.
  • the processing circuit may be the processing circuit 102 of the dedicated hardware shown in FIG. 8A, or the processor 103 that executes the program stored in the memory 104 shown in FIG. 8B.
  • the image processing device 1A includes a ground surface albedo calculation unit 18 for calculating the ground surface albedo of the area to be photographed.
  • the plurality of aerosol models in the aerosol model database 17D are associated with the surface albedo.
  • the aerosol / visibility estimation unit 13A estimates the aerosol model and visibility corresponding to the atmospheric scattering radiance detected from the surface albedo and the multispectral image IB i calculated by the surface albedo calculation unit 18. By estimating the aerosol model and visibility considering the surface albedo of the area to be photographed, the estimation accuracy of the aerosol model and visibility is improved.
  • FIG. 11 is a functional block diagram showing the configuration of the image processing device 1B according to the third embodiment.
  • the image processing device 1B includes a data storage unit 17 in which a skylight database 17B, a ground surface albedo database 17C, and an aerosol model database 17D are stored.
  • the image processing device 1B includes a radiance conversion unit 11, a scattering characteristic detection unit 12, an aerosol / visibility estimation unit 13A, an atmospheric characteristic extraction unit 14, an incident illuminance calculation unit 15, a ground surface reflectance calculation unit 16A, and a ground surface albedo calculation unit. It has 18.
  • the ground surface reflectance calculation unit 16A includes a correction unit 19.
  • Each function and each configuration of the calculation unit 18 is as described in the first embodiment and the second embodiment.
  • the surface reflectance calculation unit 16A performs the radiance L sensor, i , the atmospheric scattering radiance L scatt, i , the atmospheric transmittance ⁇ i, and the incident illuminance E i of the multispectral image IB i in the same procedure as in the first embodiment. Is used to calculate the surface reflectance ⁇ i of the area to be photographed.
  • the correction unit 19 corrects the change in the surface reflectance ⁇ i according to the interaction of light between the surface of the earth and the atmosphere.
  • FIG. 12 is a flowchart showing the image processing method according to the third embodiment, and shows the procedure of image processing by the image processing device 1B. Since the processes from step ST1b to step ST7b in FIG. 12 are the same as the processes from step ST1a to step ST7a in FIG. 10, the description of these steps will be omitted.
  • the correction unit 19 corrects the change in the ground surface reflectance ⁇ i according to the interaction of sunlight between the ground surface and the atmosphere (step ST8b).
  • FIG. 13 is a diagram showing the interaction of sunlight between the surface of the earth and the atmosphere. The magnitude of atmospheric scattered light changes according to the difference between the environment in which the surface albedo was obtained in the surface albedo database 17C and the observation environment in which the observation image MSI was taken.
  • the magnitude of atmospheric scattered light or the BRDF (bidirectional scattering distribution function) of the subject may change. This can change the scattered light from the atmosphere to the earth and change the illuminance that illuminates the surface of the earth. Therefore, the incident illuminance E i, the spectral characteristics by a combination of the direct light component and sky light component of sunlight changes.
  • Correcting unit 19 corrects the surface reflectance [rho i, the surface reflectance [rho i 'close to the original value influence is reduced in the change of the skylight component to be incident to the surface.
  • the ground surface albedo stored in the ground surface albedo database 17C is set as the threshold value ⁇ th.
  • the correction unit 19 causes the skylight component of sunlight to interact with the surface of the earth and the atmosphere to be incident on the surface of the earth. It is estimated that the skylight component has changed.
  • Each function of is realized by a processing circuit. That is, the image processing apparatus 1B includes a processing circuit for executing the processes from step ST1b to step ST8b shown in FIG.
  • the processing circuit may be the processing circuit 102 of the dedicated hardware shown in FIG. 8A, or the processor 103 that executes the program stored in the memory 104 shown in FIG. 8B.
  • the image processing apparatus 1B includes a correction unit 19 that corrects a change in the surface reflectance ⁇ i according to the interaction of light between the ground surface and the atmosphere. .. Sky light component to be incident to the surface by interaction between the surface and the atmosphere is changed, it is possible to obtain the surface reflectance [rho i 'with reduced this effect.
  • FIG. 14 is a functional block diagram showing the configuration of the image processing device 1C according to the fourth embodiment.
  • the image processing device 1C includes a data storage unit 17 in which the skylight database 17B and the aerosol model database 17E are stored. Further, the image processing device 1C includes a radiance conversion unit 11, a scattering characteristic detection unit 12, an aerosol / visibility estimation unit 13B, an atmospheric characteristic extraction unit 14, an incident illuminance calculation unit 15, and a ground surface reflectance calculation unit 16.
  • the aerosol / visibility estimation unit 13B includes a conversion unit 13C. The functions and configurations of the skylight database 17B, the radiance conversion unit 11, the scattering characteristic detection unit 12, the atmospheric characteristic extraction unit 14, the incident illuminance calculation unit 15, and the ground surface reflectance calculation unit 16 are described in the first embodiment. As explained.
  • the solar zenith angle ⁇ is small, even if the solar zenith angle ⁇ changes, the absolute value of the atmospheric scattering radiance changes, but the change in the spectral characteristics is small. found. That is, in the range where the solar zenith angle ⁇ is small, if the correction coefficient of the absolute value with respect to the solar zenith angle ⁇ is given, the atmospheric scattering radiance set in the aerosol model will be the atmosphere when the solar zenith angle ⁇ is 0 degrees. It means that even if it is replaced with the scattered radiance, it does not deviate significantly from the correspondence with the atmospheric scattered radiance before the replacement. This characteristic satisfies this relationship when the solar zenith angle ⁇ is about 50 degrees or less when the zenith angle ⁇ sat of the artificial satellite equipped with the optical sensor is about 40 degrees.
  • spectral characteristic data of atmospheric scattered radiance and atmospheric transmittance are set for each solar zenith angle ⁇ .
  • the spectral characteristic data may be associated with the atmospheric scattering radiance corresponding to the solar zenith angle ⁇ of 0 degrees in the range where the solar zenith angle ⁇ is small. Therefore, the aerosol model database 17E significantly reduces the amount of data in the range where the solar zenith angle ⁇ is smaller than that of the aerosol model database 17A.
  • a change in the atmospheric scattering radiance value with respect to the solar zenith angle ⁇ is set as a look-up table (hereinafter referred to as LUT).
  • LUT a look-up table
  • the conversion unit 13C refers to this LUT and corresponds to the atmospheric scattered radiance detected by the scattering characteristic detection unit 12 when the solar zenith angle ⁇ is 0 degrees. Convert to atmospheric scattering radiance.
  • the aerosol model / visibility estimation unit 13B uses the atmospheric scattering radiance converted by the conversion unit 13C for the aerosol model and visibility estimation.
  • FIG. 15 is a flowchart showing the image processing method according to the fourth embodiment, and shows the procedure of image processing by the image processing device 1C. Since the processes of steps ST1c, step ST2c and steps ST6c to ST8c in FIG. 15 are the same as those of steps ST1, step ST2 and steps ST4 to ST6 in FIG. 4, the description of these steps will be omitted.
  • the aerosol / visibility estimation unit 13B determines whether or not the solar zenith angle ⁇ when the image in which the atmospheric scattering radiance is detected by the scattering characteristic detection unit 12 is observed is equal to or less than the threshold value. Determine (step ST3c).
  • the solar zenith angle ⁇ is known data attached to the observed image MSI.
  • the threshold value is, for example, 60 degrees.
  • the conversion unit 13C refers to the LUT and determines the atmospheric scattering radiance detected by the scattering characteristic detection unit 12 when the solar zenith angle ⁇ is 0. It is converted to the corresponding atmospheric scattering radiance when it is a degree (step ST4c).
  • the aerosol / visibility estimation unit 13B estimates the aerosol model and visibility corresponding to the atmospheric scattering radiance converted by the conversion unit 13C from the aerosol model database 17A (step ST5c).
  • the aerosol / visibility estimation unit 13B is the atmosphere detected from the multispectral image IB i by the scattering characteristic detection unit 12 from the aerosol model database 17A.
  • the aerosol model and visibility corresponding to the scattered radiance are estimated (step ST5c).
  • the image processing device 1C includes a processing circuit for executing the processes from step ST1c to step ST8c shown in FIG.
  • the processing circuit may be the processing circuit 102 of the dedicated hardware shown in FIG. 8A, or the processor 103 that executes the program stored in the memory 104 shown in FIG. 8B.
  • the aerosol / visibility estimation unit 13B when the aerosol / visibility estimation unit 13B has a solar zenith angle of less than or equal to a threshold value, the atmospheric scattering radiance detected from the multispectral image is the solar zenith angle.
  • the value converted to the corresponding atmospheric scatter radiance when is 0 degrees is used for the aerosol model and visibility estimation.
  • the aerosol model database 17E may associate the spectral characteristic data with the atmospheric scattering radiance corresponding to the case where the solar zenith angle ⁇ is 0 degrees in the range where the solar zenith angle ⁇ is small. Therefore, it is possible to significantly reduce the amount of data corresponding to the range where the solar zenith angle ⁇ is small in the aerosol model database 17E.
  • the image processing device can be used, for example, for image processing of an observed image observed by an optical sensor mounted on an artificial satellite or an aircraft.
  • 1,1A, 1B, 1C image processing device 10 artificial satellite, 11 radiance conversion unit, 12 scattering characteristic detection unit, 13,13A, 13B aerosol / visibility estimation unit, 13C conversion unit, 14 atmospheric characteristic extraction unit, 15 incident Illumination calculation unit, 16, 16A surface reflectance calculation unit, 17 data storage unit, 17A, 17D, 17E aerosol model database, 17B skylight database, 17C surface albedo database, 18 surface albedo calculation unit, 19 correction unit, 100 input / output Interface, 101 storage device, 102 processing circuit, 103 processor, 104 memory.

Landscapes

  • Engineering & Computer Science (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Image Processing (AREA)
  • Image Analysis (AREA)

Abstract

画像処理装置（１）は、観測画像を構成するマルチスペクトル画像の画素値を放射輝度に変換する放射輝度変換部（１１）と、マルチスペクトル画像から大気散乱放射輝度を検出する散乱特性検出部１２と、エアロゾルに対応する視程と大気散乱放射輝度および大気透過率の分光特性データとが対応付けて設定された複数のエアロゾルモデルから、マルチスペクトル画像から検出された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定するエアロゾル・視程推定部（１３）と、推定されたエアロゾルモデルから、推定された視程に対応する大気散乱放射輝度および大気透過率のデータを抽出する大気特性抽出部（１４）と、被撮影領域の入射照度を算出する入射照度算出部（１５）と、マルチスペクトル画像の放射輝度値、大気散乱放射輝度、大気透過率および入射照度を用いて、地表反射率を算出する地表反射率算出部（１６）を備える。

Description

画像処理装置および画像処理方法
　本発明は、複数の観測波長帯の画像情報（以下、マルチスペクトル画像と記載する）を有した観測画像を処理する画像処理装置および画像処理方法に関する。
　マルチスペクトルセンサまたはハイパースペクトルセンサと呼ばれる光学センサは、観測対象から放射された電磁波のうち、複数の観測波長帯の電磁波を検出することにより、観測波長帯ごとの画像情報を有した観測画像を生成することが可能である。例えば、この種の光学センサは、人工衛星または航空機に搭載されて、地表の観測に使用されている。
　光学センサが高所からの地表の観測に使用される場合、光学センサによって検出される電磁波には、大気で散乱した太陽光の成分と、大気を透過してセンサに到達した電磁波の成分が含まれる。このため、光学センサによって生成された観測画像の分光特性は、地表で観測される分光特性とは異なるものとなる。従って、観測画像を用いて地表の分光特性を解析する場合、観測画像の成分から、大気の影響を受けた成分を除去する必要がある。
　例えば、非特許文献１には、ＭＯＤＩＳ（Ｍｏｄｅｒａｔｅ　Ｒｅｓｏｌｕｔｉｏｎ　Ｉｍａｇｉｎｇ　Ｓｐｅｃｔｒｏｒａｄｉｏｍｅｔｅｒ）と呼ばれる光学センサを用い、ＭＯＤＩＳによって生成された観測画像の成分から、大気の影響を受けた成分を除外する技術が記載されている。ＭＯＤＩＳとは、３６バンドの観測波長帯を有した光学センサであり、可視域から熱赤外域までの広い波長の範囲の電磁波を検出することが可能である。
　非特許文献１に記載された従来の技術は、ＭＯＤＩＳによって観測可能な３６バンドの画像のうち、可視域から熱赤外域までの２０バンドの画像を用いてセンサの観測条件（エアロゾル、水蒸気量、オゾン量および雲領域など）を推定し、推定したデータを用いて大気の影響を低減させた可視域から赤外域までの７バンドの画像の分光反射率を算出する。大気の影響を低減させた分光反射率を用いることで、地表で観測される分光特性を持つ被写体（植生または土壌、人工物、水域など）の識別および分類が容易となり、あるいは、機械学習などを行うことで、観測条件に依存しない被写体の特徴量を把握することも可能である。
　非特許文献１に記載された従来の技術は、地表面の分光反射率を推定するために用いる７バンドの画像に加えて、２０バンドの画像を大気の影響を推定するために用いるので、少なくとも２７バンドの画像を観測可能な光学センサが必要であった。ＭＯＤＩＳのように広い波長の範囲を多くのバンドで観測する光学センサは、観測条件ごとに適正なバンドを割り当てることはできる。しかしながら、観測可能なバンド数が多い分、光学センサを含む装置の規模が増大し、さらに、画像処理に使用されるバンド数が多いと、画像処理の負荷も高くなるという課題があった。
　本発明は上記課題を解決するものであり、画像処理に使用される画像数および処理負荷を低減することができる画像処理装置および画像処理方法を得ることを目的とする。
　本発明に係る画像処理装置は、地表の被撮影領域を観測して複数の観測波長帯のそれぞれに対応する複数のマルチスペクトル画像を生成する光学センサから、複数のマルチスペクトル画像を入力し、入力した複数のマルチスペクトル画像の各画素の画素値を放射輝度に変換する放射輝度変換部と、複数のマルチスペクトル画像の放射輝度値に基づいて、大気散乱光の放射輝度成分である大気散乱放射輝度をマルチスペクトル画像から検出する散乱特性検出部と、エアロゾルに対応する視程と大気散乱放射輝度および大気透過率の分光特性データとが対応付けて設定された複数のエアロゾルモデルから、マルチスペクトル画像から検出された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定するエアロゾル・視程推定部と、エアロゾル・視程推定部によって推定されたエアロゾルモデルから、推定された視程に対応する大気散乱放射輝度および大気透過率のデータを抽出する大気特性抽出部と、大気特性抽出部によって抽出されたデータを用いて、被撮影領域の入射照度を算出する入射照度算出部と、マルチスペクトル画像の放射輝度値、大気特性抽出部によって抽出されたデータ、および入射照度算出部によって算出された入射照度を用いて、被撮影領域における地表面での分光反射率を算出する地表反射率算出部とを備える。
　本発明によれば、エアロゾルに対応する視程と大気散乱放射輝度および大気透過率の分光特性データとが対応付けて設定された複数のエアロゾルモデルから、マルチスペクトル画像から検出された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定し、推定したエアロゾルモデルから抽出された視程に対応する大気散乱放射輝度および大気透過率のデータと、マルチスペクトル画像の放射輝度値および入射照度を用いて、地表面での分光反射率を算出する。エアロゾルモデルを用いて大気散乱放射輝度および大気透過率が推定されるので、地表面の分光反射率の算出に使用するマルチスペクトル画像とは別に、大気の影響を推定するためのマルチスペクトル画像とその画像処理が不要である。このため、画像処理に使用される画像数および処理負荷を低減することができる。
実施の形態１に係る画像処理装置の構成を示す機能ブロック図である。 観測画像の構成を示す図である。 人工衛星に搭載された光学センサによって地表が観測される様子を模式的に示す図である。 実施の形態１に係る画像処理方法を示すフローチャートである。 観測画像から大気散乱放射輝度を検出する処理を示す概要図である。 図６Ａは、エアロゾルに対応する視程と大気散乱放射輝度との関係を示すグラフであり、図６Ｂは、視程に対応した大気散乱放射輝度の分光特性データを示すグラフであり、図６Ｃは、エアロゾルに対応した視程と大気透過率との関係を示すグラフである。 地表に入射される太陽光の成分を示す図である。 図８Ａは、実施の形態１に係る画像処理装置の機能を実現するハードウェア構成を示す機能ブロック図であり、図８Ｂは、実施の形態１に係る画像処理装置の機能を実現するソフトウェアを実行するハードウェア構成を示す機能ブロック図である。 実施の形態２に係る画像処理装置の構成を示す機能ブロック図である。 実施の形態２に係る画像処理方法を示すフローチャートである。 実施の形態３に係る画像処理装置の構成を示す機能ブロック図である。 実施の形態３に係る画像処理方法を示すフローチャートである。 地表と大気との間での太陽光の相互作用を示す図である。 実施の形態４に係る画像処理装置の構成を示す機能ブロック図である。 実施の形態４に係る画像処理方法を示すフローチャートである。
実施の形態１．
　図１は、実施の形態１に係る画像処理装置１の構成を示す機能ブロック図である。画像処理装置１は、地表の被撮影領域を観測した光学センサから観測画像を入力する。光学センサは、マルチスペクトルセンサまたはハイパースペクトルセンサと呼ばれるセンサであり、大気を有する惑星の地表の被撮影領域を観測することができるように、人工衛星、航空機または宇宙探査機といったプラットフォームに搭載されている。画像処理装置１は、光学センサを搭載するプラットフォームと通信することによって、当該プラットフォームから観測画像のデータを受信することができる。
　図２は、観測画像ＭＳＩの構成を示す図である。観測画像ＭＳＩは、互いに異なる中心波長λを有するＮ個のバンド（観測波長帯）Ｂ１，Ｂ２，・・・，ＢＮにそれぞれ対応するＮ枚のマルチスペクトル画像ＩＢ１，ＩＢ２，・・・，ＩＢＮから構成されている。なお、図２には、バンド数Ｎが３以上である場合を示したが、画像処理装置１は、これに限定されるものではなく、バンド数Ｎは２であってもよい。マルチスペクトル画像ＩＢ１～ＩＢＮは、全て同じ被撮影領域を表している。また、マルチスペクトル画像ＩＢ１～ＩＢＮの各画素は、２次元座標値ｘ，ｙの組みで特定することができる。
　図３は、人工衛星１０に搭載された光学センサによって地表が観測される様子を模式的に示す図である。人工衛星１０に搭載された光学センサは、地表の被撮影領域の方向から入射した電磁波を観測することで、バンドＢｉの放射輝度Ｌｓｅｎｓｏｒ，ｉ（単位；Ｗ／ｓｒ／ｍ２／μｍ）を検出する。番号ｉは、１～Ｎのいずれかである。また、放射輝度Ｌｓｅｎｓｏｒ，ｉは、下記式（１）によって表現される。
Ｌｓｅｎｓｏｒ，ｉ＝τｉ×Ｌｉ＋Ｌｓｃａｔｔ，ｉ　　　・・・（１）
　上記式（１）において、τｉ×Ｌｉは、地表から放射された後に大気を透過して光学センサに到達した電磁波の観測成分である。τｉは、バンドＢｉに依存する大気透過率であり、Ｌｉは、地表の放射輝度である。上記式（１）におけるＬｓｃａｔｔ，ｉは、大気で散乱された光が光学センサに入射した成分である。
　図１に示すように、画像処理装置１は、エアロゾルモデルデータベース１７Ａおよびスカイライトデータベース１７Ｂが格納されたデータ記憶部１７を備えている。エアロゾルモデルデータベース１７Ａには、大気中のエアロゾルに対応する視程と、大気散乱放射輝度Ｌｓｃａｔｔ，ｉおよび大気透過率τｉの分光特性データとが対応付けて設定された複数のエアロゾルモデルが記憶されている。なお、エアロゾルは、大気中の浮遊微粒子である。視程は、大気中のエアロゾルによって、大気を伝搬する光が散乱または吸収されることで、その強度が一定割合（例えば、２％）の値に減少するまでの距離である。エアロゾルに対応する視程は、大気の特性を表す物理量である。
　エアロゾルモデルにおいて、バンドＢｉに依存した大気散乱放射輝度Ｌｓｃａｔｔ，ｉ（分光特性データ）およびバンドＢｉに依存した大気透過率τｉ（分光特性データ）のそれぞれの値が視程に対応している。エアロゾルモデルに設定される各データは、大気の影響を考慮した光波の伝搬環境をモデル化した大気散乱モデルおよび大気透過モデルに基づくシミュレーション計算によって得られる。例えば、ＭＯＤＴＲＡＮ（ＭＯＤｅｒａｔｅ　ｒｅｓｏｌｕｔｉｏｎ　ａｔｍｏｓｐｈｅｒｉｃ　ＴＲＡＮｓｍｉｓｓｉｏｎ）と呼ばれる光波大気伝搬計算シミュレータを用いて算出することが可能である。なお、エアロゾルモデルには、様々な環境下で実測されたデータを設定することができる。
　スカイライトデータベース１７Ｂは、エアロゾルモデルおよび視程とスカイライト成分の照度Ｅｓｋｙ，ｉの関係が設定されたデータベースである。スカイライト（天空光）は、太陽からの直達光とは別に、大気中で散乱された後に地表に到達した太陽光、または、雲で反射された後に地表に到達した太陽光などである。スカイライト成分の照度Ｅｓｋｙ，ｉは、例えば、ＭＯＤＴＲＡＮを用いたシミュレーション計算で得られた値、もしくは様々な環境下で実測された値を使用することができる。
　さらに、画像処理装置１は、放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６を備えている。放射輝度変換部１１は、光学センサからマルチスペクトル画像ＩＢｉ（ｉ＝１～Ｎ）を入力し、マルチスペクトル画像ＩＢｉの各画素の画素値を放射輝度Ｌｓｅｎｓｏｒ，ｉに変換する。散乱特性検出部１２は、マルチスペクトル画像ＩＢｉの放射輝度Ｌｓｅｎｓｏｒ，ｉの値に基づいて、マルチスペクトル画像ＩＢｉから大気散乱放射輝度Ｌ’ｓｃａｔｔ，ｉを検出する。
　エアロゾル・視程推定部１３は、エアロゾルモデルデータベース１７Ａに記憶されている複数のエアロゾルモデルの中から、散乱特性検出部１２によってマルチスペクトル画像ＩＢｉから検出された大気散乱放射輝度Ｌ’ｓｃａｔｔ，ｉに対応するエアロゾルモデルおよび視程を推定する。大気特性抽出部１４は、エアロゾル・視程推定部１３によって推定されたエアロゾルモデルから、推定された視程に対応する大気散乱放射輝度Ｌｓｃａｔｔ，ｉおよび大気透過率τｉの分光特性データを抽出する。
　入射照度算出部１５は、被撮影領域の入射照度Ｅｉを算出する。入射照度Ｅｉは、太陽光の直達光成分の照度Ｅｄｉｒｅｃｔ，ｉとスカイライト成分の照度Ｅｓｋｙ，ｉとの和であり、バンドＢｉに依存する。入射照度算出部１５は、スカイライトデータベース１７Ｂから、大気特性抽出部１４が大気散乱放射輝度Ｌｓｃａｔｔ，ｉおよび大気透過率τｉを抽出したときに用いたエアロゾルモデルおよび視程に対応するスカイライト成分の照度Ｅｓｋｙ，ｉを取得して、入射照度Ｅｉを算出する。地表反射率算出部１６は、放射輝度Ｌｓｅｎｓｏｒ，ｉ、大気散乱放射輝度Ｌｓｃａｔｔ，ｉ、大気透過率τｉ、および入射照度Ｅｉを用いて、被撮影領域における地表面での分光反射率を算出する。以下の説明では、地表面での分光反射率を、地表反射率ρｉと記載する。
　図４は、実施の形態１に係る画像処理方法を示すフローチャートであり、画像処理装置１による画像処理の手順を示している。まず、放射輝度変換部１１は、光学センサを搭載するプラットフォームと通信することで、プラットフォームに搭載された光学センサから観測画像ＭＳＩを受信する。また、放射輝度変換部１１は、プラットフォームに搭載された光学センサから観測画像ＭＳＩを入力する代わりに、観測画像ＭＳＩのデータが蓄積された地上のサーバにアクセスして、サーバから観測画像ＭＳＩを入力することもできる。
　放射輝度変換部１１は、観測画像ＭＳＩを構成するマルチスペクトル画像ＩＢ１～ＩＢＮの各画素の画素値を、放射輝度Ｌｓｅｎｓｏｒ，ｉに変換する（ステップＳＴ１）。例えば、マルチスペクトル画像ＩＢ１～ＩＢＮを構成している各画素の画素値は、ＤＮ（デジタルナンバー）と呼ばれる数値である。バンドＢｉが可視域である場合に、ＤＮは、可視域を示すデジタル値ＤＮ（ｂａｎｄ）である。
　マルチスペクトル画像ＩＢｉの放射輝度Ｌｓｅｎｓｏｒ，ｉは、下記式（２）を用いて算出することができる。α（ｂａｎｄ）は、光学センサにおける輝度値の感度に応じた放射輝度変換係数であり、β（ｂａｎｄ）は、光学センサにおける輝度値のオフセットに応じた放射輝度変換係数である。放射輝度変換部１１は、観測画像ＭＳＩに付随されたα（ｂａｎｄ）とβ（ｂａｎｄ）を取得し、取得したα（ｂａｎｄ）およびβ（ｂａｎｄ）を用いて、下記式（２）に従い、マルチスペクトル画像ＩＢｉの各画素の画素値を放射輝度Ｌｓｅｎｓｏｒ，ｉに変換する。
Ｌｓｅｎｓｏｒ，ｉ＝α（ｂａｎｄ）×ＤＮ（ｂａｎｄ）＋β（ｂａｎｄ）・・・（２）
　散乱特性検出部１２は、マルチスペクトル画像ＩＢｉの放射輝度Ｌｓｅｎｓｏｒ，ｉに基づいて、マルチスペクトル画像ＩＢｉから、大気散乱放射輝度Ｌｓｃａｔｔ，ｉを検出する（ステップＳＴ２）。観察対象の地表に存在している被写体のうち、最も反射率が低い被写体を“低反射体”と呼ぶ。低反射体としては、例えば、海面などの水域または物体の陰影領域が挙げられる。散乱特性検出部１２は、マルチスペクトル画像ＩＢｉの各画素の放射輝度Ｌｓｅｎｓｏｒ，ｉから、低反射体に対応する領域の放射輝度を検出する。
　図５は、観測画像ＭＳＩから、大気散乱放射輝度Ｌｓｃａｔｔ，ｉを検出する処理の概要を示す概要図である。散乱特性検出部１２は、マルチスペクトル画像ＩＢｉの放射輝度Ｌｓｅｎｓｏｒ，ｉをヒストグラム解析することで、図５中の右側のグラフを生成する。そして、散乱特性検出部１２は、マルチスペクトル画像ＩＢｉを構成する画素のうち、低反射体に対応する領域の放射輝度を検出する。
　低反射体が広範囲に一様かつ平坦な被写体（例えば、海面）である場合、地表面の分光反射率がほぼ０である。この場合、上記式（１）における大気透過量τｉ×Ｌｉが０になるので、マルチスペクトル画像ＩＢｉから検出された放射輝度Ｌｓｅｎｓｏｒ，ｉの値は、ほぼ大気散乱放射輝度Ｌｓｃａｔｔ，ｉに一致する。この値は、マルチスペクトル画像から検出された推定値であるため、Ｌ’ｓｃａｔｔ，ｉと記載する。
　次に、エアロゾル・視程推定部１３は、エアロゾルモデルデータベース１７Ａに記憶されている複数のエアロゾルモデルの中から、マルチスペクトル画像ＩＢｉから検出された大気散乱放射輝度Ｌ’ｓｃａｔｔ，ｉに対応するエアロゾルモデルおよび視程を推定する（ステップＳＴ３）。なお、エアロゾルモデルには、様々な環境下における、エアロゾルに対応した視程と、大気散乱放射輝度Ｌｓｃａｔｔ，ｉおよび大気透過率τｉの分光特性データとが対応付けて設定されている。
　図６Ａは、エアロゾルに対応した視程と大気散乱放射輝度Ｌｓｃａｔｔ，ｉとの関係を示すグラフである。エアロゾルモデルには、例えば、図６Ａに示すデータがバンドＢｉごとに設定されている。まず、エアロゾル・視程推定部１３は、散乱特性検出部１２によって検出された大気散乱放射輝度Ｌ’ｓｃａｔｔ，ｉを用いて、エアロゾルモデルデータベース１７Ａに記憶されている複数のエアロゾルモデルを参照する。そして、エアロゾル・視程推定部１３は、図６Ａに示す関係に基づいて、複数のエアロゾルモデルの中から、大気散乱放射輝度Ｌ’ｓｃａｔｔ，ｉに対応したエアロゾルモデルおよび大気散乱放射輝度Ｌ’ｓｃａｔｔ，ｉに対応した視程ｖ１を推定する。
　次に、大気特性抽出部１４は、エアロゾル・視程推定部１３によって推定されたエアロゾルモデルから、推定された視程ｖ１に対応する大気散乱放射輝度Ｌｓｃａｔｔ，ｉおよび大気透過率τｉの分光特性データを抽出する（ステップＳＴ４）。図６Ｂは、視程に対応した大気散乱放射輝度Ｌｓｃａｔｔ，ｉの分光特性データを示すグラフであり、バンドＢ１～Ｂ８のそれぞれにおける大気散乱放射輝度Ｌｓｃａｔｔ，ｉを示している。図６Ｃは、エアロゾルに対応した視程ｖ１と大気透過率τｉとの関係を示すグラフである。
　大気特性抽出部１４は、エアロゾル・視程推定部１３によって推定されたエアロゾルモデルにおける図６Ｂに示す関係を参照することにより、視程ｖ１に対応した大気散乱放射輝度Ｌｓｃａｔｔ，ｉの分光特性データを抽出する。ここで、データは離散的な視程の値に対して入力されるため、画像から推定した大気散乱放射輝度Ｌ’ｓｃａｔｔ，ｉとエアロゾルモデルデータベース１７Ａから抽出された大気散乱放射輝度（モデル値）Ｌｓｃａｔｔ，ｉを回帰分析し、エアロゾルモデルデータベース１７Ａから抽出された大気散乱放射輝度Ｌｓｃａｔｔ，ｉに基づいた画像の大気散乱放射輝度値を算出する。例えば、大気散乱の大きさのみ変化するものとして、全バンドに固有の係数αを考慮したα×Ｌｓｃａｔｔ，ｉを画像の大気散乱放射輝度値とする。さらに、大気特性抽出部１４は、推定されたエアロゾルモデルにおける図６Ｃに示す関係を参照することにより、視程ｖ１に対応した大気透過率τｉを抽出する。バンドＢｉごとに推定された視程を用いて大気透過率τｉの分光特性データが抽出される。
　このように、画像処理装置１は、大気の影響の推定に使用されるデータ（大気散乱放射輝度および大気透過率）をエアロゾルモデルから推定するので、大気の影響を推定するためのマルチスペクトル画像が不要である。このため、観測画像ＭＳＩから地表反射率ρｉを算出する画像処理に使用される画像数を低減することができ、これに伴い画像処理の回数が減るので、処理負荷も低減することができる。
　次に、入射照度算出部１５は、被撮影領域の入射照度Ｅｉを算出する（ステップＳＴ５）。図７は、地表に入射される太陽光の成分を示す図である。図７に示すように、地表には、太陽からの直達光に加え、スカイライトが到達する。入射照度算出部１５は、下記式（３）に従い、入射照度Ｅｉを算出する。下記式（３）において、Ｅｄｉｒｅｃｔ，ｉは、バンドＢｉに依存する直達光成分の照度であり、Ｅｓｋｙ，ｉは、バンドＢｉに依存するスカイライト成分の照度である。なお、照度Ｅｄｉｒｅｃｔ，ｉは、大気上端（Ｔｏｐ　Ｏｆ　Ａｔｍｏｓｐｈｅｒｅ）の太陽照度Ｅ（ＴＯＡ）および太陽天頂角θから、下記式（４）を用いて算出することが可能である。
Ｅｉ＝Ｅｄｉｒｅｃｔ，ｉ＋Ｅｓｋｙ，ｉ　　　　　　　　・・・（３）
Ｅｄｉｒｅｃｔ，ｉ＝τｉ×Ｅ（ＴＯＡ）ｃｏｓθ　　　　・・・（４）
　入射照度算出部１５は、スカイライトデータベース１７Ｂを参照して、エアロゾル・視程推定部１３によって推定されたエアロゾルモデルおよび視程に対応したスカイライト成分の照度Ｅｓｋｙ，ｉを取得する。入射照度算出部１５は、観測画像ＭＳＩに付随した既定値である太陽照度Ｅ（ＴＯＡ）を取得し、太陽照度Ｅ（ＴＯＡ）、スカイライト成分の照度Ｅｓｋｙ，ｉ、および大気特性抽出部１４によって抽出された大気透過率τｉを用いて、上記式（３）および（４）に従い、入射照度Ｅｉを算出する。
　地表反射率算出部１６は、放射輝度Ｌｓｅｎｓｏｒ，ｉ、大気散乱放射輝度Ｌｓｃａｔｔ，ｉ、大気透過率τｉおよび入射照度Ｅｉを用いて、下記式（５）に従い、地表反射率ρｉを算出する（ステップＳＴ６）。地表反射率ρｉは、地表で観測される本来の分光反射率に近い値になるように大気の影響が低減されている。この地表反射率ρｉを用いることで、観測画像ＭＳＩを、地表で観測される分光特性に合うように補正することが可能である。
ρｉ＝（Ｌｓｅｎｓｏｒ，ｉ－Ｌｓｃａｔｔ，ｉ）／（τｉ×Ｅｉ／π）　・・・（５）
　次に、画像処理装置１の機能を実現するハードウェア構成について説明する。
　画像処理装置１における放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６の各機能は、処理回路によって実現される。すなわち、画像処理装置１は、図４のステップＳＴ１からステップＳＴ６までの処理を実行する処理回路を備える。処理回路は、専用のハードウェアであってもよいし、メモリに記憶されたプログラムを実行するＣＰＵ（Ｃｅｎｔｒａｌ　Ｐｒｏｃｅｓｓｉｎｇ　Ｕｎｉｔ）であってもよい。
　図８Ａは、画像処理装置１の機能を実現するハードウェア構成を示す機能ブロック図であり、図８Ｂは、画像処理装置１の機能を実現するソフトウェアを実行するハードウェア構成を示す機能ブロック図である。図８Ａおよび図８Ｂにおいて、入出力インタフェース１００は、光学センサ（またはプラットフォーム）から画像処理装置１へ出力される観測画像ＭＳＩおよびこれに付随したデータを中継する。さらに、入出力インタフェース１００は、画像処理装置１から後段の装置（不図示）へ出力される地表反射率のデータを中継する。記憶装置１０１は、図１のデータ記憶部１７として機能する。
　処理回路が、図８Ａに示す専用のハードウェアの処理回路１０２である場合、処理回路１０２は、例えば、単一回路、複合回路、プログラム化したプロセッサ、並列プログラム化したプロセッサ、ＡＳＩＣ（Ａｐｐｌｉｃａｔｉｏｎ　Ｓｐｅｃｉｆｉｃ　Ｉｎｔｅｇｒａｔｅｄ　Ｃｉｒｃｕｉｔ）、ＦＰＧＡ（Ｆｉｅｌｄ－Ｐｒｏｇｒａｍｍａｂｌｅ　Ｇａｔｅ　Ａｒｒａｙ）、または、これらを組み合わせたものが該当する。画像処理装置１における放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６の各機能を別々の処理回路で実現してもよいし、これらの機能をまとめて１つの処理回路で実現してもよい。
　処理回路が、図８Ｂに示すプロセッサ１０３である場合、画像処理装置１における放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６の機能は、ソフトウェア、ファームウェアまたはソフトウェアとファームウェアの組み合わせによって実現される。ソフトウェアまたはファームウェアは、プログラムとして記述されてメモリ１０４に記憶される。
　プロセッサ１０３は、メモリ１０４に記憶されたプログラムを読み出して実行することにより、画像処理装置１における放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６の機能を実現する。例えば、画像処理装置１は、プロセッサ１０３によって実行されるときに、図４に示すフローチャートにおけるステップＳＴ１からステップＳＴ６の処理が結果的に実行されるプログラムを記憶するためのメモリ１０４を備える。これらのプログラムは、放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６の手順または方法を、コンピュータに実行させる。メモリ１０４は、コンピュータを、放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６として機能させるためのプログラムが記憶されたコンピュータ可読記憶媒体であってもよい。
　メモリ１０４は、例えば、ＲＡＭ（Ｒａｎｄｏｍ　Ａｃｃｅｓｓ　Ｍｅｍｏｒｙ）、ＲＯＭ（Ｒｅａｄ　Ｏｎｌｙ　Ｍｅｍｏｒｙ）、フラッシュメモリ、ＥＰＲＯＭ（Ｅｒａｓａｂｌｅ　Ｐｒｏｇｒａｍｍａｂｌｅ　Ｒｅａｄ　Ｏｎｌｙ　Ｍｅｍｏｒｙ）、ＥＥＰＲＯＭ（Ｅｌｅｃｔｒｉｃａｌｌｙ－ＥＰＲＯＭ）などの不揮発性または揮発性の半導体メモリ、磁気ディスク、フレキシブルディスク、光ディスク、コンパクトディスク、ミニディスク、ＤＶＤなどが該当する。
　画像処理装置１における放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６の機能の一部を専用ハードウェアで実現し、一部をソフトウェアまたはファームウェアで実現してもよい。例えば、放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３および大気特性抽出部１４は、専用のハードウェアである処理回路１０２によって機能を実現し、入射照度算出部１５および地表反射率算出部１６は、プロセッサ１０３がメモリ１０４に記憶されたプログラムを読み出し実行することで、機能を実現する。このように、処理回路は、ハードウェア、ソフトウェア、ファームウェアまたはこれらの組み合わせにより、上記機能を実現することができる。
　以上のように、実施の形態１に係る画像処理装置１は、観測画像ＭＳＩを構成するマルチスペクトル画像ＩＢｉの画素値を放射輝度に変換する放射輝度変換部１１と、マルチスペクトル画像ＩＢｉの放射輝度の値に基づいて、マルチスペクトル画像ＩＢｉから大気散乱放射輝度を検出する散乱特性検出部１２と、複数のエアロゾルモデルの中から、マルチスペクトル画像ＩＢｉから検出された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定するエアロゾル・視程推定部１３と、推定されたエアロゾルモデルから、推定された視程に対応する大気散乱放射輝度および大気透過率の分光特性データを抽出する大気特性抽出部１４と、大気透過率データを用いて、入射照度Ｅｉを算出する入射照度算出部１５と、マルチスペクトル画像ＩＢｉの放射輝度値、大気散乱放射輝度および大気透過率のデータおよび入射照度Ｅｉを用いて、地表反射率ρｉを算出する地表反射率算出部１６を備えている。エアロゾルモデルから大気散乱放射輝度および大気透過率が推定されるので、これらのデータを推定するためにマルチスペクトル画像および画像処理が不要である。これにより、画像処理に使用される画像数および処理負荷を低減することが可能である。
実施の形態２．
　図９は、実施の形態２に係る画像処理装置１Ａの構成を示す機能ブロック図である。図９に示すように、画像処理装置１Ａは、スカイライトデータベース１７Ｂ、地表アルベドデータベース１７Ｃおよびエアロゾルモデルデータベース１７Ｄが格納されたデータ記憶部１７を備えている。画像処理装置１Ａは、放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３Ａ、大気特性抽出部１４、入射照度算出部１５、地表反射率算出部１６および地表アルベド算出部１８を備えている。ここで、スカイライトデータベース１７Ｂ、放射輝度変換部１１、散乱特性検出部１２、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６の各機能および各構成は、実施の形態１で説明した通りである。
　地表アルベドデータベース１７Ｃは、地表の位置情報とその位置における地表面のアルベドとが対応付けて設定されたデータベースである。地表面のアルベドは、太陽光が地球の地表面で反射される割合を示す指標値であり、反射率に相当する。以下、地表面のアルベドを“地表アルベド”と記載する。地表アルベドが高い場合、地表面で反射して大気中で散乱される太陽光の割合が増加するため、大気での散乱光が大きくなる。例えば、地表がコンクリートまたは砂漠領域であれば、地表アルベドが高く、これに伴って大気散乱光も大きくなる。一方、海面または陰影領域では、地表アルベドが低いため、大気散乱光も小さくなる。地表アルベドデータベース１７Ｃに設定される地表アルベドには、例えば、ＭＯＤＴＲＡＮを用いたシミュレーション計算で得られた値、もしくは、様々な環境下で実測された値を使用することができる。
　エアロゾルモデルデータベース１７Ｄにおいて、複数のエアロゾルモデルのそれぞれが地表アルベドに対応付けられている。例えば、エアロゾルモデルには、エアロゾルおよび地表アルベドに対応する視程と、大気散乱放射輝度および大気透過率の分光特性データとが対応付けて設定されている。
　地表アルベド算出部１８は、地表アルベドデータベース１７Ｃに格納されている地表アルベドを用いて、被撮影領域における地表アルベドを算出する。なお、観測画像ＭＳＩに撮影された地表で反射された光は、大気中で多重散乱されている。そこで、地表アルベド算出部１８は、観測画像ＭＳＩの被撮影領域に対応する全ての地表アルベドを地表アルベドデータベース１７Ｃから取得すると、取得した地表アルベドの平均値を算出する。地表アルベド算出部１８は、地表アルベドの平均値を、被撮影領域における地表アルベドとしてエアロゾル・視程推定部１３Ａに出力する。
　エアロゾル・視程推定部１３Ａは、エアロゾルモデルデータベース１７Ｄにおける複数のエアロゾルモデルのうち、地表アルベドおよび大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定する。例えば、エアロゾル・視程推定部１３Ａは、地表アルベド算出部１８によって算出された地表アルベドを入力し、散乱特性検出部１２によってマルチスペクトル画像ＩＢｉから検出された大気散乱放射輝度を入力して、入力した地表アルベドおよび大気散乱放射輝度に対応するエアロゾルモデルおよび視程を、エアロゾルモデルデータベース１７Ｄから推定する。
　図１０は、実施の形態２に係る画像処理方法を示すフローチャートであって、画像処理装置１Ａによる画像処理の手順を示している。図１０のステップＳＴ１ａおよびステップＳＴ２ａは、図４のステップＳＴ１およびステップＳＴ２と同様の処理であり、ステップＳＴ５ａからステップＳＴ７ａまでの処理は、図４のステップＳＴ４からステップＳＴ６と同様の処理であるので、これらのステップの説明は省略する。
　ステップＳＴ２ａの処理が終了した後、地表アルベド算出部１８は、被撮影領域の地表アルベドを算出する（ステップＳＴ３ａ）。地表アルベド算出部１８によって算出される地表アルベドは、例えば、観測画像ＭＳＩの被撮影領域の位置情報に基づいて、地表アルベドデータベース１７Ｃから取得された地表アルベドの平均値である。
　エアロゾル・視程推定部１３Ａは、エアロゾルモデルデータベース１７Ｄにおける複数のエアロゾルモデルのうち、地表アルベド算出部１８によって算出された地表アルベドおよび散乱特性検出部１２によって検出された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定する（ステップＳＴ４ａ）。エアロゾル・視程推定部１３Ａによって推定されたエアロゾルモデルおよび視程は、大気特性抽出部１４に出力される。
　なお、画像処理装置１Ａにおける放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３Ａ、大気特性抽出部１４、入射照度算出部１５、地表反射率算出部１６および地表アルベド算出部１８の各機能は、処理回路によって実現される。すなわち、画像処理装置１Ａは、図１０に示したステップＳＴ１ａからステップＳＴ７ａまでの処理を実行するための処理回路を備えている。処理回路は、図８Ａに示した専用のハードウェアの処理回路１０２であってもよいし、図８Ｂに示したメモリ１０４に記憶されたプログラムを実行するプロセッサ１０３であってもよい。
　以上のように、実施の形態２に係る画像処理装置１Ａは、被撮影領域の地表アルベドを算出する地表アルベド算出部１８を備えている。エアロゾルモデルデータベース１７Ｄにおける複数のエアロゾルモデルは、地表アルベドに対応付けられている。エアロゾル・視程推定部１３Ａは、地表アルベド算出部１８によって算出された地表アルベドおよびマルチスペクトル画像ＩＢｉから検出された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定する。被撮影領域の地表アルベドが考慮されたエアロゾルモデルおよび視程を推定することで、エアロゾルモデルおよび視程の推定精度が向上する。
実施の形態３．
　図１１は、実施の形態３に係る画像処理装置１Ｂの構成を示す機能ブロック図である。図１１に示すように、画像処理装置１Ｂは、スカイライトデータベース１７Ｂ、地表アルベドデータベース１７Ｃおよびエアロゾルモデルデータベース１７Ｄが格納されたデータ記憶部１７を備えている。さらに、画像処理装置１Ｂは、放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３Ａ、大気特性抽出部１４、入射照度算出部１５、地表反射率算出部１６Ａおよび地表アルベド算出部１８を備えている。また、地表反射率算出部１６Ａは、補正部１９を備えている。
　なお、スカイライトデータベース１７Ｂ、地表アルベドデータベース１７Ｃ、エアロゾルモデルデータベース１７Ｄ、放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３Ａ、大気特性抽出部１４、入射照度算出部１５および地表アルベド算出部１８の各機能および各構成は、実施の形態１および実施の形態２で説明した通りである。
　地表反射率算出部１６Ａは、実施の形態１と同様の手順で、マルチスペクトル画像ＩＢｉの放射輝度Ｌｓｅｎｓｏｒ，ｉ、大気散乱放射輝度Ｌｓｃａｔｔ，ｉ、大気透過率τｉおよび入射照度Ｅｉを用いて、被撮影領域の地表反射率ρｉを算出する。補正部１９は、地表と大気との間での光の相互作用に応じた地表反射率ρｉの変化分を補正する。
　図１２は、実施の形態３に係る画像処理方法を示すフローチャートであって、画像処理装置１Ｂによる画像処理の手順を示している。図１２のステップＳＴ１ｂからステップＳＴ７ｂまでの処理は、図１０のステップＳＴ１ａからステップＳＴ７ａまでと同様の処理であるので、これらのステップの説明は省略する。
　補正部１９は、地表と大気との間における太陽光の相互作用に応じた地表反射率ρｉの変化分を補正する（ステップＳＴ８ｂ）。図１３は、地表と大気との間での太陽光の相互作用を示す図である。地表アルベドデータベース１７Ｃにおける地表アルベドが得られた環境と観測画像ＭＳＩが撮影された観測環境との差異に応じて大気散乱光の大きさは変化する。
　例えば、地表反射率ρｉの算出に使用された観測画像ＭＳＩの観測環境と、地表アルベドデータベース１７Ｃにおける地表アルベドの算出に使用された観測画像ＭＳＩの観測環境との間で、光学センサの視線方向または太陽高度に差異があると、大気散乱光の大きさまたは被写体のＢＲＤＦ（双方向散乱分布関数）が変化し得る。これにより、大気から地球への散乱光が変化し、地表を照らす照度が変化し得る。従って、入射照度Ｅｉは、太陽光の直達光成分とスカイライト成分の組み合わせによってその分光特性が変化する。補正部１９は、地表反射率ρｉを、地表へ入射されるスカイライト成分の変化の影響が低減された本来の値に近い地表反射率ρｉ’に補正する。
　補正部１９には、例えば、地表アルベドデータベース１７Ｃに格納された地表アルベドが閾値ρｔｈとして設定されている。補正部１９は、地表反射率ρｉが閾値ρｔｈとの間で下記式（６）の関係にある場合、地表と大気との間で太陽光のスカイライト成分が相互作用して地表へ入射されるスカイライト成分が変化したと推定する。
ρｉ≦ρｔｈ　　　　　　　　　　　　　　・・・（６）
　地表と大気との間での相互作用によるスカイライト成分の変化を、元のスカイライト成分からの割合として係数ｒｉとした場合に、スカイライトの照度成分の変化分を含む入射照度Ｅｉ’は、下記式（７）で表される。入射照度Ｅｉ’を用いて算出される地表反射率ρｉ’は、地表反射率ρｉおよびその算出に使用された入射照度Ｅｉを用いることで、下記式（８）に示す関係で表すことができる。補正部１９は、地表へ入射されるスカイライト成分が変化したと推定すると、下記式（７）に従って、地表反射率ρｉにおけるスカイライト成分に応じた変化分を低減させた地表反射率ρｉ’に補正する。これは、例えば、太陽高度が低く、建物などの影が多いような画像を想定している。逆にρｉがρｔｈ以上である場合に、スカイライトの成分がより大きくなる（ｒｉ＞１）ような補正を行ってもよい。これは、例えば、順光によって正反射成分が強く作用するような画像を想定している。
Ｅｉ’＝Ｅｄｉｒｅｃｔ，ｉ＋ｒｉ×Ｅｓｋｙ，ｉ　　　・・・（７）
ρｉ’＝ρｉ（Ｅｉ／Ｅｉ’）　　　　　　　　　　　　・・・（８）
　なお、画像処理装置１Ｂにおける放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３Ａ、大気特性抽出部１４、入射照度算出部１５、地表反射率算出部１６Ａおよび地表アルベド算出部１８の各機能は、処理回路によって実現される。すなわち、画像処理装置１Ｂは、図１２に示したステップＳＴ１ｂからステップＳＴ８ｂまでの処理を実行するための処理回路を備える。処理回路は、図８Ａに示した専用のハードウェアの処理回路１０２であってもよいし、図８Ｂに示したメモリ１０４に記憶されたプログラムを実行するプロセッサ１０３であってもよい。
　以上のように、実施の形態３に係る画像処理装置１Ｂは、地表と大気との間での光の相互作用に応じた地表反射率ρｉの変化分を補正する補正部１９を備えている。地表と大気との間での相互作用によって地表へ入射されるスカイライト成分が変化しても、この影響を低減させた地表反射率ρｉ’を得ることができる。
実施の形態４．
　図１４は、実施の形態４に係る画像処理装置１Ｃの構成を示す機能ブロック図である。図１４に示すように、画像処理装置１Ｃは、スカイライトデータベース１７Ｂおよびエアロゾルモデルデータベース１７Ｅが格納されたデータ記憶部１７を備えている。さらに、画像処理装置１Ｃは、放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３Ｂ、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６を備えている。エアロゾル・視程推定部１３Ｂは、変換部１３Ｃを備える。なお、スカイライトデータベース１７Ｂ、放射輝度変換部１１、散乱特性検出部１２、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６の各機能および各構成は、実施の形態１で説明した通りである。
　本発明の発明者による実験解析の結果、太陽天頂角θが小さい範囲では、太陽天頂角θが変化しても大気散乱放射輝度の絶対値の変化はあるものの、分光特性の変化は小さいことが判明した。すなわち、太陽天頂角θが小さい範囲では、太陽天頂角θに対する絶対値の補正係数を与えれば、エアロゾルモデルに設定されている大気散乱放射輝度を、太陽天頂角θが０度であるときの大気散乱放射輝度に置き換えても、置き換え前の大気散乱放射輝度との対応関係と大きく乖離しないことを意味する。なお、この特性は、光学センサを搭載する人工衛星の天頂角θｓａｔが約４０度であるときには、太陽天頂角θが約５０度以下でこの関係を満たす。
　実施の形態１で示したエアロゾルモデルデータベース１７Ａにおけるエアロゾルモデルには、大気散乱放射輝度および大気透過率の分光特性データが太陽天頂角θごとに設定されている。これに対して、エアロゾルモデルデータベース１７Ｅでは、太陽天頂角θが小さい範囲において、太陽天頂角θが０度であるときに相当する大気散乱放射輝度に分光特性データを対応付ければよい。このため、エアロゾルモデルデータベース１７Ｅは、エアロゾルモデルデータベース１７Ａに比べて、太陽天頂角θが小さい範囲におけるデータの量が大幅に削減される。
　エアロゾル・視程推定部１３Ｂには、例えば、太陽天頂角θに対する大気散乱放射輝度値の変化が、ルックアップテーブル（以下、ＬＵＴと記載する）として設定されている。太陽天頂角θが閾値以下である場合、変換部１３Ｃは、このＬＵＴを参照して、散乱特性検出部１２によって検出された大気散乱放射輝度を、太陽天頂角θが０度であるときに相当する大気散乱放射輝度に変換する。エアロゾルモデル・視程推定部１３Ｂは、変換部１３Ｃによって変換された大気散乱放射輝度を、エアロゾルモデルおよび視程の推定に使用する。
　図１５は、実施の形態４に係る画像処理方法を示すフローチャートであって、画像処理装置１Ｃによる画像処理の手順を示している。図１５のステップＳＴ１ｃ、ステップＳＴ２ｃおよびステップＳＴ６ｃからステップＳＴ８ｃの処理は、図４におけるステップＳＴ１、ステップＳＴ２およびステップＳＴ４からステップＳＴ６までと同様の処理であるので、これらのステップの説明は省略する。
　ステップＳＴ２ｃの処理が終了すると、エアロゾル・視程推定部１３Ｂは、散乱特性検出部１２によって大気散乱放射輝度が検出された画像が観測されたときの太陽天頂角θが閾値以下であるか否かを判定する（ステップＳＴ３ｃ）。なお、太陽天頂角θは、観測画像ＭＳＩに付随している既知データである。また、閾値は、例えば６０度とする。
　太陽天頂角θが閾値以下であった場合（ステップＳＴ３ｃ；ＹＥＳ）、変換部１３Ｃは、ＬＵＴを参照して、散乱特性検出部１２によって検出された大気散乱放射輝度を、太陽天頂角θが０度であるときに相当する大気散乱放射輝度に変換する（ステップＳＴ４ｃ）。この後、エアロゾル・視程推定部１３Ｂは、エアロゾルモデルデータベース１７Ａから、変換部１３Ｃによって変換された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定する（ステップＳＴ５ｃ）。
　一方、太陽天頂角θが閾値よりも大きかった場合（ステップＳＴ２ｃ；ＮＯ）、エアロゾル・視程推定部１３Ｂは、エアロゾルモデルデータベース１７Ａから、散乱特性検出部１２によってマルチスペクトル画像ＩＢｉから検出された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定する（ステップＳＴ５ｃ）。
　なお、画像処理装置１Ｃにおける放射輝度変換部１１、散乱特性検出部１２、エアロゾル・視程推定部１３Ｂ、大気特性抽出部１４、入射照度算出部１５および地表反射率算出部１６の各機能は、処理回路により実現される。すなわち、画像処理装置１Ｃは、図１５に示したステップＳＴ１ｃからステップＳＴ８ｃまでの処理を実行するための処理回路を備えている。処理回路は、図８Ａに示した専用のハードウェアの処理回路１０２であってもよいし、図８Ｂに示したメモリ１０４に記憶されたプログラムを実行するプロセッサ１０３であってもよい。
　以上のように、実施の形態４に係る画像処理装置１Ｃにおいて、エアロゾル・視程推定部１３Ｂが、太陽天頂角が閾値以下である場合、マルチスペクトル画像から検出された大気散乱放射輝度が太陽天頂角が０度であるときに相当する大気散乱放射輝度に変換された値を、エアロゾルモデルおよび視程の推定に使用する。エアロゾルモデルデータベース１７Ｅが、太陽天頂角θが小さい範囲において、太陽天頂角θが０度であるときに相当する大気散乱放射輝度に分光特性データを対応付ければよい。このため、エアロゾルモデルデータベース１７Ｅにおける、太陽天頂角θが小さい範囲に対応するデータの量を大幅に削減することが可能である。
　なお、本発明は上記実施の形態に限定されるものではなく、本発明の範囲内において、実施の形態のそれぞれの自由な組み合わせまたは実施の形態のそれぞれの任意の構成要素の変形もしくは実施の形態のそれぞれにおいて任意の構成要素の省略が可能である。
　本発明に係る画像処理装置は、例えば、人工衛星または航空機に搭載された光学センサによって観測された観測画像の画像処理に利用可能である。
　１，１Ａ，１Ｂ，１Ｃ　画像処理装置、１０　人工衛星、１１　放射輝度変換部、１２　散乱特性検出部、１３，１３Ａ，１３Ｂ　エアロゾル・視程推定部、１３Ｃ　変換部、１４　大気特性抽出部、１５　入射照度算出部、１６，１６Ａ　地表反射率算出部、１７　データ記憶部、１７Ａ，１７Ｄ，１７Ｅ　エアロゾルモデルデータベース、１７Ｂ　スカイライトデータベース、１７Ｃ　地表アルベドデータベース、１８　地表アルベド算出部、１９　補正部、１００　入出力インタフェース、１０１　記憶装置、１０２　処理回路、１０３　プロセッサ、１０４　メモリ。

Claims (8)

  1. 　地表の被撮影領域を観測して複数の観測波長帯のそれぞれに対応する複数のマルチスペクトル画像を生成する光学センサから、複数の前記マルチスペクトル画像を入力し、入力した複数の前記マルチスペクトル画像の各画素の画素値を放射輝度に変換する放射輝度変換部と、
    　複数の前記マルチスペクトル画像の放射輝度値に基づいて、大気散乱光の放射輝度成分である大気散乱放射輝度を前記マルチスペクトル画像から検出する散乱特性検出部と、
    　エアロゾルに対応する視程と大気散乱放射輝度および大気透過率の分光特性データとが対応付けて設定された複数のエアロゾルモデルから、前記マルチスペクトル画像から検出された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定するエアロゾル・視程推定部と、
    　前記エアロゾル・視程推定部によって推定されたエアロゾルモデルから、推定された視程に対応する大気散乱放射輝度および大気透過率のデータを抽出する大気特性抽出部と、
    　前記大気特性抽出部によって抽出されたデータを用いて、前記被撮影領域の入射照度を算出する入射照度算出部と、
    　前記マルチスペクトル画像の放射輝度値、前記大気特性抽出部によって抽出されたデータ、および前記入射照度算出部によって算出された入射照度を用いて、前記被撮影領域における地表面での分光反射率を算出する地表反射率算出部と、
    　を備えたことを特徴とする画像処理装置。
  2. 　前記被撮影領域における地表面のアルベドを算出する地表アルベド算出部を備え、
    　複数のエアロゾルモデルは、前記アルベドに対応付けられており、
    　前記エアロゾル・視程推定部は、複数のエアロゾルモデルのうち、前記地表アルベド算出部によって算出された前記アルベドおよび前記マルチスペクトル画像から検出された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定すること
    　を特徴とする請求項１記載の画像処理装置。
  3. 　地表と大気との間での光の相互作用に応じた前記分光反射率の変化分を補正する補正部を備えたこと
    　を特徴とする請求項２記載の画像処理装置。
  4. 　前記エアロゾル・視程推定部は、太陽天頂角が閾値以下である場合、前記マルチスペクトル画像から検出された大気散乱放射輝度が前記太陽天頂角が０度であるときに相当する大気散乱放射輝度に変換された値を、エアロゾルモデルおよび視程の推定に使用すること
    　を特徴とする請求項１から請求項３のいずれか１項記載の画像処理装置。
  5. 　放射輝度変換部が、地表の被撮影領域を観測して複数の観測波長帯のそれぞれに対応する複数のマルチスペクトル画像を生成する光学センサから、複数の前記マルチスペクトル画像を入力し、入力した複数の前記マルチスペクトル画像の各画素の画素値を放射輝度に変換するステップと、
    　散乱特性検出部が、複数の前記マルチスペクトル画像の放射輝度値に基づいて、大気散乱光の放射輝度成分である大気散乱放射輝度を前記マルチスペクトル画像から検出するステップと、
    　エアロゾル・視程推定部が、エアロゾルに対応する視程と大気散乱放射輝度および大気透過率の分光特性データとが対応付けて設定された複数のエアロゾルモデルから、前記マルチスペクトル画像から検出された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定するステップと、
    　大気特性抽出部が、前記エアロゾル・視程推定部によって推定されたエアロゾルモデルから、推定された視程に対応する大気散乱放射輝度および大気透過率の分光特性データを抽出するステップと、
    　入射照度算出部が、前記大気特性抽出部によって抽出されたデータを用いて、前記被撮影領域の入射照度を算出するステップと、
    　地表反射率算出部が、前記マルチスペクトル画像の放射輝度値、前記大気特性抽出部によって抽出されたデータ、および前記入射照度算出部によって算出された入射照度を用いて、前記被撮影領域における地表面での分光反射率を算出するステップと、
    　を備えたことを特徴とする画像処理方法。
  6. 　地表アルベド算出部が、前記被撮影領域における地表面のアルベドを算出するステップを備え、
    　複数のエアロゾルモデルは、前記アルベドに対応付けられており、
    　前記エアロゾル・視程推定部は、複数のエアロゾルモデルのうち、前記地表アルベド算出部によって算出された前記アルベドおよび前記マルチスペクトル画像から検出された大気散乱放射輝度に対応するエアロゾルモデルおよび視程を推定すること
    　を特徴とする請求項５記載の画像処理方法。
  7. 　補正部が、地表と大気との間での光の相互作用に応じた前記分光反射率の変化分を補正するステップを備えたこと
    　を特徴とする請求項６記載の画像処理方法。
  8. 　前記エアロゾル・視程推定部が、太陽天頂角が閾値以下である場合、前記マルチスペクトル画像から検出された大気散乱放射輝度が前記太陽天頂角が０度であるときに相当する大気散乱放射輝度に変換された値を、エアロゾルモデルおよび視程の推定に使用するステップを備えたこと
    　を特徴とする請求項５から請求項７のいずれか１項記載の画像処理方法。
PCT/JP2019/032984 2019-08-23 2019-08-23 画像処理装置および画像処理方法 WO2021038621A1 (ja)

Priority Applications (2)

Application Number Priority Date Filing Date Title
PCT/JP2019/032984 WO2021038621A1 (ja) 2019-08-23 2019-08-23 画像処理装置および画像処理方法
JP2021541757A JP6964834B2 (ja) 2019-08-23 2019-08-23 画像処理装置および画像処理方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/JP2019/032984 WO2021038621A1 (ja) 2019-08-23 2019-08-23 画像処理装置および画像処理方法

Publications (1)

Publication Number Publication Date
WO2021038621A1 true WO2021038621A1 (ja) 2021-03-04

Family

ID=74685002

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2019/032984 WO2021038621A1 (ja) 2019-08-23 2019-08-23 画像処理装置および画像処理方法

Country Status (2)

Country Link
JP (1) JP6964834B2 (ja)
WO (1) WO2021038621A1 (ja)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113092368A (zh) * 2021-03-16 2021-07-09 上海机电工程研究所 一种基于无人机的红外波段大气透过率测量方法及系统
CN113218874A (zh) * 2021-04-13 2021-08-06 中国科学院合肥物质科学研究院 一种基于遥感影像获取地表目标物反射率方法及系统
CN113533262A (zh) * 2021-03-24 2021-10-22 北京航空航天大学 一种大气气溶胶红外散射透过率确定方法
WO2024087464A1 (zh) * 2022-10-24 2024-05-02 河北先河环保科技股份有限公司 离水辐射光谱计算方法、装置、终端及存储介质

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100008595A1 (en) * 2008-07-08 2010-01-14 Harris Corporation Automated atmospheric characterization of remotely sensed multi-spectral imagery
JP2016126566A (ja) * 2015-01-05 2016-07-11 三菱電機株式会社 画像処理装置及び画像処理方法
JP2017068456A (ja) * 2015-09-29 2017-04-06 三菱電機株式会社 画像処理装置及び画像処理方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20100008595A1 (en) * 2008-07-08 2010-01-14 Harris Corporation Automated atmospheric characterization of remotely sensed multi-spectral imagery
JP2016126566A (ja) * 2015-01-05 2016-07-11 三菱電機株式会社 画像処理装置及び画像処理方法
JP2017068456A (ja) * 2015-09-29 2017-04-06 三菱電機株式会社 画像処理装置及び画像処理方法

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113092368A (zh) * 2021-03-16 2021-07-09 上海机电工程研究所 一种基于无人机的红外波段大气透过率测量方法及系统
CN113533262A (zh) * 2021-03-24 2021-10-22 北京航空航天大学 一种大气气溶胶红外散射透过率确定方法
CN113533262B (zh) * 2021-03-24 2022-11-04 北京航空航天大学 一种大气气溶胶红外散射透过率确定方法
CN113218874A (zh) * 2021-04-13 2021-08-06 中国科学院合肥物质科学研究院 一种基于遥感影像获取地表目标物反射率方法及系统
WO2024087464A1 (zh) * 2022-10-24 2024-05-02 河北先河环保科技股份有限公司 离水辐射光谱计算方法、装置、终端及存储介质

Also Published As

Publication number Publication date
JPWO2021038621A1 (ja) 2021-10-21
JP6964834B2 (ja) 2021-11-10

Similar Documents

Publication Publication Date Title
JP6964834B2 (ja) 画像処理装置および画像処理方法
US11461994B2 (en) Methods for in-scene shadow compensation using sunlit and skylit illumination factors
US10832390B2 (en) Atmospheric compensation in satellite imagery
US8073279B2 (en) Automated atmospheric characterization of remotely sensed multi-spectral imagery
US9530188B2 (en) Image processing method, image processing system, and image processing program
US11039076B2 (en) Information processing apparatus, information processing method, and storage medium
JP2015032205A (ja) 画像処理装置及び画像処理方法
US11164297B2 (en) Image processing device, image processing method, and computer-readable recording medium for enhancing quality of an image after correction
Schläpfer et al. Correction of shadowing in imaging spectroscopy data by quantification of the proportion of diffuse illumination
JP5921311B2 (ja) 画像処理装置及び画像処理方法
JP2016126566A (ja) 画像処理装置及び画像処理方法
CA3077923C (en) Methods for in-scene compensation using water vapor content
JP6747436B2 (ja) 画像処理装置、画像処理システム、画像処理方法およびコンピュータプログラム
JP6856066B2 (ja) 情報処理装置、情報処理システム、情報処理方法およびコンピュータプログラム
JP6556409B2 (ja) 画像処理装置
Ma et al. Adjacency effect estimation by ground spectra measurement and satellite optical sensor synchronous observation data
Yang et al. Radiometric calibration for MWIR cameras
JPWO2018016555A1 (ja) 画像処理装置および画像処理方法

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

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2021541757

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 19943675

Country of ref document: EP

Kind code of ref document: A1