2006189

SAR를 이용한 파랑분석 원천기술 개발

강 원 대 학 교

한 국 해 양 연 구 원

□□□□□□□□□□□□□□□□□□ □ 주 의 □ □ (편집순서8) □ □ □□□□□□□□□□□□□□□□□□ (16 포인트 고딕체) ↑ 7cm ↓	2 0 0 6 1 8 9 S A R 를 이 용 한 파 랑 분 석 원 천 기 술 개 발 강 원 대 학 교	표지

주        의

1. 이 보고서는 한국해양연구원에서 수행한 기본연구사업의 연구결과보고서 입니다.

2. 이 보고서 내용을 발표할 때에는 반드시 한국해양연구원에서 수행한 기본연구사업의 연구결과임을 밝혀야 합니다.

3. 국가과학기술 기밀유지에 필요한 내용은 대외적으로 발표 또는 공개하여서는 안됩니다.

제   출   문

한국해양연구원장 귀하

본 보고서를 “SAR를 이용한 파랑분석 원천기술 개발”과제의 최종보고서로 제출합니다.

2006.   11.

위탁연구기관명 : 강원대학교

위탁연구책임자 : 이 훈 열

참 여 연 구 원 : 강 문 경

박 용 욱

한 승 하

김 준 수

요   약   문

Ⅰ. 개     요

위탁 과제명: SAR를 이용한 파랑분석 원천기술 개발

본 과제명: 기상재난파에 의한 해양구조물의 안전성 평가기술 개발

연구 기간: 2006년 2월 1일 - 2006년 11월 30일

위탁 연구비: 20,000,000원

Ⅱ. 연구개발의 목적 및 필요성

최근 위성원격탐사는 해양학 분야에서 경제적이고 필수적인 기술로 부각되고 있다. 특히 해양과 관련된 원격탐사 활용분야에서는 대기나 기상상태에 따른 제약점이 없고, 낮/밤의 영향에 무관하게 수시로 변화하는 해상의 상태를 고해상도의 영상 자료로 제공할 수 있는 SAR 영상의 활용이 활발히 이루어지고 있다. 본 연구에서는 해양의 해풍(wind), 파랑(wave), 해류(current)와 같은 다양한 해양 현상에 대하여 인공위성 SAR 영상으로부터 풍속의 크기와 방향, 해파의 파장과 진행방향, 유속의 크기와 방향과 같은 해양물리학적 정보를 추출할 수 있는 종합적인 프로세서의 개발과 적용을 목표로 하고 있다.

Ⅲ. 연구개발의 내용 및 범위

이 연구는 해양의 해풍, 파랑, 해류 관련 파라미터를 추출하기 위해 일반적으로 적용되고 있는 주요 이론 및 연구 기법들을 토대로 종합적 프로세서인 SOP (SAR Ocean Processor)을 개발하여 SAR 자료를 이용하여 자료처리를 수행하고, 추출된 결과를 기상자료나 해양자료와 같은 현장자료와 비교 검토하였다.

1. SAR 영상 자료로부터 후방산란계수(σ0)를 계산하고, 해상풍의 정보를 추출하기 위해 CMOD4와 CMOD_IFR2 경험 모델과 편광 비(polarization ratio, PR)를 이용한 경험식을 적용하여 풍속을 추출하였다.

2. Wave-SAR 역산 방법을 통해 파랑의 스펙트럼을 추출하고, 파장과 진행방향을 산출하였으며, 이때 발생되는 방향에 대한 180° 모호성 문제를 해결하기 위해 inter-look cross-spectra 기법을 적용하였다.

3. 해류의 속도 정보는 단일 산란체에 대하여 SAR 영상이 얻어지는 시간 동안에 해상의 산란체가 움직이면서 발생되는 Doppler shift를 이용한 기법으로 산출하였다.

4. 개발된 SOP 프로세서를 이용하여 자료처리 결과 SAR 영상으로부터 얻어진 해양학적 파라미터를 기상자료나 해양자료와 같은 현장자료와 비교 검토하였다.

Ⅳ. 연구개발결과

본 연구 목표인 해양물리학적 파라미터를 추출하기 위한 종합적인 프로세서 개발 및 적용 결과는 다음과 같다.

1. 개발된 SOP 프로그램은 Windows 환경에서 Linux 환경을 모사하는 Cygwin이라는 프로그램에서 컴파일되고 수행되며, ANSI-C 언어로 이루어져 있어 사용자의 편의에 따라 다양한 컴파일러를 사용할 수 있다. 명령어 사용 환경에서 운영되기 때문에 사용자는 먼저 프로그램에 필요한 간단한 변수와 실행 제어 변수를 입력하여 파일로 만들고, 이를 통해 프로그램을 제어한다.

2. 입력파일로는 원시자료(raw data)로부터 SAR focusing 과정을 거친 CEOS 포멧의 SLC (single look complex) 자료를 사용한다. SOP 프로세서 수행결과 출력되는 파일은 해풍 정보인 CMOD4, CMOD_IFR2 자료, 파랑 정보를 갖고 있는 파랑 스펙트럼(wave spectrum), inter-look cross-spectra 자료이며, 해류의 속도 정보를 산출할 수 있는 Doppler shift 결과 자료이다. 모든 결과는 원격탐사 및 GIS 전용 상용 프로그램으로 손쉽게 열어볼 수 있도록 설계하였다.

3. SOP 프로세서를 한반도 연안 지역에 적용하여 Radarsat-1 영상자료로부터 해풍, 파랑, 해류의 물리적 정보를 추출하였으며, 이를 현장자료와 비교 검토한 결과 신뢰할만한 결과를 얻을 수 있었다.

Ⅴ. 연구개발결과의 활용계획

본 연구에서 개발된 SOP 프로세서를 수행하여 SAR 영상 자료로부터 추출된 해풍, 파랑, 해류와 관련된 물리적 정보는 현장관련 자료와 비교 검토한 결과 신뢰할만한 결과를 보였다. 이를 토대로 추후 해양에서 발생되는 다양한 현상과 관련된 영상자료 및 현장 실측자료들이 확보되고, 정확하고 정밀한 정보 추출을 위한 알고리즘 보완 및 수정 단계를 거치게 되면 더욱 향상된 정보 추출이 가능할 것으로 판단된다. 또한 다양하고 역동적인 해상의 상태를 파악하고 모니터링하며, 이러한 현상과 관련된 기작을 밝히고 연구하는데 활용 가능성이 높을 것으로 판단된다.

목      차

제 1 장 서 론

제 2 장 국내외 기술개발 현황

제 3 장 연구개발수행 내용 및 결과

제 1 절 이론적 배경

1. 해풍(wind)

가. CMOD4 모델

나. CMOD_IFR2 모델

다. 편광 비례식

2. 파랑(wave)

가. SAR Wave 스펙트럼

나. Inter-look cross-spectra 기법

3. 해류(current): Doppler shift 기법

제 2 절 SAR Ocean Processor(SOP) 개발 및 자료처리

1. RADARSAT-1 영상 자료

2. SLC Header Off 영상 추출

3. Multilook 영상 추출

4. 해풍(wind) 정보 추출

5. 파랑(wave) 정보 추출

가. 파랑 스펙트럼(SAR wave spectrum) 추출

나. Inter-look cross-spectra 추출

6. 해류(current) 정보 추출

가. Doppler 영상 추출

나. 속도 자료 추출

제 3 절 SOP 적용 결과 및 현장자료와의 비교

제 4 장 연구개발목표 달성도 및 대외기여도

제 5 장 연구개발결과의 활용계획

제 6 장 참고문헌

부록 A  SAR Ocean Processor 원시 프로그램

제 1 장 서 론

최근 위성원격탐사(satellite remote sensing)는 해양학 분야에서 경제적이고 필수적인 도구로 부각되고 있다. 해양과 관련된 원격탐사 활용분야에서는 대기나 기상 상태에 따른 제약점이 없고, 낮/밤에 무관하게 수시로 변화하는 해상의 상태를 고해상도의 영상 정보로 제공하는 SAR 영상의 활용이 필수적이라 할 수 있다(Curlander and McDonough, 1991).

합성구경레이더(synthetic aperture radar, SAR)는 영상레이더(imaging radar)의 일종으로서, 해양에서 중력이 주요 복원력이 되는 중력파(gravity wave)나 파장이 1.73 ㎝ 미만의 표면장력파(capillary wave) 등 파랑의 분포상태를 2차원적인 영상으로 관측할 수 있는 원격탐사 시스템이다. 해상에 나타나는 표면파(surface waves), 내부파(internal waves), 해류(currents), 바람(wind cells), 해저 지형에 따른 특성(bathymetric features), 선박 항적(ship wakes), 기름 유출(oil spills) 등과 같은 해양학적 현상들은 후방산란계수(backscattering coefficient)나 공간적 스펙트럼 분포(spectral distribution)에 따라 2차원의 SAR 영상자료로 가시화된다(Elachi, 1988; Mouchot and Garello, 1998).

원격탐사에서 사용되는 레이더 시스템의 종류는 크게 영상레이더, 고도계(altimeter), 산란계(scatterometer)로 구별되며, 이 중 영상레이더는 고해상도의 영상자료를 취득할 수 있는 시스템으로 주로 지질학적 구조, 토지이용, 해상표면파, 극빙등과 같은 지표 특성을 영상화하여 관측하는데 이용되고 있다. SAR는 1 ㎝ ～ 1 m 사이의 마이크로파 영역에서 측정되는 원격탐사 시스템으로 일반적으로 대상체나 대상지역의 경사, 지형, 거칠기(roughness) 등과 같은 물리적 성질과 유전상수(dielectric constant), 흡수도(absorption), 전도도(conductivity) 등과 같은 전기적 특성에 따라 영향을 받게 되며, 이를 통하여 대상체나 대상지역의 성질과 특성을 영상에 반영하게 된다.

시․공간적으로 변동하는 해황을 가능한 신속히 파악하고 모니터링 하는 것은 효율적인 어장탐색을 비롯한 합리적․경제적 어업경영과 연계되어 있을 뿐 아니라 지진해일이나 기상재난파(freak wave)와 같이 자연재해와 관련된 해상의 급격한 변화에 따른 상황을 분석하고 모니터링하여 신속한 대처방안을 모색하는데 필수적으로 선결되어야 하는 중요한 과정이라 할 수 있다. 따라서 국내․외적으로 SAR 영상자료를 이용하여 해양의 대표적인 물리량인 주변지역의 해류 속도, 파랑, 바람, 내부파 등에 대한 정성적, 정량적 정보를 추출하기 위한 다양한 연구들이 진행되고 있는 상황이다.

1978년 최초 민간용 마이크로파 인공위성 SAR 센서인 SEASAT(USA)이 발사된 이후 SAR 위성 영상자료를 이용한 해양의 역학적 현상과 특징들을 수집하고 모니터링하는 연구들이 활발하게 진행되어 왔으며, 전 지구적인 육상, 해상, 기상의 정보를 추출하기 위한 다양한 인공위성 SAR 센서들이 개발되어 왔다. 이러한 인공위성 센서들은 ERS-1 (EU, 1991-2000), ERS-2 (EU, 1995-현재), JERS-1 (Japan, 1992-1998), RADARSAT-1 (Canada, 1995-현재), ENVISAT (EU, 2002-현재), ALOS PALSAR (Japan, 2006), RADARSAT-2 (Canada, 2006 예정) 등이 발사되어 지표면의 다양한 자료를 제공하고 있다. 최근 개발되고 있는 SAR 시스템들은 다양한 주파수(multi-frequency) 범위 내에서 측정 가능하고, 수평, 수직 편광 상태에 따른 모든 편광 조합이 가능한 다중편광(multi-polarization) 영상자료 취득이 가능하도록 개발되는 추세이며, 고해상도의 영상레이더와 반복 촬영 주기가 짧은 산란계의 특성을 조합한 방향으로도 발전하고 있다(예, ESA의 Cryosat-2).

또한, 우리나라도 2008년 다목적 실용위성 5호(KOMPSAT-5)와 향후 다목적 실용위성 7호(KOMPSAT-7)에 SAR 센서를 탑재한 관측위성이 발사될 계획이며, 이를 위한 인공위성 SAR 시스템의 도입 및 활용기술에 관한 기초적인 연구가 진행되고 있는 상황이다. 따라서, 향후 우리나라 주변지역에 대한 다양한 SAR 영상자료 취득이 용이할 것으로 기대되며, 특히 취득된 SAR 영상자료를 기반으로 파랑 정보를 비롯한 해상의 역학적 변화 및 특징과 관련된 해양물리학적 정보 추출을 위한 연구가 더욱 활발히 진행되어질 것으로 기대된다.

이를 대비한 기초 원천 기술의 확보가 필수적일 것으로 여겨지며, 해상의 파랑 특성과 변화 양상을 파악하는 작업은 해양의 기본적인 정보를 획득하는 의미뿐 아니라 해양과 관련되어 발생되는 모든 생산 활동과 그 안에서 일어나는 현상들의 예측 및 재난 방재 등과 같은 응용적인 측면에서의 활용을 위해 선행되어져야 할 기본적인 과정이라 할 수 있다.

본 연구는 해양의 바람(wind), 파랑(wave), 해류(current)와 같은 다양한 해양 현상에 대하여 인공위성 SAR 영상으로부터 풍속의 크기와 방향, 해파의 파장과 진행방향, 유속의 크기와 방향과 같은 해양물리학적 정보를 추출할 수 있는 종합적인 프로세서의 개발과 적용을 목표로 하였다.

이 연구에 적용한 연구 기법은 SAR 영상의 후방산란계수로부터 해상풍의 정보를 추출하기 위해 CMOD4 (Stoffelen and Anderson, 1997a; 1997b)와 CMOD_IFR2 (IFREMER-CERSAT, 1999; Quilfen et al., 1998) 경험 모델과 편광 비(polarization ratio)를 이용한 비례식(Horstmann et al., 2000)을 적용하여 풍속을 얻었다. 또한 Wave-SAR 역산 방법을 통해 파랑의 스펙트럼을 추출하고 파장과 파의 진행방향을 산출하였으며, 이 때 발생되는 파랑 방향에 대한 180° 모호성(ambiguity)의 문제를 해결하기 위하여 inter-look cross-spectra 기법(Dowd et al., 2001; Bao and Alpers, 1998; Engen and Johnsen, 1995)을 적용하였다. 해류의 속도 정보는 단위 산란체에 대하여 SAR 영상이 얻어지는 시간 동안에 해상의 산란체가 움직이면서 발생되는 Doppler shift를 이용한 기법(Chapron et al., 2005)으로 산출하였다. 또한 개발된 프로세서를 수행하여 SAR 영상으로부터 추출된 해풍, 파랑, 해류의 해양학적 파라미터를 기상자료나 해양자료 등과 같은 현장자료와 비교하여 추출된 결과를 검토하였다.

제 2 장 국내․외 기술개발 현황

1978년 미국의 SEASAT 위성이 발사된 이후 SAR를 이용한 해양 부분의 초기 연구는 해양과 관측 센서간의 기하학적 관계 등 기본원리에 관한 연구(Elachi and Brown, 1977; Tomiyasu, 1978; Jain, 1981; Ivanov, 1982; Moore, 1985) 내용이 주를 이루었으며, 이후 해양의 역학적 움직임과 현상을 이해하고 파악하기 위한 현상학적 해양 파라미터의 추정 및 검증을 위한 연구(Alpers and Rufenach, 1979; Alpers et al., 1981; Brüning et al., 1991)들이 시도되었다. 또한 선형적 거동을 보이는 파의 움직임뿐 아니라 준선형적(quasi-linear) 및 비선형적(nonlinear) 거동 특성을 보이는 파에 관한 연구(Swift and Wilson, 1979; Raney, 1980; Raney, 1981)들이 이루어져 왔다. Hasselmann and Hasselmann (1991)은 역산방법을 이용한 비선형파의 맵핑에 관한 연구에서 파의 진행방향 및 스펙트럼 특성을 분석하여 이후 다른 연구들에 도움을 주었다.

최근에는 선박, 구조물, 기상재해 등 해양활동과 관련된 다양한 활용 기술적 응용 분야와 알고리즘의 수정보완 및 개발(Dowd et al., 2001; Bao and Alpers, 1998; Engen and Johnsen, 1995; Chapron et al., 2005)에 관한 연구들이 활발히 진행되고 있다. 또한 interferometry 기술이나 polarimetry 기술을 적용한 응용 연구들과 dual-beam ATI (Along-Track Interferometry)와 같은 새로운 시스템 개발과 관련된 연구(Toporkov et al., 2005; Farquharson et al., 2004; Frasier and Camps, 2001)들이 시도되고 있다.

이(2006)는 SAR 관련 국내․외 연구들을 통계 분석하여 SAR 시스템, 기술, 활용분야별로 고찰하였는데 SAR를 이용한 해양 활용 분야를 wind, surface wave, internal wave, current, bathymetry, slick, ship, intertidal/coastal 나누어 분류하였다. 또한 최근 국내에서 발표된 해양관련 연구들을 살펴보면, 양(2006)은 SAR를 이용한 선박의 안전운항에 중요한 요소인 파랑 방향스펙트럼 추정에 관한 연구에서 SAR 영상자료를 활용하여 파랑을 계측하고 추정하는 기법을 시도하였고, 김(2001; 2005)은 ERS와 RADARSAT 영상자료로부터 경험모델을 이용하여 해상풍의 방향벡터추출 등 기존의 경험모델의 적용 및 수정에 관한 연구, 파랑, 해류, 내부파에 적용 가능성과 알고리즘의 수정보완을 위한 연구, 그리고 해양표면에 존재하는 파랑에 대한 준선형(quasi-linear) 변환을 시도한 바 있다.

제 3 장 연구개발수행 내용 및 결과

제 1 절 이론적 배경

SAR 영상에서 얻어질 수 있는 후방산란크기, 2차원적 스펙트럼, 도플러 등 측정할 수 있는 물리량을 통하여 해풍, 파랑, 해류의 특성을 추출해 내는 연구가 그동안 많은 진전을 보이고 있으나, 아직까지 각각의 방법에 있어서 한계성을 보이고 있다. Table 1은 본 연구에서 적용한 해풍, 파랑, 해류별로 각각의 연구기법과 추출 가능한 해양특성의 물리량과 한계점을 정리한 표이다.

측정 대상이 해풍인 경우 CMOD4와 CMOD_IFR2 모델은 10～20 m의 영상 해상도 내의 풍속(wind speed) 정보를 얻을 수 있으나 방향에 대한 정보는 현장자료나 관련 자료를 참조하여 해결해야 하는 한계점이 있다. 파랑 정보 추출 기법 중 wave-SAR 변환 기법은 방향에 대한 180° 모호성을 해결하지 못하는 한계가 있으며, 반면 inter-look cross-spectra 기법으로 이를 해결할 수 있다. 그러나 파랑의 파장이 SAR 영상의 해상도와 비교하여 작은 경우에는 파랑을 추출할 수 없는 한계가 있다. 해류와 같이 해상 표면의 움직임에 대한 이동속도를 구하기 위한 방법으로 Doppler shift 기법이 응용되고 있는데 이는 시야방향인 LOS (Line of Sight) 방향에 대한 정보만을 얻을 수 있는 한계가 있다. Dual-Beam ATI의 경우는 이러한 한계점을 해결할 수 있는 방법이나 최근 개발되고 있는 새로운 SAR 시스템을 필요로 한다.

본 연구에서는 SAR 영상에서 해풍, 파랑, 해류의 특성과 관련된 물리량을 종합적으로 추출하기 위한 프로그램을 개발하기 위해 이와 관련된 연구 기법들을 적용하였으며, 주요 내용들은 다음과 같다.

Table 1 SAR를 이용한 해양 연구
Target	Wind	Wave		Surface Velocity
Methods Values	Wind Retrieval Models	Wave-SAR Transform	inter-look Cross-Spectra	Doppler Shift	Dual-Beam ATI
Magnitude	○	○	○	△ LOS only	○
Direction	×	△ 180° ambiguity	○	△ LOS only	○
Spatial Resolution	10 m SAR Resolution	2 ㎞ SAR imagette	2 ㎞ SAR imagette	2 ㎞ SAR imagette	10 m SAR Resolution
Limitation	A priori wind direction required	Azimuth cutoff, strongly nonlinear		wind+current decomposition	New system

1. 해풍(wind)

가. CMOD4 모델

일반적으로 해상풍의 정보를 추출하는 알고리즘으로 사용되는 CMOD4 (Stoffelen and Anderson, 1997a) 경험 모델에서 후방산란계수 [eqn: sigma ^{0}]와 해상풍의 속도 [eqn: V]의 관계식은 다음 식과 같다.

[eqn: sigma ^{0} =b _{0} (1+b _{1} cos phi +b _{3} tanhb _{2} cos2 phi ) ^{1.6}]	(1)

여기서, [eqn: b _{0} =b _{r} TIMES 10 ^{alpha + gamma f _{1} (V+ beta )}]이며,

