Himawari画像から表層流速を推定する：MCC/PIVとOptical Flow

最初の画像 field1 を作ります。乱数ノイズだけではPIVが難しいため、複数のガウス状パッチ、背景勾配、波状の濃淡を重ねます。ここで作るのは海洋そのものではなく、SSTやCHLAに似た「追跡できる模様」を持つ教材用画像です。

• xg, yg：格子座標。
• field1：時刻1の疑似トレーサー画像。
• robust_normalize：外れ値の影響を抑え、相関計算に使いやすくする関数。

真値の変位場 uTrue, vTrue を作ります。ここでは単純な一様移動ではなく、渦のような回転成分と弱い背景流を含めます。学生が「真値」と「推定値」を比べられるようにするための重要ブロックです。

• uTrue：x方向の真の変位。pixel単位です。
• vTrue：y方向の真の変位。pixel単位です。
• この段階では物理単位 m/s ではなく、画像格子上の移動量として扱います。

field1 を真値変位場で移動させ、時刻2の画像 field2 を作ります。ここでは「時刻2の各格子点に来た水塊は、時刻1ではどこにいたか」を補間で求めるため、interp2 を使います。

入力画像と真値ベクトルを図にします。ここで、模様がどの方向へ流されているかを目で確認します。PIVの前に、画像に追跡可能な構造があるかを確認するためのブロックです。

PIVのパラメータを設定します。ここでの設定が結果の見た目と信頼性を大きく変えます。

このブロックがPIVの本体です。robust_normalize で画像を標準化し、mcc_piv_core で各窓の最大相互相関を探します。その後、basic QCとneighbor QCを順にかけます。

I1 = robust_normalize(field1, 3.0);
I2 = robust_normalize(field2, 3.0);

P = mcc_piv_core(I1, I2, win, step, maxShift, ...
                 minValidFraction, peakExcludeRadius);

basic_qc = raw_mask & ...
           P.valid_fraction >= minValidFraction & ...
           P.rmax >= minCorrelation & ...
           P.peak_gap >= minPeakGap & ...
           P.texture >= minTextureStd & ...
           ~P.edge_flag;

neighbor_qc = neighbor_consistency_qc(P.dx, P.dy, basic_qc, maxNeighborDiff);
final_qc = basic_qc & neighbor_qc;

basic_qc は各ベクトル単体の品質、neighbor_qc は周囲との整合性です。どちらか一方だけでは不十分です。

真値ベクトルと推定ベクトルを並べます。ここでは「推定できている場所」と「推定が怪しい場所」を目で確認します。真値があるrun06aでしかできない重要な診断です。

dx, dy, rmax, peak_gap などを地図状に表示します。ベクトル図だけでは、なぜ落ちたのか・なぜ残ったのかが分かりにくいため、QC診断図を必ず確認します。

結果をMATファイルに保存します。後から比較・図の再作成・誤差計算をするため、推定ベクトルだけでなく、真値、QCマスク、相関診断量も保存します。

スクリプト末尾には、計算を支えるローカル関数があります。授業では中身を全部書く必要はありませんが、役割は理解します。

run06aと同じです。雲を入れる前の「本当なら追跡できる画像」を作ります。ここをrun06aと同じ構造にすることで、雲の影響だけを比較できます。

人工的な雲マスク cloudMask1, cloudMask2 を作り、雲の場所をNaNにします。雲は「値が見えない領域」です。PIVではNaNを相関計算から除外しますが、それだけでは雲縁や中心点の問題が残ります。

• cloudMask1：画像1の雲。
• cloudMask2：画像2の雲。
• field1_cloudy, field2_cloudy：雲部分をNaNにした画像。

雲なし画像と雲あり画像を比較します。雲によって相関窓がどの程度欠けるかを視覚的に確認します。

run06aの設定に加えて、雲QC用の cloudBufferPix を設定します。

