SpacetySliceData
/
Manual-Labeling-Tool
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity
Manual-Labeling-Tool/Manual-Labeling-Client/WindInverstLib/WindInverserToolGPU.cu

139 lines
4.7 KiB
Plaintext
Raw Permalink Normal View History

增加了风场cmod5反演
2025-11-20 09:32:46 +00:00
#include "WindInverserToolGPU.cuh"
#include "GPUTool.cuh"
// ʹ<><CAB9> __device__ <20>޶<EFBFBD><DEB6><EFBFBD><EFBFBD><EFBFBD>ʹ<EFBFBD><CAB9><EFBFBD><EFBFBD><EFBFBD><EFBFBD>GPU<50><55>ִ<EFBFBD><D6B4><EFBFBD><EFBFBD>ֻ<EFBFBD>ܱ<EFBFBD><DCB1><EFBFBD><EFBFBD><EFBFBD>GPU<50><55><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
__device__ float windspeedcmd5_device(float angle, float incidenceangle, float sigma0) {
// CMOD5ģ<35><C4A3>ϵ<EFBFBD><CFB5><EFBFBD><EFBFBD><EFBFBD><EFBFBD> - <20><><EFBFBD><EFBFBD><EFBFBD>ڴ<EFBFBD>ͨ<EFBFBD><CDA8><EFBFBD><EFBFBD><EFBFBD>ʺϴ<CABA><CFB4><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
const float cmod5_c[30] = {
1.0f, -0.688f, -0.793f, 0.338f, -0.173f, 0.0f, 0.004f,
0.111f, 0.0162f, 6.34f, 2.57f, -2.18f, 0.4f,
-0.6f, 0.045f, 0.007f, 0.33f, 0.012f, 22.0f,
1.95f, 3.0f, 8.39f, -3.44f, 1.36f, 5.35f, 1.99f, 0.29f, 3.80f, 1.53f
};
// ȷ<><C8B7><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Dz<EFBFBD>С<EFBFBD><D0A1>16<31><36>
if (incidenceangle + 0.5f <= 16.0f) {
incidenceangle = 16.0f;
}
// <20><>0-26m/s<><73>Χ<EFBFBD><CEA7><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
for (float U = 0.0f; U <= 26.0f; U += 0.01f) {
// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD>һ<EFBFBD><D2BB><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
float y0 = cmod5_c[19]; // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>C<EFBFBD><43><EFBFBD><EFBFBD><EFBFBD><EFBFBD>0<EFBFBD><30>ʼ
float n = cmod5_c[20];
float a = y0 - (y0 - 1.0f) / n;
float b = 1.0f / (n * powf(y0 - 1.0f, n - 1.0f)); // ʹ<><CAB9>powf<77><66><EFBFBD><EFBFBD>^
// <20><>һ<EFBFBD><D2BB><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
float x = (incidenceangle - 40.0f) / 25.0f;
// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɢ<EFBFBD><C9A2>ϵ<EFBFBD><CFB5>
float a0 = cmod5_c[1] + cmod5_c[2] * x + cmod5_c[3] * x * x + cmod5_c[4] * x * x * x;
float a1 = cmod5_c[5] + cmod5_c[6] * x;
float a2 = cmod5_c[7] + cmod5_c[8] * x;
float gamma = cmod5_c[9] + cmod5_c[10] * x + cmod5_c[11] * x * x;
float s0 = cmod5_c[12] + cmod5_c[13] * x;
float s = a2 * U;
// <20><><EFBFBD>㱥<EFBFBD>ͺ<EFBFBD><CDBA><EFBFBD>
float f;
if (s < s0) {
float g = 1.0f / (1.0f + expf(-s0)); // ʹ<><CAB9>expf
float alpha = s0 * (1.0f - g);
f = powf(s / s0, alpha) * g; // ʹ<><CAB9>powf
}
else {
f = 1.0f / (1.0f + expf(-s));
}
// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ͬ<EFBFBD>Բ<EFBFBD><D4B2><EFBFBD>
float b0 = powf(f, gamma) * powf(10.0f, a0 + a1 * U); // ʹ<><CAB9>powf
// <20><><EFBFBD><EFBFBD>һ<EFBFBD>׵<EFBFBD><D7B5><EFBFBD>ϵ<EFBFBD><CFB5>
float b1_term = cmod5_c[15] * U * (0.5f + x - tanhf(4.0f * (x + cmod5_c[16] + cmod5_c[17] * U)));
float b1 = (cmod5_c[14] * (1.0f + x) - b1_term) / (expf(0.34f * (U - cmod5_c[18])) + 1.0f);
// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>׵<EFBFBD><D7B5><EFBFBD>ϵ<EFBFBD><CFB5>
float v0 = cmod5_c[21] + cmod5_c[22] * x + cmod5_c[23] * x * x;
float d1 = cmod5_c[24] + cmod5_c[25] * x + cmod5_c[26] * x * x;
float d2 = cmod5_c[27] + cmod5_c[28] * x;
float y = (U / v0 + 1.0f);
float v2;
if (y < y0) {
v2 = a + b * powf(y - 1.0f, n); // ʹ<><CAB9>powf
}
else {
v2 = y;
}
float b2 = (-d1 + d2 * v2) * expf(-v2);
// <20><><EFBFBD><EFBFBD>ģ<EFBFBD><C4A3><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ɢ<EFBFBD><C9A2>ϵ<EFBFBD><CFB5><EFBFBD><EFBFBD><EFBFBD><EFBFBD>ʵ<EFBFBD><CAB5>ֵ<EFBFBD>Ƚ<EFBFBD>
float angle_rad = angle * 3.14159265f / 180.0f;
float sigma0_model = 10.0f * log10f(b0 * powf(1.0f + b1 * cosf(angle_rad) + b2 * cosf(2.0f * angle_rad), 1.6f));
if (sigma0_model - sigma0 >= 0.0f) {
return U;
}
}
return 26.0f; // <20><><EFBFBD><EFBFBD>δ<EFBFBD>ҵ<EFBFBD><D2B5><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
}
__global__ void CUDA_WindSpeedInversered_Kernel(float* in_wind_dir_U_arr_cuda,
float* in_wind_dir_V_arr_cuda,
float* in_inc_angle_arr_cuda,
float* in_sar_sigma0_arr_cuda,
float* out_wind_speed_arr_cuda, int64_t pixelCount) {
int64_t idx = blockIdx.x * blockDim.x + threadIdx.x;
if (idx < pixelCount) {
float sar_sigma0 = in_sar_sigma0_arr_cuda[idx];
if (isinf(sar_sigma0) || isnan(sar_sigma0)) {
out_wind_speed_arr_cuda[idx] = -1;
}
else {
float U10 = in_wind_dir_U_arr_cuda[idx];
float V10 = in_wind_dir_V_arr_cuda[idx];
float inc_angle = in_inc_angle_arr_cuda[idx];
float wind_dir_angle = 270 - atan2(V10, U10) * 180.0 / PI;
wind_dir_angle = fmodf(wind_dir_angle, 360.0f);
float out_wind_speed = windspeedcmd5_device(wind_dir_angle, inc_angle, sar_sigma0);
//printf("windspeedcmd5(%f, %f, %f)-%f \n", wind_dir_angle, inc_angle, sar_sigma0, out_wind_speed);
out_wind_speed_arr_cuda[idx] = out_wind_speed;
}
}
}
void WindSpeedInverseredFunc(float* in_wind_dir_U_arr_cuda, float* in_wind_dir_V_arr_cuda, float* in_inc_angle_arr_cuda, float* in_sar_sigma0_arr_cuda, float* out_wind_speed_arr_cuda, int64_t pixelCount)
{
// <20><><EFBFBD><EFBFBD> CUDA <20>˺<EFBFBD><CBBA><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>Ϳ<EFBFBD><CDBF>ijߴ<C4B3>
int64_t blockSize = 256; // ÿ<><C3BF><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>߳<EFBFBD><DFB3><EFBFBD>
int64_t numBlocks = (pixelCount + blockSize - 1) / blockSize; // <20><><EFBFBD><EFBFBD> pixelcount <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>С
// <20><><EFBFBD><EFBFBD> CUDA <20>˺<EFBFBD><CBBA><EFBFBD>
CUDA_WindSpeedInversered_Kernel << <numBlocks, blockSize >> > (in_wind_dir_U_arr_cuda, in_wind_dir_V_arr_cuda, in_inc_angle_arr_cuda, in_sar_sigma0_arr_cuda, out_wind_speed_arr_cuda, pixelCount);
#ifdef __CUDADEBUG__
cudaError_t err = cudaGetLastError();
if (err != cudaSuccess) {
printf("CUDAdistanceAB CUDA Error: %s\n", cudaGetErrorString(err));
exit(2);
}
#endif // __CUDADEBUG__
cudaDeviceSynchronize();
}