CN104121884B - Satellite image pixel view zenith angle and azimuthal computational methods - Google Patents
Satellite image pixel view zenith angle and azimuthal computational methods Download PDFInfo
- Publication number
- CN104121884B CN104121884B CN201410352963.8A CN201410352963A CN104121884B CN 104121884 B CN104121884 B CN 104121884B CN 201410352963 A CN201410352963 A CN 201410352963A CN 104121884 B CN104121884 B CN 104121884B
- Authority
- CN
- China
- Prior art keywords
- azimuth
- observation
- angle
- image
- point
- Prior art date
- Legal status (The legal status 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 status listed.)
- Expired - Fee Related
Links
- 238000000205 computational method Methods 0.000 title description 3
- 238000003384 imaging method Methods 0.000 claims abstract description 18
- 238000000034 method Methods 0.000 claims abstract description 6
- 238000004364 calculation method Methods 0.000 claims description 4
- 230000009466 transformation Effects 0.000 abstract description 4
- 238000012937 correction Methods 0.000 abstract description 3
- 238000012892 rational function Methods 0.000 abstract description 3
- 238000000844 transformation Methods 0.000 abstract description 2
- 238000013519 translation Methods 0.000 abstract description 2
- 238000010586 diagram Methods 0.000 description 4
- 230000000007 visual effect Effects 0.000 description 4
- 230000000694 effects Effects 0.000 description 3
- 230000002457 bidirectional effect Effects 0.000 description 2
- 238000005315 distribution function Methods 0.000 description 2
- 238000013507 mapping Methods 0.000 description 2
- 230000005855 radiation Effects 0.000 description 2
- 230000005540 biological transmission Effects 0.000 description 1
- 230000007123 defense Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C1/00—Measuring angles
Landscapes
- Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Image Analysis (AREA)
- Image Processing (AREA)
Abstract
It is a kind of that pixel observation zenith angle and azimuthal method are calculated according to satellite image geometry imaging model：The imaging model for only needing to know satellite image (can be rigorous geometry model, can also be the general three-dimensional imaging model of rational function model etc), view zenith angle and the azimuth of each pixel in image can just be calculated, when the geometric transformations such as geometric correction, re-projection are carried out to image, the observation angle of each pixel can together be entered line translation, therefore the observation angle information of image will not be caused to lose because of geometric transformation.
Description
Technical field
The present invention relates to satellite-remote-sensing image view zenith angle and azimuthal calculating, can apply in agricultural, forestry, gas
As remote sensing departments such as, ecological environment and national defense and military.
Background technology
Observation angle (zenith angle and azimuth) information has important effect for the treatment and application of remotely-sensed data.It is first
First, the difference of moonscope angle can cause the difference of air path in radiation transmission, therefore observation angle is typically remote sensing shadow
As |input paramete critically important in atmospheric correction treatment.Secondly, under the conditions of given incidence, the reflection received by sensor
Radiation intensity is also influenceed by observation angle, therefore generally utilizes bidirectional reflectance distribution function (Bidirectional
Reflectance Distribution Function, BRDF) describe how light is reflected in body surface.Additionally,
When with hypsography, observation angle can also produce large effect to the geometrical property of image.
For the less satellite image of the angle of visual field, the change of observation angle is smaller in a scape image capturing range, metadata provider
Generally a view zenith angle and azimuth (such as Landsat-5) are provided to a scape image.The image larger for the angle of visual field, then
Need to provide more observation angles (such as MODIS, HJ-1A/B) for different regions.In fact, if necessary to accurate land productivity
With the observation angle information of satellite, even for the less satellite image of Landsat-5 etc angles of visual field, a scape gives one group of sight
Measuring angle is also inadequate, especially in high-latitude area.(such as GF-1 wide cuts shadow in addition, many larger remote sensing images of the angle of visual field
Picture) one group of observation angle only also is given to a scape image, it is thus more inaccurate.
In view of important effect of observation angle (the zenith angle and azimuth) information for treatment and the application of remotely-sensed data,
And observation angle information that existing satellite image product is provided is generally imperfect, therefore research satellite image each pixel is seen
The computational methods of measuring angle have important practical value.
The content of the invention
It is an object of the invention to solve the deficiencies in the prior art, it is proposed that one kind is according to satellite image geometry imaging model
Calculate pixel observation zenith angle and azimuthal method.Only need to know the imaging model of satellite image, it is possible to calculate shadow
The view zenith angle of each pixel and azimuth as in, when the geometric transformations such as geometric correction, re-projection are carried out to image, can
The observation angle of each pixel is entered into line translation together, therefore the observation angle information of image will not be caused to lose because of geometric transformation
Lose.
Technical scheme is as follows：
A kind of pixel observation zenith angle and azimuth calculation method based on satellite image geometry imaging model, methods described
Middle geometry imaging model can be the rigorous geometry model, or general three-dimensional imaging model described with ephemeris parameter
(such as rational function model)；Comprise the following steps,
Step 1, recovers observation sight line：For any point p on image (pixel coordinate be (x, y)), take two it is different
Height value h1And h2(h1＜ h2), obtain the corresponding two spaces point P of p points using the imaging model of satellite image1(φ1, λ1, h1)
And P2(φ2, λ2, h2), recover the observation sight line of image plane point p
Step 2, calculating observation zenith angle：Connection the earth's core O and P1Point, calculates OP1With P2P1Angle, as image plane point p
View zenith angle θzenith；
Step 3, calculating observation azimuth initial value：Calculate P2P1With P excessively1The angle of the meridian plane of point, as image plane point p
The initial value θ of observed azimuthazimuth；
Step 4, north-south amendment is carried out to observed azimuth：If φ2＜ φ1, then θ is takenazimuth=180 ° of-θazimuth；
Otherwise, step 5 is gone to；
Step 5, East and West direction amendment is carried out to observed azimuth：If λ2＜ λ1, then θ is takenazimuth=-θazimuth；Otherwise,
Terminate to calculate.
Brief description of the drawings
Fig. 1 is the flow chart of the embodiment of the present invention one；
Fig. 2 is moonscope zenith angle schematic diagram；
Fig. 3 is OP '1P2Place great circle floor map；
Fig. 4 is satellite observation direction angle schematic diagram；
Fig. 5 is plane O ' IP2Schematic diagram；
Fig. 6 is meridian plane OIP '1Schematic diagram.
Specific embodiment
The present invention proposes a kind of according to satellite image geometry imaging model calculating pixel observation zenith angle and azimuthal
Method.Technical solution of the present invention is described in detail below in conjunction with one embodiment and accompanying drawing 1~6.
Example is provided and calculates pixel observation zenith angle and azimuthal method according to satellite image geometry imaging model, referring to
Accompanying drawing 1：
Step 1, recovers observation sight line：For any point p on image (pixel coordinate be (x, y)), take two it is different
Height value h1And h2(h1＜ h2), obtain the corresponding two spaces point P of p points using the positive calculation model of satellite image1(φ1, λ1, h1)
And P2(φ2, λ2, h2), recover the observation sight line of image plane point pReference can be made to accompanying drawing 2；
The imaging model of satellite image is used for describing pixel coordinate (ranks number) and geodesic latitude and longitude coordinates (longitude, latitude
And elevation) between mapping relations, generally represented from pixel coordinate and given elevation to geodesic latitude and longitude coordinates with positive calculation model
Mapping, such as shown in (1),
(φ, λ)=F (x, y, h) (1)
X represents the row coordinate of image plane point in formula (1), and y represents the row coordinate of image plane point, and φ represents latitude, and λ represents warp
Degree, h represents elevation, and F () represents positive and calculates model.
Imaging model is not limited to certain specific model, both can be the rigorous geometry model described with ephemeris parameter,
It can be general three-dimensional imaging model (such as rational function model).
Here h1Can value be 0 meter, h2Can be taken as 1000 meters.
Step 2, calculating observation zenith angle：Connection the earth's core O and P1Point, calculates OP1With P2P1Angle, as image plane point p
View zenith angle θzenith；
Extension OP1 hand over P2 points where elevation face in P ' 1, then θ zenith=∠ P2P1P ' 1, the latitude and longitude coordinates that 1 point of P ' are P ' 1
(φ 1, λ 1, h2), for the sake of directly perceived, fall into a trap in the big disk where OP ' 1P2 calculate ∠ P2P1P ' below, referring to accompanying drawing 3.Cross P2 points
Make the vertical OP ' 1 of P2M in M, make ∠ P2OP ' 1=γ.By " Haversine formula " [H.B.Goodwin, The
Haversine in nautical astronomy, Naval Institute Proceedings, vol.36, no.3
(1910), pp.735-746：Evidently if a Table of Haversines is employed we shall be
saved in the first instance the trouble of dividing the sum of the logarithms
By two, and in the second place of multiplying the angle taken from the tables
by the same number.This is the special advantage of the form of table first
Introduced by Professor Inman, of the Portsmouth Royal Navy College, nearly a
Century ago.] γ can be calculated obtain
If earth radius is R ≈ 6.4 × 106Rice, then OP1=R+h1, OP2=R+h2, then
Therefore can obtain view zenith angle is：
Step 3, calculating observation azimuth initial value：Calculate P2P1With P excessively1The angle of the meridian plane of point, as image plane point p
The initial value θ of observed azimuthazimuth；
If I points are P2Parallel and P excessively where putting1The intersection point of the meridian plane of point, crosses P2Make plane OIP1Vertical line (intersection point
It is designated as M), then ∠ P2P1M is the observed azimuth (or supplementary angle) of satellite, referring to accompanying drawing 4.From solid geometry knowledge, vertical line
P2M must be in P2In weft plane where point, and P2M ⊥ O ' I, wherein O ' are P2Weft plane where point and earth's axis
Intersection point.Then
P is calculated separately below2M and P1M。
For the sake of directly perceived, first in plane O ' IP2Middle calculating P2The length of M, referring to accompanying drawing 5.
The radius of circle O ' is by P2The latitude and elevation of point determine there are O ' P2=O ' I=(R+h2)cosφ2, and ∠ P2O ' I are
P2Point and P1The difference of longitude of point, i.e. ∠ P2O ' I=| λ2-λ1|, can then obtain
O ' M=(R+h2)cosφ2cos|λ2-λ1|---(5)
P2M=(R+h2)cosφ2sin|λ2-λ1|
Identical with step 2, R represents the radius of the earth here.
Then in meridian plane OIP '1Middle calculating P1The length of M, referring to accompanying drawing 6：
In accompanying drawing 6, OE is the intersection of meridian plane and equatorial plane.Cross P1Point makees vertical line and meets at A with OE and O ' I respectively to OE
And B, then
Wherein, P1B=AP1- AB=AP1- OO '=(R+h1)sinφ1-(R+h2)sinφ2,
BM=O ' M-O ' B=O ' M-OA=(R+h2)cosφ2cos(λ2-λ1)-(R+h1)cosφ1。
P can be utilized below2M and P1The azimuthal initial value of M calculating observations：
Due to R ＞ ＞ h1And R ＞ ＞ h2, it is believed thatThen (7) formula can be reduced to：
Step 4, north-south amendment is carried out to observed azimuth：If φ2＜ φ1, then θ is takenazimuth=180 ° of-θazimuth；
Otherwise, step 5 is gone to；
Step 5, East and West direction amendment is carried out to observed azimuth：If λ2＜ λ1, then θ is takenazimuth=-θazimuth；Otherwise,
Terminate to calculate.
Claims (1)
1. a kind of pixel observation zenith angle and azimuth calculation method based on satellite image geometry imaging model, in methods described
Geometry imaging model includes the rigorous geometry model and general three-dimensional imaging model that are described with ephemeris parameter；It is characterized in that：Bag
Include following steps,
Step 1, recovers observation sight line：For any one pixel p (x, y) on image, two different height value h are taken1=0 He
h2, and cause h1＜ h2, the corresponding two spaces point P of p points is obtained using the imaging model of satellite image1(φ1, λ1, h1) and P2
(φ2, λ2, h2), recover the observation sight line of image plane point p
Step 2, calculating observation zenith angle：Connection the earth's core O and P1Point, calculates OP1With P2P1Angle, the as sight of image plane point p
Observation vertex angle thetazenith；
Step 3, calculating observation azimuth initial value：Calculate P2P1With P excessively1The angle of the meridian plane of point, as image plane point p observations
Azimuthal initial value θazimuth；
Step 4, north-south amendment is carried out to observed azimuth：If φ2＜ φ1, then θ is takenazimuth=180 ° of-θazimuth；Otherwise,
Go to step 5；
Step 5, East and West direction amendment is carried out to observed azimuth：If λ2＜ λ1, then θ is takenazimuth=-θazimuth；Otherwise, meter is terminated
Calculate.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410352963.8A CN104121884B (en) | 2014-07-24 | 2014-07-24 | Satellite image pixel view zenith angle and azimuthal computational methods |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201410352963.8A CN104121884B (en) | 2014-07-24 | 2014-07-24 | Satellite image pixel view zenith angle and azimuthal computational methods |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104121884A CN104121884A (en) | 2014-10-29 |
| CN104121884B true CN104121884B (en) | 2017-06-06 |
Family
ID=51767398
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201410352963.8A Expired - Fee Related CN104121884B (en) | 2014-07-24 | 2014-07-24 | Satellite image pixel view zenith angle and azimuthal computational methods |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104121884B (en) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN104951656B (en) * | 2015-06-23 | 2018-01-12 | 中国科学院遥感与数字地球研究所 | Wide ken satellite image Reflectivity for Growing Season inversion method |
| CN107389617A (en) * | 2017-06-27 | 2017-11-24 | 环境保护部卫星环境应用中心 | The inversion method and equipment of aerosol optical depth based on No. four satellites of high score |
| CN109932341B (en) * | 2019-03-11 | 2021-03-23 | 北京环境特性研究所 | Bidirectional reflection distribution function measuring method of typical target in field environment |
| CN111060078A (en) * | 2019-12-20 | 2020-04-24 | 彭耿 | Positioning method based on satellite observation angle error estimation |
Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102338871A (en) * | 2010-07-22 | 2012-02-01 | 曹春香 | Method and device for calculating reflectivity of earth surface |
-
2014
- 2014-07-24 CN CN201410352963.8A patent/CN104121884B/en not_active Expired - Fee Related
Patent Citations (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102338871A (en) * | 2010-07-22 | 2012-02-01 | 曹春香 | Method and device for calculating reflectivity of earth surface |
Non-Patent Citations (2)
| Title |
|---|
| 大幅面可见光遥感图像典型目标识别关键技术研究;韩现伟;《中国博士学位论文全文数据库信息科技辑》;20140115;全文 * |
| 缺少控制点的高分辨率卫星遥感影像几何纠正;张过;《中国优秀博硕士学位论文全文数据库（博士）基础科学辑》;20060515;全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104121884A (en) | 2014-10-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US11635481B2 (en) | Stellar atmospheric refraction measurement correction method based on collinearity of refraction surfaces | |
| CN104535979B (en) | A kind of remote sensing inversion method and system of land cloud optical thickness | |
| CN102928861B (en) | Target positioning method and device for airborne equipment | |
| CN103345737B (en) | A kind of UAV high resolution image geometric correction method based on error compensation | |
| US10670689B2 (en) | System and method for determining geo location of a target using a cone coordinate system | |
| US10467726B2 (en) | Post capture imagery processing and deployment systems | |
| CN104913780B (en) | The high-precision deviation of plumb line method for fast measuring of integrated GNSS and CCD zenith telescopes | |
| CN104121884B (en) | Satellite image pixel view zenith angle and azimuthal computational methods | |
| CN106767714A (en) | Improve the equivalent mismatch model multistage Calibration Method of satellite image positioning precision | |
| CN102778675B (en) | Atmospheric correction method and atmospheric correction module for satellite remote-sensing image | |
| CN104180808A (en) | Aerial autonomous refueling circular taper sleeve vision position and attitude resolving method | |
| CN102589544B (en) | Three-dimensional attitude acquisition method based on space characteristics of atmospheric polarization mode | |
| CN107490364A (en) | A kind of wide-angle tilt is imaged aerial camera object positioning method | |
| CN105444778B (en) | A kind of star sensor based on imaging geometry inverting is in-orbit to determine appearance error acquisition methods | |
| CN102654576A (en) | Image registration method based on synthetic aperture radar (SAR) image and digital elevation model (DEM) data | |
| CN104406583B (en) | Combined defining method for carrier attitude of double-star sensor | |
| CN105527620B (en) | The automatic calibration method that a kind of aerosol thickness postpones with laser radar range | |
| CN110487266A (en) | A kind of airborne photoelectric passive high-precision localization method suitable for sea-surface target | |
| CN106407560A (en) | A building method for a troposphere mapping function model representing atmospheric anisotropy | |
| CN104239678B (en) | A kind of method and apparatus for realizing interferometer direction finding positioning | |
| CN108733711B (en) | Distribution line space distance obtaining method based on three-dimensional GIS technology | |
| CN104613956A (en) | Atmospheric polarization neutral point-based navigation orientation method | |
| CN110081905A (en) | A kind of light wave Atmospheric Refraction Error calculation method based on single station electro-optic theodolite | |
| CN103791889A (en) | Cross structure light assisted monocular vision pose measurement method | |
| KR20230040921A (en) | Atmospheric refraction positioning error correction method for optical remote sensing satellite image in qinghai-tibet plateau region |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20170606 |