[eqn: f _{1} (y)= cases{pile{{-10}#{logy}}&pile{{if`y LEQ 10 ^{-10}}#{````````````if`10 ^{-10} <y LEQ 5}}#sqrt {y} /3.2&```````````if`y>5}] 이다.

[eqn: alpha ,` beta ,` gamma ,`b _{1} ,`b _{2} ,`b _{3}]는 Table 2에 나타난 18개의 계수(coefficients)로 전개된다.

[eqn: alpha =c _{1} P _{0} +c _{2} P _{1} +c _{3} P _{2}# gamma =c _{4} P _{0} +c _{5} P _{1} +c _{6} P _{2}# beta =c _{7} P _{0} +c _{8} P _{1} +c _{9} P _{2}# b _{1} =c _{10} P _{0} +c _{11} V+(c _{12} P _{0} +c _{13} V)f _{2} (x)# b _{2} =c _{14} P _{0} +c _{15} (1+P _{1} )V# b _{3} =0.42(1+c _{16} (c _{17} +x)(c _{18} +V))# b _{r} =LUT( theta )# f _{2} (x)=tanh{+2.5(x+0.35)}-0.61(x+0.35)]	(2)

여기서 [eqn: b _{r}] 은 residual correction factor이며, [eqn: phi]는 풍향(relative wind direction)이고 [eqn: theta]는 입사각(incidence angel)이다.  [eqn: b _{r}]과 [eqn: theta]의 값은 Table 3과 같이 주어진다. Legendre polynomials 식은 아래와 같다.

[eqn: P _{0} =1``````````P _{1} =x``````````P _{2} =(3x ^{2} -1)/2# with```x=( theta -40)/25]	(3)

Table 2 CMOD4 coefficients
Model		CMOD4
[eqn: alpha]	[eqn: c _{1}]	-2.301523
	[eqn: c _{2}]	-1.632686
	[eqn: c _{3}]	0.761210
[eqn: gamma]	[eqn: c _{4}]	1.156619
	[eqn: c _{5}]	0.595955
	[eqn: c _{6}]	-0.293819
[eqn: beta]	[eqn: c _{7}]	-1.015244
	[eqn: c _{8}]	0.342175
	[eqn: c _{9}]	-0.500786
[eqn: b _{1}]	[eqn: c _{10}]	0.014430
	[eqn: c _{11}]	0.002484
	[eqn: c _{12}]	0.074450
	[eqn: c _{13}]	0.004023
[eqn: b _{2}]	[eqn: c _{14}]	0.148810
	[eqn: c _{15}]	0.089286
[eqn: b _{3}]	[eqn: c _{16}]	-0.006667
	[eqn: c _{17}]	3.000000
	[eqn: c _{18}]	-10.00000

Table 3 Residual factors for CMOD4
[eqn: theta ^{0}]	[eqn: b _{r}]	[eqn: theta ^{0}]	[eqn: b _{r}]
16	1.075	39	0.988
17	1.075	40	0.998
18	1.075	41	1.009
19	1.072	42	1.021
20	1.069	43	1.033
21	1.066	44	1.042
22	1.056	45	1.050
23	1.030	46	1.054
24	1.004	47	1.053
25	0.979	48	1.052
26	0.967	49	1.047
27	0.958	50	1.038
28	0.949	51	1.028
29	0.941	52	1.016
30	0.934	53	1.002
31	0.927	54	0.989
32	0.923	55	0.965
33	0.930	56	0.941
34	0.937	57	0.929
35	0.944	58	0.929
36	0.955	59	0.929
37	0.967	60	0.929
38	0.978

나. CMOD_IFR2 모델

SAR 영상에서 해상풍을 추출하는 또 다른 경험식인 CMOD_IFR2 (IFREMER-CERSAT, 1999; Quilfen et al., 1998) 경험모델은 아래와 같이 주어진다.

[eqn: sigma ^{0} =10 ^{alpha + beta sqrt {V}} TIMES (1+b _{1} TIMES cos PHI +tanhb _{2} TIMES cos2 PHI )]	(4)

이때, [eqn: alpha ,` beta]는 다음과 같다.

[eqn: alpha =c _{1} +c _{2} P _{1} +c _{3} P _{2} +c _{4} P _{3}# beta =c _{5} +c _{6} P _{1} +c _{7} P _{2}]	(5)

여기서 [eqn: P _{1} ,`P _{2} ,`P _{3}]는 다음과 같다.

[eqn: P _{1} =x``````````P _{2} =(3x ^{2} -1)/2``````````P _{3} =x(5x ^{2} -3)/2# ``````````x=( theta -36)/19]	(6)

그리고 [eqn: b _{1} ,`b _{2}]는 아래 식으로 산출된다.

[eqn: b _{1} =c _{8} +c _{9} V _{1} +c _{10} P _{1} +c _{11} P _{1} V _{1} +c _{12} P _{2} +c _{13} P _{2} V _{1}# b _{2} =c _{14} +c _{15} P _{1} +c _{16} P _{2} +(c _{17} +c _{18} P _{1} +c _{19} P _{2} )V _{1}# ```````````````+(c _{20} +c _{21} P _{1} +c _{22} P _{2} )V _{2} +(c _{23} +c _{24} P _{1} +c _{25} P _{2} )V _{3}]	(7)

여기서, [eqn: P _{1}], [eqn: P _{2}], [eqn: V _{1}], [eqn: V _{2}], [eqn: V _{3}] 는 다음 식과 같으며, 파라미터 값(parameter values)은 Table 4에 정리하였다.

[eqn: P _{1} =x``````````P _{2} =2x ^{2} -1# ``````````x=(2 theta -76)/40]	(8)

[eqn: V _{1} =(2V-28)/22``````````V _{2} =2V _{1} ^{2} -1``````````V _{3} =(2V _{2} -1)V _{1}]	(9)

CMOD_IFR2 모델을 통하여 해상풍 속도 [eqn: V]를 추출한 이후 다음과 같은 보정 과정을 거치게 된다.

[eqn: V _{corrected} =V+V _{bias}]

여기서,

[eqn: V _{bias} =0`````if```V LEQ 10`m/s# V _{bias} =0.0831 TIMES V-0.0173 TIMES V ^{2} +0.0009 TIMES V ^{3} ```if`10<V LEQ 22m/s# V _{bias} =arctan(V-22)+3.0382`````if`V>22m/s]	(10)

Table 4 CMOD_IFR2 parameters
Parameter	Value	Parameter	Value
[eqn: c _{1}]	-2.437597	[eqn: c _{14}]	0.412754
[eqn: c _{2}]	-1.567031	[eqn: c _{15}]	0.121785
[eqn: c _{3}]	0.370824	[eqn: c _{16}]	-0.024333
[eqn: c _{4}]	-0.040590	[eqn: c _{17}]	0.072163
[eqn: c _{5}]	0.404678	[eqn: c _{18}]	-0.062954
[eqn: c _{6}]	0.188397	[eqn: c _{19}]	0.015958
[eqn: c _{7}]	-0.027262	[eqn: c _{20}]	-0.069514
[eqn: c _{8}]	0.064650	[eqn: c _{21}]	-0.062945
[eqn: c _{9}]	0.054500	[eqn: c _{22}]	0.035538
[eqn: c _{10}]	0.086350	[eqn: c _{23}]	0.023049
[eqn: c _{11}]	0.055100	[eqn: c _{24}]	0.074654
[eqn: c _{12}]	-0.058450	[eqn: c _{25}]	-0.014713
[eqn: c _{13}]	-0.096100

다. 편광 비례식

위에서 보듯이 식(1)과 식(4)는 C-band VV 편광인 경우에만 적용되는 경험식이며, CMOD4와 CMOD_IFR2 경험식에 적용되는 매개변수는 각각 다른 값을 가진다. 본 연구에서 사용된 RADARSAT-1 영상자료와 같이 C-band HH 편광 자료인 경우는 다음과 같은 편광 비(polarization ratio, PR)에 대한 비례식을 적용하게 된다.

[eqn: sigma _{0} ^{HH} = {(1+ alpha tan ^{2} theta ) ^{2}} over {(1+2tan ^{2} theta ) ^{2}} sigma _{0} ^{VV} (u,` phi ,` theta )]	(11)

[eqn: sigma _{0} ^{HH} = {(1+2sin ^{2} theta ) ^{2}} over {(1+2tan ^{2} theta ) ^{2}} sigma _{0} ^{VV} (u,` phi ,` theta )]	(12)

식(10)에서 [eqn: alpha]값에 따라 Bragg ([eqn: alpha]=0), Krichhoff ([eqn: alpha]=1), Unal ([eqn: alpha]=0.6) 산란 모델로 구분하며,  Elfouhaily 산란 모델은 식(11)로부터 산출된다. 실제 SAR 영상에 적용할 때에는 해양의 상태에 따라 적절한 경험식과 비례식을 선택하여 사용하게 되며, 이 때 현장 관측 자료와의 상호 비교 및 검증이 필수적이다.

2. 파랑(wave)

가. SAR Wave 스펙트럼

Hasselmann et al. (1996)은 파랑의 스펙트럼과 SAR의 2차원 스펙트럼과의 관계식을 다음과 같이 나타내었다.

[eqn: P( vec{k} )=exp(-k _{x} ^{2} xi ^{prime 2} ) sum _{n=1} ^{INF } sum _{m=2n-2} ^{2n} (k _{x} beta ) ^{m} P _{nm} ( vec{k} )]	(13)

[eqn: P( vec{k} )]: a SAR image spectrum

[eqn: k _{x}]: the wavenumber component of the long waves in the azimuthal direction

[eqn: beta]: the ratio of the slant range to the platform velocity (R/V)

[eqn: P _{nm}]: the spectral factors

[eqn: n]: the nonlinearity order with respect to the input wave spectrum

[eqn: m]: the order with respect to the velocity bunching parameter [eqn: beta]

식(13)에서 [eqn: xi ^{prime 2}]에 관한 식은 아래 식과 같다.

[eqn: xi ^{prime 2} = beta ^{2} LEFT < u _{r} ^{2} RIGHT > = beta ^{2} int _{} ^{} {LEFT | T _{k} ^{v} RIGHT | ^{2} F( vec{k} )d vec{k}}]	(14)

[eqn: xi ^{prime 2}]: the mean square azimuthal displacement

[eqn: u _{r}]: the radial component of the orbital velocity of the ocean waves

[eqn: F( vec{k} )]: an ocean wave spectrum

또한 range-velocity 변환함수 [eqn: T _{k} ^{v}]는 다음 식으로 산출된다.

[eqn: T _{k} ^{v} =- omega (sin theta {k _{y}} over {LEFT | k RIGHT |} +icos theta )]	(15)

[eqn: T _{k} ^{v}]: the range-velocity transfer function

[eqn: omega]: the radian frequency of the long waves

([eqn: omega =± sqrt {gk}] )

[eqn: theta]: the incidence angle

[eqn: k _{y}]: the wavenumber component of the long waves in the range direction

신호 처리 과정은 SAR 영상을 2차원적 Fourier 변환을 실시한 후 wave-SAR 변환을 적용한 결과에서 파랑의 주요 스펙트럼을 추출함으로써 대표적 파랑의 파장과 진행 방향을 결정하게 된다. 이 연구에서 사용된 FFT (Fast Fourier Transform, 고속 푸리에 변환)기법은 함수의  근사 값을 계산하는 알고리즘으로써 푸리에변환에 근거하여 근사공식을 이용한 이산푸리에변환(discrete Fourier transform)을 계산할 때 연산횟수를 줄일 수 있도록 고안된 알고리즘이다.

일반적으로 일정한 파수벡터 [eqn: vec{k _{0}}]와 각주파수 [eqn: omega _{0}], 그리고 진폭 A를 가지는 파동은 공간 [eqn: vec{r}] 과 시간 t의 함수인 [eqn: f]로 다음과 같이 나타낼 수 있다.

[eqn: f( vec{r} ,t)=Ae ^{j( vec{k _{0}} ․ vec{r} - omega _{0} t)}]	(16)

여기서 파랑은 2차원적인 파동이므로,

[eqn: vec{r} =(x,y)=x hat{x} +y hat{y}]

[eqn: vec{k _{0}} =(k _{x} ,k _{y} )=k _{x} hat{x} +k _{y} hat{y}]	(17)

와 같이 나타낼 수 있다.

이러한 선형의 파랑을 SAR 위성이 지나가면서 시간 [eqn: t=t _{0} `]에서 순간적으로 촬영하므로, SAR　영상에서는 [eqn: f]는 [eqn: vec{r}]만의 함수가 된다.

[eqn: f( vec{r} )=Ae ^{j( vec{k _{0}} ․ vec{r} - omega _{0} t _{0} )}]	(18)

이렇게 얻어진 2차원적 영상을 2D Fourier 변환을 시켜주면 다음과 같다.

[eqn: F( vec{k} )= int _{- INF } ^{INF } { int _{- INF } ^{INF } {f( vec{r} )e ^{-j vec{k} ․ vec{r}} dxdy}}# ```````````````````# ```````````````````= int _{- INF } ^{INF } { int _{- INF } ^{INF } {Ae ^{j( vec{k _{0}} ․ vec{r} - omega _{0} t _{0} )} e ^{-j( vec{k} ․ vec{r} )} dxdy}}]	(19)

[eqn: =Ae ^{-j omega _{0} t _{0}} int _{- INF } ^{INF } { int _{- INF } ^{INF } {e ^{-j( vec{k} - vec{k _{0}} )․ vec{r}} dxdy}}# =Ae ^{-j omega _{0} t _{0}} delta ( vec{k} - vec{k _{0}} )]

이 때  파수 [eqn: vec{k}]는 공간상에서 [eqn: vec{k _{0}}]에 해당하는 곳에서만 값을 가지고 다른 곳에서는 0의 값을 가지는 2차원 delta 함수의 형태가 된다. 이 함수의 절대 값을 취하면 다음과 같다.

[eqn: LEFT | F( vec{k} ) RIGHT | =A delta ( vec{k} - vec{k _{0}} )]	(20)

하지만 실제 SAR 영상에서는 real 값만이 존재하므로 다음 식과 같이 표현된다.

[eqn: f( vec{r} )=Acos( vec{k _{0}} ․ vec{r} - omega _{0} t _{0} )]	(21)

이를 Fourier 변환하면

[eqn: F( vec{k} )= int _{- INF } ^{INF } { int _{- INF } ^{INF } {Acos( vec{k _{0}} ․ vec{r} - omega _{0} t _{0} )e ^{-j( vec{k} ․ vec{r} )} dxdy}}# ```````````````````=A int _{- INF } ^{INF } { int _{- INF } ^{INF } {e ^{j( vec{k _{0}} ․ vec{r} - omega _{0} t _{0} )} +e ^{-j( vec{k _{0}} ․ vec{r} - omega _{0} t _{0} )} e ^{j( vec{k} ․ vec{r} )} dxdy}}]

[eqn: = {A} over {2} e ^{-j omega _{0} t _{0}} int _{- INF } ^{INF } { int _{- INF } ^{INF } {e ^{-j( vec{k} - vec{k _{0}} )․ vec{r}} dxdy}}# `````````+ {A} over {2} e ^{j omega _{0} t _{0}} int _{- INF } ^{INF } { int _{- INF } ^{INF } {e ^{-j( vec{k} + vec{k _{0}} )․ vec{r}} dxdy}}]

[eqn: = {A} over {2} e ^{-j omega _{0} t _{0}} delta ( vec{k} - vec{k _{0}} )+ {A} over {2} e ^{j omega _{0} t _{0}} delta ( vec{k} + vec{k _{0}} )]	(22)

가 된다. 이 때, [eqn: cos theta =(e ^{i theta } +e ^{-i theta } )/2] 관계를 사용하였다.

따라서[eqn: F( vec{k} )]에 절대 값을 취하여 스펙트럼(spectrum)을 구하면,

[eqn: LEFT | F( vec{k} ) RIGHT | = {A} over {2} delta ( vec{k} - vec{k _{0}} )+ {A} over {2} delta ( vec{k} + vec{k _{0}} )]	(23)

가 된다. 즉, [eqn: vec{k _{0}}]와 -[eqn: vec{k _{0}}] 위치에서만 값을 가지며, 나머지 영역에서는 값을 가지지 않는다. 또한 파랑의 방향성 판단에서 [eqn: vec{k _{0}}] 방향인지 -[eqn: vec{k _{0}}] 방향인지 알 수 없는 180° 모호성(ambiguity)이 발생되는 원인이 된다. Fig. 1과 같이 선형적으로 일정한 패턴을 보이는 파랑의 경우 real 자료를 2D Fourier 변환하면 Fig. 2a와 같이 원점대칭(symmetric)의 형태를 띠게 된다.

영상에서 볼 때 [eqn: vec{k _{0}}]은 [eqn: k _{x}]와 [eqn: k _{y}] 성분으로 나눌 수 있으며(Fig. 2b), [eqn: k _{x}]와 [eqn: k _{y}] 은 각각 다음 식으로 구할 수 있다.

[eqn: k _{x} =x축`픽셀`수 TIMES DELTA k _{x}] ,  [eqn: k _{y} =y축`픽셀`수 TIMES DELTA k _{y}]	(24)

이때 [eqn: DELTA k _{x}]와 [eqn: DELTA k _{y}]는 파수 영역에서 한 픽셀(pixel)당 파수의 값으로서 FFT의 성질에 따라 다음과 같이 나타낼 수 있다.

[eqn: DELTA k _{x} = {2 pi } over {X}]  (X : FFT를 한 imagette 영역의 x축 거리, meter)	(25)

[eqn: DELTA k _{y} = {2 pi } over {Y}]  (Y : FFT를 한 imagette 영역의 y축 거리, meter)	(26)

여기서 [eqn: vec{k _{0}}]는 2D FFT 결과로 나온 원점대칭(symmetric)의 영상에서 원점에서 중심 값까지의 벡터 값으로 이것은 곧 해양에서의 파랑에 적용 시 현재 진행하고 있는 파랑(wave) 혹은 너울(swell)의 파장 [eqn: lambda]값을 아는데 실마리를 제공한다.

위에서 구한 [eqn: k _{x}]와 [eqn: k _{y}]는 다음과 같이 파수(wavenumber)로 바꿀 수 있다.

[eqn: k= sqrt {k _{x} ^{2} +k _{y} ^{2}}]	(27)

이렇게 산출된 [eqn: k]값으로부터 아래와 같이 [eqn: lambda](파랑의 파장)를 구할 수 있게 된다.

[eqn: lambda = {2 pi } over {k}]	(28)

Fig. 1. 선형적 일정한 파랑의 모식도(사선은 파정을 나타냄).

(a) 2D FFT 결과 영상(원점대칭)	(b) 2D FFT 결과 [eqn: vec{k _{0}}]의 성분

Fig. 2. 2D FFT result of linear wave

나. Inter-look cross-spectra 기법

위에서 기술한 바와 같이 wave-SAR 변환식으로 추출된 파랑의 방향은 180° 모호성(ambiguity)을 가지며, 이러한 모호성을 해결하기 위해 inter-look cross-spectra 연구기법(Dowd et al., 2001; Bao and Alpers, 1998; Engen and Johnsen, 1995)을 적용하였다. 단일 지점에 대하여 SAR 영상을 얻을 때 약 0.5～1초 사이의 coherent integration time 동안 영상화가 이루어지는데, 이 시간 동안 파랑이 진행하면서 위치가 변하게 된다. 따라서 SAR 영상화 전처리 단계에서 coherent integration time을 나누어 서로 다른 시간에 얻어진 영상을 2개 이상 생성하고, 이 영상간의 부정합 관계를 이용하여 파랑의 이동 방향을 추출해 낸다. Cross-spectra 기법을 적용 시 파랑의 방향은 imaginary 스펙트럼 자료로부터 얻을 수 있다. Real 스펙트럼 결과는 대칭(symmetric) 형태를 보이고, imaginary 스펙트럼 결과는 비대칭(antisymmetric) 형태를 보이는데 imaginary 스펙트럼 결과에서 (+)쪽이나 (-)쪽으로 파랑의 진행방향이 결정된다. 본 연구에서는 (+)쪽이 파랑의 진행방향이다.

3. 해류(Current): Doppler shift 기법

SAR 영상이 얻어질 당시 산란체의 움직임은 Doppler shift로 기록된다. 해류의 이동속도(speed, [eqn: U _{D}])와 Doppler shift([eqn: f _{D}])와의 관계식은 다음과 같다(Chapron et al., 2005).

[eqn: U _{D} `=`-f _{D} `/`k _{e} sin theta _{I} `]	(29)

[eqn: U _{D}] : target speed towards radar (m/s)

[eqn: f _{D}] : Doppler shift

([eqn: f _{D} `=` DELTA f TIMES DELTA "p"ixel], [eqn: DELTA f=prf/"imagtte"`"size"])

[eqn: k _{e}] : the electromagnetic wavenumber ([eqn: 2 pi / lambda])

[eqn: theta _{I}] : the angle of incidence of the radar beam

[eqn: f _{D}]는 [eqn: DELTA f]와 nominal centroid와 예측된 Doppler centroid 간의 픽셀 차([eqn: DELTA "pixel"])로부터 계산된다. [eqn: DELTA f]는 PRF(Pulse Repetition Frequency) 값을 영상처리 시 일정한 영역에 해당하는 imagette size로 나눈 값이다. 해류 이동 속도의 부호는 산란체(해류)가 안테나 쪽을 향하여 움직이는 경우는 (+)값을, 멀어져가는 경우는 (-)값을 나타낸다. 일반적으로 Doppler shift는 SAR 영상화 전처리 단계(focusing)에서 Range 방향으로 처리가 이루어진 후 Azimuth 방향으로 1D Fourier 변환한 상태에서 이루어진다. SAR SLC 영상이 입력 자료인 경우에는 적당한 크기의 imagette을 잘라 Azimuth 방향으로 1D Fourier 변환된 상태에서 Doppler shift가 산출된다.

제 2 절 SAR Ocean Processor (SOP) 개발 및 자료처리

SOP (SAR Ocean Processor)는 Radarsat-1 SLC (single look complex) SAR 영상자료를 이용하여 해풍, 파랑, 해류 정보 추출 알고리즘을 기반으로 일련의 해역 특성 정보를 추출하도록 개발되었다. 이 프로그램은 Windows 환경에서 Linux 환경을 모사하는 Cygwin (www.cygwin.com)이라는 프로그램에서 컴파일되고 수행되지만, ANSI-C 언어로 이루어져 있기 때문에 사용자의 편의에 따라 다양한 컴파일러를 사용할 수 있다. 명령어 사용 환경에서 운영되기 때문에 사용자는 먼저 프로그램에 필요한 간단한 변수와 실행 제어 변수를 입력하여 파일로 만들고, 이를 통해 프로그램을 제어한다.

입력파일로는 Radarsat-1 원시자료(raw data)로부터 SAR focusing 과정을 거친  CEOS 포멧의 SLC (single look complex) 자료를 사용한다. 프로세서 수행결과 출력되는 파일은 해풍 정보인 CMOD4, CMOD_IFR2 자료, 파랑 정보를 갖고 있는 파랑 스펙트럼(wave spectrum), inter-look cross-spectra 자료이며, 해류의 속도 정보를 산출할 수 있는 Doppler shift 결과 자료이다(Table 5). 모든 결과는 원격탐사 및 GIS 전용 상용 프로그램으로 손쉽게 열어볼 수 있도록 설계하였다. 부록 A에서는 SOP 프로그램의 원시 코드를 수록하였다.

Table 5 SOP 프로세서 입출력

Input: RADARSAT-1 SLC (CEOS Format)

Output:

1. SLC Header Off

2. Multilook Image (by average)

3. Wind:

- Backscattering Coefficient (Sigma_naught, dB)

- CMOD_4 Wind Speed (m/s)

- CMOD_IFR2 Wind Speed (m/s)

4. Wave:

- Wave Spectrum

- Interlook Cross Spectrum

- Multilook (from interlook processing)

5. Current:

- Doppler Image

- Doppler shift estimation(vector file)

- Velocity estimation (text file)

- Velocity vectors (vector file)

Fig. 3. SOP (SAR Ocean Processor) 흐름도.

1. RADARSAT-1 영상 자료

RADARSAT-1 시스템은 캐나다항공우주국(Canadian Space Agency)에서 발사한 SAR 시스템으로 5.3 ㎓ 주파수의 C-band HH 편광을 사용하며, Fine, Standard, Wide, ScanSAR narrow, ScanSAR wide, Extended high, Extended low 등 다양한 영상모드(imaging mode)에서 관측이 가능하다. Table 6은 RADARSAT-1 SAR의 기계적 특성을 정리한 표이며, Table 7과 Fig. 4는 영상모드별로 특징을 정리한 표와 그림이다.

Table 6 RADARSAT-1 SAR 특성
Frequency/Wavelength	5.3 ㎓/C-band 5.6 ㎝
Radio frequency bandwidth	11.6, 17.3 or 30.0 ㎒
Transmitter power (peak)	5 ㎾
Transmitter power (average)	300 W
Maximum data rate	85 Mb/s (recorded) - 105 Mb/s (R/T)
Antenna size	15 m × 15 m
Antenna polarization	HH

Table 7 영상모드(imaging modes)
Mode	Nominal Resolution (m)	No. of Positions / Beams	Swath Width (㎞)	Incidence Angles (degrees)
Fine	8	15	45	37-47
Standard	30	7	100	20-49
Wide	30	3	150	20-45
ScanSAR narrow	50	2	300	20-49
ScanSAR wide	100	2	500	20-49
Extended high	18-27	3	75	52-58
Extended low	30	1	170	10-22

Fig. 4. 영상모드(imaging mode)

이 연구에서 사용된 RADARSAT-1 영상 자료는 standard 모드에서 관측된 헤더정보가 포함된 SLC (single look complex, CEOS format) 자료이며, 관측날짜, 시각, 관측모드, 관측지역을 정리하면 Table 8과 같다. 제주도 근해를 포함한 1999년 11월 15일, 11월 25일, 12월 19일 3개 영상 자료는 해풍과 파랑 정보 추출 과정에서 사용되었으며, 서해안 일부 해역을 포함한 2004년 11월 8일 1개 영상 자료는 해류의 정보 추출 과정에서 사용되었다.

Table 8 RADARSAT-1 영상 자료
	관측날짜	관측시각	관측모드	영상지역
1	1999.11.15	06:34	Descending	제주도
2	1999.11.25	06:42	Descending	제주도
3	1999.12.19	06:42	Descending	제주도
4	2004.11.08	18:32	Ascending	서해안

2. SLC Header Off 영상 추출

SOP 프로세서를 수행하면 먼저 입력 영상 자료로부터 헤더정보가 제거된 SLC 영상자료가 만들어지며, 이후 이 SLC 영상 자료로부터 일련의 해풍, 파랑, 해류 정보 추출 과정이 수행된다. Fig. 5에서 Fig. 8은 Table 8에 정리된 사용 영상 자료들의 SLC amplitude 영상들이다. SLC 영상의 경우 전형적으로 SAR 영상에서 나타나는 스펙클 잡음들이 산재하여 나타나는데 이는 다음 단계에서 수행되는 multi-look 과정을 통해 제거할 수 있다.

Fig. 5. SLC header off image (1999. 11. 15).

Fig. 6. SLC header off image (1999. 11. 25).

Fig. 7. SLC header off image (1999. 12. 19).

Fig. 8. SLC header off image (2004. 11. 8).

3. Multilook 영상 추출

SAR SLC 영상에는 스펙클 잡음이 혼재되어 있어 영상의 밝기를 추정하는데 많은 에러를 줄 수 있다. 이러한 스펙클을 제거하는 방법은 보통 두 가지가 있는데, 첫째는 SAR 원시자료를 focusing하는 단계에서 coherent integration time을 수 개로 나누어 각각 영상을 얻고 이를 합산하는 방법을 쓴다. SLC 영상을 얻은 경우에는 SAR focusing 단계로 거슬러 올라가서 coherent integration time을 나누어 여러 영상을 만든 후 합산하는 방식을 사용할 수 있다. 두 번째는 SLC 영상 자료를 azimuth 방향으로 averaging filter 처리하면 위의 방법과 매우 흡사한 스펙클 제거 효과를 볼 수 있어 ML (multilook) 영상 자료를  얻을 수 있다. SOP 프로그램에서는 위의 두가지 방법을 모두 다 적용하였는데, 효율적인 자료 처리를 위하여 inter-look cross spectrum 방법에만 SAR focusing 개념의 multilook 방식으로 자료처리를 하였고, wind나 wave 추출시에는 SLC에 average filter를 한 multilook 영상을 사용하였다. 본 연구에서 사용된 4개의 영상 자료의 ML 영상 결과는 Fig. 9에서 Fig. 12와 같다. ML 영상 결과에서 보듯이 SLC 영상 내에 존재하였던 스펙클 잡음을 어느 정도 제거할 수 있다. SAR 영상은 표면 거칠기와 유전율에 의해 영상의 밝기 값이 결정되는데 영상 내에서 해양은 해상이 잔잔한 경우는 어둡게 나타나며, 바람에 의해 해상의 표면 거칠기가 증가할수록 밝게 나타난다.

4. 해풍(wind) 정보 추출

SOP 수행결과 해풍과 관련된 결과는 후방산란계수(backscattering coefficient, σ0)와 CMOD4, CMOD_IFR2 모델 기법을 적용하여 산출된 풍속 결과로 도출된다. 후방산란계수의 단위는 [㏈]이며, CMOD4와 CMOD_IFR2 모델로 추출된 풍속의 단위는 [㎧]이다.

바람에 의한 영향은 영상 내에 밝고 어두운 지역으로 살펴볼 수 있는데 바람이 센 지역은 밝고, 바람이 약한 지역은 다소 어둡게 나타난다. 해상의 표면이 매우 잔잔한 경우는 아주 어두운 검정색으로 보인다. 어두운 검정색으로 나타나는 지역(잔잔한 해상)은 -20 ㏈ 이하의 값을 보이며, 표면이 거친 지역(바람이 센 지역)으로 예상되는 밝은 지역은 -5 ㏈ 정도의 값을 나타내고 있다. 회색 지역으로 나타나는 중간영역은 바람의 정도에 따라 -20 ㏈에서 -5 ㏈ 사이 값을 보이고 있다.

CMOD4와 CMOD_IFR2 경험식은 C-band VV 편광인 경우의 모델로 개발되었으므로 이를 HH 편광모드에 적용하기 위해서 앞서 기술한 PR 비례식인 식(11), 식(12)를 사용하게 되는데 이때 산란모델에 따라 Bragg, Unal, Krichhoff, Elfouhaily로 나누어서 결과를 도출하였다. 일반적으로 태풍과 같이 강한 바람이 부는 경우는 CMOD_IFR2 모델보다 CMOD4 모델이 적합한 결과를 도출하는 것으로 알려져 있으며, 특히 SAR 영상의 입사각이 큰 경우와 바람이 센 경우는 경험모델과 비례식 선택 시 주의해야 한다.

Fig. 9. SAR Multilook image (1999. 11. 15).

Fig. 10. SAR Multilook image (1999. 11. 25).

Fig. 11. SAR Multilook image (1999. 12. 19).

Fig. 12. SAR Multilook image (2004. 11. 8).

5. 파랑(wave) 정보 추출

가. 파랑 스펙트럼(SAR wave spectrum) 추출

파랑 특성을 추출하기 위한 wave SAR 변환이나 inter-look cross-spectra 과정은 일정한 크기(imagette)로 나누어 수행되어지는데 Fig. 13에서 Fig. 16은 1999년 12월 19일 영상자료의 같은 영역에 대해 imagette 크기에 따라 64×64, 128×128, 256×256, 512×512의 추출 결과를 도시한 것이다. imagette의 크기가 작을수록 해상도는 향상되나 imagette의 크기가 너무 작거나 너무 큰 경우는 영상 내의 스펙트럼 특성을 추출하는데 부적합하게 된다. 따라서 이러한 문제를 해결하기 위해 너울의 특성이 잘 발달된 1999년 12월 19일 영상 자료를 이용하여 영상의 스펙트럼 특성에 적합한 imagette 크기를 결정하고자 하였다. 너울과 같은 발달된 파랑의 유무를 고려할 때 전반적으로 256×256과 128×128 크기일 때 영상의 특성을 추출하는데 적합한 것으로 판단되었다.

Fig. 17에서 Fig. 19는 제주도 근해 지역의 1999년 11월 15일, 11월 25일, 12월 19일 3개 영상 자료를 2D FFT 기법을 적용하여 추출된 파랑 스펙트럼 결과이다. 영상 전체 지역을 도시한 그림이며, imagette 크기는 256×256이다. 너울의 특성을 뚜렷하게 보이는 12월 19일 파랑 스펙트럼 추출 영상에서는 스펙트럼 형태와 방향이 뚜렷하게 나타나는 반면 11월 15일과 11월 25일 영상에서는 스펙트럼의 형태가 불확실하거나 분산된 형태로 나타나는 경우가 많다. 이는 12월 19일 영상과 달리 너울과 같이 충분히 발달된 파랑의 특성이 나타나지 않았기 때문으로 부분적으로 뚜렷한 파랑 스펙트럼 형태를 보이는 영역이 있으나 전반적으로 스펙트럼 형태가 분명하게 나타나지 않는다.

Fig. 13. SAR wave spectrum 결과(1999년 12월 19일, 64×64 imagette).

Fig. 14. SAR wave spectrum 결과(1999년 12월 19일, 128×128 imagette).

Fig. 15. SAR wave spectrum 결과(1999년 12월 19일, 256×256 imagette).

Fig. 16. SAR wave spectrum 결과(1999년 12월 19일, 512×512 imagette).

Fig. 17. SAR wave spectrum (1999년 11월 15일, 256×256 imagette). GCP 보정 전 영상.

Fig. 18. SAR wave spectrum (1999년 11월 25일, 256×256 imagette). GCP 보정 전 영상.

Fig. 19. SAR wave spectrum (1999년 12월 19일, 256×256 imagette). GCP 보정 전 영상.

나. Inter-look cross-spectra 추출

Fig. 20은 inter-look cross-spectra 과정을 통해 추출된 real, imaginary 스펙트럼 결과를 도시한 것이다. 1999년 12월 19일 영상 자료를 사용하였으며, 같은 영상 영역을 64×64, 128×128, 256×256 imagette 크기에 따라 추출한 결과이다. Fig. 20 (e)와 (f)에서 보듯이 real 스펙트럼은 대칭(symmetry) 형태를 보이며, imaginary 스펙트럼은 비대칭(antisymmetry) 형태를 나타내는데 이때 파랑의 진행 방향은 imaginary 스펙트럼의 (+)쪽(빨간색 영역)이다.

Fig. 20. Inter-look cross-spectra 추출.

6. 해류(current) 정보 추출

SOP에서 해류 관련 파라미터를 입력 후 처리하여 도출되는 결과는 Doppler 영상과 벡터 파일 형식의 Doppler shift 예측 정보, 벡터와 텍스트 파일 형식의 속도 정보로 나눌 수 있다. 이 연구에서는 조류의 흐름이 예상되는 서해안 해역의 2004년 11월 8일 영상 자료를 사용하였으며, 영상처리 결과를 검토하기 위해 국립해양조사원(www.nori.go.kr)에서 제공하는 같은 날짜와 시각에 대한 예측조류도 자료를 활용하였다.

가. Doppler 영상 추출

Doppler shift 기법을 적용하여 Fig. 21과 같이 도플러 영상을 추출하였으며, 추출된 도플러 영상에 이론적 도플러 중심(nominal Doppler centroid)과 추출된 도플러 중심(Doppler centroid)을 도시할 수 있는 Doppler shift estimation 벡터 파일을 중첩시켜 Fig. 22에 도시하였다. Fig. 23은 Fig. 22의 빨간색 상자에 해당하는 영역의 확대 영상이며, 이때 파란색 점은 이론적 도플러 중심, 빨간색 점은 Doppler shift 기법을 적용하여 추출된 도플러 중심이다. Fig. 23 영상에서 보듯이 빨간색 점이 위 또는 아래로 이동되어 있는데 이는 조석의 흐름과 같이 해상에서 발생되는 해류의 영향으로 발생되며, 이러한 차이로부터 레이더 파의 진행방향인 LOS에 대한 벡터 정보를 추출할 수 있다. 해류의 이동방향이 시야방향에 가까워지는 경우와 멀어지는 경우로 나누어지며, 다른 방향의 벡터 성분은 산출하지 못하는 한계가 있다.

Fig. 21. Doppler image (2004년 11월 8일 서해안 인근 해역).

Fig. 22. Doppler shift 추출 결과(2004년 11월 8일 서해안 인근 해역). 파란점은 이론적 도플러 중심(nominal Doppler center), 빨간점은 해류에 의한 Doppler shift.

Fig. 23. Doppler shift 추출 결과(덕적도 부근, 2004년 11월 8일). Fig. 22의 빨간색 상자 영역의 확대 영상.

나. 속도 자료 추출

Fig. 24와 Fig. 25는 ML 영상위에 추출된 속도 벡터(velocity vector) 추출 결과를 중첩한 그림이다. 그림에서 보듯이 덕적도 부근 남쪽(A)과 동쪽(B) 해역에 해당하는 영상지역의 벡터 추출 결과는 부분적으로 왼쪽과 오른쪽 방향 모두 나타나지만 주로 왼쪽 방향이 우세하게 나타나고 있다. 또한 국립해양조사원에서 제공하는 조류예측도 자료(Fig. 32 참조)와 비교하여 볼 때 양호한 결과를 보였다.

Doppler shift 기법으로 추출된 이동속도는 텍스트 파일(velocity estimation text file)로도 출력되며, range 방향(x), azimuth 방향(y), Doppler shift(㎐), 속도(㎧) 순으로 나열되어 있다. Fig. 26은 2004년 11월 8일 서해안 인근 해역의 영상 처리 결과 산출된 속도 텍스트 파일을 워드패드에서 열어본 그림이다. 영상 전체적으로 추출된 해류 이동속도는 수십 ㎝에서 4 ㎧ 범위의 값을 나타냈다.

Fig. 24. Doppler shift 벡터 추출 결과(2004년 11월 8일, 덕적도 부근 남쪽 해역(A)).

Fig. 25. Doppler shift 벡터 추출 결과(2004년 11월 8일, 덕적도 부근 동쪽 해역(B)).

Fig. 26. Velocity estimation text file.

제 3 절 SOP 적용 결과 및 현장자료와의 비교

본 절에서 사용된 자료는 RADARSAT-1 영상자료로 제주도 인근 해역을 관측한 1999년 11월 15일, 1999년 11월 25일, 1999년 12월 19일 3개 영상과 2004년 11월 8일 관측된 서해안 인근 해역의 1개 영상이다. SOP 프로세서 수행 결과 도출된 자료 분석 결과를 비교 검토하기 위한 현장자료로는 우도에서 관측된 AWS (Automatic Weather System) 기상자료와 국립해양조사원(www.nori.go.kr)에서 제공하는 인천 인근의 서해안 해역의 조류예측도 자료를 참조하였다.

CMOD4와 CMOD_IFR2 모델을 적용하여 산출된 풍속과 후방산란계수(σ0)의 관계를 Fig. 27에 도시하였다. Fig. 27에 도시된 그래프는 1999년 11월 15일 제주 인근 해역을 관측한 영상자료로부터 추출된 값으로 대략 -20～-5 dB 범위의 후방산란계수 값을 나타내고 있으며, 각각의 산란모델에 따라 풍속은 1～11 m/s의 값을 보이고 있다. 이는 영상의 일부지역에서 측정된 결과이며, 다른 연구(김, 2005)에서 보고되어진 후방산란계수와 풍속의 관계와 부합하는 결과를 보이고 있다.

Fig. 27. CMOD4와 CMOD-IFR2 모델별 후방산란계수(σ0)와 풍속(wind speed) 관계.

Fig. 28은 1999년 12월 19일 제주도 인근 해역의 RADARSAT-1 영상자료로부터 추출된 파랑 스펙트럼과 cross-spectra 과정을 통해 도출된 amplitude, real, imaginary 스펙트럼 결과를 도시한 것이다. Wave-SAR 변환식을 적용하여 추출된 파랑 스펙트럼과 amplitude의 Cross-Spectrum 결과는 [eqn: vec{k _{0}}], [eqn: - vec{k _{0}}] 위치에서 값을 보이는 유사한 대칭 형태를 보이며, Cross-Spectrum 결과에서 real과 달리 imaginary 스펙트럼 결과는 비대칭 형태를 나타내고 있다. 이에 추출된 imaginary 스펙트럼 결과에서 (+)쪽을 파랑의 진행방향이 결정되며, 이는 현장에서 관측된 AWS 자료와도 일치된다(Table 9 참조).

Fig. 29는 wave-SAR 변환식을 적용한 256×256 imagette 결과 영상이며, 파랑의 크기인 파장(wavelength)을 구하기 위해 추출된 파랑 스펙트럼으로부터 파수벡터 [eqn: vec{k _{0}}]를 구하고, 앞서 기술된 식(28)로부터 파랑의 파장([eqn: lambda])을 산출하여 Table 9에 정리하였다.

Fig. 28. SAR wave spectrum과 inter-look cross-spectra 기법 적용 결과(128×128 imagette). GCP 보정(ground control point correction) 전 영상자료 임.

Fig. 29. SAR wave spectrum (2D FFT) 결과(256×256 imagette). GCP 보정 후 영상자료: (A) 1999년 11월 15일, (B) 1999년 11월 25일, (C) 1999년 12월 19일.

Fig. 29에 도시된 영상지역은 제주도 인근 해역이며, 1999년 11월 15일, 1999년 11월 25일, 1999년 12월 19일에 관측된 자료들이다. Table 9에 정리된 바와 같이 12월 19일(C) 영상은 평균풍향 280°, 평균풍속 7 ㎧의 바람이 7일 이상 지속되어 영상 내에서 너울(swell)의 형태를 육안으로도 살펴볼 수 있다. 상대적으로 11월 15일(A)과 11월 25일(B) 영상에서는 너울의 형태를 육안으로 식별하기는 어려우며, 부분적인 영역에서 파랑의 특성을 보이고 있다. 또한 파랑이 너울과 같이 충분히 발달된 경우는 육안식별로도 영상 특징을 구별할 수 있을 뿐 아니라 파랑 스펙트럼 분석 결과도 신뢰할만한 값을 보이고 있다. 일반적으로 너울과 같이 충분히 발달된 파로 발전하기 위해서는 3일 이상의 지속시간이 필요함을 고려해볼 때 11월 15일과 11월 25일 영상은 뚜렷한 파랑의 특성을 나타내기엔 부족한 기상상태였음을 알 수 있다. 이는 영상 스펙트럼 결과에서도 나타나고 있으며, 산출된 측정파장의 오차범위 또한 상대적으로 크게 나타남을 알 수 있다.

Table 9 RADARSAT-1 영상자료의 wave 결과와 AWS 자료
영상취득날짜 측정값		1999년 11월 15일 (A)	1999년 11월 25일 (B)	1999년 12월 19일 (C)
파랑자료 (위성영상 분석)	측정파향	281° (± 4°)	286° (± 10°)	281° (± 1°)
	측정파장	94.5 m (± 14 m)	99.5 m (± 28 m)	156 m (± 8 m)
바람자료 (AWS)	평균풍향	285°	277°	280°
	지속시간	2일 정도	1일 이하	7일 이상
	평균풍속	9 ㎧	11 ㎧	7 ㎧

Fig. 30과 Fig. 31은 Doppler shift 기법을 적용하여 추출된 속도 벡터 자료를 ML 영상 위에 중첩하여 도시한 그림이다. 속도 벡터 정보가 잘 보이도록 충분히 확대한 영상이며, 덕적도 부근의 서쪽, 동쪽, 남쪽 해역에 해당한다. 이때 추출된 이동속도(current speed)는 수십 ㎝에서 4 ㎧ 정도의 값을 보였다. Doppler shift 기법은 위성체의 진행방향에 수직하며, 레이더 파의 진행방향인 시야방향(LOS)에 대한 정보만을 추출할 수 있는 한계가 있다. 따라서 해류의 이동방향에 대한 정보는 시야방향으로 가까워지는 경우와 멀어지는 경우로 나누어지며, 다른 방향의 벡터 성분은 산출하지 못한다. 앞서 도시된 Fig. 24와 Fig. 25 결과에서 추출된 벡터 방향은 주로 영상의 왼쪽 방향을 향하고 있음을 살펴보았다. Fig. 30과 Fig. 31에서도 빨간색 화살표의 방향은 영상 내 왼쪽과 오른쪽 방향 모두 나타나고 있지만 주된 방향은 왼쪽 방향임을 알 수 있다. 따라서 덕적도 인근 해역의 주된 흐름은 영상에서 왼쪽 방향인 서쪽 방향으로 이동하였을 것으로 예측되며, 이를 검토하기 위해 국립해양조사원에서 제공하고 있는 예측조류도 자료를 활용하였다.

Fig. 30. 속도 추출 결과(덕적도 부근 서쪽 해역, 2004년 11월 8일). 빨간색 화살표는 속도 벡터 방향.

Fig. 31. 속도 추출 결과(덕적도 부근 동쪽과 남쪽 해역, 2004년 11월 8일). 빨간색 화살표는 속도 벡터 방향.

Fig. 32는 국립해양조사원에서 제공하는 예측조류도 자료이다. 사용된 RADARSAT-1 영상자료가 관측된 2004년 11월 8일 오후 6시 30분 자료를 검색한 결과이며, 실측자료가 아닌 예측자료이므로 오차가 있을 수 있으나 Fig. 30과 Fig. 31 추출 결과를 비교하여 보면 전반적으로 이동방향이 일치하는 것을 볼 수 있다. Fig. 32에서 보듯이 덕적도 주변의 유속은 근처 다른 영역에 비해 높은 값을 보이고 있으며, 추출된 이동속도 값과 비교하였을 때 유사한 값을 나타냄을 알 수 있다. 이동 방향은 대략 270～290° 범위로 섬주변의 북서쪽과 서쪽 방향으로 나타나고 있다. Doppler shift 기법을 적용하여 추출된 속도 벡터 자료의 이동 방향은 서쪽 방향으로 단순히 시야방향에 대한 정보의 한계성이 있으나 조류의 예측 방향과 부합되는 결과를 보이는 것으로 판단된다.

Fig. 32. 예측조류도 (2004년 11월 8일 18시 30분, 국립해양조사원).

제 4 장 연구개발목표 달성도 및 대외기여도

1. 본 연구의 목표인 해양물리학적 정보 추출을 위한 종합적인 프로세서 개발은 해상에서 발생되는 해풍, 파랑, 해류의 특성과 관련 물리량을 산출하기 위한 것이며, 개발되어진 SOP (SAR Ocean Processor) 수행결과와 현장관련 자료를 비교 검토하여 신뢰할 만할 결과를 얻을 수 있었다.

2. 개발된 SOP 프로세서는 향후 해양에서 발생되는 다양한 현상을 관측한 영상자료 및 동시 관측 가능한 현장 실측 자료들이 확보되고, 정확하고 정밀한 정보 추출을 위한 알고리즘 보완 및 수정 단계를 거치게 되면 더욱 향상된 해양물리학적 정보 추출이 가능할 것으로 판단된다.

3. SAR 영상 자료로부터 해양 파라미터를 추출할 수 있는 종합적인 프로세서의 필요성을 고려해볼 때 본 연구에서 개발된 SOP은 아직 초기단계이나 이러한 기대에 부흥하고 있으며, 향후 발전 가능성이 매우 높다고 할 수 있다.

4. 개발된 SOP 프로세서를 수행하여 다양한 영상 자료와 현장 실측 자료를 기초로 추출된 자료들은 정보의 신뢰성 면에서 정확도 및 정밀도를 향상시킬 수 있을 뿐 아니라 해양에서 발생되는 역동적이고 다양한 현상들을 규명하기 위한 정보와 해법을 제공할 수 있을 것으로 판단된다.

제 5 장 연구개발결과의 활용계획

1. 연구개발 된 SOP 프로세서는 해양 현상과 관련하여 해풍, 파랑, 해류 정보를 제공할 수 있는 종합 프로그램으로 활용 가능할 것으로 기대된다.

2. SOP 프로세서는 연구 또는 교육 목적에 따라 발전 가능하므로 해양 연구와 관련하여 기초 자료로도 활용 가능할 것으로 기대된다.

3. 해양학적 현상과 특징을 분석하는 작업은 해양의 기작을 규명하는 것과 같은 순수 목적뿐 아니라 해상의 상태가 수시로 변화하는 역학적 특성으로 인해 어업, 해상운송, 기상예측 등 다양한 분야에 종사하는 사람들로부터 가능한 실시간 정보(near-real time information) 제공을 요구받고 있다. 따라서 인공위성을 통해 해상 특성과 관련된 정보를 추출하고 분석하는 작업은 매우 중요한 과정들이며, SOP 처리 결과 추출된 자료는 제공 가능한 범위의 해양물리학적 정보를 제공할 수 있을 것으로 기대된다.

제 6 장 참고문헌

[1] 김덕진 (2001), 인공위성 SAR 자료를 이용한 해상풍 벡터 추출, 석사학위논문, 서울대학교.

[2] 김덕진 (2005), SAR를 이용한 해양의 해상풍, 해류, 파랑 및 내부파의 정량적 원격탐사, 박사학위논문, 서울대학교.

[3] 양찬수 (2006), SAR에 의한 파랑 방향스펙트럼 추정에 관한 연구, 2006 춘계학술대회 논문집, 대한원격탐사학회, pp. 71-74.

[4] 이훈열 (2006), SAR 관련 논문 통계 분석에 의한 SAR 시스템, 기술, 활용분야 고찰, 대한원격탐사, 22권, 2호, pp. 153-174.

[5] Alpers, W. R. and Brüning, C. (1986), On the Relative Importance of Motion-Related Contributions to the SAR Imaging Mechanism of Ocean Surface Waves, IEEE Transaction on Geoscience and Remote Sensing, vol. GE-24, no. 6, pp. 873-885.

[6] Alpers, W. R. and Rufenach, C. L. (1979), The Effect of Orbital Motions on Synthetic Aperture Radar Imagery of Ocean Waves, IEEE Transaction on antennas and propagation, vol. AP-27, no. 5, pp. 685-690.

[7] Alpers, W. R., Ross, D. B. and Rufenach, C. L. (1981), On the Detect Ability of Ocean Surface Waves by Real and Synthetic Aperture Radar, Journal of geophysical research, vol. 86, no. C7, pp. 6481-6498.

[8] Bao, M. and Alpers, W. (1998), On the Cross Spectrum Between Individual-Look Synthetic Aperture Radar Images of Ocean Waves, IEEE Transactions on Geoscience and Remote Sensing, vol. 36, no. 3, pp. 922-932.

[9] Brüning C., Alpers, W. R. and Schroeter, J. G. (1991), On the Focusing Issue of Synthetic Aperture Radar Imaging of Ocean Waves, IEEE Transaction of Geoscience and Remote Sensing, vol. GE-29, no. 1, pp. 120-128.

[10] Chapron, B., Collard, F., and Ardhuin, F. (2005), Direct Measurements of Ocean Surface Velocity from Space: Interpretation and Validation, J. of Geophysical Research, vol. 110, pp. 1-17.

[11] Curlander, J. C. and McDonough, R. M. (1991), Synthetic Aperture Radar Systems and Signal Processing, John Wiley & Sons, Inc., New York.

[12] Dowd, M., Vachon, P. W., Dobson, F. W., and Olsen R. B. (2001), Ocean Wave Extraction from RADARSAT Synthetic Aperture Radar Inter-Look Image Cross-Spectra, IEEE Transaction on Geoscience and Remote Sensing, vol. 39, no. 1, pp. 21-37.

[13] Elachi, C. (1988), Spaceborne Radar Remote Sensing: Applications and Techniques, IEEE Press, NY., pp. 38-50.

[14] Elachi, C. and Brown, W. E. (1977), Models of Radar Imaging of the Ocean Surface Waves, IEEE transaction on antennas and propagation, vol. AP-25, no. 1, pp. 84-94.

[15] Engen, G. and Johnsen, H. (1995), SAR-Ocean Wave Inversion Using Image Cross Spectra, IEEE transaction on antennas and propagation, vol. 33, no. 4, pp. 1047-1056.

[16] Farquharson, G., Junek, W. N., Ramanathan, A., Frasier, S. J., Tessier, R., McLaughlin, D. J., Sletten, M. A., and Toporkov, J. V. (2004), IEEE Geoscience and Remote Sensing Letters, vol. 1, no. 2, pp. 62-65.

[17] Frasier, S. J. and Camps, A. J. (2001), Dual-Beam Interferometry for Ocean Surface Current Vector Mapping, IEEE Transaction on Geoscience and Remote Sensing, vol. 39, no. 2, pp. 401-414.

[18] Hasselmann, K. and Hasselmann, S. (1991), On the Nonlinear Mapping of an Ocean Wave Spectrum into a Synthetic Aperture Radar Imang Spectrum and Its Inversion, J. Geophysical Research, vol. 96, pp. 10713-10729.

[19] Hasselmann, S., Brüning, C., Hasselmann K. and Heimbach P. (1996), An Improved Algorithm for the Retrieval of Ocean Wave Spectra from Synthetic Aperture Radar Image Spectra, J. of Geophysical Research, vol. 101, no. C7, pp. 16615-16629.

[20] Horstmann, J., Koch, W., Lehner, S., and Tonboe, R. (2000), Wind Retrieval over the Ocean using Synthetic Aperture Radar with C-band HH Polarization, IEEE Transactions on Geoscience and Remote Sensing, vol. 38, no. 5, pp. 2122-2131.

[21] IFREMER-CERSAT (1999), Off-Line Wind Scatterometer ERS Products: User Manual, Technical Report C2-MUT-W-01-IF, IFREMER-CERSAT.

[22] Ivanov, A. V. (1982), On the Synthetic Aperture Radar Imaging of Ocean Surface Waves, IEEE Journal of oceanic engineering, vol. OE-7, no. 2, pp. 96-103.

[23] Jain, A. (1981), SAR Imaging of Ocean Waves: Theory, IEEE Journal of oceanic engineering, vol. OE-6, no. 4, pp. 130-139.

[24] Moore, R. K. (1985), Radar Sensing of the Ocean, IEEE Journal of oceanic engineering, vol. OE-10, no. 2, pp. 84-113.

[25] Mouchot, M. -C., Garello, R. (1998), SAR for Oceanography in Henderson, F. M. and Lewis, A. J., Ed., Principles & Application of Imaging Radar, John Wiley & Sons, NY., pp. 631-669.

[26] Quilfen, Y., Chapron, B., Elfouhaily, T., Katsaros, K., and Tournadre, J. (1998), Observation of Tropical Cyclones by High-Resolution Scatterometry, J. of Geophysical Research, vol. 103, no. C4, pp. 7767-7787.

[27] Raney, R. K. (1980), SAR Response to Partially Coherent Phenomena, IEEE transaction on antennas and propagation, vol. AP-28, no. 6, pp. 777-787.

[28] Raney, R. K. (1981), Wave Orbital Velocity, Fade and SAR Response to Azimuth Waves, Journal of oceanic engineering, vol. OE-6, no. 4, pp. 140-146.

[29] Stoffelen, A. and Anderson, D. (1997a), Scatterometer Data Interpretation: Estimation and Validation of the Transfer Function CMOD4, J. of Geophysical Research, vol. 102, no. C3, pp. 5767-5780.

[30] Stoffelen, A. and Anderson, D. (1997b), Scatterometer Data Interpretation: Measurement Space and Inversion, J. of Atmospheric and Oceanic Technology, vol. 14, pp. 1298-1313.

[31] Swift, C. T. and Wilson, L. R. (1979), Synthetic Aperture Radar Imaging of Moving Ocean Waves, IEEE transaction on antennas and propagation, vol. AP-27, no. 6, pp. 725-729.

[32] Tomiyasu, K. (1978), Tutorial Review of Synthetic Aperture radar (SAR) with Applications to Imaging of the Ocean Surface, Proceedings of IEEE, vol. 66, no. 5, pp. 563-583.

[33] Toporkov, J. V., Perkovic, D., Farquharson, G., Sletten, M., and Frasier, S. J. (2005), Sea Surface Velocity Vector Retrieval Using Dual-Beam Interferometry: First Demonstration, IEEE Transaction on Geoscience and Remote Sensing, vol. 41, no. 11, pp. 2494-2502.

부록 A. SAR Ocean Processor 원시프로그램

1. Program Main: sop.c

/*******************************************************

sop.c

SAR Ocean Processor

to extract Ocean Wind, Current and Wave from SAR

by

Prof. Hoonyol Lee

Department of Geophysics

Kangwon National Univeristy

email: hoonyol@kangwon.ac.kr

Input: Radarsat-1 SLC

Output:

1. SLC Header Off (ERMapper format)

2. Multilook Image (by average)

3. Wind:

- Backscattering Coefficient (Sigma_naught, dB)

- CMOD_4 Wind Speed (m/s)

- CMOD_IFR2 Wind Speed (m/s)

4. Wave:

- Wave Spectrum

- Interlook Cross Spectrum

- Multilook (from interlook processing)

5. Current:

- Doppler Image

- Doppler shift (vector file)

- Velocity (vectr file)

- Velocity (text file)

Versions:

Initial: 17 July 2006

0.1 (18 July 2006): SLC extraction

0.2 (24 July 2006): Multilook and SAR wave

0.3 (25 July 2006): Doppler

0.4 (12 August 2006): Cross Spectrum

0.9 (10 November 2006): Wind, Wave and Current

1.0 (20 November 2006): true incidence angle

1.1 (22 November 2006): split wave and current

**********************************************************/

#define VERSION "1.1"

#define NBLOCK 1024

#include <stdio.h>

#include <stdlib.h>

#include <math.h>

#include <string.h>

#include <time.h>

#include "para.h"

#include "ma.c"

#include "swapbytes.c"

#include "gen_ers.c"

#include "gen_ers_cmod.c"

#include "gen_ers_crs.c"

#include "load_para.c"

#include "bip2bil.c"

#include "bip2biln.c"

#include "load_slc.c"

#include "multilook_resmp.c"

#include "cfft.c"

#include "cshift.c"

#include "maverage.c"

#include "generv.c"

#include "wind.c"

#include "wave.c"

#include "wave_crs.c"

#include "current.c"

void usage();

void intro();

main(int argc, char *argv[])

{

time_t clock_start, clock_end, clock_duration;

intro();

if(argc != 2){usage(); exit(0);}

clock_start=time(NULL);

/* Load Parameters and Data */

printf("Load Parameters and Data\n");

load_para(argv[1]);

if(c_load_slc) load_slc();

/* Multilook (average and resample SLC)  */

if(c_multilook_resmp) multilook_resmp(fname_slc, fname_ml,

nazimuth, nrange, nx, ny, nx_res, ny_res);

/** Wind **/

if(c_wind) wind();

/** Wave **/

if(c_wave) wave();

/** Wave - Interlook Cross Spectrum **/

if(c_wave_crsflag>0) wave_crs();

/* Current */

if(c_current) current();

clock_end=time(NULL);

clock_duration= clock_end - clock_start;

printf("\nProgram finished successfully (%d sec). \n",	clock_duration);

}

void usage()

{

printf("\n     Usage: sop parameter_file\n");

}

void intro()

{

printf("\n");

printf(" *************************************************** \n");

printf("*	SAR Ocean Processor (SOP)          *\n");

printf("*	v%s	*\n", VERSION);

printf("*     All right reserved by Dr. Hoonyol Lee,        *\n");

printf("* Dept. of Geophysics, Kangwon National University. *\n");

printf("* Email:hoonyol@kangwon.ac.kr, Tel:+82-33-250-8587  *\n");

printf(" *************************************************** \n");

}

2. Parameter Definition: para.h

/* Parameter Definition for eSAR */

/*   Input Parameters       */

char fname_in[200], fname_leader[200], fhead[80];

int c_load_slc, c_multilook_resmp;

int nx, ny, nx_res, ny_res;

int c_wind;

double wind_direction, vv2hh_alpha;

int c_wave, nptxw, nptyw;

int c_wave_crsflag;

int c_current, nptyc, dc_xstep, dc_ystep, dc_avx, dc_avy;

double vel_scale;

/* Data parameter */

int nrange_whole, nazimuth_whole;

int asc_des;

int nbyte_header_record, nbyte_data;

int nbyte_per_record, nskip_record, nskip_byte_per_record;

double earth_radius_major, earth_radius_minor;

int nrange, nazimuth;

int nrange_ml, nazimuth_ml; /* multilook */

double nominal_incidence_angle, wavelength;

double prf, azimuth_bw;

double Dc0, Dc1, Dc2;	/* cross track doppler centroid, Hz,	Hz/pixel, Hz/pixel/pixel */

double range_pixel_size, azimuth_pixel_size;

double ground_range_pixel_size;

double drange, dazimuth; /* nominal pixel size for ERMapper	output */

double drange_ml, dazimuth_ml; /* nominal multilook pixel	size for ERMapper output */

/* wind */

int inc_step=30; /* LUT every inc_step range pixels */

double sigma0_min, sigma0_max;

double **cmod4_lut, **cmodifr2_lut;

double sig_min_cmod4, sig_max_cmod4;

double sig_min_cmodifr2, sig_max_cmodifr2;

/* SAR Leader File - Radiometric Data Record */

int samp_inc, n_data_pixel;

double rad[512], plat_lat, plat_long, plat_head;

double eph_orb_data, srgr_coef[6];

double *inc_ang, *inc_angr, inc_ang_min, inc_ang_max;

/* Filename */

char fname_slc[200], fname_ml[200], fname_wave[200], fname_crs[200];

char fname_cur[200], fname_sigma0[200], fname_cmod[200];

char fname_crs_ml[200], fname_cur_erv[200], fname_cur_txt[200];

char fname_cur_vec[200];

/* General Paramters */

double cspeed, pi;

3. Load Parameter File: load_para.c

/*  load_para.c

1. Load parameter files

2. Derive relevant paramters

3. Allocate memory

*/

void load_para(char *fname_para)

{

FILE *fp, *fp1, *fp2;

char vbuf[80], *buf;

int i, j, malloc_total, irad;

double emin, emax, tp2, er, rs, h;

printf(" load_para:\n");

/* Input parameter */

printf("  <Input Parameters>\n");

if((fp=fopen(fname_para, "r")) == NULL)

{

printf("  (Error) Cannot open %s.\n\a", fname_para);

exit(0);

}

fscanf(fp, "%s %s", vbuf, fname_in);

printf("   %s=%s\n", vbuf, fname_in);

fscanf(fp, "%s %s", vbuf, fname_leader);

printf("   %s=%s\n", vbuf, fname_leader);

fscanf(fp, "%s %s", vbuf, fhead);

printf("   %s=%s\n", vbuf, fhead);

printf("\n  <Control Paramters>\n");

fscanf(fp, "%s %d", vbuf, &c_load_slc);

printf("   %s= %d \n", vbuf, c_load_slc);

fscanf(fp, "%s %d", vbuf, &c_multilook_resmp);

printf("   %s= %d \n", vbuf, c_multilook_resmp);

fscanf(fp, "%s %d %d", vbuf, &nx, &ny);

printf("   multilook %s= %d, %d \n", vbuf, nx, ny);

fscanf(fp, "%s %d %d", vbuf, &nx_res, &ny_res);

printf("   multilook %s= %d, %d\n", vbuf, nx_res, ny_res);

fscanf(fp, "%s %d", vbuf, &c_wind);

printf("   %s= %d \n", vbuf, c_wind);

fscanf(fp, "%s %lf", vbuf, &wind_direction);

printf("   %s= %g \n", vbuf, wind_direction);

fscanf(fp, "%s %lf", vbuf, &vv2hh_alpha);

printf("   %s= %g \n", vbuf, vv2hh_alpha);

if(vv2hh_alpha>=0 && vv2hh_alpha<=1)

printf("    Thompson Model\n");

else if(vv2hh_alpha==2) printf("     Elfouhaily Model\n");

else

{

printf(" Error.. Wrong VV to HH conversion model.\n");

exit(0);

}

fscanf(fp, "%s %d", vbuf, &c_wave);

printf("   %s= %d \n", vbuf, c_wave);

fscanf(fp, "%s %d %d", vbuf, &nptxw, &nptyw);

printf("    %s= %d, %d\n", vbuf, nptxw, nptyw);

fscanf(fp, "%s %d", vbuf, &c_wave_crsflag);

printf("   %s= %d \n", vbuf, c_wave_crsflag);

if(c_wave_crsflag==2) printf("    IQ Cross Spectrum\n");

else if(c_wave_crsflag==1) printf("    AP Cross Spectrum\n");

fscanf(fp, "%s %d", vbuf, &c_current);

printf("   %s= %d \n", vbuf, c_current);

fscanf(fp, "%s %d", vbuf,  &nptyc);

printf("    %s= %d\n", vbuf, nptyc);

fscanf(fp, "%s %d %d", vbuf, &dc_xstep, &dc_ystep);

printf("    %s= %d, %d\n", vbuf, dc_xstep, dc_ystep);

fscanf(fp, "%s %d %d", vbuf, &dc_avx, &dc_avy);

printf("    %s= %d, %d\n", vbuf, dc_avx, dc_avy);

fscanf(fp, "%s %lf", vbuf, &vel_scale);

printf("   %s= %lf\n", vbuf, vel_scale);

printf("\n");

fclose(fp);

/* Data Parameter - from

Image Options File Descriptor Record (SAR Data file) */

printf("  <Data File Parameters>\n");

if((fp1=fopen(fname_in, "r")) == NULL)

{

printf("  (Error) Cannot open %s.\n\a", fname_in);

exit(0);

}

fseek(fp1, 5, SEEK_SET);

fread(&nskip_byte_per_record, 1, 1, fp1);

printf("   nskip_byte_per_record=%d\n", nskip_byte_per_record);

fseek(fp1, 8, SEEK_SET);

fread(&nbyte_header_record, 4, 1, fp1);

byteswap(&nbyte_header_record, 4);

printf("   nbyte_header_record=%d\n", nbyte_header_record);

buf=cvector(100);

fseek(fp1, 180, SEEK_SET);

fread(buf, 6, 1, fp1);

nazimuth_whole=atoi(buf);

printf("   nazimuth_whole=%d\n", nazimuth_whole);

free(buf);

buf=cvector(100);

fseek(fp1, 186, SEEK_SET);

fread(buf, 6, 1, fp1);

nbyte_per_record=atoi(buf);

printf("   nbyte_per_record=%d\n", nbyte_per_record);

free(buf);

buf=cvector(100);

fseek(fp1, 216, SEEK_SET);

fread(buf, 4, 1, fp1);

nbyte_data=atoi(buf)/8;

printf("   nbyte_data=%d\n", nbyte_data);

free(buf);

buf=cvector(100);

fseek(fp1, 248, SEEK_SET);

fread(buf, 8, 1, fp1);

nrange_whole=atoi(buf);

printf("   nrange_whole=%d\n", nrange_whole);

free(buf);

printf("\n");

fclose(fp1);

/* Data Parameter - from SAR Leader File */

printf("  <Leader File Parameters>\n");

if((fp2=fopen(fname_leader, "r")) == NULL)

{

printf("  (Error) Cannot open %s.\n\a", fname_leader);

exit(0);

}

buf=cvector(100);

fseek(fp1, 720+100, SEEK_SET);

fread(buf, 16, 1, fp1);

if(strncmp(buf, "DESCENDING", 3)==0) asc_des= -1;

else if(strncmp(buf, "ASCENDING", 3)==0) asc_des=1;

else

{

printf("  Error Ascending/Descending unknown\n");

exit(0);

}

printf("   Orbit orientation(asc_des)=%s(%d)\n", buf, asc_des);

free(buf);

buf=cvector(100);

fseek(fp2, 720+180, SEEK_SET);

fread(buf, 16, 1, fp2);

earth_radius_major = atof(buf);

printf("   earth_radius_major=%lf\n", earth_radius_major);

free(buf);

buf=cvector(100);

fseek(fp2, 720+196, SEEK_SET);

fread(buf, 16, 1, fp2);

earth_radius_minor = atof(buf);

printf("   earth_radius_minor=%lf\n", earth_radius_minor);

free(buf);

buf=cvector(100);

fseek(fp2, 720+452, SEEK_SET);

fread(buf, 8, 1, fp2);

plat_lat = atof(buf);

printf("   plat_lat=%lf(deg)\n", plat_lat);

free(buf);

buf=cvector(100);

fseek(fp2, 720+460, SEEK_SET);

fread(buf, 8, 1, fp2);

plat_long = atof(buf);

printf("   plat_long=%lf(deg)\n", plat_long);

free(buf);

buf=cvector(100);

fseek(fp2, 720+468, SEEK_SET);

fread(buf, 8, 1, fp2);

plat_head = atof(buf);

printf("   plat_head=%lf(deg)\n", plat_head);

free(buf);

buf=cvector(100);

fseek(fp2, 720+484, SEEK_SET);

fread(buf, 8, 1, fp2);

nominal_incidence_angle = atof(buf);

printf("   nominal_incidence_angle=%lf\n", nominal_incidence_angle);

free(buf);

buf=cvector(100);

fseek(fp2, 720+500, SEEK_SET);

fread(buf, 16, 1, fp2);

wavelength = atof(buf);

printf("   wavelength=%lf\n", wavelength);

free(buf);

buf=cvector(100);

fseek(fp2, 720+934, SEEK_SET);

fread(buf, 16, 1, fp2);

prf = atof(buf);

printf("   prf=%lf\n", prf);

free(buf);

buf=cvector(100);

fseek(fp2, 720+1206, SEEK_SET);

fread(buf, 16, 1, fp2);

azimuth_bw = atof(buf);

printf("   azimuth_bw=%lf\n", azimuth_bw);

free(buf);

buf=cvector(100);

fseek(fp2, 720+1478, SEEK_SET);

fread(buf, 16, 1, fp2);

Dc0 = atof(buf);

printf("   Dc0=%e\n", Dc0);

free(buf);

buf=cvector(100);

fread(buf, 16, 1, fp2);

Dc1 = atof(buf);

printf("   Dc1=%e\n", Dc1);

free(buf);

buf=cvector(100);

fread(buf, 16, 1, fp2);

Dc2 = atof(buf);

printf("   Dc2=%e\n", Dc2);

free(buf);

buf=cvector(100);

fseek(fp2, 720+1686, SEEK_SET);

fread(buf, 16, 1, fp2);

azimuth_pixel_size = atof(buf);

printf("   azimuth_pixel_size=%lf\n", azimuth_pixel_size);

free(buf);

buf=cvector(100);

fseek(fp2, 720+1702, SEEK_SET);

fread(buf, 16, 1, fp2);

range_pixel_size = atof(buf);

printf("   range_pixel_size=%lf\n", range_pixel_size);

free(buf);

/* Detailed Processing Parameter Record */

irad=720+4096+1620+16920*2; /* begin of record*/

fseek(fp2, irad+4648, SEEK_SET);

buf=cvector(100);

fread(buf, 16, 1, fp2);

eph_orb_data=atof(buf);

printf("   eph_orb_data=%g\n", eph_orb_data);

free(buf);

fseek(fp2, irad+4907, SEEK_SET);

for(i=0; i<6; i++)

{

buf=cvector(100);

fread(buf, 16, 1, fp2);

srgr_coef[i]=atof(buf);

free(buf);

}

/* Radiometric Data Record */

irad=720+4096+1620+16920*2+7726+8960+8960; /* begin of radio. cal */

buf=cvector(100);

fseek(fp2, irad+84, SEEK_SET);

fread(buf, 4, 1, fp2);

samp_inc=atoi(buf);

printf("   Radiometric samp_inc =%d\n", samp_inc);

free(buf);

/*printf(" Radiometric Para number\n");*/

for(i=0; i<512; i++)

{

buf=cvector(100);

fread(buf, 16, 1, fp2);

rad[i]=atof(buf);

/*	printf(" %d, %g\n", i, rad[i]);

*/	free(buf);

}

printf("\n");

fclose(fp2);

/* Generic Parameters */

cspeed=2.997925e8;

pi = 4.*atan(1.);

printf("\n  <Generic Parameters>\n");

printf("   cspeed= %lf m/s\n", cspeed);

printf("   pi= %lf\n", pi);

printf("\n");

/* Other Parameters */

printf("  <Other Parameters>\n");

nskip_record=1;

printf("   nskip_record=%d\n", nskip_record);

nrange=nrange_whole;

nazimuth=nazimuth_whole;

ground_range_pixel_size = range_pixel_size

/ sin(nominal_incidence_angle*pi/180.);

printf("   ground_range_pixel_size=%lf\n", ground_range_pixel_size);

drange = ground_range_pixel_size;

dazimuth = azimuth_pixel_size;

printf("   Output drange=%lf, dazimuth=%lf\n", drange, dazimuth);

drange_ml = drange*nx_res;

dazimuth_ml = dazimuth*ny_res;

nrange_ml = nrange/nx_res;

nazimuth_ml=nazimuth/ny_res;

printf("   Multilook drange_ml=%lf, dazimuth_ml=%lf\n", drange_ml, dazimuth_ml);

printf("   Multilook nrange_ml=%d, nazimuth_ml=%d\n", nrange_ml, nazimuth_ml);

/* incidence angles */

emin=earth_radius_minor;

emax=earth_radius_major;

tp2=tan(plat_lat*pi/180.)*tan(plat_lat*pi/180.);

er=emin*sqrt((1.+tp2)/(emin*emin/emax/emax+tp2))*1000.;

h=eph_orb_data-er;

printf("   Earth radius = %g(m)\n", er);

printf("   Orbit Altitude = %g(m)\n", h);

/*	printf("      range(pix), incidence angle(deg)\n");

*/

inc_ang=dvector(nrange_whole);

inc_angr=dvector(nrange_whole);

for(i=0; i<nrange_whole; i++)

{

if(asc_des==-1)  /* descending: far range first */

rs=srgr_coef[0]+drange_ml*(nrange_whole-i);

else /* ascending: near range first */

rs=srgr_coef[0]+drange_ml*i;

inc_angr[i]=acos( (h*h-rs*rs+2.*er*h)/(2.*rs*er) );

inc_ang[i]=inc_angr[i]*180./pi;

/*printf("    %d  %g\n", i, inc_ang[i]); */

}

if(asc_des==-1)	/* descending: far range first */

{

inc_ang_max=inc_ang[0];

inc_ang_min=inc_ang[nrange_whole-1];

}

else

{

inc_ang_min=inc_ang[0];

inc_ang_max=inc_ang[nrange_whole-1];

}

printf("   Incidence Angles \n");

printf("     first range %d,  %g(deg)\n", 0, inc_ang[0]);

printf("     mid range %d, %g(deg)\n", nrange_whole/2, inc_ang[nrange_whole/2]);

printf("     far range %d, %g(deg)\n", nrange_whole, inc_ang[nrange_whole-1]);

/* output filenames */

strcpy(fname_slc, "01_");

strcat(fname_slc, fhead);

strcat(fname_slc, "_SLC");

printf("   fname_slc=%s\n", fname_slc);

strcpy(fname_ml, "02_");

strcat(fname_ml, fhead);

strcat(fname_ml, "_ML");

printf("   fname_ml=%s\n", fname_ml);

strcpy(fname_sigma0, "03_");

strcat(fname_sigma0, fhead);

strcat(fname_sigma0, "_Wind_Simga0");

printf("   fname_sigma0=%s\n", fname_sigma0);

strcpy(fname_cmod, "03_");

strcat(fname_cmod, fhead);

strcat(fname_cmod, "_Wind_CMOD");

printf("   fname_cmod=%s\n", fname_cmod);

strcpy(fname_wave, "04_");

strcat(fname_wave, fhead);

strcat(fname_wave, "_Wave");

printf("   fname_wave=%s\n", fname_wave);

strcpy(fname_crs, "04_");

strcat(fname_crs, fhead);

strcat(fname_crs, "_Wave_ILC");

printf("   fname_crs=%s\n", fname_crs);

strcpy(fname_crs_ml, "04_");

strcat(fname_crs_ml, fhead);

strcat(fname_crs_ml, "_Wave_ML");

printf("   fname_crs_ml=%s\n", fname_crs_ml);

strcpy(fname_cur, "05_");

strcat(fname_cur, fhead);

strcat(fname_cur, "_Current");

printf("   fname_cur=%s\n", fname_cur);

strcpy(fname_cur_erv, "05_");

strcat(fname_cur_erv, fhead);

strcat(fname_cur_erv, "_Current_Doppler");

printf("   fname_cur_erv=%s\n", fname_cur_erv);

strcpy(fname_cur_vec, "05_");

strcat(fname_cur_vec, fhead);

strcat(fname_cur_vec, "_Current_Vector");

printf("   fname_cur_vec=%s\n", fname_cur_vec);

strcpy(fname_cur_txt, "05_");

strcat(fname_cur_txt, fhead);

strcat(fname_cur_txt, "_Current_Vector.txt");

printf("   fname_cur_txt=%s\n", fname_cur_txt);

}

4. Load SAR SLC Data: load_slc.c

/*****************************************************************

load_slc

1. load SLC in 2-byte CEOS format

2. write output in 2-byte RAW format

by Dr. Hoonyol Lee

Dept. of Geophysics

Kangwon National Univeristy

15 June 2004

12 July 2004: oflag added

18 July 2006: single file

*****************************************************************/

void load_slc()

{

FILE *fp_in, *fp_out;

int i;

short *sbuf;

int nremnant_byte_per_record;

printf("\n load_slc:\n");

printf("  Extract SAR SLC data. \n");

sbuf=svector(nrange*2);

nremnant_byte_per_record = nbyte_per_record

- nskip_byte_per_record

- nrange*nbyte_data*2;

/*** READ AND WRITE CEOS SLC ***/

if((fp_in=fopen(fname_in, "r"))==NULL)

{

printf(" (Error) Cannot open %s.\n\a", fname_in);

exit(0);

}

if((fp_out=fopen(fname_slc, "w"))==NULL)

{	printf(" (Error) Cannot open %s.\n\a", fname_slc);

exit(0);

}

printf("  SAR SLC data output %s.\n", fname_slc);

/*header + azimuth record skip*/

fseek(fp_in, nbyte_header_record*nskip_record, SEEK_SET);

/* Data input and output */

for(i=0; i<nazimuth; i++)

{

/* skip first bytes per record*/

fseek(fp_in, nskip_byte_per_record, SEEK_CUR);

/* read range data */

fread(sbuf, 2, nrange*2, fp_in);

/* write range data */

swapbytes(sbuf, nrange*2, nbyte_data);

bip2bil(sbuf, nrange, 2, nbyte_data);

fwrite(sbuf, 2, nrange*2, fp_out);

/* skip rest of the range per record */

fseek(fp_in, nremnant_byte_per_record, SEEK_CUR);

}

fclose(fp_in);

fclose(fp_out);

gen_ers(fname_slc, nazimuth, nrange,

dazimuth, drange, 2, "120");

free(sbuf);

}

5. Multilook by Resampling: multilook_resmp.c

/*****************************************************

multilook_resmp.c

Input: 2-byte short integer complex in BIL file format

Process: calculate magnitude and do average

Output: RESAMPLED multilook magnitude (4-byte Real) file

void multilook(fname_in, fname_out, nrow, ncol, nx, ny, row_block)

fname_in: input file name

fname_out: output file name

nrow: number of row

ncol: number of column

nx: averaging factor in column

ny: averaging factor in row

nx_res: resampling factor in x

ny_res: resampling factor in y

by Hoonyol Lee

Department of Geophysics, Kangwon National University

23 July 2006

******************************************************/

void multilook_resmp(fname_in, fname_out, nrow, ncol, nx, ny, nx_res, ny_res)

char *fname_in,  *fname_out;

int nrow, ncol, nx, ny, nx_res, ny_res;

{

FILE *fp_in, *fp_out;

int i, j, k, l, kx, write_size_col, write_size_row, pstep, nstep;

short *sbuf;

double tmp1, tmp2, tmp3;

float *fbuf;

unsigned seekcount, readcount, nyr, ncount;

printf("\n Multilook by Average and Resampling \n");

write_size_row = nrow/ny_res;

write_size_col=ncol/nx_res;

sbuf=svector(ny*ncol*2); /* input */

fbuf=fvector(write_size_col); /* output */

if((fp_in=fopen(fname_in, "r"))==NULL)

{

printf(" (Error) Cannot open %s.\n\a", fname_in);

exit(0);

}

if((fp_out=fopen(fname_out, "w"))==NULL)

{	printf(" (Error) Cannot open %s.\n\a", fname_out);

exit(0);

}

nstep=0;

pstep=0;

for(i=0; i<nrow; i+= ny_res)

{

if(i>pstep && nstep<10)

{

nstep ++;

printf("   %d/10\n", nstep);

pstep += nrow/10;

}

seekcount=(i-ny/2)*ncol*2*2;

nyr=ny;

readcount=nyr*ncol*2;

if((i-ny/2)<0) /*first block */

{

seekcount=0;

nyr=i+ny/2;

readcount=nyr*ncol*2;

}

if((i+ny/2)>=nrow) /* last block */

{

nyr -= (i+ny/2-nrow);

readcount -= (i+ny/2-nrow)*ncol*2;

}

fseek(fp_in, seekcount, SEEK_SET);

fread(sbuf, 2, readcount, fp_in);

for(j=0; j<ncol; j+=nx_res)

{

tmp3=0.;

ncount=0;

for(k=0; k<nyr; k++)

{

kx=k*ncol*2;

for(l=(j-nx/2); l<(j+nx/2+1); l++)

{

if(l>=0 && l<ncol)

{

ncount++;

tmp1=(double)sbuf[kx+l];

tmp2=(double)sbuf[kx+l+ncol];

tmp3+=sqrt(tmp1*tmp1+tmp2*tmp2);

}

}

}

tmp3 /= ncount;

fbuf[j/nx_res]=(float)tmp3;

}

fwrite(fbuf, 4, write_size_col, fp_out);

} /* end of i */

gen_ers(fname_out, write_size_row, write_size_col,

dazimuth*ny_res, drange*nx_res, 1, "141");

free(sbuf);

free(fbuf);

}

6. Wind by Sigma0, CMOD4, CMOD_IFR2: wind.c

/***********************************

wind.c

Sigma0, CMOD4 and CMOD_IFR2

by

Hoonyol Lee

Department of Geophyscis

Kangwon National University

8 November 2006

***********************************/

double sigma0();

double sigma0_vv2hh();

double cmod4_v2sigma();

double cmodifr2_v2sigma();

void cmod4_LUT();

void cmodifr2_LUT();

double cmod4();

double cmodifr2();

void wind()

{

FILE *fp_ml, *fp_sig, *fp_cmod;

int i, ix, j, jx;

float *fbuf, *fbufout;

double tmp, tmp1, tmp2;

printf("\n");

printf(" Sigma0, CMOD4 and CMOD-IFR2 for Wind Extraction\n");

fbuf=fvector(nrange_ml);

fbufout=fvector(nrange_ml*2);

if((fp_ml=fopen(fname_ml, "r"))==NULL)

{

printf(" (Error) Cannot open %s.\n\a", fname_ml);

exit(0);

}

if((fp_sig=fopen(fname_sigma0, "w"))==NULL)

{	printf(" (Error) Cannot open %s.\n\a", fname_sigma0);

exit(0);

}

if((fp_cmod=fopen(fname_cmod, "w"))==NULL)

{	printf(" (Error) Cannot open %s.\n\a", fname_cmod);

exit(0);

}

/* Sigma Naught from Multilook SAR */

printf("   Sigma0:\n");

sigma0_min=1e15;

sigma0_max=0.;

for(i=0; i<nazimuth_ml; i++)

{

fread(fbuf, 4, nrange_ml, fp_ml);

for(j=0; j<nrange_ml; j++)

{

tmp=(double)fbuf[j];

tmp=sigma0(tmp, j);

if(tmp>sigma0_max) sigma0_max=tmp;

if(tmp<sigma0_min) sigma0_min=tmp;

if(tmp<-100.) printf("%d,%d,%g  \n", i, j, tmp);

fbufout[j]=(float)tmp;

}

fwrite(fbufout, 4, nrange_ml, fp_sig);

}

printf("     min=%g, max=%g\n", sigma0_min, sigma0_max);

gen_ers(fname_sigma0, nazimuth_ml, nrange_ml, dazimuth_ml,

drange_ml, 1, "141");

fclose(fp_ml);

fclose(fp_sig);

/* LUT for CMOD4 and CMOD_IFR2 */

cmod4_LUT();

cmodifr2_LUT();

/* Sigma Naught, CMOD4, and CMOD_IFR2 */

if((fp_sig=fopen(fname_sigma0, "r"))==NULL)

{	printf(" (Error) Cannot open %s.\n\a", fname_sigma0);

exit(0);

}

for(i=0; i<nazimuth_ml; i++)

{

fread(fbuf, 4, nrange_ml, fp_sig);

for(j=0; j<nrange_ml; j++)

{

tmp=(double)fbuf[j];

fbufout[j]=(float)cmod4(tmp, j);

fbufout[j+nrange_ml]=(float)cmodifr2(tmp, j);

}

fwrite(fbufout, 4, nrange_ml*2, fp_cmod);

}

gen_ers_cmod(fname_cmod, nazimuth_ml, nrange_ml, dazimuth_ml, drange_ml, 2, "141");

free(fbuf);

free(fbufout);

fclose(fp_sig);

fclose(fp_cmod);

}

double sigma0(double dn, int irange)

{

int j1, j2;

double x, a, s, s0;

j1=irange/samp_inc;

x=fmod(irange, samp_inc)/samp_inc;

j2=j1+1;

a=rad[j1]*(1.-x)+ rad[j2]*x;

s=20.*log10(dn/a);

s0= s+10.*log10(sin(inc_angr[irange]));

return s0;

}

double sigma0_vv2hh(double sigma0, double thetar, double alpha)

{

double tmp;

if(alpha>=0 && alpha<=1.) /* Thompson Model */

{

tmp=(1.+alpha*tan(thetar)*tan(thetar))

/(1.+2.*tan(thetar)*tan(thetar));

tmp=pow(tmp, 2);

sigma0 *= tmp;

}

else if(alpha==2.) /* Elfouhaily Model */

{

tmp=(1.+2.*sin(thetar)*sin(thetar))

/(1.+2.*tan(thetar)*tan(thetar));

tmp=pow(tmp, 2);

sigma0 *= tmp;

}

else

{

printf("  Wrong input for VV2HH\n");

exit(0);

}

return(sigma0);

}

double cmod4_v2sigma(double w_speed, double w_direction, double incidence)

{

double alpha, gamma, beta, b0, b1, b2, b3;

double c1, c2, c3, c4, c5, c6, c7, c8, c9, c10;

double c11, c12, c13, c14, c15, c16, c17, c18;

double br[61], p0, p1, p2;

double theta, x, f1, f2, y;

int itheta;

double v, phi, phir, tmp, sigma0, thetar;

v=w_speed;

phi=w_direction;

theta=incidence;

phir=phi*pi/180.;

thetar=theta*pi/180.;

br[16]=1.075; br[17]=1.075; br[18]=1.075; br[19]=1.072;

br[20]=1.069; br[21]=1.066; br[22]=1.056; br[23]=1.030;

br[24]=1.004; br[25]=0.979; br[26]=0.967; br[27]=0.958;

br[28]=0.949; br[29]=0.941; br[30]=0.934; br[31]=0.927;

br[32]=0.923; br[33]=0.930; br[34]=0.937; br[35]=0.944;

br[36]=0.955; br[37]=0.967; br[38]=0.978; br[39]=0.988;

br[40]=0.998; br[41]=1.009; br[42]=1.021; br[43]=1.033;

br[44]=1.042; br[45]=1.050; br[46]=1.054; br[47]=1.053;

br[48]=1.052; br[49]=1.047; br[50]=1.038; br[51]=1.028;

br[52]=1.016; br[53]=1.002; br[54]=0.989; br[55]=0.965;

br[56]=0.941; br[57]=0.929; br[58]=0.929; br[59]=0.929;

br[60]=0.929;

c1=-2.301523; c2=-1.632686; c3=0.761210;

c4=1.156619; c5=0.595955; c6=-0.293819;

c7=-1.015244; c8=0.342175; c9=-0.500786;

c10=0.014430; c11=0.002484; c12=0.074450; c13=0.004023;

c14=0.148810; c15=0.089286;

c16=-0.006667; c17=3.; c18=-10.;

x=(theta-40.)/25.;

p0=1.;

p1=x;

p2=(3.*x*x-1.)/2.;

alpha=c1*p0+c2*p1+c3*p2;

gamma=c4*p0+c5*p1+c6*p2;

beta=c7*p0+c8*p1+c9*p2;

f2=atan(2.5*(x+0.35))*180./pi - 0.61*(x+0.35);

b1=c10*p0+c11*v +(c12*p0+c13*v)*f2;

b2=c14*p0+c15*(1.+p1)*v;

b3=0.42*(1.+c16*(c17+x)*(c18+v));

y=v+beta;

if(y<=1.e-10) f1=-10.;

else if(y>5.) f1=sqrt(y)/3.2;

else f1=log10(y);

itheta=(int)(theta+0.5);

b0=br[itheta]*pow(10., alpha+gamma*f1);

tmp=1.+b1*cos(phir)+b3* (atan(b2)*180./pi) * cos(2*phir);

tmp=fabs(tmp);

tmp=pow(tmp, 1.6);

sigma0=b0*tmp;

/* conversion from VV to HH */

sigma0=sigma0_vv2hh(sigma0, thetar, vv2hh_alpha);

sigma0 = 10.*log10(sigma0);

return(sigma0);

}

double cmodifr2_v2sigma(double w_speed, double w_direction, double incidence)

{

double alpha, beta, b1, b2, v, phi;

double c1, c2, c3, c4, c5, c6, c7, c8, c9, c10;

double c11, c12, c13, c14, c15, c16, c17, c18;

double c19, c20, c21, c22, c23, c24, c25;

double p1, p2, p3;

double v1, v2, v3;

double theta, x;

double phir, tmp, sigma0, thetar;

v=w_speed;

phi=w_direction;

theta=incidence;

phir=phi*pi/180.;

thetar=theta*pi/180.;

c1=-2.437597; c2=-1.567031; c3=0.370824; c4=-0.04059;

c5=0.404678; c6=0.188397; c7=-0.027262; c8=0.06465;

c9=0.0545; c10=0.08635; c11=0.0551; c12=-0.05845;

c13=-0.0961; c14=0.412754; c15=0.121785; c16=-0.024333;

c17=0.072163; c18=-0.062954; c19=0.015958; c20=-0.069514;

c21=-0.062945; c22=0.035538; c23=0.023049; c24=0.074654;

c25=-0.014713;

x=(theta-36.)/19.;

p1=x;

p2=(3.*x*x-1)/2.;

p3=x*(5.*x*x-3)/2.;

alpha=c1+c2*p1+c3*p2+c4*p3;

beta=c5+c6*p1+c7*p2;

x=(2.*theta-76)/40.;

p1=x;

p2=2.*x*x-1.;

v1=(2.*v-28.)/22.;

v2=2.*v1*v1-1.;

v3=(2.*v2-1.)*v1;

b1=c8+c9*v1+c10*p1+c11*p1*v1+c12*p2+c13*p2*v1;

b2=c14+c15*p1+c16*p2+ (c17+c18*p1+c19*p2)*v1 +

(c20+c21*p1+c22*p2)*v2+(c23+c24*p1+c25*p2)*v3;

tmp=pow(10., alpha+beta*sqrt(v));

sigma0=tmp*fabs(1.+b1*cos(phir)+(atan(b2)*180./pi)*cos(2.*phir));

/* conversion from VV to HH */

sigma0=sigma0_vv2hh(sigma0, thetar, vv2hh_alpha);

sigma0 = 10.*log10(sigma0);

return(sigma0);

}

/* creat LUT for cmod4_lut[sigma0 in 1/100 dB]= wind speed in cm/s */

void cmod4_LUT()

{

double vmin, vmax, v;

double phi, tmp;

int i, ix, j, nlut_sig, nlut_ran, isigma;

printf("   CMOD4 LUT: \n");

phi=wind_direction;

/* find min and max wind speed (vmax, vmin) */

vmin=100.;

vmax=0.;

tmp=0.;

do

{

vmax += 1.;

tmp=cmod4_v2sigma(vmax, phi, inc_ang_max);

/*	printf("v=%g, phi=%g, theta=%g, sigma0=%g\n",

vmax, phi, inc_ang_max, tmp); */

}while(tmp<sigma0_max);

do

{

vmin -= 1.;

tmp=cmod4_v2sigma(vmin, phi, inc_ang_min);

/*	printf("v=%g, phi=%g, theta=%g, sigma0=%g\n",

vmin, phi, inc_ang_min, tmp); */

}while(tmp>sigma0_min && vmin>0);

/*printf("     vmin=%g. vmax=%g\n", vmin, vmax);*/

sig_max_cmod4=cmod4_v2sigma(vmax, phi, inc_ang_min);

sig_min_cmod4=cmod4_v2sigma(vmin, phi, inc_ang_max);

nlut_ran=nrange_ml/inc_step+1;

nlut_sig=(int)((sig_max_cmod4 - sig_min_cmod4)*10. + 1);

printf("     LUT[range, sigma0]=[%d, %d]\n",nlut_ran, nlut_sig);

printf("     Incidence angle min, max(deg) = %g, %g\n",

inc_ang_min, inc_ang_max);

if(asc_des==-1) printf("       far range first\n");

else printf("       near range first\n");

printf("     LUT sigma0 min, max(dB)= %g, %g\n",

sig_min_cmod4, sig_max_cmod4);

printf("     LUT wind min, max(m/s)= %g, %g\n", vmin, vmax);

/*	printf("     LUT range, index, sigma0(dB), windspeed(m/s)\n");

*/

cmod4_lut=dmatrix(nlut_ran, nlut_sig);

for(i=0; i<nlut_ran; i++)

{

ix=i*inc_step;

j=0;

for(v=vmin; v<vmax; v+=0.1)

{

tmp=cmod4_v2sigma(v, phi, inc_ang[ix]);

isigma=(int)((tmp-sig_min_cmod4)*10);

/*printf("isigma=%d\n", isigma);*/

while(j<=isigma)

{

cmod4_lut[i][j]=v;

/*printf(" %d, %d, %g, %g\n",

ix, j, j/10.+sig_min_cmod4, v);

*/

j++;

}

}

}

}

/* creat LUT for cmodifr2_lut[sigma0 in 1/100 dB]= wind speed in cm/s */

void cmodifr2_LUT()

{

double vmin, vmax, v;

double phi, theta, tmp;

int i, ix, j, nlut_ran, nlut_sig, isigma;

printf("   CMOD_IFR2 LUT: \n");

phi=wind_direction;

/* find max and min wind speed (vmax, vmin) */

vmin=100.;

vmax=0.;

tmp=0.;

do

{

vmax += 1.;

tmp=cmodifr2_v2sigma(vmax, phi, inc_ang_max);

/*	printf("v=%g, phi=%g, theta=%g, sigma0=%g\n",

vmax, phi, inc_ang_max, tmp); */

}while(tmp<sigma0_max);

do

{

vmin -= 1.;

tmp=cmodifr2_v2sigma(vmin, phi, inc_ang_min);

/*	printf("v=%g, phi=%g, theta=%g, sigma0=%g\n",

vmin, phi, inc_ang_min, tmp); */

}while(tmp>sigma0_min && vmin>0);

/*printf("    vmin=%g. vmax=%g\n", vmin, vmax); */

sig_max_cmodifr2=cmodifr2_v2sigma(vmax, phi, inc_ang_min);

sig_min_cmodifr2=cmodifr2_v2sigma(vmin, phi, inc_ang_max);

nlut_ran=nrange_ml/inc_step+1;

nlut_sig=(int)((sig_max_cmodifr2-sig_min_cmodifr2)*10.+ 1);

printf("     LUT[range, sigma0]=[%d, %d]\n",nlut_ran, nlut_sig);

printf("     Incidence angle min, max(deg) = %g, %g\n",

inc_ang_min, inc_ang_max);

if(asc_des==-1) printf("       far range first\n");

else printf("       near range first\n");

printf("     LUT sigma0 min, max(dB)= %g, %g\n",

sig_min_cmod4, sig_max_cmod4);

printf("     LUT wind min, max(m/s)= %g, %g\n", vmin, vmax);

/*	printf("     LUT range, index, sigma0(dB), windspeed(m/s)\n");

*/

cmodifr2_lut=dmatrix(nlut_ran, nlut_sig);

for(i=0; i<nlut_ran; i++)

{

ix=i*inc_step;

j=0;

for(v=vmin; v<vmax; v+=0.1)

{

tmp=cmodifr2_v2sigma(v, phi, inc_ang[ix]);

isigma=(int)((tmp-sig_min_cmodifr2)*10.);

while(j<=isigma)

{

cmodifr2_lut[i][j]=v;

/*printf(" %d, %d, %g, %g\n",

ix, j, j/10.+sig_min_cmodifr2, v);

*/

j++;

}

}

}

}

double cmod4(double sigma0, int j)

{

int isigma;

isigma=(int)((sigma0-sig_min_cmod4)*10.);

if(isigma<0) return(0);

else return(cmod4_lut[j/inc_step][isigma]);

}

double cmodifr2(double sigma0, int j)

{

int isigma, jx;

double v, vbias;

isigma=(int)((sigma0-sig_min_cmodifr2)*10);

if(isigma<0) v=0.;

else v=cmodifr2_lut[j/inc_step][isigma];

/* High Wind Bias Correction */

if(v<=10.) vbias=0.;

else if(v>22.) vbias=atan(v-22.)*180./pi + 3.0382;

else vbias=0.0831*v-0.0173*v*v+0.0009*v*v*v;

v += vbias;

return(v);

}

7. Wave by SAR-Wave Spectrum: wave.c

/*****************************************************

wave.c

SAR Wave Spectrum

Input: SAR Multilook

Output: SAR Wave Spectrum

by Hoonyol Lee

Department of Geophysics, Kangwon National University

23 July 2006

22 Nov 2006: Split wave and current

******************************************************/

void wave()

{

FILE *fp_ml, *fp_wave;

int i, j, k, l, kx, kxx, kx1, m, ntotal;

float *fbuf, *fbuf1;

double *dbuf;

double tmp, tmp1, tmp2, tmp3;

double del_f;

printf("\n SAR Wave Spectrum\n");

del_f=prf/nptyw;

ntotal=nptxw*nptyw;

fbuf=fvector(nptyw*nrange_ml);

dbuf=dvector(nptxw*nptyw*2); /* complex */

if((fp_ml=fopen(fname_ml, "r"))==NULL)

{

printf(" (Error) Cannot open %s.\n\a", fname_ml);

exit(0);

}

if((fp_wave=fopen(fname_wave, "w"))==NULL)

{	printf(" (Error) Cannot open %s.\n\a", fname_wave);

exit(0);

}

/* Sarwave processing */

fbuf1=fvector(nptyw*nrange_ml);

for(i=0; i<(nazimuth_ml-nptyw+1); i+= nptyw)

{

printf("   %d/%d block\n", i/nptyw+1, nazimuth_ml/nptyw+1);

fseek(fp_ml, i*nrange_ml*4, SEEK_SET);

fread(fbuf, 4, nptyw*nrange_ml, fp_ml);

for(j=0; j<(nrange_ml-nptxw+1); j+=nptxw)

{

/* read multilook data */

for(k=0; k<nptyw; k++)

{

kx=k*nrange_ml+j;

kxx=k*nptxw*2;

for(l=0; l<nptxw; l++)

{

dbuf[kxx+l*2]= (double)fbuf[kx+l];

dbuf[kxx+l*2+1]=0.;

}

}

/*printf("      read ml ok\n");*/

/* 2D FFT and shift of Multilook Data*/

cfft2d(dbuf, nptyw, nptxw, -1);

cshift2d(dbuf, nptyw, nptxw);

/*printf("      2D FFT and shift ok\n");*/

/* calculate magnitude of 2D-FFTed ML */

for(k=0; k<nptyw; k++)

{

kx=k*nrange_ml+j;

kxx=k*nptxw*2;

for(l=0; l<nptxw; l++)

{

tmp1 = dbuf[kxx+2*l];

tmp2 = dbuf[kxx+2*l+1];

tmp3 = sqrt(tmp1*tmp1 + tmp2*tmp2);

fbuf1[kx+l]= (float)tmp3;

fbuf1[kx+l] /= ntotal;

}

}

/*printf("      mag ok\n");*/

} /* end of j */

fwrite(fbuf1, 4, nptyw*nrange_ml, fp_wave);

} /* end of i */

free(fbuf);

free(fbuf1);

/* final null output */

printf("   final null block\n");

fbuf1=fvector(nrange_ml);

for(i=((nazimuth_ml-nptyw)/nptyw)*nptyw; i<nazimuth_ml; i++)

fwrite(fbuf1, 4, nrange_ml, fp_wave);

free(fbuf1);

/*printf("      final null ok\n");*/

fclose(fp_ml);

fclose(fp_wave);

gen_ers(fname_wave, nazimuth_ml, nrange_ml, dazimuth_ml, drange_ml, 1, "141");

free(dbuf);

}

8. Wind from Interlook Cross Spectrum: wave_crs.c

/*****************************************************

wave_crs.c

Input: 2-byte short integer SLC in BIL file format

Process: Multilook and Cross Spectrum

Output: RESAMPLED multilook magnitude (2-byte unsigned integer) and Cross Spectrum files

by Hoonyol Lee

Department of Geophysics, Kangwon National University

8 August 2006

******************************************************/

void wave_crs()

{

FILE *fp_in, *fp_out, *fp_crs;

int i, j, k, m, n, jx, kx, kxx, mx, mxx, nblock, block_byte, nblock_fft;

short *sbuf, *usbuf;

double *dbuf, *dbuf_s1, *dbuf_s2, *dbuf_s3, *dbuf_s4, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6;

double *dbuf_s1x, *dbuf_s2x, *dbuf_s3x, *dbuf_s4x;

double dc, sub_bw,  del_f;

int sub_bw_px, sub_bw_px_slap, dc_px, dc_pxs;

int seek_offset, read_count, write_start, write_size_az, write_count;

int nlook;

float *fbuf;

double *taper;

printf("\n Inter-Look Cross Spectrum Process\n");

nlook=4; /* 4-look fixed*/

nblock=nptyw*nlook;

if(nblock<NBLOCK) nblock=NBLOCK;

printf("  nblock=%d \n", nblock);

sub_bw=azimuth_bw/nlook;

nblock_fft=nblock*2;

del_f=prf/nblock_fft;

sub_bw_px= (int)(sub_bw/del_f);  /* sub-aperture bandwidth pixel */

tmp1=(double)(sub_bw_px)*0.1;

sub_bw_px_slap=(int)tmp1; /* half of band side-lap in pixel */

block_byte=nblock_fft*nrange*2;

if((fp_in=fopen(fname_slc, "r"))==NULL)

{

printf(" (Error) Cannot open %s.\n\a", fname_slc);

exit(0);

}

if((fp_out=fopen(fname_crs_ml, "w"))==NULL)

{	printf(" (Error) Cannot open %s.\n\a", fname_crs_ml);

exit(0);

}

if((fp_crs=fopen(fname_crs, "w"))==NULL)

{	printf(" (Error) Cannot open %s.\n\a", fname_crs);

exit(0);

}

/* Hann window taper design */

taper=dvector(nblock_fft/nlook);

for(i=0; i<sub_bw_px_slap; i++)

taper[i]=0.5-0.5*cos(pi*i/sub_bw_px_slap);

for(i=sub_bw_px_slap; i<(sub_bw_px_slap+sub_bw_px); i++)

taper[i]=1.;

for(i<(sub_bw_px_slap+sub_bw_px), j=0; i<(sub_bw_px_slap*2+sub_bw_px);

i++, j++) taper[i]=0.5+0.5*cos(pi*j/sub_bw_px_slap);

/* Block Processing */

for(i=0; i<nazimuth; i+=nblock)

{

printf("   %d/%d block\n", i/nblock+1, nazimuth/nblock+1);

if(i==0) /* first block */

{

printf("    first block \n");

seek_offset = 0;

read_count=(nblock+nblock/2)*nrange*2;

write_start=0;

write_size_az=nblock/nlook;

write_count=nblock/nlook*nrange;

}

else if((i+nblock+nblock/2)>=nazimuth) /* last block */

{

printf("    last block\n");

seek_offset = (i - nblock/2)*nrange*2*2;

read_count=(nblock/2+nazimuth-i)*nrange*2;

write_start=nblock/2/nlook;

write_size_az=(nazimuth-i)/nlook;

write_count=(nazimuth-i)/nlook*nrange;

}

else /* mid block */

{

seek_offset = (i - nblock/2)*nrange*2*2;

read_count= nblock_fft*nrange*2;

write_start=nblock/2/nlook;

write_size_az=nblock/nlook;

write_count=nblock/nlook*nrange;

}

sbuf=svector(block_byte);  /* input */

dbuf=dvector(block_byte); /* complex */

dbuf_s1=dvector(block_byte/nlook); /* look 1*/

dbuf_s2=dvector(block_byte/nlook); /* look 2 */

dbuf_s3=dvector(block_byte/nlook); /* look 3 */

dbuf_s4=dvector(block_byte/nlook); /* look 4 */

fbuf=fvector(block_byte/nlook); /* output */

/* read slc */

fseek(fp_in, seek_offset, SEEK_SET);

fread(sbuf, 2, read_count, fp_in);

for(j=0; j<read_count; j++) dbuf[j]=(double)sbuf[j];

bil2bipn(dbuf, nblock_fft, nrange, 2);

printf("      read slc ok\n");

/* Forward FFT azimuth */

cfft1d(dbuf, nblock_fft, nrange, 1, -1);

printf("      Az FFT ok\n");

/* Split Aperture */

for(j=0; j<nrange; j++)

{

if(asc_des) m=j; /* ascending, right looking */

else m=nrange-j-1; /* descending, right looking */

dc=Dc0+Dc1*m+Dc2*m*m;

dc=fmod(dc, prf);

if(asc_des==-1) dc += prf/2.;

if(dc>=(0.5*prf)) dc -= prf;

if(dc<(-0.5*prf)) dc += prf;

dc_px=(int)dc/del_f;

if(dc_px<=0) dc_px *= -1;

else dc_px=nblock_fft-dc_px;

for(k=0; k<(sub_bw_px+sub_bw_px_slap*2); k++)

{

/* look 1 */

dc_pxs=dc_px - sub_bw_px*2 - sub_bw_px_slap + k;

while(dc_pxs<0) dc_pxs += nblock_fft;

while(dc_pxs>=nblock_fft) dc_pxs -= nblock_fft;

kx=k*nrange*2+j*2;

kxx=dc_pxs*nrange*2+j*2;

dbuf_s1[kx]=dbuf[kxx]*taper[k];

dbuf_s1[kx+1]=dbuf[kxx+1]*taper[k];

/* look 2 */

dc_pxs=dc_px - sub_bw_px - sub_bw_px_slap + k;

while(dc_pxs<0) dc_pxs += nblock_fft;

while(dc_pxs>=nblock_fft) dc_pxs -= nblock_fft;

kxx=dc_pxs*nrange*2+j*2;

dbuf_s2[kx]=dbuf[kxx]*taper[k];

dbuf_s2[kx+1]=dbuf[kxx+1]*taper[k];

/* look 3 */

dc_pxs=dc_px - sub_bw_px_slap + k;

while(dc_pxs<0) dc_pxs += nblock_fft;

while(dc_pxs>=nblock_fft) dc_pxs -= nblock_fft;

kxx=dc_pxs*nrange*2+j*2;

dbuf_s3[kx]=dbuf[kxx]*taper[k];

dbuf_s3[kx+1]=dbuf[kxx+1]*taper[k];

/* look 4 */

dc_pxs=dc_px + sub_bw_px - sub_bw_px_slap + k;

while(dc_pxs<0) dc_pxs += nblock_fft;

while(dc_pxs>=nblock_fft) dc_pxs -= nblock_fft;

kxx=dc_pxs*nrange*2+j*2;

dbuf_s4[kx]=dbuf[kxx]*taper[k];

dbuf_s4[kx+1]=dbuf[kxx+1]*taper[k];

}/* end of k */

} /* end of j */

printf("      split aperture ok\n");

/* Inverse FFT azimuth */

cfft1d(dbuf_s1, nblock_fft/nlook, nrange, 1, 1);

cfft1d(dbuf_s2, nblock_fft/nlook, nrange, 1, 1);

cfft1d(dbuf_s3, nblock_fft/nlook, nrange, 1, 1);

cfft1d(dbuf_s4, nblock_fft/nlook, nrange, 1, 1);

printf("      Az iFFT ok\n");

/* detect sub-aperture images */

for(j=write_start*nrange, k=0; j<(write_start+write_size_az)*nrange; j++, k++)

{

tmp4=0.;

tmp1= dbuf_s1[j*2];

tmp2= dbuf_s1[j*2+1];

tmp3= sqrt(tmp1*tmp1+tmp2*tmp2);

dbuf_s1[k*2]=tmp3;

dbuf_s1[k*2+1]=0.;

tmp4 += tmp3;

tmp1= dbuf_s2[j*2];

tmp2= dbuf_s2[j*2+1];

tmp3= sqrt(tmp1*tmp1+tmp2*tmp2);

dbuf_s2[k*2]=tmp3;

dbuf_s2[k*2+1]=0.;

tmp4 += tmp3;

tmp1= dbuf_s3[j*2];

tmp2= dbuf_s3[j*2+1];

tmp3= sqrt(tmp1*tmp1+tmp2*tmp2);

dbuf_s3[k*2]=tmp3;

dbuf_s3[k*2+1]=0.;

tmp4 += tmp3;

tmp1= dbuf_s4[j*2];

tmp2= dbuf_s4[j*2+1];

tmp3= sqrt(tmp1*tmp1+tmp2*tmp2);

dbuf_s4[k*2]=tmp3;

dbuf_s4[k*2+1]=0.;

tmp4 += tmp3;

fbuf[k]=(float)(tmp4/nlook);

}

printf("      detect multilook ok\n");

/* write multilook image */

fwrite(fbuf, 4, write_count, fp_out);

printf("      write multilook ok\n");

/* Cross Spectrum */

for(j=0; j<write_size_az; j+= nptyw)

{

for(k=0; k<(nrange-nptxw+1); k+=nptxw)

{

dbuf_s1x=dvector(nptxw*nptyw*2);

dbuf_s2x=dvector(nptxw*nptyw*2);

dbuf_s3x=dvector(nptxw*nptyw*2);

dbuf_s4x=dvector(nptxw*nptyw*2);

for(m=0; m<nptyw; m++)

{

mx=m*nptxw*2;

mxx=(j+m)*nrange*2+k*2;

for(n=0; n<nptxw*2; n++)

{

dbuf_s1x[mx+n]=dbuf_s1[mxx+n];

dbuf_s2x[mx+n]=dbuf_s2[mxx+n];

dbuf_s3x[mx+n]=dbuf_s3[mxx+n];

dbuf_s4x[mx+n]=dbuf_s4[mxx+n];

}

}

/* Forward 2D FFT and shift */

cfft2d(dbuf_s1x, nptyw, nptxw, -1);

cfft2d(dbuf_s2x, nptyw, nptxw, -1);

cfft2d(dbuf_s3x, nptyw, nptxw, -1);

cfft2d(dbuf_s4x, nptyw, nptxw, -1);

cshift2d(dbuf_s1x, nptyw, nptxw);

cshift2d(dbuf_s2x, nptyw, nptxw);

cshift2d(dbuf_s3x, nptyw, nptxw);

cshift2d(dbuf_s4x, nptyw, nptxw);

for(m=0; m<nptyw; m++)

{

mx=m*nptyw*2;

for(n=0; n<nptxw; n++)

{

/* look1 and look 3 */

tmp1=dbuf_s1x[mx+n*2];

tmp2=dbuf_s1x[mx+n*2+1];

tmp3=dbuf_s3x[mx+n*2];

tmp4=dbuf_s3x[mx+n*2+1];

tmp5=(tmp1*tmp3-tmp2*tmp4)/2.;

tmp6=(tmp1*tmp4+tmp3*tmp2)/2.;

/* look 2 and look 4 */

tmp1=dbuf_s2x[mx+n*2];

tmp2=dbuf_s2x[mx+n*2+1];

tmp3=dbuf_s4x[mx+n*2];

tmp4=dbuf_s4x[mx+n*2+1];

tmp5 += (tmp1*tmp3-tmp2*tmp4)/2.;

tmp6 += (tmp1*tmp4+tmp3*tmp2)/2.;

mxx=(j+m)*nrange*2+k*2;

if(c_wave_crsflag==2)

{

/* I and Q */

dbuf_s1[mxx+n*2]= tmp5;

dbuf_s1[mxx+n*2+1]=tmp6;

dbuf_s1[mxx+n*2] /= nptxw*nptyw;

dbuf_s1[mxx+n*2+1] /= nptxw*nptyw;

}

else if(c_wave_crsflag==1)

{

/* amplitude and phase */

dbuf_s1[mxx+n*2]=sqrt(tmp5*tmp5+tmp6*tmp6);

dbuf_s1[mxx+n*2+1]=atan2(tmp6, tmp5);

dbuf_s1[mxx+n*2] /= nptxw*nptyw;

dbuf_s1[mxx+n*2+1] /= nptxw*nptyw;

}

}

}

free(dbuf_s1x);

free(dbuf_s2x);

free(dbuf_s3x);

free(dbuf_s4x);

} /* end of k */

} /* end of j */

printf("      cross spec ok\n");

/* fill the null right */

for(j=0; j<write_size_az; j++)

{

jx=j*nrange*2;

for(k=(nrange/nptxw*nptxw*2); k<(nrange*2); k++)

dbuf_s1[jx+k]=0.;

}

bip2biln(dbuf_s1, write_size_az, nrange, 2);

for(j=0; j<write_count*2; j++) fbuf[j]=(float)dbuf_s1[j];

/* write cross spectrum */

if(c_wave_crsflag) fwrite(fbuf, 4, write_count*2, fp_crs);

printf("      write crs ok\n");

free(sbuf);

free(dbuf);

free(dbuf_s1);

free(dbuf_s2);

free(dbuf_s3);

free(dbuf_s4);

} /* end of i */

/* fill-in the null bottom */

/**	usbuf=usvector(nrange*2);

for(i=(nazimuth/nblock*nblock/nlook); i<(nazimuth/nlook); i++)

{

fwrite(usbuf, 2, nrange, fp_out);

fwrite(usbuf, 2, nrange*2, fp_crs);

}

free(usbuf);

**/

free(fbuf);

free(taper);

fclose(fp_in);

fclose(fp_out);

fclose(fp_crs);

gen_ers(fname_crs_ml, nazimuth/nlook, nrange, dazimuth*nlook, drange*nx_res, 1, "141");

gen_ers_crs(fname_crs, nazimuth/nlook, nrange, dazimuth*nlook, drange*nx_res, 2, "141");

}

9. Current by Doppler Shift Analysis: current.c

/*****************************************************

current.c

Doppler Shift for Current

Input: SAR SLC

Output: Doppler Shift (AP or IQ)

by Hoonyol Lee

Department of Geophysics, Kangwon National University

23 July 2006

22 Nov 2006: Split from sarwave

******************************************************/

void current()

{

FILE *fp_slc, *fp_cur, *fp_erv, *fp_vec, *fp_txt;

char ervname[200];

int i, j, jx, jx1, k, m;

short *sbuf;

float *fbuf, *fbuf1;

double *dbuf_slc;

double tmp, tmp1, tmp2, tmp3;

double dc, del_f, ampmax, dc_cur, del_dc, u;

int dc_px, dc_cur_px;

double tlx, tly, brx, bry;

double x, y, x1, y1, y2;

int ntmp, icurend;

del_f=prf/nptyc;

tlx=nrange_ml;

tly=nazimuth_ml;

brx=0.;

bry=0.;

printf("\n Current from Doppler Shift\n");

fbuf=fvector(nptyc*nrange_ml);

fbuf1=fvector(nrange_ml);

sbuf=svector(nptyc*ny_res*nrange*2);

dbuf_slc=dvector(nptyc*ny_res*nrange*2); /* complex */

if((fp_slc=fopen(fname_slc, "r"))==NULL)

{

printf(" (Error) Cannot open %s.\n\a", fname_slc);

exit(0);

}

if((fp_cur=fopen(fname_cur, "w"))==NULL)

{	printf(" (Error) Cannot open %s.\n\a", fname_cur);

exit(0);

}

if((fp_erv=fopen(fname_cur_erv, "w"))==NULL)

{	printf(" (Error) Cannot open %s.\n\a", fname_cur_erv);

exit(0);

}

if((fp_vec=fopen(fname_cur_vec, "w"))==NULL)

{	printf(" (Error) Cannot open %s.\n\a", fname_cur_vec);

exit(0);

}

if((fp_txt=fopen(fname_cur_txt, "w"))==NULL)

{	printf(" (Error) Cannot open %s.\n\a", fname_cur_txt);

exit(0);

}

/* Doppler processing */

fprintf(fp_txt,

"x(range) y(azimuth) Doppler_Shift(Hz) Velocity(m/s)\n");

icurend=0;

for(i=0; i<(nazimuth_ml-nptyc + 1); i+= dc_ystep)

{

printf("   %d/%d blocks\n", i/dc_ystep+1, nazimuth_ml/dc_ystep+1);

/* read slc */

ntmp=nptyc*ny_res*nrange*2;

fseek(fp_slc, i*ny_res*nrange*2*2, SEEK_SET);

fread(sbuf, 2, ntmp, fp_slc);

for(j=0; j<ntmp; j++) dbuf_slc[j]=(double)sbuf[j];

bil2bipn(dbuf_slc, nptyc*ny_res, nrange, 2);

printf("      read slc ok\n");

/* 1D Azimuth FFT and shift of SLC Data */

cfft1d(dbuf_slc, nptyc*ny_res, nrange, 1, -1);

cshift_y(dbuf_slc, nptyc*ny_res, nrange);

printf("      Az FFT ok\n");

/* calculate magnitude of 1D Azimuth-FFTed SLC */

ntmp=nptyc*ny_res*nrange;

for(j=0; j<ntmp; j++)

{

tmp1 = dbuf_slc[j*2];

tmp2 = dbuf_slc[j*2+1];

dbuf_slc[j] = sqrt(tmp1*tmp1+tmp2*tmp2);

}

printf("      detect ok\n");

/* multilook 1D-FFTed magnitude */

maverage(dbuf_slc, nptyc*ny_res, nrange, 1, nx, ny);

printf("      mlook ok\n");

/* resample the 1D-FFTed magnitude */

for(j=0; j<nptyc; j++)

{

jx=j*nrange;

jx1=jx*ny_res;

for(k=0; k<nrange; k++)

dbuf_slc[jx+k]=dbuf_slc[jx1+k];

}

printf("      resample ok\n");

/* average of 1D-FFTed magnitude */

maverage(dbuf_slc, nptyc, nrange, 1, dc_avx, dc_avy);

printf("      avg ok\n");

/* Doppler Shift Estimation for Current */

for(j=0; j<nrange; j+= dc_xstep)

{

/* Nominal Doppler Centroid */

if(asc_des) m=j; /* ascending, right looking */

else m=nrange-j-1;  /* descending, right looking */

dc=Dc0+Dc1*m+Dc2*m*m;

dc=fmod(dc, prf);

if(asc_des==-1) dc += prf/2.;

if(dc>=(0.5*prf)) dc -= prf;

if(dc<(-0.5*prf)) dc += prf;

dc_px= nptyc/2 - (int)(dc/del_f);

/*printf("Nominal Doppler Centroid pixel=%d\n",

dc_px);

*/

/* Doppler Centroid Estimation */

ampmax=0.;

for(k=0; k<nptyc; k++)

{

tmp=dbuf_slc[k*nrange+j];

if(tmp>ampmax)

{

ampmax=tmp;

dc_cur_px=k;

}

}

dc_cur= (nptyc/2-dc_cur_px)*del_f;

del_dc=dc_cur - dc;

if(del_dc>(prf/2.)) del_dc -= prf;

if(del_dc<(prf/(-2.))) del_dc += prf;

u=del_dc*wavelength/2./pi/sin(inc_angr[m]);

/*printf("%d %d %d %d %lf %lf\n",

j, i+nptyc/2, dc_px, dc_cur_px, del_dc, u);

*/

/* output txt */

fprintf(fp_txt, "%d %d %lf %lf\n",

j, i+nptyc/2, del_dc, u);

/* output vector files */

x=j*drange_ml;

y=(i+nptyc/2)*dazimuth_ml;

x1=x-u*vel_scale*asc_des;

if(i>=icurend)

{

y1=(i+dc_px)*dazimuth_ml;

y2=(i+dc_cur_px)*dazimuth_ml;

fprintf(fp_erv,

"point(,%lf,%lf,0,0,0,255,0).\n",

x, y1);

fprintf(fp_erv,

"point(,%lf,%lf,0,255,0,0,0).\n",

x, y2);

if(y1<tly) tly=y1;

if(y1>bry) bry=y1;

if(y2<tly) tly=y2;

if(y2>bry) bry=y2;

}

fprintf(fp_vec,

"poly(,2,[%lf,%lf,%lf,%lf],1,1,0,0,0,0,255,0,0,0).\n", x, y, x1, y);

if(x<tlx) tlx=x;

if(x>brx) brx=x;

if(y<tly) tly=y;

if(y>bry) bry=y;

if(x1<tly) tly=x1;

if(x1>bry) bry=x1;

} /* end of j */

/* output doppler image */

if(i>=icurend)

{

for(j=0; j<(i-icurend); j++) fwrite(fbuf1, 4, nrange_ml, fp_cur);

for(j=0; j<nptyc*nrange_ml; j++) fbuf[j]=(float)dbuf_slc[j];

fwrite(fbuf, 4, nptyc*nrange_ml, fp_cur);

icurend = i+nptyc;

}

printf("      doppler ok\n");

} /* end of i */

free(fbuf);

free(fbuf1);

/* final null output */

fbuf=fvector(nrange_ml);

for(i=icurend; i<nazimuth_ml; i++)

fwrite(fbuf, 4, nrange_ml, fp_cur);

free(fbuf);

fclose(fp_slc);

fclose(fp_cur);

fclose(fp_erv);

fclose(fp_vec);

fclose(fp_txt);

gen_ers(fname_cur, nazimuth_ml, nrange_ml, dazimuth_ml, drange_ml, 1, "141");

strcpy(ervname, fname_cur_erv);

strcat(ervname, ".erv");

generv(ervname, tlx, tly, brx, bry);

strcpy(ervname, fname_cur_vec);

strcat(ervname, ".erv");

generv(ervname, tlx, tly, brx, bry);

free(sbuf);

free(dbuf_slc);

}

10. BIP from/to BIL Data Conversion: bip2biln.c

/*****************************************************************

bip2biln

convert double BIP to BIL

bil2bipn

convert double BIL to BIP

data: double data

nrow: number of row

ncol: number of column

nband: number of band

nbyte: number of byte

by Dr. Hoonyol Lee

Dept. of Geophysics

Kangwon National Univeristy

29 June 2004

12 August 2006: bip2bil, bil2bip

*****************************************************************/

/*#include "ma.c" */

void bip2biln(data, nrow, ncol, nband)

double *data;

int nrow, ncol, nband;

{

int i, j, k, ix;

double *tmp;

tmp=dvector(ncol*nband);

for(i=0; i<nrow; i++)

{

ix=i*ncol*nband;

for(j=0; j<ncol; j++)

{

for(k=0; k<nband; k++)

tmp[j+k*ncol]=data[ix+j*nband+k];

}

for(j=0; j<ncol*nband; j++) data[ix+j]=tmp[j];

}

free(tmp);

}

void bil2bipn(data, nrow, ncol, nband)

double *data;

int nrow, ncol, nband;

{

int i, j, k, ix;

double *tmp;

tmp=dvector(ncol*nband);

for(i=0; i<nrow; i++)

{

ix=i*ncol*nband;

for(j=0; j<ncol; j++)

{

for(k=0; k<nband; k++)

tmp[j*nband+k]=data[ix+j+k*ncol];

}

for(j=0; j<ncol*nband; j++) data[ix+j]=tmp[j];

}

free(tmp);

}

11. Conversion from BIN to INT: btoi.c

/********************************************

btoi.c

convert nbyte binary unsinged char to 4-byte int

1. The bytes larger than sizeof(int) will be lost.

2. Byte order (MSB- or LSB-first) independent

3. Note the difference from atoi() in <stdlib.h>

which convert from asc to int.

by Dr. Hoonyol Lee

16 June 2004

*********************************************/

int btoi(unsigned char *data, int nbyte)

{

int i, idata, itmp;

idata=0;

for(i=0; i<nbyte; i++)

{

idata <<= 8;

idata |= *(data++);

}

return(idata);

}

12. Byte Swap: byteswap.c

/* byteswap.c */

void byteswap(void *buf, unsigned nbyte)

{

unsigned i;

char tmp;

for(i=0; i<nbyte/2; i++)

{

tmp=*((char *)buf+i);

*((char *)buf+i)=*((char *)buf+nbyte-i-1);

*((char *)buf+nbyte-i-1)=tmp;

}

}

13. 1D and 2D FFT: cfft.c

/**********************************

1D- and 2D-FFT of 2D Complex Data

H. Lee, 17/7/01

modified usgin ma.c, 7 May 2004

Normalized.

data: double complex in BIP

nrow, ncol:

iflag:  0 --- x directional FFT

1 --- y directional FFT

isign:  -1 --- Forward FFT

1 --- Inverse FFT

**********************************/

/*#include "ma.c"*/

#include "twon.c"

#include "dfour1.c"

void cfft1d(data, nrow, ncol, iflag, isign)

double *data;

int nrow, ncol, iflag, isign;

{

unsigned long nn;

int i,j, linen;

double *tmp, scale;

linen = ncol*2;

/* x-directional FFT */

if(iflag == 0)

{

tmp = dvector(ncol*2+1);

nn = twon(ncol);

for(i=0; i<nrow; i++){

for(j=0; j<ncol; j++){

tmp[2*j+1] = data[linen*i + 2*j];

tmp[2*j+2] = data[linen*i + 2*j + 1];

}

dfour1(tmp, nn, isign);

for(j=0; j<ncol; j++){

data[linen*i+2*j] = tmp[2*j+1];

data[linen*i+2*j+1] = tmp[2*j+2];

}

}

/* scaling at iFFT */

if(isign==1) for(i=0; i<nrow*ncol*2; i++) data[i] /= ncol;

}

/* y-directional FFT */

if(iflag == 1)

{

tmp = dvector(nrow*2+1);

nn = twon(nrow);

for(i=0; i< ncol; i++){

for(j=0; j<nrow; j++){

tmp[2*j+1] = data[linen*j+2*i];

tmp[2*j+2] = data[linen*j+2*i+1];

}

dfour1(tmp, nn, isign);

for(j=0; j<nrow; j++){

data[linen*j+2*i] = tmp[2*j+1];

data[linen*j+2*i+1] = tmp[2*j+2];

}

}

/* scaling at iFFT*/

if(isign==1) for(i=0; i<nrow*ncol*2; i++) data[i] /= nrow;

}

free(tmp);

}

void cfft2d(data, nrow, ncol, isign)

double *data;

int nrow, ncol, isign;

{

cfft1d(data, nrow, ncol, 0, isign);

cfft1d(data, nrow, ncol, 1, isign);

}

void cfft1d_ri(realdata, imagdata, nrow, ncol, iflag, isign)

double *realdata, *imagdata;

int nrow, ncol, iflag, isign;

{

int i;

double *data;

data=dvector(nrow*ncol*2);

for(i=0; i<nrow*ncol; i++)

{

data[2*i]= realdata[i];

data[2*i+1]=imagdata[i];

}

cfft1d(data, nrow, ncol, iflag, isign);

for(i=0; i<nrow*ncol; i++)

{

realdata[i]=data[2*i];

imagdata[i]=data[2*i+1];

}

free(data);

}

void cfft2d_ri(realdata, imagdata, nrow, ncol, isign)

double *realdata, *imagdata;

int nrow, ncol, isign;

{

cfft1d_ri(realdata, imagdata, nrow, ncol, 0, isign);

cfft1d_ri(realdata, imagdata, nrow, ncol, 1, isign);

}

14. Cyclic Data Shift: cnshift.c

/****************************************

cnshift.c

: cyclic shift of complex double data in BIP

by integer n

27 Oct 2006: cnshift

24 July 2006: cshift

by Hoonyol Lee

***************************************/

/*#include "ma.c"*/

void cnshift_x(double *data, int nrow, int ncol, int nshift)

{

int i, j, ix, jj, ncol2, nshift2;

double *tmp;

tmp=dvector(ncol*2);

ncol2=ncol*2;

nshift2=nshift*2;

for(i=0; i<nrow; i++)

{

ix=i*ncol2;

for(j=0; j<ncol2; j++)

{

jj=j+nshift2;

while(jj>=ncol2) jj-=ncol2;

while(jj<0) jj+=ncol2;

tmp[jj]=data[ix+j];

}

for(j=0; j<ncol2; j++) data[ix+j]=tmp[j];

}

free(tmp);

}

void cnshift_y(double *data, int nrow, int ncol, int nshift)

{

int i, j, jj, ncol2;

double *tmp;

tmp=dvector(nrow);

ncol2=ncol*2;

for(i=0; i<ncol2; i++)

{

for(j=0; j<nrow; j++)

{

jj=j+nshift;

while(jj>=nrow) jj-=nrow;

while(jj<0) jj+=nrow;

tmp[jj]=data[j*ncol2+i];

}

for(j=0; j<nrow; j++) data[j*ncol2+i]=tmp[j];

}

free(tmp);

}

void cnshift_xy(double *data, int nrow, int ncol, int xshift, int yshift)

{

cnshift_x(data, nrow, ncol, xshift);

cnshift_y(data, nrow, ncol, yshift);

}

15. Cyclic Data Half-Shift: cshift.c

/****************************************

cshift.c

: cyclic shift of complex double data in BIP

: only shift by half of nrow and ncol

24 July 2006

by Hoonyol Lee

***************************************/

void cshift_x(double *data, int nrow, int ncol)

{

int i, j, ix;

double tmp1, tmp2;

for(i=0; i<nrow; i++)

{

ix=i*ncol*2;

for(j=0; j<ncol/2; j++)

{

tmp1=data[ix+j*2];

tmp2=data[ix+j*2+1];

data[ix+j*2]=data[ix+j*2+ncol];

data[ix+j*2+1]=data[ix+j*2+ncol+1];

data[ix+j*2+ncol]=tmp1;

data[ix+j*2+ncol+1]=tmp2;

}

}

}

void cshift_y(double *data, int nrow, int ncol)

{

int i, j, jx;

double tmp1, tmp2;

for(i=0; i<ncol; i++)

{

for(j=0; j<nrow/2; j++)

{

jx=j*ncol*2;

tmp1=data[jx+i*2];

tmp2=data[jx+i*2+1];

data[jx+i*2]=data[jx+nrow*ncol+i*2];

data[jx+i*2+1]=data[jx+nrow*ncol+i*2+1];

data[jx+nrow*ncol+i*2]=tmp1;

data[jx+nrow*ncol+i*2+1]=tmp2;

}

}

}

void cshift2d(double *data, int nrow, int ncol)

{

cshift_x(data, nrow, ncol);

cshift_y(data, nrow, ncol);

}

16. FFT Subroutine: dfour1.c

/*#include <math.h>*/

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

void dfour1(double data[], unsigned long nn, int isign)

{

unsigned long n,mmax,m,j,istep,i;

double wtemp,wr,wpr,wpi,wi,theta;

double tempr,tempi;

n=nn << 1;

j=1;

for (i=1;i<n;i+=2) {

if (j > i) {

SWAP(data[j],data[i]);

SWAP(data[j+1],data[i+1]);

}

m=n >> 1;

while (m >= 2 && j > m) {

j -= m;

m >>= 1;

}

j += m;

}

mmax=2;

while (n > mmax) {

istep=mmax << 1;

theta=isign*(6.28318530717959/mmax);

wtemp=sin(0.5*theta);

wpr = -2.0*wtemp*wtemp;

wpi=sin(theta);

wr=1.0;

wi=0.0;

for (m=1;m<mmax;m+=2) {

for (i=m;i<=n;i+=istep) {

j=i+mmax;

tempr=wr*data[j]-wi*data[j+1];

tempi=wr*data[j+1]+wi*data[j];

data[j]=data[i]-tempr;

data[j+1]=data[i+1]-tempi;

data[i] += tempr;

data[i+1] += tempi;

}

wr=(wtemp=wr)*wpr-wi*wpi+wr;

wi=wi*wpr+wtemp*wpi+wi;

}

mmax=istep;

}

}

#undef SWAP

/* (C) Copr. 1986-92 Numerical Recipes Software */

17. Generate ERMapper Header: gen_ers.c

void gen_ers (fname, nrow, ncol, drow, dcol, nband, bflag)

char fname[], bflag[];

int  ncol, nrow, nband;

double drow, dcol;

{

char fnamers[80];

FILE *fp;

int i, nbyte;

strcpy(fnamers, fname);

strcat(fnamers, ".ers");

if((fp = fopen(fnamers, "w"))==NULL)

{

printf("Error... Can not open %s.\n", fnamers);

exit(0);

}

fprintf(fp,"DatasetHeader Begin\n");

fprintf(fp,"        Version         = \"6.0\"\n");

fprintf(fp,"        Name	= \"%s\"\n", fnamers);

fprintf(fp,"        Description	= \"Generated by Hoonyol Lee\"\n");

fprintf(fp,"        LastUpdated     = Thu Jun 12 00:00:00 GMT 2004\n");

fprintf(fp,"        SensorName      = \"Unknown\"\n");

fprintf(fp,"        SenseDate       = Thu Jun 12 00:00:00 GMT 2004\n");

fprintf(fp,"        DataSetType     = ERStorage\n");

fprintf(fp,"        DataType        = Raster\n");

fprintf(fp,"        ByteOrder       = LSBFirst\n");

fprintf(fp,"        CoordinateSpace Begin\n");

fprintf(fp,"                Datum           = \"RAW\"\n");

fprintf(fp,"                Projection      = \"RAW\"\n");

fprintf(fp,"                Projection      = \"RAW\"\n");

fprintf(fp,"                CoordinateType  = RAW\n");

fprintf(fp,"                Rotation        = 0:0:0.0\n");

fprintf(fp,"        CoordinateSpace End\n");

fprintf(fp,"        RasterInfo Begin\n");

if(strncmp(bflag, "010", 3) == 0)

fprintf(fp,"                CellType        = Unsigned8BitInteger\n");

else if(strncmp(bflag, "110",3) == 0)

fprintf(fp,"                CellType        = Signed8BitInteger\n");

else if(strncmp(bflag, "020",3) == 0)

fprintf(fp,"                CellType        = Unsigned16BitInteger\n");

else if(strncmp(bflag, "120",3) == 0)

fprintf(fp,"                CellType        = Signed16BitInteger \n");

else if(strncmp(bflag, "040",3) == 0)

fprintf(fp,"                CellType        = Unsigned32BitInteger\n");

else if(strncmp(bflag, "140",3) == 0)

fprintf(fp,"                CellType        = Signed32BitInteger \n");

else if(strncmp(bflag, "141",3) == 0)

fprintf(fp,"                CellType        = IEEE4ByteReal\n");

else if(strncmp(bflag, "181",3) == 0)

fprintf(fp,"                CellType        = IEEE8ByteReal\n");

else

{	nbyte = atoi(bflag);

fprintf(fp,"                CellType        = Unsigned%dBitInteger\n", nbyte);

}

fprintf(fp,"                CellInfo Begin\n");

fprintf(fp,"                        Xdimension      = %lf\n", dcol);

fprintf(fp,"                        Ydimension      = %lf\n", drow);

fprintf(fp,"                CellInfo End\n");

fprintf(fp,"                NrOfLines       = %d\n", nrow);

fprintf(fp,"                NrOfCellsPerLine        = %d\n", ncol);

fprintf(fp,"                NrOfBands       = %d\n", nband);

for(i=0; i<nband; i++){

fprintf(fp,"                BandId Begin\n");

fprintf(fp,"                        Value           = \"Band%d\"\n", i+1);

fprintf(fp,"                BandId End\n");

}

fprintf(fp,"        RasterInfo End\n");

fprintf(fp,"DatasetHeader End\n");

fclose(fp);

}

void ers_usage()

{

printf(" Usage: gen_ers datafile nrow ncol nband bit_flag(or nbyte)\n ");

printf("\n");

printf("        where datafile : input file name\n");

printf("                  nrow : number of row\n");

printf("                  ncol : number of column\n");

printf("                 nband : number of band\n\n");

printf("                 nbyte : number of byte per cell\n\n");

printf("              bit_flag : bit flag\n\n");

printf("                         number of byte, or \n");

printf("                         010   Unsigned  8 bit Integer\n");

printf("                         110     Signed  8 bit Integer\n");

printf("                         020   Unsigned 16 bit Integer\n");

printf("                         120     Signed 16 bit Integer\n");

printf("                         040   Unsigned 32 bit Integer\n");

printf("                         140     Signed 32 bit Integer\n");

printf("                         141       IEEE  4 Byte Real\n");

printf("                         181       IEEE  8 Byte Real\n");

printf("\n");

}

18. Generate ERMapper Header for CMOD: gen_ers_cmod.c

void gen_ers_cmod (fname, nrow, ncol, drow, dcol, nband, bflag)

char fname[], bflag[];

int  ncol, nrow, nband;

double drow, dcol;

{

char fnamers[80];

FILE *fp;

int i, nbyte;

strcpy(fnamers, fname);

strcat(fnamers, ".ers");

if((fp = fopen(fnamers, "w"))==NULL)

{

printf("Error... Can not open %s.\n", fnamers);

exit(0);

}

fprintf(fp,"DatasetHeader Begin\n");

fprintf(fp,"        Version         = \"6.0\"\n");

fprintf(fp,"        Name	= \"%s\"\n", fnamers);

fprintf(fp,"        Description	= \"Generated by Hoonyol Lee\"\n");

fprintf(fp,"        LastUpdated     = Thu Jun 12 00:00:00 GMT 2004\n");

fprintf(fp,"        SensorName      = \"Unknown\"\n");

fprintf(fp,"        SenseDate       = Thu Jun 12 00:00:00 GMT 2004\n");

fprintf(fp,"        DataSetType     = ERStorage\n");

fprintf(fp,"        DataType        = Raster\n");

fprintf(fp,"        ByteOrder       = LSBFirst\n");

fprintf(fp,"        CoordinateSpace Begin\n");

fprintf(fp,"                Datum           = \"RAW\"\n");

fprintf(fp,"                Projection      = \"RAW\"\n");

fprintf(fp,"                Projection      = \"RAW\"\n");

fprintf(fp,"                CoordinateType  = RAW\n");

fprintf(fp,"                Rotation        = 0:0:0.0\n");

fprintf(fp,"        CoordinateSpace End\n");

fprintf(fp,"        RasterInfo Begin\n");

if(strncmp(bflag, "010", 3) == 0)

fprintf(fp,"                CellType        = Unsigned8BitInteger\n");

else if(strncmp(bflag, "110",3) == 0)

fprintf(fp,"                CellType        = Signed8BitInteger\n");

else if(strncmp(bflag, "020",3) == 0)

fprintf(fp,"                CellType        = Unsigned16BitInteger\n");

else if(strncmp(bflag, "120",3) == 0)

fprintf(fp,"                CellType        = Signed16BitInteger \n");

else if(strncmp(bflag, "040",3) == 0)

fprintf(fp,"                CellType        = Unsigned32BitInteger\n");

else if(strncmp(bflag, "140",3) == 0)

fprintf(fp,"                CellType        = Signed32BitInteger \n");

else if(strncmp(bflag, "141",3) == 0)

fprintf(fp,"                CellType        = IEEE4ByteReal\n");

else if(strncmp(bflag, "181",3) == 0)

fprintf(fp,"                CellType        = IEEE8ByteReal\n");

else

{	nbyte = atoi(bflag);

fprintf(fp,"                CellType        = Unsigned%dBitInteger\n", nbyte);

}

fprintf(fp,"                CellInfo Begin\n");

fprintf(fp,"                        Xdimension      = %lf\n", dcol);

fprintf(fp,"                        Ydimension      = %lf\n", drow);

fprintf(fp,"                CellInfo End\n");

fprintf(fp,"                NrOfLines       = %d\n", nrow);

fprintf(fp,"                NrOfCellsPerLine        = %d\n", ncol);

fprintf(fp,"                NrOfBands       = %d\n", nband);

fprintf(fp,"                BandId Begin\n");

fprintf(fp,"                        Value           = \"CMOD4\"\n");

fprintf(fp,"                BandId End\n");

fprintf(fp,"                BandId Begin\n");

fprintf(fp,"                        Value           = \"CMOD_IFR2\"\n");

fprintf(fp,"                BandId End\n");

fprintf(fp,"        RasterInfo End\n");

fprintf(fp,"DatasetHeader End\n");

fclose(fp);

}

19. Generate ERMapper Header for wave_crs: gen_ers_crs.c

void gen_ers_crs(fname, nrow, ncol, drow, dcol, nband, bflag)

char fname[], bflag[];

int  ncol, nrow, nband;

double drow, dcol;

{

char fnamers[80];

FILE *fp;

int i, nbyte;

strcpy(fnamers, fname);

strcat(fnamers, ".ers");

if((fp = fopen(fnamers, "w"))==NULL)

{

printf("Error... Can not open %s.\n", fnamers);

exit(0);

}

fprintf(fp,"DatasetHeader Begin\n");

fprintf(fp,"        Version         = \"6.0\"\n");

fprintf(fp,"        Name	= \"%s\"\n", fnamers);

fprintf(fp,"        Description	= \"Generated by Hoonyol Lee\"\n");

fprintf(fp,"        LastUpdated     = Thu Jun 12 00:00:00 GMT 2004\n");

fprintf(fp,"        SensorName      = \"Unknown\"\n");

fprintf(fp,"        SenseDate       = Thu Jun 12 00:00:00 GMT 2004\n");

fprintf(fp,"        DataSetType     = ERStorage\n");

fprintf(fp,"        DataType        = Raster\n");

fprintf(fp,"        ByteOrder       = LSBFirst\n");

fprintf(fp,"        CoordinateSpace Begin\n");

fprintf(fp,"                Datum           = \"RAW\"\n");

fprintf(fp,"                Projection      = \"RAW\"\n");

fprintf(fp,"                Projection      = \"RAW\"\n");

fprintf(fp,"                CoordinateType  = RAW\n");

fprintf(fp,"                Rotation        = 0:0:0.0\n");

fprintf(fp,"        CoordinateSpace End\n");

fprintf(fp,"        RasterInfo Begin\n");

if(strncmp(bflag, "010", 3) == 0)

fprintf(fp,"                CellType        = Unsigned8BitInteger\n");

else if(strncmp(bflag, "110",3) == 0)

fprintf(fp,"                CellType        = Signed8BitInteger\n");

else if(strncmp(bflag, "020",3) == 0)

fprintf(fp,"                CellType        = Unsigned16BitInteger\n");

else if(strncmp(bflag, "120",3) == 0)

fprintf(fp,"                CellType        = Signed16BitInteger \n");

else if(strncmp(bflag, "040",3) == 0)

fprintf(fp,"                CellType        = Unsigned32BitInteger\n");

else if(strncmp(bflag, "140",3) == 0)

fprintf(fp,"                CellType        = Signed32BitInteger \n");

else if(strncmp(bflag, "141",3) == 0)

fprintf(fp,"                CellType        = IEEE4ByteReal\n");

else if(strncmp(bflag, "181",3) == 0)

fprintf(fp,"                CellType        = IEEE8ByteReal\n");

else

{	nbyte = atoi(bflag);

fprintf(fp,"                CellType        = Unsigned%dBitInteger\n", nbyte);

}

fprintf(fp,"                CellInfo Begin\n");

fprintf(fp,"                        Xdimension      = %lf\n", dcol);

fprintf(fp,"                        Ydimension      = %lf\n", drow);

fprintf(fp,"                CellInfo End\n");

fprintf(fp,"                NrOfLines       = %d\n", nrow);

fprintf(fp,"                NrOfCellsPerLine        = %d\n", ncol);

fprintf(fp,"                NrOfBands       = %d\n", nband);

if(c_wave_crsflag==1)

{

fprintf(fp,"                BandId Begin\n");

fprintf(fp,"                        Value           = \"Amplitude\"\n");

fprintf(fp,"                BandId End\n");

fprintf(fp,"                BandId Begin\n");

fprintf(fp,"                        Value           = \"Phase\"\n");

fprintf(fp,"                BandId End\n");

}

else  /* c_crsflag==2 */

{

fprintf(fp,"                BandId Begin\n");

fprintf(fp,"                        Value           = \"Real\"\n");

fprintf(fp,"                BandId End\n");

fprintf(fp,"                BandId Begin\n");

fprintf(fp,"                        Value           = \"Imaginary\"\n");

fprintf(fp,"                BandId End\n");

}

fprintf(fp,"        RasterInfo End\n");

fprintf(fp,"DatasetHeader End\n");

fclose(fp);

}

20. Generate ERMapper Vector Header: generv.c

void generv(filename, tlx, tly, brx, bry)

char filename[];

double tlx, tly, brx, bry;

{

FILE *fp1;

char anow[64];

strcpy(anow,"Thu Jul 08 16:56:11 GMT 2004");

if((fp1 = fopen(filename, "w")) == NULL){

printf("ERROR... can not creat %s.\n", filename);

exit(0);

}

fprintf(fp1, "DatasetHeader Begin\n");

fprintf(fp1, "	Version         = \"6.0\"\n");

fprintf(fp1, "	Name            = \"%s\"\n", filename);

fprintf(fp1, "	LastUpdated     = %s\n", anow);

fprintf(fp1, "	DataSetType     = None\n");

fprintf(fp1, "	DataType        = Vector\n");

fprintf(fp1, "	ByteOrder       = LSBFirst\n");

fprintf(fp1, "	CoordinateSpace Begin\n");

fprintf(fp1, "        Datum           = \"RAW\"\n");

fprintf(fp1, "        Projection      = \"RAW\"\n");

fprintf(fp1, "        CoordinateType  = RAW\n");

fprintf(fp1, "        Rotation        = 0:0:0.0\n");

fprintf(fp1, "	CoordinateSpace End\n");

fprintf(fp1, "	VectorInfo Begin\n");

fprintf(fp1, "        Type            = ERVEC\n");

fprintf(fp1, "        FileFormat      = ASCII\n");

fprintf(fp1, "        Extents Begin\n");

fprintf(fp1, "                TopLeftCorner Begin\n");

fprintf(fp1, "                        MetersX         = %lf\n", tlx);

fprintf(fp1, "                        MetersY         = %lf\n", tly);

fprintf(fp1, "                TopLeftCorner End\n");

fprintf(fp1, "                BottomRightCorner Begin\n");

fprintf(fp1, "                        MetersX         = %lf\n", brx);

fprintf(fp1, "                        MetersY         = %lf\n", bry);

fprintf(fp1, "                BottomRightCorner End\n");

fprintf(fp1, "        Extents End\n");

fprintf(fp1, "	VectorInfo End\n");

fprintf(fp1, "DatasetHeader End\n");

fclose(fp1);

}

21. Memory Allocation: ma.c

/* memory allocation for vectors and matrices */

/* modified from nrutil.c of Numerical Recipes in C */

/* index notation follows C convention */

/*#include <stdlib.h>*/

unsigned char *ucvector(int n)

{

unsigned char *v;

v = (unsigned char *)calloc((size_t) n, (size_t) sizeof( unsigned char));

if(!v)

{       printf("malloc failure in unsigned char vector\n");

exit(0);

}

return v;

}

char *cvector(int n)

{

char *v;

v = (char *)calloc((size_t) n, (size_t) sizeof(char));

if(!v)

{       printf("malloc failure in char vector\n");

exit(0);

}

return v;

}

unsigned short *usvector(int n)

{

unsigned short *v;

v = (unsigned short *)calloc((size_t) n, (size_t) sizeof(unsigned short));

if(!v)

{       printf("malloc failure in unsigned short vector\n");

exit(0);

}

return v;

}

short *svector(int n)

{

short *v;

v = (short *)calloc((size_t) n, (size_t) sizeof(short));

if(!v)

{       printf("malloc failure in short vector\n");

exit(0);

}

return v;

}

unsigned int *uivector(int n)

{

unsigned int *v;

v = (unsigned int *)calloc((size_t) n, (size_t) sizeof(unsigned int));

if(!v)

{       printf("malloc failure in unsigned int vector\n");

exit(0);

}

return v;

}

int *ivector(int n)

{

int *v;

v = (int *)calloc((size_t) n, (size_t) sizeof(int));

if(!v)

{       printf("malloc failure in int vector\n");

exit(0);

}

return v;

}

float *fvector(int n)

{

float *v;

v = (float *)calloc((size_t) n, (size_t) sizeof(float));

if(!v)

{	printf("malloc failure in float vector\n");

exit(0);

}

return v;

}

double *dvector(int n)

{

double *v;

v = (double *)calloc((size_t) n, (size_t) sizeof(double));

if(!v)

{       printf("malloc failure in double vector\n");

exit(0);

}

return v;

}

int **imatrix(int nrow, int ncol)

{

int i;

int **m;

/* allocate pointers to rows */

m=(int **)calloc((size_t) nrow, (size_t) sizeof(int*));

if(!m)

{       printf("malloc failure 1 in int matrix\n");

exit(0);

}

/* allocate rows and set pointers to them */

m[0]=(int *)calloc((size_t) nrow*ncol, (size_t)sizeof(int));

if(!m[0])

{       printf("malloc failure 2 in int matrix\n");

exit(0);

}

for(i=1; i< nrow; i++) m[i] = m[i-1] + ncol;

return m;

}

float **fmatrix(int nrow, int ncol)

{

int i;

float **m;

/* allocate pointers to rows */

m=(float **)calloc((size_t) nrow, (size_t) sizeof(float*));

if(!m)

{       printf("malloc failure 1 in float matrix\n");

exit(0);

}

/* allocate rows and set pointers to them */

m[0]=(float *)calloc((size_t) nrow*ncol, (size_t) sizeof(float));

if(!m[0])

{       printf("malloc failure 2 in float matrix\n");

exit(0);

}

for(i=1; i< nrow; i++) m[i] = m[i-1] + ncol;

return m;

}

double **dmatrix(int nrow, int ncol)

{

int i;

double **m;

/* allocate pointers to rows */

m=(double **)calloc((size_t) nrow, (size_t) sizeof(double*));

if(!m)

{       printf("malloc failure 1 in double matrix\n");

exit(0);

}

/* allocate rows and set pointers to them */

m[0]=(double *)calloc((size_t) nrow*ncol, (size_t) sizeof(double));

if(!m[0])

{       printf("malloc failure 2 in double matrix\n");

exit(0);

}

for(i=1; i< nrow; i++) m[i] = m[i-1] + ncol;

return m;

}

22. Data Average Filter: maverage.c

/*********************************************

maverage.c

multi-band average (replace)

maverage : double

maverage_f : float

input: double BIP

output: double BIP averaged (replace)

note nx, ny should be odd numbers

programmed by Dr. Hoonyol Lee

5 March 2002

modifiend to double

19 April 2002

memory allocation

extreme extra memory saving:

[nrow*ncol*nband] -> [nrow_buf*ncol*nband]

21 May 2004

**********************************************/

/*

#include <math.h>

#include "ma.c"

*/

void maverage(data, nrow, ncol, nband, nx, ny)

double *data;

int ncol, nrow, nband, nx, ny;

{

double *buf, *tmp;

int i, j, k, l, m, jx, kx, lx, ncount;

int nrow_buf, ibuf, ibufx, ri, rix;

/*printf(" maverage: %dx%d\n", nx, ny);*/

if(nx==(nx/2*2)) nx++;

if(ny==(ny/2*2)) ny++;

nrow_buf= ny/2+1; /* number of buffer row */

buf = dvector(nrow_buf*ncol*nband);

tmp = dvector(nband);

for(i=0; i< (nrow+nrow_buf); i++)

{

/* buffer row */

ibuf = (int)fmod(i, nrow_buf);

ibufx = ibuf*ncol*nband;

/* replace data with this buffer row */

ri = i - nrow_buf;

if(ri>=0)

{

rix=ri*ncol*nband;

for(j=0; j<ncol*nband; j++)

data[rix+j]=buf[ibufx+j];

}

/* make new averaged row */

if(i<nrow)

{

for(j=0; j<ncol; j++)

{

for(m=0; m<nband; m++) tmp[m]=0.;

ncount = 0;

for(k=i-ny/2; k<=i+ny/2; k++)

{

kx=k*ncol*nband;

for(l=j-nx/2; l<=j+nx/2; l++)

if(k>=0 && k<nrow && l>=0 && l<ncol)

{

lx=l*nband;

for(m=0; m<nband; m++)

tmp[m] += data[kx+lx+m];

ncount ++;

}

}

jx=j*nband;

for(m=0; m<nband; m++)

buf[ibufx+jx+m] = tmp[m]/ncount;

}

}

}

free(buf);

}

void maverage_f(data, nrow, ncol, nband, nx, ny)

float *data;

int ncol, nrow, nband, nx, ny;

{

float *buf, *tmp;

int i, j, k, l, m, jx, kx, lx, ncount;

int nrow_buf, ibuf, ibufx, ri, rix;

if(nx==(nx/2*2)) nx++;

if(ny==(ny/2*2)) ny++;

nrow_buf= ny/2+1; /* number of buffer row */

buf = fvector(nrow_buf*ncol*nband);

tmp = fvector(nband);

for(i=0; i< (nrow+nrow_buf); i++)

{

/* buffer row */

ibuf = (int)fmod(i, nrow_buf);

ibufx = ibuf*ncol*nband;

/* replace data with this buffer row */

ri = i - nrow_buf;

if(ri>=0)

{

rix=ri*ncol*nband;

for(j=0; j<ncol*nband; j++)

data[rix+j]=buf[ibufx+j];

}

/* make new averaged row */

if(i<nrow)

{

for(j=0; j<ncol; j++)

{

for(m=0; m<nband; m++) tmp[m]=0.;

ncount = 0;

for(k=i-ny/2; k<=i+ny/2; k++)

{

kx=k*ncol*nband;

for(l=j-nx/2; l<=j+nx/2; l++)

if(k>=0 && k<nrow && l>=0 && l<ncol)

{

lx=l*nband;

for(m=0; m<nband; m++)

tmp[m] += data[kx+lx+m];

ncount ++;

}

}

jx=j*nband;

for(m=0; m<nband; m++)

buf[ibufx+jx+m] = tmp[m]/ncount;

}

}

}

free(buf);

free(tmp);

}

23. Input Parameter File: sop.pam

fname_in	x:/RS/Radarsat1/jeju/R1_19991219_21528D_SLC/SCENE01/DAT_01.001

fname_leader	x:/RS/Radarsat1/jeju/R1_19991219_21528D_SLC/SCENE01/LEA_01.001

fhead	jeju

c_load_slc	1

c_multilook_resmp	1

nx,ny	1  5

nx_res,ny_res	1  4

c_wind	1

wind_direction	280.

vv2hh_alpha	0.6

c_wave	1

nptxw,nptyw	256 256

c_wave_crsflag	1

c_current	1

nptyc	256

dc_xstep,dc_ystep	128 128

dc_avx,dc_avy	3  11

vel_scale	100.

<< Parameter Description >>

<File Control>

fname_in: location of SAR SLC data file in CEOS format

fname_leader: location of SAR SLC leader file

fhead: file header name for new file products

<Process Controls>

c_load_slc : Load SLC Image

0	Skip

1	Run

c_multilook_resmp : Produce Multilook Image by average and Resample

0	Skip

1	Run

nx, ny: averaging factor for multilook by slc amplitude averaging

nx_res, ny_res: resampling factor for multilook

c_wind: SAR Wind by Sigma0, CMOD4 and CMOD_IFR2

0	Skip

1	Run

wind_direction:

(degree) [0,360)

vv2hh_alpha : VV to HH conversion model selection

0.	Bragg Scattering

0.6	Unal model

1.	Kirchhoff Scattering

2	Elfouhaily Model

c_wave: SAR Wave and Doppler Shift for Current

0	Skip

1	Run

nptxw, nptyw: 2D-FFT size in range and azimuth

c_wave_crsflag : Wave direction (solve 180deg ambiguity) by Interlook Cross Spectrum

0	Skip

1	AP Cross Spectrum

2	IQ  Cross Spectrum

c_current: Current estimation via Doppler shift

0	Skip

1	Run

nptyc: 1D azimuth FFT size

dc_xstep,dc_ystep: current estimation steps

dc_avx,dc_avy: averaging factor during doppler center estimation

vel_scale: velocity vector scale (velocity = vector_length/vel_scale m/s )

24. Convert String to Integer: stoi.c

/********************************************

stoi.c

:converts nbyte string to 4-byte int

1. The bytes larger than sizeof(int) will be lost.

2. Byte order dependent

3. Note the difference from atoi() in <stdlib.h>

which convert from asc to int.

by Dr. Hoonyol Lee

17 July 2006

*********************************************/

int stoi(unsigned char *data, int nbyte)

{

int i, idata, *itmp;

itmp=ivector(nbyte);

idata=0;

for(i=0; i<nbyte; i++)

{

itmp=0;

itmp |= *(data++);

idata |= itmp;

}

return(idata);

}

25. Swap Byte of an Array: swapbytes.c

/***********************************

swapbytes.c

swap byte order of an array

by Dr. Hoonyol Lee

Dept. of Geophysics

Kangwon National University

29 June 2004

***********************************/

void byteswap(void *buf, unsigned nbyte)

{

unsigned i;

char tmp;

for(i=0; i<nbyte/2; i++)

{

tmp=*((char *)buf+i);

*((char *)buf+i)=*((char *)buf+nbyte-i-1);

*((char *)buf+nbyte-i-1)=tmp;

}

}

void swapbytes(buf, narray, nbyte)

void *buf;

unsigned narray, nbyte;

{

int i;

for(i=0; i<narray; i++)

byteswap(buf+i*nbyte, nbyte);

}

26. Get N Power of Two: town.c

/* program to get nn of 2*nn */

/* H.Lee 17/7/01 */

unsigned long twon(int numb)

{

unsigned long nn = 1;

int n;

n = numb - 1;

do{

n >>= 1;

nn <<= 1;

}while(n != 0);

return(nn);

}

27. Program Run Log: log

***************************************************

*	SAR Ocean Processor (SOP)          *

*	v1.1	*

*     All right reserved by Dr. Hoonyol Lee,        *

* Dept. of Geophysics, Kangwon National University. *

* Email:hoonyol@kangwon.ac.kr, Tel:+82-33-250-8587  *

***************************************************

Load Parameters and Data

load_para:

<Input Parameters>

fname_in=x:/RS/Radarsat1/jeju/R1_19991219_21528D_SLC/SCENE01/DAT_01.001

fname_leader=x:/RS/Radarsat1/jeju/R1_19991219_21528D_SLC/SCENE01/LEA_01.001

fhead=jejuD

<Control Paramters>

c_load_slc= 1

c_multilook_resmp= 1

multilook nx,ny= 1, 5

multilook nx_res,ny_res= 1, 4

c_wind= 1

wind_direction= 280

vv2hh_alpha= 0.6

Thompson Model

c_wave= 1

nptxw,nptyw= 256, 256

c_wave_crsflag= 1

AP Cross Spectrum

c_current= 1

nptyc= 256

dc_xstep,dc_ystep= 128, 128

dc_avx,dc_avy= 3, 11

vel_scale= 100.000000

<Data File Parameters>

nskip_byte_per_record=192

nbyte_header_record=16252

nazimuth_whole=18077

nbyte_per_record=22364

nbyte_data=2

nrange_whole=5543

<Leader File Parameters>

Orbit orientation(asc_des)=DESCENDING      (-1)

earth_radius_major=6378.137000

earth_radius_minor=6356.752314

plat_lat=33.091000(deg)

plat_long=130.187000(deg)

plat_head=193.519000(deg)

nominal_incidence_angle=23.552000

wavelength=0.056565

prf=1306.373779

azimuth_bw=900.000000

Dc0=5.544929e+03

Dc1=2.460548e-01

Dc2=-4.636476e-06

azimuth_pixel_size=5.112507

range_pixel_size=8.117024

eph_orb_data=7.16709e+06

Radiometric samp_inc =11

<Generic Parameters>

cspeed= 299792500.000000 m/s

pi= 3.141593

<Other Parameters>

nskip_record=1

ground_range_pixel_size=20.313827

Output drange=20.313827, dazimuth=5.112507

Multilook drange_ml=20.313827, dazimuth_ml=20.450028

Multilook nrange_ml=5543, nazimuth_ml=4519

Earth radius = 6.37174e+06(m)

Orbit Altitude = 795347(m)

Incidence Angles

first range 0,  35.3185(deg)

mid range 2771, 28.7752(deg)

far range 5543, 19.1547(deg)

fname_slc=01_jejuD_SLC

fname_ml=02_jejuD_ML

fname_sigma0=03_jejuD_Wind_Simga0

fname_cmod=03_jejuD_Wind_CMOD

fname_wave=04_jejuD_Wave

fname_crs=04_jejuD_Wave_ILC

fname_crs_ml=04_jejuD_Wave_ML

fname_cur=05_jejuD_Current

fname_cur_erv=05_jejuD_Current_Doppler

fname_cur_vec=05_jejuD_Current_Vector

fname_cur_txt=05_jejuD_Current_Vector.txt

load_slc:

Extract SAR SLC data.

SAR SLC data output 01_jejuD_SLC.

Multilook by Average and Resampling

1/10

2/10

3/10

4/10

5/10

6/10

7/10

8/10

9/10

10/10

Sigma0, CMOD4 and CMOD-IFR2 for Wind Extraction

Sigma0:

min=-31.1184, max=12.4245

CMOD4 LUT:

LUT[range, sigma0]=[185, 450]

Incidence angle min, max(deg) = 19.1547, 35.3185

far range first

LUT sigma0 min, max(dB)= -26.8595, 18.0841

LUT wind min, max(m/s)= 1, 67

CMOD_IFR2 LUT:

LUT[range, sigma0]=[185, 532]

Incidence angle min, max(deg) = 19.1547, 35.3185

far range first

LUT sigma0 min, max(dB)= -26.8595, 18.0841

LUT wind min, max(m/s)= 0, 41

SAR Wave Spectrum

1/18 block

2/18 block

3/18 block

4/18 block

5/18 block

6/18 block

7/18 block

8/18 block

9/18 block

10/18 block

11/18 block

12/18 block

13/18 block

14/18 block

15/18 block

16/18 block

17/18 block

final null block

Inter-Look Cross Spectrum Process

nblock=1024

1/18 block

first block

read slc ok

Az FFT ok

split aperture ok

Az iFFT ok

detect multilook ok

write multilook ok

cross spec ok

write crs ok

2/18 block

read slc ok

Az FFT ok

split aperture ok

Az iFFT ok

detect multilook ok

write multilook ok

cross spec ok

write crs ok

... 중략 ...

18/18 block

last block

read slc ok

Az FFT ok

split aperture ok

Az iFFT ok

detect multilook ok

write multilook ok

cross spec ok

write crs ok

Current from Doppler Shift

1/36 blocks

read slc ok

Az FFT ok

detect ok

mlook ok

resample ok

avg ok

doppler ok

2/36 blocks

read slc ok

Az FFT ok

detect ok

mlook ok

resample ok

avg ok

doppler ok

... 중략 ...

34/36 blocks

read slc ok

Az FFT ok

detect ok

mlook ok

resample ok

avg ok

doppler ok

Program finished successfully (1002 sec).