雲あり画像をrobust normalizeしてMCC/PIVを実行します。相関計算では isfinite(A) & isfinite(B) により、NaN画素を除外します。しかし、窓の一部が雲で欠けていても相関は出るため、追加でcloud QCを行います。

cloudMask1_buffer = buffer_logical_mask(cloudMask1, cloudBufferPix);
cloudMask2_buffer = buffer_logical_mask(cloudMask2, cloudBufferPix);

origin_clear = sample_logical_mask(~cloudMask1_buffer, P.ix, P.iy);
target_clear = sample_logical_mask(~cloudMask2_buffer, P.ix + P.dx, P.iy + P.dy);

cloud_qc = raw_mask & origin_clear & target_clear;

origin_clear はベクトル始点、target_clear は推定された移動先を調べます。これにより、中心点が雲上なのに窓の周辺だけで相関が出てしまう問題を防ぎます。

QC前のベクトルをあえて表示します。雲や雲縁の上にもベクトルが出てしまうことを確認するためです。この図は「なぜQCが必要か」を理解するために重要です。

valid_fraction, rmax, peak_gap, cloud_qc, basic_qc, neighbor_qc, final_qc を段階的に表示します。どのQCでどのベクトルが落ちたかを確認します。

最終QC後のベクトルを表示します。雲ありでは推定可能領域が減りますが、雲のない場所ではPIVが成立することを確認します。雲縁に少し残るベクトルは、実データでも起こり得る注意点です。

最終QC後の dx, dy, 誤差、相関を表示します。雲あり条件で、残ったベクトルがどの程度真値に合っているかを確認します。

雲マスク、cloud QC、最終QC、真値、推定値を保存します。保存しておくことで、雲バッファや閾値を変えたときの比較ができます。

入力ファイル、出力フォルダ、PIVモード、ROIを設定します。ここを変えるだけで、広域・fine・ROI集中解析を切り替えられます。

SSTのMATファイル2つを読み込みます。想定される変数は lon, lat, SST です。緯度が降順になっている場合は、図化・計算しやすいように昇順へ直します。

• ensure_lat_ascending：緯度配列とデータを南から北へ並べ直します。
• crop_domain：ROIモードでは指定範囲だけ切り出します。
• SSTの非現実的な値や欠測はNaNとして扱います。

格子間隔、時間差、PIV窓、QC閾値を決めます。SSTはCHLAより滑らかなので、テクスチャ閾値は厳しすぎるとベクトルが残りません。

SST画像をMCC/PIV用画像 I1, I2 に変換します。ここでは画像全体から広域場を引きません。raw SSTを robust_normalize し、相関計算の各小窓内で平均を引きます。

% SST matching images: no broad high-pass subtraction
I1 = robust_normalize(SST1, 3.0);
I2 = robust_normalize(SST2, 3.0);

SSTは南北勾配や水塊境界が強いため、模様が滑らかな場所では移動先が曖昧になります。そのため rmax だけでなく peak_gap と neighbor_qc を必ず確認します。

mcc_piv_core で変位 dx, dy を求め、格子間隔と時間差から U, V, speed に変換します。ここで初めて pixel変位が m/s になります。

P = mcc_piv_core(I1, I2, lon, lat, WINDOW_SIZE, STEP, MAX_SHIFT, ...
                 MIN_VALID_FRACTION, PEAK_EXCLUDE_RADIUS);

[P.U, P.V, P.speed] = grid_displacement_to_velocity(P.dx, P.dy, ...
    P.lat, dlon, dlat, dt_seconds);

basic_qc = P.valid_fraction >= MIN_VALID_FRACTION & ...
           P.rmax >= MIN_CORRELATION & ...
           P.peak_gap >= MIN_PEAK_GAP & ...
           P.texture >= MIN_TEXTURE_STD & ...
           P.speed <= MAX_SPEED_MS & ...
           ~P.edge_flag & ...
           isfinite(P.dx) & isfinite(P.dy);

入力画像、QC前後、最終ベクトル、QC診断、ROIズームを作ります。PIVは数値だけでは信頼性を判断しにくいため、必ずQC診断図を一緒に見ます。

図、MATファイル、CSVを保存します。MATファイルには、lon_piv, lat_piv, U_piv, V_piv, speed_piv, mask_final などを保存します。run06eはこの出力を読みます。

入力CHLAファイル、出力フォルダ、PIVモード、ROIを設定します。roi_eddy では0.01°のCHLAを活かすこともできますが、通常のfine modeではSSTと比較しやすいように0.02°程度へ粗視化します。

CHLAのMATファイル2つを読み込みます。CHLAは海色アルゴリズム由来の値なので、SSTより欠測・ノイズ・沿岸影響が多くなります。緯度方向の並びをそろえ、必要ならROIを切り出します。

CHLA用のPIVパラメータを設定します。CHLAは細かいパッチやフィラメントが多く、SSTよりベクトルが出やすい一方、ノイズや大気補正の影響も拾いやすいです。

CHLAはまず正値だけを残し、log10(CHLA) に変換します。その後、必要に応じて粗視化し、robust normalizeします。SSTと同様、画像全体のhigh-passは行いません。

CHLA1_native(CHLA1_native <= 0) = NaN;
CHLA2_native(CHLA2_native <= 0) = NaN;

LCHLA1_native = log10(CHLA1_native);
LCHLA2_native = log10(CHLA2_native);

[lon, lat, LCHLA1] = block_nanmean_2d(lon0, lat0, LCHLA1_native, COARSEN_FACTOR);
[~,   ~,   LCHLA2] = block_nanmean_2d(lon0, lat0, LCHLA2_native, COARSEN_FACTOR);

I1 = robust_normalize(LCHLA1, 3.0);
I2 = robust_normalize(LCHLA2, 3.0);

run06cと同じMCC/PIV関数を使います。SSTとCHLAで同じ思想の相関計算を使うことで、run06eで比較しやすくします。違うのは、CHLAでは log10 と粗視化を先に行う点です。

入力画像、QC前後、最終ベクトル、QC診断、ROIズームを保存します。CHLAはベクトルが多く出るので、図としては密になります。これは必ずしも全部が正しいという意味ではなく、追跡できる濃淡パターンが多いという意味です。

MATとCSVを保存します。run06eは、SSTと同じ形式の piv 構造体を読み、SST-PIVとCHLA-PIVを比較します。

run06cとほぼ同じ関数群に加え、CHLA粗視化のための block_nanmean_2d が重要です。

run06c/dの出力MATファイル名、出力フォルダ、ROI、比較半径を設定します。run06eは穴埋めではありません。比較の考え方を読むための完全版です。

run06c_SST_PIV_vectors.mat と run06d_CHLA_PIV_vectors.mat を読みます。これらには、ベクトル位置、U/V成分、速度、QCマスクが入っています。

全域でSST-PIVを黒、CHLA-PIVを赤で重ねます。これは全体像を見る図です。CHLAの方がベクトル数が多いため、見た目だけで一致・不一致を判断するのは危険です。

148E, 38.5N付近など、注目領域だけを拡大して重ねます。広域図では見えない局所的な一致・不一致を確認します。

各SSTベクトルに対して、近傍のCHLAベクトルを探し、速度差と方向差を計算します。現行版では、SSTベクトル位置から一定距離内のCHLAベクトルを対応させる比較を行います。

% Concept:
% For each SST-PIV vector, find a nearby CHLA-PIV vector.
% Then compare:
%   speed difference    = CHLA speed - SST speed
%   direction difference = angle between SST and CHLA vectors

速度散布図、速度差ヒストグラム、方向差ヒストグラムを作ります。速度差が0付近にピークを持つなら、速度スケールは大きく破綻していません。方向差が0°側に山を持つなら、少なくとも一部では同じ移流場を反映していると考えられます。

SSTとCHLAの対応ベクトル、速度差、方向差をMATファイルやCSVに保存します。後で条件を変えて再解析するために、図だけでなく数値表も残します。